Pressure Driven Nozzle Flow with Shock – rhoCentralFoam

In this post I will go over the set up and solution of a pressure driven nozzle flow with a shock located in the diverging section. This refers to the type of flow problem described by region b in my page covering stationary normal shock-waves.

Here is a summary of the contents of this post:

  1. Nozzle Design
  2. Mesh Generation
  3. Boundary and Initial Conditions
  4. Constants and Schemes
  5. Running the Simulation/Results
  6. Conclusion

The OpenFOAM file that is ready to download and run can be downloaded from here:

Nozzle with Shock Tutorial File Download

Simulation summary:

Table

Nozzle Design Calculation

We want to simulate the flow in a converging-diverging nozzle, with a shock at a specified location in the diverging section. Based on the nozzle geometry we create we must calculate the required back pressure. This boundary condition will generate the shock in our nozzle.

Here is the nozzle geometry we will work with:

  • Known nozzle height at through (A*=A_{throat})
  • Known nozzle height at the outlet (A_e/A_{throat}=2, the outlet height is twice the throat height)
  • Shock located when flow reaches M=1.88 (corresponding to A/A*=1.53)
  • Stagnation Properties, Po=3 Bar, To=298K
  • Fluid is air (γ=1.4)

To know how to set the back pressure, we must calculate the flow through the nozzle. This can be easily calculated by hand, using tables, or by using the Matlab app available from curiosityFluids.

The shock is to occur when the Mach number is 1.88. Therefore, we can calculate the preshock and post shock properties easily.

P_1=\frac{P_1}{P_o}P_o = 0.1539 (3000000 \textnormal{ Pa}) = 0.462  \textnormal{ MPa}

P_2=\frac{P_2}{P_1}P_1 = 3.9568 (461700  \textnormal{ Pa}) = 1.83 \textnormal{ MPa}

After the shock, the (now) subsonic flow is expanded out of the diverging section of the nozzle. We need to calculate the Mach number and pressure at the exit.

Recall that A_{throat}=A_1^*

\frac{A_e}{A_2^*}=\left(\frac{A_1^*}{A_2^*}\right)\left(\frac{A_e}{A_{throat}}\right)= (0.7765)(2) = 1.55

therefore

M_e=0.41, \frac{P_e}{P_{o2}}=0.89

and

P_e=\left(\frac{P_e}{P_{o2}}\right)\left(\frac{P_{o2}}{P_{o1}}\right)P_{o1} =\left(0.89\right)\left(0.7765\right) 3\times 10^6 \textnormal{ Pa} = 2.0674 \textnormal { MPa}

This is the approximate pressure at the outlet! This will be used for our downstream boundary condition!

Mesh Generation

To build the mesh we will use the native openFoam meshing utility blockMesh. For a detailed overview of the blockMesh utility it is best to start at with the OpenFoam user guide. Block mesh uses definitions of vertices, blocks, edges, and faces to define the mesh.

Vertices

Vertices define the shape of the boundary. I find it helpful to number the vertices in the file, and draw a diagram:

vertices
Fig: Vertex Numbering
vertices
(
 (-0.6 0 -0.005) //0
 (-0.2121321 0 -0.005) //1
(0 0 -0.005) //2
 (2 0 -0.005) //3
(-0.6 0.45 -0.005) //4
 (-0.3 0.45 -0.005) //5
(0 0.15 -0.005) //6
 (2 0.3 -0.005) //7
(-0.6 0 0.005) //8
 (-0.2121321 0 0.005) //9
(0 0 0.005) //10
 (2 0 0.005) //11
(-0.6 0.45 0.005) //12
 (-0.3 0.45 0.005) //13
 (0 0.15 0.005) //14
 (2 0.3 0.005) //15 
 (-0.2121321 0.237867965 -0.005) //16
 (-0.2121321 0.237867965 0.005) //17
(-0.6 0.237867965 -0.005) //18
 (-0.6 0.237867965 0.005) //19
 (-4 0 -0.005) //20
 (-4 0.237867965 -0.005) //21
 (-4 0.45 -0.005) //22
 (-4 3.85 -0.005) //23
 (-0.6 3.85 -0.005) //24
 (-0.3 3.85 -0.005) //25
(-4 0 0.005) //26
 (-4 0.237867965 0.005) //27
 (-4 0.45 0.005) //28
 (-4 3.85 0.005) //29
 (-0.6 3.85 0.005) //30
 (-0.3 3.85 0.005) //31
);

Blocks

The blocks are used to define the different regions within the mesh. The blocks are defined here in a way to allow mesh inflation into the reservoir where the interesting flow is not occurring.

blocks
Fig: Block Number

The blocks are defined below. The blocks with “simpleGrading” defined are the blocks that have an inflation applied to them.

blocks
(
 hex (9 8 19 17 1 0 18 16) (20 20 1) simpleGrading (1 1 1)
 hex (17 19 12 13 16 18 4 5) (20 20 1) simpleGrading (1 1 1)
 hex (10 9 17 14 2 1 16 6) (20 20 1) simpleGrading (1 1 1)
 hex (11 10 14 15 3 2 6 7) (150 20 1) simpleGrading (1 1 1)
 hex (12 28 29 30 4 22 23 24) (40 40 1) simpleGrading (10 10 1)
 hex (8 26 27 19 0 20 21 18) (40 20 1) simpleGrading (10 1 1)
 hex (19 27 28 12 18 21 22 4) (40 20 1) simpleGrading (10 1 1)
 hex (13 12 30 31 5 4 24 25) (20 40 1) simpleGrading (1 10 1)
);

Edges

We need to use an edge definition here in order to make a smooth subsonic inlet into the nozzle.

Edges
Fig: Edges

Below is the edge definition from the blockMeshDict:

edges
(
 arc 5 16 (-0.295442325 0.397905546 -0.005)
 arc 16 6 (-0.052094453 0.154557674 -0.005)
 arc 13 17 (-0.295442325 0.397905546 0.005)
 arc 17 14 (-0.052094453 0.154557674 0.005)
);

Boundaries

The face definition is very important as it defines the locations of all the pertinent boundaries, names them, and sets their boundary type. The faces are shown in the figure below. Notice that there is no definition for the front and back faces. This is because if you leave them undefined blockMesh adds them to a group of faces called defaultFaces which we can then later define as empty.

faces
Fig: Faces

The boundary definitions from the blockMeshDict are shown below:

boundary
(
 inlet1
 {
 type patch;
 faces
 (
 (20 26 27 21)
 (21 27 28 22)
 (22 28 29 23)
 );
 }
 top
 {
 type wall;
 faces
 (
 (23 29 30 24)
 (24 30 31 25)
 );
 }
outlet
 {
 type patch;
 faces
 (
 (11 3 7 15)
 );
 }
 bottom
 {
 type symmetryPlane;
 faces
 (
 (8 0 1 9)
 (9 1 2 10)
 (10 2 3 11)
 (26 20 0 8)
 );
 }
 nozzle
 {
 type wall;
 faces
 (
 (5 13 17 16)
 (16 17 14 6)
 (6 14 15 7)
 (13 5 25 31)
 );
 }
);

Build Mesh

To build the mesh simply type blockMesh into the command terminal from inside you case file. To ensure the mesh built properly, then run checkMesh .

Boundary and Initial Conditions

Once the mesh is built, it is time to set up the case. The boundary and initial conditions are defined in the “0” folder. The solver for this case is rhoCentralFoam which requires P, T, U definitions in the 0 folder. A summary of the boundary conditions and the mesh resulting from the previous section is shown in the next picture:

setup
Fig: Boundary Condition Set-up and Mesh

Initial Conditions

Here I/we have initialized the internalField the simulation with a small x-velocity (7 m/s), pressure equal to the stagnation pressure, and temperature equal to 298 K.

Boundary Conditions

There are some key aspects of this setup that are essential to setting up this problem. The first is the totalPressure  definition for the inlet pressure. This boundary condition maintains the inlet at a constant stagnation pressure. This allows this boundary to behave like an infinite reservoir. To complement this boundary condition we use the zeroGradient boundary conditions for velocity. By setting velocity to zeroGradient, we are allowing flow to be induced by the nozzle. If you set a fixedValue for the velocity, you would be constraining the amount of mass flow into the system, and therefore would require some extra set-up calculations. Additionally, you would be making the use of the totalPressure boundary condition obsolete. For the temperature we set a fixed value at 298 K. This maintains the flow induced from the nozzle at a constant temperature.

For the outlet we want to ensure that the correct back pressure is maintained. This will make sure that we have a shock forming in the nozzle. Temperature and velocity are set to zeroGradient and therefore will be determined by the nozzle flow.

The symmetry plane is set by setting p, T, and U to be type symmetry. 

The nozzle boundary should be set to type slip for velocity. This ensures that flow tangency is maintained for velocity (but we are not enforcing the no-slip condition as this is an inviscid simulation). For pressure, and temperature zeroGradient should be used. This makes the wall adiabatic (i.e. no heat flux) and makes the pressure calculated by the flow.

Below is a summary of p, T, and U files that should be in your “0” folder:

p

dimensions [1 -1 -2 0 0 0 0];
internalField uniform 3000000;
boundaryField
{
 inlet1
 {
 type totalPressure;
 p0 uniform 3000000;
 gamma 1.4;
 value $internalField;
 }
outlet
 {
 type fixedValue;
 value uniform 2067400;
 }
bottom
 {
 type symmetryPlane;
 }
 top
 {
 type zeroGradient;
 }
nozzle
 {
 type zeroGradient;
 }
defaultFaces
 {
 type empty;
 }
}

T

dimensions [0 0 0 1 0 0 0];
internalField uniform 298;
boundaryField
{
 inlet1
 {
 type fixedValue;
 value uniform 298;
 }
outlet
 {
 type zeroGradient;
 }
 top
 {
 type slip;
 }
 bottom
 {
 type symmetryPlane;
 }
nozzle
 {
 type zeroGradient;
 }
defaultFaces
{
 type empty;
 }
}

U

dimensions [0 1 -1 0 0 0 0];
internalField uniform (7 0 0);
boundaryField
{
 inlet1
 {
 type zeroGradient;
 }
outlet
 {
 type zeroGradient;
 }
bottom
 {
 type symmetryPlane;
 }
 top
 {
 type slip;
 }
nozzle
 {
 type slip;
 }
defaultFaces
 {
 type empty;
 }
}

Constants and Schemes

thermophysicalProperties

As the name states, in the thermophysicalProperties file we will specify the thermophysicalProperties of the gas.

thermoType
{
 type hePsiThermo;
 mixture pureMixture;
 transport const;
 thermo hConst;
 equationOfState perfectGas;
 specie specie;
 energy sensibleInternalEnergy;
}
mixture
{
 specie
 {
 nMoles 1;
 molWeight 29;
 }
 thermodynamics
 {
 Cp 1005;
 Hf 0;
 }
 transport
 {
 mu 0;
 Pr 1;
 }
}

In the file shown above, we have defined the simulation fluid is an ideal gas with Cp=1005, and a molar mass of approximately 29. This gives the approximate properties of air. We set the viscosity mu=0 because the simulation is inviscid.

fvSchemes

Usually, the fvSchemes from the rhoCentralFoam tutorials are good enough. However, I always change the flux limiter from vanLeer to vanAlbada. This is because, in my experience, it works much better at removing spurious oscillations from the flow-field. There are other options that I change occasionally, but for this wedge simulation that is the only parameter I have altered. The limiter option is found in the interpolationSchemes section:

interpolationSchemes
{
 default linear;
 reconstruct(rho) vanAlbada;
 reconstruct(U) vanAlbadaV;
 reconstruct(T) vanAlbada;
}

Running the Simulation/Results

Running rhoCentralFoam

Once the case is set up, it is time to run the simulation. To run the simulation you simply type rhoCentralFoam into the terminal from inside the case file.

Establishing convergence in time

In order to reach the final answer, the simulation must be run for a long time. Why? This is a limitation due to how we have chosen to approach this problem. Here we are using rhoCentralFoam which is an unsteady solver. The problem we are examining here is in fact a steady problem. So in order to achieve the answer we have to run from our arbitrary (but not insignificant) initial conditions, until the simulation stops changing in time. In steady-state solvers, the \frac{\partial}{\partial t} terms of the governing equations are manually fixed to zero. The remaining equation is then solved iteratively until convergence. Here, we solve forward in time until the \frac{\partial}{\partial t} terms approach zero eventually.

Recall that because of our boundary conditions, the inlet velocity is an output of the simulation. Here, I have chosen to use this as the variable for which I will test convergence in time. The following graph shows the inlet velocity as a function of time at a node along the inlet face.

timeConvergence
Fig: Simulation Time Convergence

As you can see from the picture, the inlet velocity fluctuates significantly at the beginning of the simulation. However by the end of the simulation time it has converged on a constant value.

NozzleShock
Fig: Shock-wave moving around as convergence achieved

This is a downside of the current approach. However, there are ways to speed up convergence such as the application of a relaxation factor in time, or the use of a non-reflective boundary condition. I will not go into detail regarding those here.

Results and Analysis

Below are contour plots showing the solution results.

Are the results what we expected?

At first glance it appears that the simulation is good. As designed, we have a stationary normal shock standing in the diverging section of the nozzle. To analyze this further I have taken the Mach number along the centerline of the diverging section and calculated the corresponding A/A*.

Analysis
Fig: A/A* and nozzle geometry (A/Athroat) in diverging section

We designed the nozzle to produce a shock at a specific location (A/A*) in the nozzle. From the figure above, we can see that the A/A* where the shock occurs is slightly higher than the design point. Similarly the  A/A* at the nozzle exit is higher than the design point. There are several reasons that this could be. For one, this simulation is two-dimensional and the flow properties are not constant in the y-direction. If we were to be rigorously comparing to the 1D calculation, we would probably need to either (a) average the Mach number at each length-wise station or (b) extend the nozzle to be significantly longer. Another possibility is (again since the simulation is 2D) that the Prandtl-Meyer waves emanating from the throat upset the accuracy of the simulation. These waves are visible in the figure above. These waves occur because of the finite change in angle between the throat and the straight diverging section. This could be avoid in future simulations by: (a) modelling the transition from throat to diverging section smoothly, or by (b) extending the length of the diverging section

What are some of the observed limitations? Improvements?

The biggest and most obvious limitation of this simulation (in my opinion) is the finite angle change that was used to create the geometry. This produced a Prandtl-Meyer expansion fan that hurt the quality of the simulation. If we were to improve this simulation, I would extend the diverging section to be much longer to reduce this finite angle change.

Another limitation of this simulation is that it is inviscid. While it is interesting and useful to compare to the analytical 1D solutions, in real life viscosity has a very real and important affect on the flow result.

Some improvements could be made to the simulation numerically as well. First of all, the grid used in the simulation was made coarse so that the simulation time would be small. A finer grid will produce a higher resolution image of the shock wave and an overall better result. Similarly, the size of the timestep could be decrease to lower the simulation CFL number. Some spurious oscillations can be seen in the post-shock section of the nozzle. These will be reduced by lowering the CFL.

Conclusions

