Skip to content # 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

Simulation summary: ## 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
(
(-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.

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.

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.

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:

### 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
{
}
nozzle
{
}
defaultFaces
{
type empty;
}
}

T

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

U

dimensions [0 1 -1 0 0 0 0];
internalField uniform (7 0 0);
boundaryField
{
inlet1
{
}
outlet
{
}
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;
}

## 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.

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.

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*. 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 Categories

Tags ## 20 thoughts on “Pressure Driven Nozzle Flow with Shock – rhoCentralFoam” Leave a comment ›

1. Shehab Salem says:

thanks

• Shehab Salem says:

could you plzz give me access to download the files, i am waiting long time, but still you did not reply my email

2. cey says:

Hello,
I tried this case but with rhoCentralFoam I always always always get negative T0 error. How did you avoid it? I am using OpenFoam5x.

3. hussein says:

Hi..thanks for the nice sharing!
In the 0 folder of your case file, there is a quantity Ma. I would appreciate it if you could explain how to define this quantity, particularly in the nozzle boundary.

• curiosityFluidsAdmin1 says:

The Ma parameter is actually a post-processing parameter. You do not need it in the 0 folder for this case

• hussein says:

4. Yvette says:

Hi There
I have downloaded the files and ran blockMesh but got a few errors. Appreciate if there should be a correction.

Best regards
Yvette

Create time

Not deleting polyMesh directory
“/mnt/c/Users/brigh/tutorial/NozzleShock/constant/polyMesh”
because it contains blockMeshDict
Creating block mesh from
“/mnt/c/Users/brigh/tutorial/NozzleShock/constant/polyMesh/blockMeshDict”
Creating block edges
No non-planar block faces defined
Creating topology blocks
Creating topology patches

Creating block mesh topology
–> FOAM Warning :
From function Foam::polyMesh::polyMesh(const Foam::IOobject&, Foam::pointField&&, const cellShapeList&, const faceListList&, const wordList&, const Foam::PtrList&, const Foam::word&, const Foam::word&, bool)
in file meshes/polyMesh/polyMeshFromShapeMesh.C at line 873
Found 16 undefined faces in mesh; adding to default patch.

Check topology
Curved edge 21 23
does not correspond to a block edge.
Curved edge 27 29
does not correspond to a block edge.

Basic statistics
Number of internal faces : 9
Number of boundary faces : 30
Number of defined boundary faces : 30
Number of undefined boundary faces : 0
Checking patch -> block consistency

–> FOAM FATAL ERROR:
Block mesh topology incorrect, stopping mesh generation!

From function void Foam::blockMesh::check(const Foam::polyMesh&, const Foam::dictionary&) const
in file blockMesh/blockMeshCheck.C at line 228.

FOAM exiting

• curiosityFluidsAdmin1 says:

Which version of OF are you using? Did you change anything in the blockMeshDict I uploaded?

• curiosityFluidsAdmin1 says:

I have updated the location and content of the blockMeshDict and it should now work

• delta says:

Hi，Yvette!
I attempt to download the file but I failed since it told me that permission denied.Did you encounter such a situation？
Best wishes!

5. mar hernan says:

Thank you so much for this exercise so explicity. Could you share one of converging-diverging nozzle?

6. masoud says:

hi
i simulate external flow on body by rhocentralfoam,but in paraview i see temperature decrease low than 273 k
what’s is your opinion?
thanks

7. Gerardo says:

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!

8. Ansgar says:

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

• curiosityFluidsAdmin1 says:

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.

9. Ansgar says:

Thank you so much for this nice example!
Am I mistaken, or is the pressure at the inlet set to 30 bar, not 3 bar?
All the best,
Ansgar

10. Vlad says:

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.

• curiosityFluidsAdmin1 says:

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.

• kaize says:

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?

• kaize says: