Mach 1.5 flow over 23 degree wedge – rhoCentralFoam

In this post I will simulate inviscid Mach 1.5 flow over a 23 degree wedge of a finite height. A finite height is required because at Mach 1.5 the leading edge shock-wave is detached for a deflection angle of 23 degrees. Without a finite height the shock wave standoff distance would be indeterminate.

Here is a summary of the contents of this post:

1. Calculation of Inlet Conditions
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:

Mach 1.5 Wedge Flow Tutorial File

Calculation of inlet conditions

We want to simulate flow at a Mach number of 1.5. However, in rhoCentralFoam, you cannot specify Mach number explicitly as a boundary condition. It must be implied by assigning the correct conditions in the p, T, and U files within the “zero” directory. The Mach number is fixed entirely by the temperature and velocity files.

To start this simulation I simply chose a arbitrary inlet temperature of 145.77 K (actually approximately equal to room temperature flow expanded to M=2). This gives a speed of sound (assuming γ=1.4 and R=287 J/kgK) of 242.01 m/s. Therefore, for Mach 1.5 flow given our inlet temperature, the inlet velocity will be 363.02 m/s.

Since the pressure in the inviscid case does not matter, I just set it to approximately 81.1 kPa. NOTE: In a viscous simulation, the pressure would be the parameter used to control the Reynolds number given a known viscosity.

Mesh Generation

For this simulation, the mesh is built using 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

The vertices are used to define the geometry of the problem. The following is a diagram of the geometric design of this simulation.

vertices
(
 (-0.0041 0 -0.005) //0
 ( 0 0 -0.005) //1
 ( 0.0041 0.001740346746 -0.005) //2
 (-0.0041 0.03 -0.005) //3
 ( 0 0.03 -0.005) //4
 (0.0041 0.03 -0.005) //5 
 (0.0082 0.001740346746 -0.005) //6
 (0.0082 0.03 -0.005) //7
 (-0.0041 0 0.005) //8
 ( 0 0 0.005) //9
 ( 0.0041 0.001740346746 0.005) //10
 (-0.0041 0.03 0.005) //11
 ( 0 0.03 0.005) //12
 ( 0.0041 0.03 0.005) //13
 (0.0082 0.001740346746 0.005) //14
 (0.0082 0.03 0.005) //15
);

Blocks

The blocks are used to define the different regions within the mesh. Here I have split up the grid into probably the most simple blocking pattern possible. Because of the simplicity of the grid design, the resulting grid would probably be unsuitable for more complicated viscous simulations without very high grid refinement.

blocks
(
 hex (0 1 4 3 8 9 12 11 ) (80 200 1) simpleGrading (1 8 1)
 hex (1 2 5 4 9 10 13 12) (80 200 1) simpleGrading (1 8 1)
 hex (2 6 7 5 10 14 15 13) (20 200 1) simpleGrading (5 8 1)
);

Boundaries

The boundaries section of the blockMeshDict file sets the faces and the boundary types for the simulation. Thus, it is very important to the simulation set up.

The 23 degree wedge obstacle is defined as type wall. The inlet, outlet and top boundary faces are set to type patch. To reduce the simulation time, the bottom boundary before the 23 degree wedge obstacle is set to type symmetryPlane.  The front and back faces are left undefined. By leaving them as undefined they get assigned to the group defaultFaces which we can later ensure is set to empty.

The boundary definitions from the blockMeshDict are shown below:

boundary
(
 inlet
 {
 type patch;
 faces
 (
 (0 8 11 3)
 );
 }
 outlet
 {
 type patch;
 faces
 (
 (6 7 15 14)
 );
 }
 bottom
 {
 type symmetryPlane;
 faces
 (
 (0 1 9 8)
 );
 }
 top
 {
 type patch;
 faces
 (
 (3 11 12 4)
 (4 12 13 5)
 (5 13 15 7)
 );
 }
 obstacle
 {
 type wall;
 faces
 (
 (1 2 10 9)
 (2 6 14 10)
 );
 }
);

Building the 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:

For the inlet, p, T, and U fully define the flow. By setting the outlet as zeroGradient, this enforces that the value on the face will be assumed to be the same as value of the first cell within the domain. For an outlet in supersonic flow this is usually an appropriate boundary condition. However, there are times when this will mess up the simulation results. I will probably write a post about this topic at some point. The obstacle enforces flow tangency, but not the no-slip condition (i.e. the flow is inviscid). Since the flow is symmetrical, a symmetry boundary condition is used for p, T and U on the bottom boundary.

For the initial conditions (the flow-field at time=0), I set the field to a pressure lower than the inlet pressure, the temperature equal to the inlet temperature, and set the velocity to 0. These conditions are essentially unimportant because they will not affect the end steady state solution. However, if these initial conditions are wild, it could cause the simulation to take longer to converge, or even crash.

