# Equations for Steady 1D Isentropic Flow

The equations used to describe steady 1D isentropic flow are derived from conservation of mass, momentum, and energy, as well as an equation of state (typically the ideal gas law).

These equations are typically described as ratios between the local static properties (p, T, $\rho$) and their stagnation property as a function of Mach number and the ratio of specific heats, $\gamma$. Recall that Mach number is the ratio between the velocity and the speed of sound, a.

These ratios are given here:

Temperature: $T_o/T = \left(1+\frac{\gamma -1}{2} M^2\right)$

Pressure: $P_o/P = \left(1+\frac{\gamma -1}{2} M^2\right)^{\frac{\gamma}{\gamma-1}}$

Density: $\rho_o/\rho = \left(1+\frac{\gamma -1}{2} M^2\right)^{\frac{1}{\gamma-1}}$

In addition to the relationships between static and stagnation properties, 1D nozzle flow offers an equation regarding the choked cross-sectional flow area (recall that the flow is choked when M=1.)

$A/A^* = \frac{1}{M}\left(\left(\frac{2}{\gamma+1}\right)\left(1+\frac{\gamma -1}{2} M^2\right)\right)^{\frac{\gamma+1}{2\left(\gamma-1\right)}}$

Some excellent references for these equations are:

• Gas Dynamics Vol. I – Zucrow and Hoffman – 1976
• Gas Dynamics – John and Keith – 2nd Ed. – 2006

# Turbulent Zero Pressure Gradient Flat Plate – simpleFOAM, komegaSST

An excellent test case, and case to familiarize yourself with some of the turbulence models available in OpenFOAM is a 2D flat plate with zero pressure gradient. I will solve this problem using the solver simpleFoam and the komegaSST model.

Here are the sections of this post:

1. Quick Overview: kω-SST (komegaSST) Boundary Conditions
2. Case set-up and mesh
3. Results
• u+ vs. y+
• Coefficient of friction
4. Conclusions and useful resources

## Quick Overview: kω-SST (komegaSST) Boundary Conditions

In this section I will describe the boundary set-up for komegaSST  where no wall functions are implemented. This requires that the y+ along the wall is less than or equal to one. For the komegaSST turbulence model the boundary conditions are as follows:

At the wall:

• ω (omega) – specific dissipation rate
• BC type: fixedValue
• BC value: $\omega_{wall}=10\frac{6\nu_\infty}{\beta_1\left(\Delta y_{wall}\right)^2}$
• k – turbulent kinetic energy
• BC type: fixedValue
• BC value: 0
• nut – turbulent viscosity
• BC type: fixedValue
• BC value: 0

In the free-stream:

• ω (omega) – specific dissipation rate
• BC type: fixedValue
• BC value: $\frac{U_\infty}{L} < \omega_{\infty} < \frac{10 U_\infty}{L}$
• k – turbulent kinetic energy
• BC type: fixedValue
• BC value:$\frac{10^{-5}U_\infty^2}{Re_L} < k_{\infty} < \frac{0.1U_\infty^2}{Re_L}$
• nut – turbulent viscosity
• BC type: calculated
• BC value: 0 (this is just an initial value)

where $\beta_1=0.075$, and $\Delta y_{wall}$ is the normal distance from the wall to the first cell center.

## Case set-up and mesh

### Free-Stream Properties

For this case I have followed a similar set up to the 2D flat plate case used on the NASA turbulence modelling resource website. By doing this it gives me something to compare to! However, I am not going to use the exact set-up on that website. Since we are using simpleFoam, I am going to set up this case the way I prefer which is to use a velocity (U) of 1 m/s and scale all other properties accordingly. The simulation properties that I used are :

• U∞=1 m/s
• ν=4e-7 m2/s
• L=2 m

These correspond to a Reynolds number at L=2m of 5 million .

### Grid Generation

The grid used was generated in blockMesh. High inflation was used in the boundary layer region in order to achieve the desired y+ value of less than 1. For more details on grid generation using blockMesh see the OpenFOAM manual!