In this post, we designed and simulated a converging diverging nozzle connected to a reservoir using the OpenFOAM solver rhoCentralFoam. The simulation produced the expected results. However, the finite angle change used to model the nozzle geometry produced expansion waves which hurt the quality of the simulation.

Please provide comments and suggestions on improvements!

curiosityFluids

Logo1

8 thoughts on “Pressure Driven Nozzle Flow with Shock – rhoCentralFoam”

  1. Hey guys!! nice website and axamples… i always enjoy with someone like you is supporting the openFoam community.

    during next months we are going to upload to our website the presentation of our work on noozle calculations with rhoCentral, that we presented in a conference. I will keep you posted!

    Like

  2. Thanks for this nice example!
    The pressure in your model is set to 30 bar, not 3 bar, right?
    Does it make sense to model such a nozzle without turbulence? / Is it possible to include turbulence in the rhoCentralFoam application?
    All the best,
    Ansgar

    Like

    1. You are right! It is 30 bar, not 3 bar. I need to fix the page.

      Definitely, if you were simulating such a nozzle for any real application, turbulence (or at least viscosity) should be included. In this case, I have simulated it without any viscosity (Euler equations). This was just to demonstrate the simple set up, and the calculations for the isentropic flow equations.

      If you simulated with turbulence you might need to alter the simulation set up. This was pointed out by a previous commenter.

      Like

  3. Hello,
    thank you for this fine example of a nozzle flow with the rhoCentralFoam. I think its very nice tutorial for beginners to start with.
    Regarding the results and analytics comparison, deviation in the shock position evaluation is expected mainly due to physics reasons and not numerics. As mentioned, compressible equations listed for this case are for an ideal shockless nozzle, contrary here you have a conic nozzle that doesn’t allow ideal expansion and thus incorporates the weak compression waves that are observed. Secondly as also mentioned, this is an inviscid calculation which doesn’t oppose the analytic equations, which are inviscid as well, but reality case where viscous boundary layer plays a major role for flow wall detachment and shock formation.
    Finally in numerical approach, it is not good to impose outlet conditions directly on the nozzle exit section as it is done here but to make an external buffer farfield zone for the outlet. And this is the major concern if the results mismatch while other physical conditions are satisfied. Actually it is very difficult to keep convective outflows at constant section while changing the first neighbor cells inside the nozzle. Thus, it should be changed, passing to full 3D and ameliorating the mesh resolution are secondary to this.

    Like

    1. Thanks for the input! I am always looking for feedback.

      The point of this case was to show a simple set up for this case in rhoCentral and also to provide an opportunity to show the calculations of such a case using the analytical solutions. Obviously in reality there is friction and other effects that are ignored by the ideal equations (unless you incorporate fanno flow etc).

      As for the outlet condition. You are very correct. In a non-ideal case there is certainly going to be non-uniform results at the outlet and setting a uniform outlet pressure would be incorrect. But in this case where achieving close to the ideal conditions is in fact the goal this should not be a problem.

      Please look at my other post (https://curiosityfluids.com/2016/0d7/14/converging-diverging-nozzle-v2-rhocentralfoam/), where I have done the NASA benchmarking case for a CD nozzle. Because the nozzle geometry is properly designed there is no problems with the expansion waves as seen here. As well enforcing the static pressure outlet boundary condition is used to achieve the analytical solution.

      Like

    2. Hello Vlad,

      in regards to your comment:
      ‘Finally in numerical approach, it is not good to impose outlet conditions directly on the nozzle exit section as it is done here but to make an external buffer farfield zone for the outlet’

      what would be realistic conditions to set for the outlet then?

      Like

    3. Hi Vlad,
      in regards to your comment ‘Finally in numerical approach, it is not good to impose outlet conditions directly on the nozzle exit section as it is done here but to make an external buffer farfield zone for the outlet.’

      what would be appropriate boundary conditions at the outlet then?

      thanks

      Like

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s