Analysis of effect of Sub-Grid Scale on temperature in Prandtl Meyer shock flow

AIM:

The aim of this project is to simulate the Prandtl Meyer shock flow in Converge Studio.

OBJECTIVES OF THE PROJECT:

• Observe the mesh generation in ParaView.
• Effect of SGS temperature value on cell count and the shock location
• Create an animation of parameters varying through the geometry

THEORY:

What are shock Waves?

Imagine a situation where an object is moving with some speed like 5m/s while emitting some energy in one of its forms like sound, heat, etc. Now, this energy propagates with the speed of sound, and the air molecules far ahead of the object get adjusted according to the disturbance created by the object. Now if the same object were to travel with the speed of sound or greater than that, the air molecules have no way to know that the object is coming to them. When this happens the air molecules might stack up, Pressure builds up and the air loses its capacity to hold the water vapor and forms the mist near the object. The “near the object” mentioned in the previous statement depends on the shape of the object. To understand this, we look at the Prandtl Meyer expansion.

(Pic credits: CFD Basics with the application by John D Anderson Jr.)

In the above picture, a fluid is passing flowing from left to right at a supersonic speed. The geometry contains a sharp turn, and the flow expands around this sharp corner. This expansion wave contains an infinite number of weak Mach waves that fan out from the corner. In simple words, Mach waves are waves just like sound waves and along with this geometry, their intensity reduces due to expansion. Before the expansion, we can see that the streamlines are parallel to the geometry and are 1 dimensional after the expansion also, the streamlines are parallel to the geometry and are one-dimensional. But at the corner, the streamlines are bent and are two-dimensional. In short wherever there is a sudden change in geometry expansion happens and shock waves occur.

The shock wave problems can be solved using different methods. There are analytical methods and there are numerical methods like FDM which involve the Prandtl Meyer function. Analytical methods involve using the Prandtl Meyer function to find the parameters after the expansion. Numerical methods involve solving similar equations but by using the McCormack method. For this type of approach, Initial and Boundary conditions are required.

There are a variety of boundary conditions used in CFD that are differentiated based on their mathematical form. The two main types are

1. Dirichlet Boundary condition

In this type of Boundary conditions, the values such as the flow parameters are directly described. E.g., u = 300, T= 300 k.

2. Neuman Boundary condition

In this type of Boundary conditions, the values such as the flow parameters or domain are described as gradient or normal derivative.

The other types are a mix of these two boundary conditions.

Apart from these, there are other types of Boundary conditions that most of the commercial software use such as

Wall

This type of Boundary condition is used when there is no relative motion between a particular surface and the fluid flowing. There are Fixed Wall – when the surface is stationary, moving wall – When the surface is moving and No-slip – This is an ideal situation where the molecules at the wall are stationery and slip – This is when very few molecules at the wall are stationary

Inflow/Outflow

This type of boundary condition is used to denote the starting and ending point of the flow. The solver will understand that the points beyond these BCs need not be solved.

Symmetry

This type of Boundary condition applied to a particular surface means that there is no flow across that surface or the gradient across that surface is zero.

PROCEDURE:

SETTING UP THE CASE

Cleaning the geometry

We start with creating geometry in Solidworks. In this case, the geometry was already created.

• As we open a new project, we import the geometry by dragging it from the file explorer and dropping it on the converge screen or by importing it from the file menu.
• Once the geometry appears, we need to check the geometry units. We can do this by going to the Left-side geometry panel > Options > check the boundary box.
• If the boundary size is in the range of 100meters we can scale it by geometry panel > Transform > Scale > give a scaling factor of 0.001. And the geometry will come to a range of 0.1 m
• In the “Boundary” options select flag and press “+” select 4 boundaries and press ok. We keep the geometry in standard orientation.
• 4 boundaries will be shown, we rename them as

Boundary 1 > “Inlet” (Select the 2 triangles on the left.)

Boundary 2 > “Outlet” (select the 2 triangles on the right)

Boundary 3 > “Top and Bottom” (select the top 2 and the bottom 4 triangles)

Boundary 4 > “Front2d” (select the 2 triangles on the front)

Boundary 5 > “Back2d” (select the 2 triangles on the rear)

• Now left-click on a boundary name and select the triangles on the geometry to assign those triangles to the boundary. Repeat this process for all the 5 boundaries created.
• We can make sure that all the boundaries are assigned by checking the Not Assigned (Boundary ID 0), If the triangles in this boundary are none then all the triangles are assigned to the created boundaries.

Since the geometry in converge is made of triangles, the above steps will assign triangles to boundaries.

• On the left side of the geometry tab, there are some icons. Click on the icon with a blue triangle with a yellow normal pointing outward. This will show the normal of each boundary on the geometry created. If these normal are facing outward they must point inward.
• To do so on the tool on the left there is another tab called Transform > Normal.
• In this tab make sure the propagate change from a single triangle is checked and click on any triangle on the geometry and click Apply. This will change the normal direction to inward.

The normal in geometry should point towards the fluid domain. The above steps will invert the normal towards the inside where the fluid is flowing.

• On the right side of the GUI, there is another toolbox. At the bottom of this toolbox, there are three tabs. Go to the diagnosis tab and with the default options click on find. This will generate a report if all the boxes shown above are checked green select clear.
• Go to the case setup tab and click Begin Case Setup.
• A new window will pop up (If it does not, click on application type) select Time based, and click on apply, and then done.
• Now on the right side, we can see many sections which is where we can setup the simulation.
• Click on the Materials tab and select the predefined mixture as Air, Uncheck the reaction mechanism and check the species option, and press ok.
• A new message will pop up asking Do you really want to overwrite gas.dat, mech.dat, and therm.dat by loading predefined mixture as air? Click on yes.

We selected the type of fluid flowing inside the geometry to air in the above steps

• In the Gas simulation option with the default settings click apply and done
• Same for the global transport parameters with the default settings click on apply and done.
• On the species mechanism click on add and type O2 and N2 and add them from the suggestions/manual entry.
• In the run parameters in the solver: select the steady-state solver to make sure the simulation mode is set to full hydrodynamic Gas flow solver: Compressible.
• In the miscellaneous tab uncheck the use shared memory option and check the steady-state monitor option. And click on ok.

The type of solver will be set to transient in the transient condition and a steady-state monitor is not necessary for the transient case.

This selects the type of solver to be used for our case.

• In the Simulation, time parameters enter the following values.

Steady-state

Start time: 0

End time: 25000 cycles

Initial time step: 1e-7

Minimum time step: 1e-7

Maximum time step: 1

Maximum convective CFL: 0.8 click on ok.

• In the solver parameters section select the Density-based option. And in the equations tab on the left side click pressure and on the right side select the preconditioner to none. In the steady-state solver, control tab selects the Mach number to 0.8 and click OK.
• Now we go to Steady-state Monitor and Add config > Minimum number of cycles to be executed for steady-state solver: 5000, Monitor Units: cycles Monitor start delay, cycles: 200, Sample size: 100, In the left window double click the Monitor variables and select AVG_TEMP, Target type double click and select Boundary, Target as Outlet.
• Add two more variables and select AVG_VELCOITY, TOTAL_PRESSURE, MASSFLOWRATE select the target type Boundary, Target as Outlet. Press OK.
• Now before we go into the boundaries option click on regions and initialization. Click on add new region. Right-click on region 0 and rename it to the volumetric region (user-defined). On the left side near the species click on +Air gives the velocity in x-direction as 680m/s and pressure as 50000 pa and Temperature 300K click ok.
• Now in the boundaries tab assign the following to each boundary.

Inlet > inflow

Outlet > Outflow

Top and Bottom > wall

Front2d > Two_D

Back2d > Two_D

• In the inlet boundary select the pressure to “specified value” 101325 pa. Select the velocity to “specified value” and input 680m/s. Select the temperature to a specified value of 286.1 K and click on +Air in the species boundary condition.
• In the Outlet, boundary selects the pressure to “zero normal gradients” and select the velocity to “zero normal gradients” select the temperature backflow to a specified value of 300K and click on +Air in the Species Backflow.
• In the Top and Bottom boundary select the Velocity boundary condition to the option next to UDF as slip. Select the temperature condition “Zero normal gradients”. The Turbulent Kinetic energy (tke) boundary condition as Zero normal gradient. Also, the Turbulent dissipation (eps) boundary condition to wall model.
• Both the Front2d and Back2d are set to TWO_D
• Now set the regions in the left panel next to boundaries as Volumetric region User-defined.

The above steps are to be followed to set the boundary conditions.

• In the physical models > Turbulence modeling select the RANS option and select RNG-kƐ and press ok.

We select the type of flow (laminar/turbulent in the above steps)

These steps will tell converge to observe the given values, changes in them,

and note them for post-processing.

• In the Grid control select the Adaptive mesh refining option and press ok. A new option appears in grid control “Adaptive Mesh Refining”
• In the Grid Control > Base Grid > set the x, y, z sizes to 0.8m and click ok.
• In the Adaptive mesh refining option. Select the minimum number of cells to 1 and the maximum number of cells to 2,00,000. And add the volumetric region to the active region. In the bottom tab select the Sub-grid Scale and select the max embedding level: 2; Sub-grid criterion: 0.05.
• Select the timing control type: Sequential and select the start 5000, select the end time: 999999999 and click ok. (These values can be played with for a little experimentation)
• In the Output/Post-processing section,> post variable selection > In the species tab check the mass fraction tab and add N2, O2, and press ok.
• In the Output/Post-processing section > Output files > writing time intervals tab > set the time interval for writing 3d output files to 100cycles this will make sure that the files are written for every 100 cycles thereby we have more data and more clear curves.
• The time interval for writing the restarting output file. The restart file can be used to restart the simulation from a particular moment without having to start the simulation from the starting in case of any unexpected events like power loss, Out of memory, etc. But writing it every time step takes time and energy so for our convenience, we write it every 1000s. Press ok.

Save the case setup, upon pressing ctrl+S we save the case setup file which has a .cvg extension.

EXPORTING THE FILES

Now since the case is set, we need to run the simulation. Save this case file on the local computer and then go to the “taskbar on the top left” > file > export > export input files. We should see a total of 12files. The usual location is set to the same directory where the case file is saved. Set the location and export the files.

RUNNING THE SIMULATION

Open the Cygwin terminal as administrator and navigate to the location where the case files are stored using the cd command. Then we execute the following command.

mpiexec.exe -n 4 converge.exe restricted </dev/null> logfile.txt &

mpiexec.exe is the message passing interface executive file.

-n 4 represents the number of cores to be used for running the simulation.

converge.exe is the converge file that runs the simulation.

restricted limits the total number of cells to 5,00,000

</dev/null> logfile.txt writes the output to a logfile.txt

& allows to use terminal even with the process running

Now we can see many lines being printed on the Cygwin terminal. We can stop this process by pressing Ctrl + C in the Cygwin terminal. We go to the “inputs” file in the case directory, In the outputs control > log level > 0. Now run the above command again in the terminal. This should reduce the amount of output being printed in the Cygwin terminal.

This is done to know the control of the level of output being printed on the Cygwin terminal. Once the simulation is finished, we need to post-process the files.

POST-PROCESSING

We open converge > on the bottom left there are some options like postprocessing, lines, etc. We select the post-processing and that will take us to a new window. Here we select the location of the output file which is created in the same directory in which we stored the input files. Once the output file is selected a list appears in the window below the file location option > click on select all and the below the cell variables will be listed > click on select all. Above the file location option, there is space for naming the post-processed file we give a proper name (preferably without spaces/special characters). And then select the output format as “Para-view vtk inline binary”. Once we do this the click convert. And the files will be post-processed.

Once the post-processing is done open Para-view > open file select the output file location > select the converted file _..vtm – group file. This will load the results in Para-view. Click on apply to view the case we setup.

In the top panel select the surface with edges to view the mesh and we might see that with the mesh there are some diagonal lines visible on the mesh. This is because the triangle converge created. So, we create a slice of this geometry. This slice will have the normal in the Z direction. For this in the top panel in the bottom-most line of the top panel there is an option called slice we select that and then click on Z-normal and apply. This will create a slice of the geometry. In the top toolbar, there is an option called surface with edges select it. This will show the mesh

The mesh on the slice would look like this

V = 100m/s; SGS = 0.05 @ last timestep V = 680m/s; SGS = 0.05 @ last timestep

V = 680m/s; SGS = 0.01 @ last timestep V = 680m/s; SGS = 0.05 @ last timestep V = 680m/s; SGS = 0.1 @ last timestep

Now for the slice to the left of the surface option, there is another option which might be something like vtkblockcolors or idreg or solid-color. Here, select the required parameter and this will show the parameter's contour at each timestep. Here we select temperature and these are the contours.

Temperature plots

V = 100m/s; SGS = 0.05 @ last timestep V = 680m/s; SGS = 0.05 @ last timestep

V = 680m/s; SGS = 0.01 @ last timestep V = 680m/s; SGS = 0.05 @ last timestep V = 680m/s; SGS = 0.1 @ last timestep

Upon clicking the play button on the top the contours of the parameter at all the timesteps will be displayed one after the other at 1sec interval. The animation of all the time steps can be saved using the save animation option.

Temperature animation@

V = 100m/s; SGS = 0.05

V = 680m/s; SGS = 0.01

V = 680m/s; SGS = 0.05

V = 680m/s; SGS = 0.1

We can also use some of the converge features to compare results. In converge on the bottom left beside the post-processing, there is another tab called line plotting. Upon clicking this it will take us into another window. Here we select the same location where we saved the case setup file. This is because the converge takes the .out files as the input upon selecting the file location some variables are displayed in the window below. We select the file as cellcountranks.out. The x-axis will be cycles and on the y axis, since we are interested to see the variation of cell count with the subgrid cells we select cell count.

These are the graphs we generated using the converge line plotting module.

Total cell count

CONCLUSION:

From the above results, it is observed that

• With a decrease in SGS value the cell count increases. This can also be observed from the animation that more cells are used to capture the shock wave.
• The span of the wave increased with a decrease in SGS value
• At a speed, less than the sonic speed, the change in temperature is very little.