### Boundary Conditions

For the incompressible solver simpleFoam, the minimum boundary conditions required for a simulation are p and U. However, if the simulation is a RANS simulation additional boundary conditions are required. For the kω-SST model we need to have a boundary condition on k and ω as well. The boundary conditions I defined in the zero (0) folder can be found in the attached tutorial file.

The only boundary condition that really needs any comment is omega. We calculate omega using $\omega_{wall}=10\frac{6\nu_\infty}{\beta_1\left(\Delta y_{wall}\right)^2}$. In our case the wall distance to the first cell center is 5.48316E-06. Using our free stream viscosity of 4E-7 this gives a value of omega at the wall of 10643595.39.

### Tip for fvSolution

If you find that the results you are getting are wrong, it could be that the residuals for the different properties are too high! Certain properties converge before others and therefore you need to ensure that they all converge to a sufficiently low value!

## Results

First we compare the coefficient of friction to the .dat file available from the .DAT available from the NASA Turbulence modelling resource. NOTE: We simulated at Reynolds number of 5 million whereas the NASA setup is at 10 million. So the x coordinate in the following is plots is rescaled and in fact we are only using half of the data from the .DAT file!

We can see from the figure that the coefficient of friction from our simulation matches the expected data closely! Hurray!

Next let’s compare the u+ vs. y+ profile to the  universal profile for turbulent boundary layers. Recall that u+ is the velocity normalized by the friction velocity ($u_*$), and y+ is given by the following equation:

$y^+ = \frac{u_* y}{\nu}$

The u+ is vs y+ is plotted here:

We can see from the figure that our solution is pretty good! The y+ value of the first node is located around a y+ of approximately 0.5, the viscous sublayer matches very closely and the log law layer is not significantly off!

## Conclusions

In this post we simulated a zero pressure gradient flat plate at a Reynolds number of 5 million. We compared the results for shear stress to the NASA turbulence modelling resource expected results and showed close agreement. Then the u+ vs y+ profile was compared to the universal law of the wall and again the results were okay!

## Some useful references

The NASA Turbulence Modelling Resource, http://turbmodels.larc.nasa.gov/

# Compressible Aerodynamics Calculator – Matlab

Here I have written a Matlab App for calculating compressible aerodynamics using the basic gas dynamics equations. The calculations available are:

• Isentropic Flow (including Prandtl-Meyer flow)
• Normal Shocks (See article)
• Oblique Shocks
• Fanno Flow
• Rayleigh Flow

To install simply download the file from the Matlab File Exchange here:

http://www.mathworks.com/matlabcentral/fileexchange/56152-compressiblecalccf?s_tid=srchtitle

Please let me know if you find it useful! Or (and especially) if you find any mistakes or bugs etc!

# Newton-Raphson Method

I was recently asked by a class to go over the Newton-Raphson method for solving non-linear equations. Specifically in this case it was to solve 1D gas dynamics equations. Here I will just do a brief overview of the method, and how its used. I will solve two cases, one where the derivative of the function is known, and one where the derivative of the function should be approximated.

## Contents

Here are the contents you’ll find in this post:

1. Introduction to the Newton-Raphson Method
2. Using the Method
1. Example 1 (derivative of the function is known)
• Calculating Mach number from critical area ratio (M from A/A*)
2. Example 2 (derivative of the function is unknown… or to annoying to derive)
• Calculating incident shock pressure ratio from diaphragm pressure ratio
3. Conclusions and Useful References

## Introduction to the Newton-Raphson Method

What is the Newton-Raphson method? Basically it is an iterative approach for solving the roots of functions. There are tons of these. You probably don’t need to know all of them (just pick a few that work for you!) Typically I stick to the Newton-Raphson method and the bisection method and I rarely have problems.

The Newton-Raphson method is used when you have some function f(x) and you want to find the value of the dependent variable (x) when the function equals zero. AKA you want to find the roots of the equation. If you have an initial guess at some point, $x_i$, the tangent can be extended to some point that crosses 0 at an easily calculable point $x_{i+1}$. This point gives an improved estimation of the root. The basics of the method come from the fact that the first derivative is equivalent to the slope and therefore if you know it, you can calculate the tangent intersection:

$f'(x_i)=\frac{f(x_i)-0}{x_i-x_{i+1}}$  (1)

Which we can rearrange to get the “Newton-Raphson Formula”:

$x_{i+1}=x_{i}-\frac{f(x_{i})}{\frac{df}{dx}(x_{i})}$  (2)

The method is depicted in the figure below:

Equation 2 pretty much sums up the method. If you start with a known derivative and function value $x_i$ you can calculate a new prediction $x_{i+1}$. By repeating this process (perform iterations) better and better approximations of the value for x are obtained.

## Using the Method

### Example 1 (derivative of the function is known):

Lets say for this example that we know a critical area ratio (let’s say A/A*=3.3, with $\gamma=1.4$) at a point in the flow field, and we want to use it calculate Mach number. For this we will solve the following equation:

$\frac{A}{A^*}=\frac{1}{M}\left[\left( \frac{2}{\gamma+1}\right) \left(1 +\frac{\gamma -1}{2} M^2\right)\right]^{\frac{\gamma+1}{2\left(\gamma-1\right)}}$  (3)

This equation is equation for the critical area ratio for a given Mach number. Sometimes however the A/A* is known and the Mach number is desired. Let’s say A/A* =3.3 what is the Mach number? To do this, this equation can be easily solved using the Newton-Raphson approach.

Essentially the Newton-Raphson method is a root finding method. This means we are finding where an equation crosses zero. So how do we apply this to the problem above? We aren’t interested in where the equation above crosses zero, we are interested in where it crosses 3.3! Well it’s easy to change the above equation to give:

$f(M)= \left[\left( \frac{2}{\gamma+1}\right) \left(1 +\frac{\gamma -1}{2} M^2\right)\right]^{\frac{\gamma+1}{2\left(\gamma-1\right)}} - M\frac{A}{A^*} = 0$  (4)

So now the problem is in a form that we can use the Newton-Raphson approach. In general, the approach is written as:

$x_{new}=x_{old}-\frac{f(x_{old})}{\frac{df}{dx}(x_{old})}$  (5)

Here, f(x) is the equation that is to be solved (in our case Equation 4), x is the variable to be solved for (in our case Mach number) and df/dx is the derivative of the equation with respect to the unknown variable (in our case the derivative of Equation 4 with respect to Mach number).

So at this point we know f(M), but we must calculate the derivative of f(M). The result is the following:

$\frac{df}{dM}(M)= \frac{A}{A^*}-(\frac{2}{\gamma+1})^{\frac{\gamma+1}{2\left(\gamma-1\right)}-1}M\left(1+\frac{(\gamma-1)M^2}{2}\right)^{\frac{\gamma+1}{2\left(\gamma-1\right)}-1}$  (6)

The next step is to apply the iterative algorithm. Essentially, you have to make an intelligent initial guess for the unknown variable (Mach number). In our case let’s make our first guess M=3. In our first step this is $M_{old}$. Using this value we can get a guess for $M_{new}$. Then at our next iteration $M_{new}$ becomes $M_{old}$ and this process is repeated until convergence ($M_{old} \approx M_{new}$).

The following is a table showing the convergence of this problem:

So the Mach number for an A/A* of 3.3 is 2.738! But what if we had chosen a different starting point? Say… 0.5? Well then the answer looks like…

Wait why is the answer different?? Dont forget! There are two solutions to this problem, a subsonic and a supersonic solution. When we started at M=3 we converged on the supersonic solution (M=2.73805). When we start at 0.5 we converged on the subsonic solution (M=0.17875).

### Example 2 (derivative df/dx is unknown… or too annoying to derive):

