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

If you downloaded it prior to March 28th, 2016 please reinstall the app! I had mislabelled some properties.

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

Logo1

 

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:

NewtonRaphson
Fig 1) Depiction of Newton-Raphson Method

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:

NewtonRaphson

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…

NewtonRaphson2

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:

NewtonRaphson_CD

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.

Also the wikipedia page on this is pretty good… haha

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.

Nozzle_noshock
Fig 1) Nozzle flow schematic

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.

Nozzle_shocks
Fig 2) Nozzle flow shematic (with shocks)

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.

movie
Nozzle flow with shock in diverging section (simulated in rhoCentralFoam – See Related Blog Post Under OpenFOAM)

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.

CV
Fig 3) Normal Shock Control Volume Layout

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.

 

 

Solving the cavity flow problem using the streamfunction-vorticity formulation

In this post I am going to write a (hopefully) simple code in matlab to solve the cavity flow problem using the vorticity stream function formulation. I will the compare the result to the result calculated by the OpenFOAM solver, icoFoam.  Writing your own solvers is fun, rewarding, and is a practice that really cements some of the fundamental knowledge of CFD.

Contents

  1. Introduction to the Cavity Flow Problem
  2. Governing Equations
  3. Finite Differences and Solver Algorithm
    • Understanding and setting the boundary conditions
    • Semi-discretization of the governing equations
    • Solving the vorticity equation in time
    • Solving for the vorticity at the boundaries
    • Solving for the stream function
    • Putting it all together
  4. Results/Comparison to icoFoam
    • Re=10
    • Re=100
  5. Conclusion and some References

Introduction to the Cavity Flow Problem

The cavity flow problem is described in the following figure. Basically, there is a constant velocity across the top of the cavity which creates a circulating flow inside. To simulate this there is a constant velocity boundary condition applied to the lid, while the other three walls obey the no slip condition. Different Reynolds numbers give different results, so in this post I’ll simulate Re=10, and Re=800. At high Reynolds numbers we expect to see a more interesting result with secondary circulation zones forming in the corners of the cavity.

cavityflowbasic

Continue reading Solving the cavity flow problem using the streamfunction-vorticity formulation