Here is a summary of the p, T, and U files in the 0 directory:

p

dimensions [1 -1 -2 0 0 0 0];

internalField uniform 5000;

boundaryField
{
 inlet
 {
 type fixedValue;
 value uniform 81134.8794470031;
 }

outlet
 {
 type zeroGradient;
 }

bottom
 {
 type symmetryPlane;
 }

top
 {
 type zeroGradient;
 }

obstacle
 {
 type zeroGradient;
 }

defaultFaces
 {
 type empty;
 }
}

T

dimensions [0 0 0 1 0 0 0];

internalField uniform 145.768159204;

boundaryField
{
 inlet
 {
 type fixedValue;
 value uniform 145.768159204;
 }

outlet
 {
 type zeroGradient;
 }

bottom
 {
 type symmetryPlane;
 }

top
 {
 type zeroGradient;
 }

obstacle
 {
 type zeroGradient;
 }

defaultFaces
 {
 type empty;
 }
}

U

dimensions [0 1 -1 0 0 0 0];

internalField uniform (0 0 0);

boundaryField
{
 inlet
 {
 type fixedValue;
 value uniform (363.02 0 0);
 }

outlet
 {
 type zeroGradient;
 }

bottom
 {
 type symmetryPlane;
 }

top
 {
 type zeroGradient;
 }

obstacle
 {
 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

Typically, the fvSchemes from the rhoCentralFoam tutorials are good enough. However, there are a few parameters that I tend to (always) change. Specifically, 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;
}

controlDict

In the controlDict all of the simulation start and end details are set. The OpenFOAM manual should be consulted for any difficulty regarding these. For the present simulation I have clamped the CFL number at the low value of 0.1. By keeping a small CFL number it ensures stability of the results, and helps ensure smooth changes across the shock. To do this, adjustTimeStep needs to be set to yes, and maxCo must be set to the desired value.

adjustTimeStep yes;

maxCo 0.1;

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

I will probably have to repeat this in every post using rhoCentralFoam but, rhoCentralFoam is an unsteady solver. Therefore in order to simulate a steady problem like the one in this case, you must establish convergence in time. While a steady state solver sets the manually, an unsteady solver must be simulated until changes in the flowfield disappear.

For this case we will measure convergence in time using centerline post-shock pressure. To do this I picked a point approximately half way between the inlet and the front of the obstacle and plotted the pressure over time:

The first increase in pressure is the point being hit by an initial wave from the inlet. The second increase is when the shock moves across the point I selected. Here is a video of the convergence of the simulation:

Results and Analysis

Below are contour plots showing the steady state solution achieved by rhoCentralFoam. The simulation results can be qualitatively examined from the pictures above.

Are the results what we expected? 

It is clear that rhoCentralFoam is adequately capturing the desired physics in this simulation. We have a curved bow shock that is steady at a finite stand-off distance. The shock is smooth and there are no obvious spurious oscillations evident in the flow. This is probably in part due to the low CFL number (0.1) and the vanAlbada and vanAbadaV flux limiters employed.

At the centerline the bow-shock is a normal shock. Thus, the centerline pressure can be compared to the normal shock equations. From the normal shock equations, if M1 =1.5, p2/p1 across the shock should be 2.458. From the plot of center-line pressure distribution the simulation predicts p2/p1 of 2.4633 at the first node after the shock. This is a % error of approximately 0.22% indicating a good capture of the pressure jump across the shock.

What are some of the observed limitations?

 Some of the limitations of the current set-up can be seen from the images. For example, further away from the obstacle, where the mesh resolution is poor, the shock-wave appears significantly diffuse. This is especially apparent in the synthetic Schlieren image. If the objective of the study using this simulation was to observe shock-wave shape it would probably be of interest to increase the grid resolution in this area.

Improvements?

The quality of the simulation could be improved in several obvious ways. The first is to increase the grid resolution throughout. This will improve the overall quality of the simulation and specifically the shock capturing. Thermo-physical modelling could be employed to account for temperature dependent gas properties throughout the flow-field. Viscosity could be added to the simulation to account for the boundary layer on the obstacle wall. However, this is unlikely to affect the stand-off shock-wave.

Conclusions

In this post, we successfully simulated inviscid Mach 1.5 flow over a 23 degree wedge. The simulation produced the expected result of a curved shock-wave at a finite stand-off distance. We compared the pressure jump across the shock wave to the analytical results and saw good agreement. I suggested some improvements to the simulation that I may undertake in the future.

As always let me know if you have any suggestions/comments or find any errors etc.

curiosityFluids