For example 2 we are again going to solve a simple 1D gas dynamics problem using the Newton-Raphson method. For this problem, let’s say that we are given a diaphragm pressure ratio for a shock tube ($P_4/P_1=40$). We want to solve for the pressure ratio that will exist  across the normal shock that is created when the diaphragm ruptures. For this problem see John and Keith (2006) for an excellent chapter on the shock-tube problem.

Basically this problem boils down to solving the following equation

$\frac{P_4}{P_1}=\frac{P_2}{P_1}\left[1-\frac{(\gamma_4-1)(\frac{a_1}{a_4})(\frac{P_2}{P_1}-1)}{\sqrt{2\gamma_1}\sqrt{2\gamma_1+(\gamma_1+1)(\frac{P_2}{P_1}-1)}}\right]^{-\frac{2\gamma_4}{\gamma_4-1}}$          (7)

for a $\frac{P_2}{P_1}$ with a known shock tube pressure ratio $\frac{P_4}{P_1}$. Equation 5 is not the type of equation that is very attractive to differentiate, so in this example we will use a simple central different to approximate the derivative. Essentially what this means is that instead of directly calculating the derivative at each of our iterations, we will use a point slightly to the right and slightly to the left in order to approximate the slope.

First of all we rearrange (5) to give us an equation that equals zero (so that we can find its roots).  Demonstrating this is easier if we set our unknown $p=P_2/P_1$.

$f(p)=\frac{P_4}{P_1}\left[1-\frac{(\gamma_4-1)(\frac{a_1}{a_4})(p-1)}{\sqrt{2\gamma_1}\sqrt{2\gamma_1+(\gamma_1+1)(p-1)}}\right]^{\frac{2\gamma_4}{\gamma_4-1}}-p=0$ (8)

For the Newton-Raphson the above equation would be solved using:

$p_{new}=p_{old}-\frac{f(p_{old})}{\frac{df}{dp}(p_{old})}$  (9)

Since we do not know the derivative of f approximate it with a central difference (centered approximation of the local slope at that point):

$\frac{df}{dp}(p_{old})=\frac{f(p_{old}+\Delta p) -f(p_{old}-\Delta p)} {2 \Delta p}$ (10)

so:

$p_{new}=p_{old}-2\Delta p \frac{f(p_{old})}{f(p_{old}+\Delta p) -f(p_{old}-\Delta p)}$  (11)

With the above iterative equation, the problem is solved in the same fashion as Example 1. See the following table:

Using Newton-Raphson we have now calculated that for a diaphragm pressure ratio of 40, the pressure ratio across the shock ($P_2/P_1$) should be 4.7726.

## Conclusion and Useful References

Hopefully this was helpful. The Newton-Raphson method is extremely useful. However be careful! There are some cases where this method doesn’t work very well. For example if an equation has multiple roots, your initial guess must be fairly close to the answer you are looking for or you could get the completely wrong root. It often helps to plot your function beforehand or after so that you know what it looks like, and it also helps to never blindly trust the answer. If it isn’t physical… it probably isn’t the one you were looking for.

Here are some useful references

[1] https://en.wikipedia.org/wiki/Newton%27s_method

[2] John, J. E. A., & Keith, T, G. (2006) Gas Dynamics, 3rd Edition, Pearson Prentice Hall

[3] Chapra, S. C., Canale, R. P. (2010) Numerical Methods for Engineers (6th Ed.), McGraw-Hill

# Stationary Normal Shock Waves

Shock-waves are a common topic covered in undergraduate thermodynamics, fluid mechanics, and gas dynamics courses. Because shock-waves are awesome, I figured I would do one of my first posts on their basic theory. I’m going to try to cover the topic by first providing some context and a clear picture of what a shock wave is (and why it forms)… followed by how we analyze it. I try to include general features of shockwaves and the things that are particularly “important to know” … selon moi.

## Contents

Here are the sections you’ll find in this post:

1. Introduction to the normal shock
• Simple example: shock-wave in diverging nozzle
2. Analysis of the normal shock
3. Important properties of a normal shock
4. Conclusions and Useful References

## Introduction to the stationary normal shock

So first, what is a shock wave? Well, in order to understand this, you must first recall that in supersonic flow the flow is moving faster than the speed of sound. Also, recall that the speed of sound is the speed at which information is transmitted throughout the continuum. Thus, if you are a small person riding along a fluid particle that is supersonic, neither you nor the fluid particle you are hitching a ride on can see what is coming up ahead of you. There is no mechanism for you to see what is ahead and gradually slow down or correct course. You are now committed to a sudden, discontinuous change resulting in an increase in pressure and temperature, and a loss of speed. So, basically the idea of shockwaves is this: they are a discontinuity that forms in order for the flow to meet some downstream condition. If a supersonic flow never encounters something downstream, like an obstacle or a back pressure, it will never shock. And likewise, a subsonic flow can never form a shock wave (unless it is first accelerated to supersonic somewhere within the flow field). In subsonic flows, downstream obstacles and pressures are communicated upstream at a speed that allows for the curving of streamlines and/or deceleration. You can observe this by picturing the streamlines around an airfoil that curve around its shape but never contact it.

#### Simple example of a shock-wave

Let’s look at (probably) the most common learning example: a shockwave forming in a supersonic nozzle.

In Figure 1 we see the classic nozzle example. Hopefully this example is not extremely new to you. The main concept is this: When the back pressure is high enough, the flow through the nozzle never goes supersonic (including at the throat) and is simply compressed to the throat and then expanded in the diverging section (and in these cases Pe=Pb). These are the cases shown by Region a in the figure. When the back pressure (Pb) is low enough that the flow in the nozzle throat is choked (Mach 1 is achieved) there are now two possible solutions. The first solution is the subsonic solution. Here, the flow decelerates from Mach 1 to the exit and flow exits at M<1. The second solution is the supersonic solution. Here, the flow continues to accelerate and exits at M>1. So, which solution is the right solution? From the graph in Figure 1, you can probably guess… the back pressure! When Pb is equal the supersonic back pressure, the flow will obey the supersonic solution. And similarly, when Pb is equal to the subsonic back pressure the flow will remain subsonic.

But, what if the back pressure is somewhere in-between? The back pressure is too high for the supersonic solution… but too low for the subsonic solution? This situation is what exists in Region b. This is where the concept of the shockwave is really apparent. When the back pressure is in this region, the flow starts out by accelerating to supersonic after the nozzle throat. Once the flow is supersonic, it cannot receive any communication from the flow ahead of it. Why? Because (again) communication between fluid particles is at the speed of sound! So here we have a fluid particle being accelerated in a diverging section, with a downstream boundary condition that it cannot meet (because it can’t see it coming). So, what happens? A normal shock is produced at some point in the diverging section of the nozzle. The resulting scenarios are shown in the next figure.

Fig 2 shows what happens in Region b. As the back pressure is increased, the shock-wave that is formed moves upstream and is formed earlier in the nozzle. This continues until Pe=Pb(subsonic) and the shock disappears. Can you have a shock in the converging section?  No. Not unless you had supersonic flow in the converging section in the first place.

Because I love CFD… To illustrate this example I have simulated it using OpenFOAM and the solver rhoCentralFoam (hopefully I’ll write post about it’s setup at some point). As you can see the flow accelerates from the reservoir (through the converging-diverging nozzle) but due to the back pressure, a shock-wave forms in the diverging section.

## Analysis of Normal Shock-waves

So how do we analyze a shockwave? Well, normal shock waves are a one-dimensional, adiabatic, discontinuous phenomena that are governed by the equations of fluid mechanics (conservation of mass, momentum and energy). Thus, it makes sense that we simply analyze them with a 1D control volume analysis. We will have an inlet (side 1) and an outlet (side 2) and in between is the (practically) discontinuous normal shock wave.

The governing equations for the control volume above are:

Conservation of mass:       $\rho_2 U_2 =\rho_2 U_2$

Conservation of momentum:       $p_2+\rho_2 U_2^2 =p_1+\rho_1 U_1^2$

Conservation of energy:       $h_2+\frac{U_2^2}{2}=h_1+\frac{U_1^2}{2}$

Perfect gas equation of state assuming constant specific heats:

$h=\int_0^T c_p dT \, = \left(\frac{\gamma R}{\gamma-1}\right)T$

and       $p = \rho RT$

With the above equations the problem is fully defined. It is then just a matter of an exercise in algebra. From various combinations of the above equations (its worth it to derive them yourself at some point… especially if you are studying for a PhD candidacy)  the normal shock relations are obtained.

The normal shock relations are defined here:

• Temperature Ratio:     $\frac{T_2}{T_1}=\frac{\left(1+\frac{\gamma-1}{2}M_1^2\right)\left(\frac{2\gamma}{\gamma-1}M_1^2-1\right)}{\left[\frac{\left(\gamma+1\right)^2}{2\left(\gamma-1\right)}M_1^2\right]}$
• Pressure Ratio:        $\frac{p_2}{p_1}=\frac{2\gamma M_1^2}{\gamma+1}-\frac{\gamma-1}{\gamma+1}$
• Density and velocity ratio:       $\frac{\rho_2}{\rho_1}=\frac{U_1}{U_2}=\frac{\left(\gamma+1\right)M_1^2}{\left(\gamma-1\right)M_1^2+2}$
• Stagnation Pressure Ratio:        $\frac{p_{o2}}{p_{o1}}=\left[\frac{\left(\gamma+1\right)M_1^2}{\left(\gamma-1\right)M_1^2+2}\right]^\frac{\gamma}{\gamma-1}\left[\frac{\gamma+1}{2\gamma M_1^2-\left(\gamma-1\right)}\right]^\frac{1}{\gamma-1}$
• Post-shock Mach number:       $M_2^2=\frac{M_1^2+\frac{2}{\gamma-1}}{\frac{2\gamma}{\gamma-1}M_1^2-1}$

## Properties of the Normal Shock

There is lots to know about shock-waves, but here is a summary of what I think is  important to know and what sticks out in my mind. Hopefully I’m not missing anything… but I can always add stuff later!

• Normal shock-waves obey conservation of mass, momentum, and energy, AND the 2nd law of thermodynamics. It may seem obvious but it should be re-stated. A shock-wave obeys the usual laws of fluid mechanics.
• Shock-waves are adiabatic. By this I mean to say that when we analyze the control volume shown in Figure 3, Energy in=Energy out.
• When you have a perfect gas with constant specific heats stagnation temperature does not change across the shock. Sometimes students forget this and it’s quite important. This is really an extension of the last point. Note: As stated, when the gas does not have constant specific heats, energy is still conserved but the stagnation temperature across the shock will change.
• Stagnation pressure always decreases across a shock.
• The Mach number always decreases across a shock. Additionally, for a normal shock the post-shock Mach number is always subsonic. However as you will see later if I make a post on oblique shocks, the post-shock Mach number can remain supersonic.

## Conclusions and Useful References

Well I hope you found this post helpful! … And not confusing. As usual if someone spots a mistake or something I missed then I would very much appreciate a comment and I will fix it. I am not perfect obviously, and this is both an exercise in sharing knowledge, but an exercise in strengthening and confirming knowledge through the act of sharing… or something like that.

There is no shortage of books that cover normal shocks. But I have found that not all are created equal. These are the books that I are always sitting on my desk and have great sections covering this topic:

[1] Anderson, J. D. (1990). Modern compressible flow: with historical perspective (Vol. 12). McGraw-Hill.

[2] John, J. E. A., & Keith, T, G. (2006) Gas Dynamics, 3rd Edition, Pearson Prentice Hall

[3] Çengel, Y. A., & Boles, M. A. (2015). Thermodynamics: an engineering approach. M. Kanoğlu (Ed.). McGraw-Hill Education.

[4] Zucrow, M. J., & Hoffman, J. D. (1976). Gas dynamics. New York: Wiley, 1976, 1.