Turbulence Modeling in STARCCM+

Aim:

In this project, we aim to study the turbulence modelling by simulating a flow over a backward step with varying Reynolds numbers in STARCCM+. The concept of Y+ was also studied.

Calculations:

The simulation was done for air at a temperature of 25 degree Celcius.

Reynolds Number (Re) is given as:

Re = (ρ * v * D) / μ

The equation rearranges to give velocity as:

v = (Re * μ) / (ρ * D)

Case 1:

Re = 100

ρ = 1.184 kg/m^3                                      (Density of air at 25 C)

μ = 1.849*10^-5 kg-ms                            (Dynamic viscosity of air at 25 C)

D = (4*Area)/(Perimeter) = 0.046 m  

This gives, v = 0.03 m/s

Case 2:

Re = 1e5

v = 312.3 m/s

Case 3:

Re = 1e7

v = 5205 m/s

for Re = 1e6,

     v = 520.5 m/s

The flows were simulated for Reynolds number 100, 1e5 and 1e6 instead of 1e7. The reasons and related theory will be described as we progress with the solution.

Procedure:

The initial steps from creating the geometry to setting up a volume mesh is same for all the cases and hence wil be described only once at the beginning.

Step -1: Creating the Geometry:

• After opening a new simulation file in STARCCM+, right click on 3D CAD Models and select New. This will open the modeling window. 
• Right click on XY plane and select Sketch.
• In the sketch plane that will appear, draw the 2D sketch of the geometry as provided in the question.

                                    2D Sketch of the backward step geometry

• After creating the sketch, click OK. By default, the sketch will be saved as Sketch 1.
• Right click on this Sketch 1 and select Extrude. In the extrude tab, select the extrusion length as 0.1m and click OK.

                                               Completed Geometry

• Once the geometry is created, select the front face and rename it as 'Inlet'. Select the rear face and name it as 'Outlet' (refer figure below). The other faces will be defined as Default which can be later changed to Walls for ease of understanding.

• Once the geometry is created, right click on the geometry part and select New Geometry Part. This will assign the geometry to the parts section where the next step will be to check for surface defects. Here, we can rename the surfaces (Default to Walls) if we want.
• To check for surface defects, rigth click on the geometry part under Parts and select Repair Surface. In the window that will appear, make sure that all the 6 default error types are selected and hit Execute. This will show the number of errors for the respective error types.

                                                Surface Repair Window

• Two errors each of Face Quality and Face Proximity were found. Since these two error types are taken care of during the automatic surface repair step, we do not need to use the surface wrapper.

STEP -2: Meshing

• To create the mesh, right click on Operations --> New --> Mesh --> Automated Mesh.
• In the dialog box, select the following models:

       Surface Remesher

       Automatic Surface Repair

       Polyhedral Mesher (mostly used for internal flow simulations)

       Prism Layer Mesher (to capture effects at the boundary)

                                                 Volume Mesh

• To create a plane section, right click on Derived Parts --> New Part --> Section --> Plane. Make sure that the section is perpendicular to the Z direction.
• The detailed mesh and prism layers can be seen on this section. Go to the Mesh scene and hide the other surfaces except the derived plane.

                                             2D View of the Mesh

STEP -3: Creating Regions and giving BCs:

• Right click on the geometry part under Parts and select Assign Parts to Regions. In the dialog box that will appear, make sure to select 'Create a boundary for each part surface'.
• This will create a separate region for each of the surface we defined earlier while creating the geometry. 
• To give the boundary conditions, magnify the Regions tab, go to Boundaries. There, we can see the three regions - Inlet, Outlet and Walls. 
• Select the Inlet region and set the Boundary Condition as 'Velocity Inlet'. Similarly, set 'Pressure Outlet' and ' Wall' for Outlet and Walls respectively.

--------------------------------------------------------------------------------------------------

Significance of Prism Layers:

Prism Layers are prism-shaped cells that are generated at the boundaries to capture the Boundary Layer phenomenon. A boundary layer is a layer of fluid generated in the immediate vicinity of a bounding surface where the effects of viscosity are significant.

The boundary layer is a complex phenomenon which is very thin and the solution gradients are very high. To capture these steep gradients, we could use a fine mesh (of any type of cell) throughout the domain but it would be a costly and time-consuming affair.

The solution to this is to use a Prism Layer near the boundaries and the core mesher in the central/free flow region. Prism layers allow high aspect ratio cells which can compute the steep gradients. 

The effect of a prism layer is to minimize the meshing time and simulation costs incurred by a fine mesh by providing thin cells only at the boundaries to capture the boundary layer gradients. Thus, prism layers facilitate accurate boundary layer resolution without increasing the cost and time of simulation.

Prism layers can be seen in the 2D view of mesh.

--------------------------------------------------------------------------------------------------------

The next step will be to set up the Physics conditions. To do so, go to Continua, right click on Physics 1 and select 'Select Models'.

The steps till here will be the same for all cases.

Case 1:

A] For the Laminar flow case, the following models were selected:

Three Dimensional, Steady, Gas, Segregated Flow, Constant Density, Laminar. Constant density was used as the density variation of the fluid was assumed negligible as the velocity was very low and temperature effects will be negligible on the density variation.

The stopping criteria was set to 250 iterations and the simulation was run.

                                                Pressure Contour

                                                Velocity Contour

                                                    Residuals Plot

As can be seen from the above plots, the solution converges after about 60 iterations and the developed flow profile can be observed in the contour plots.

B] The same simulation for Re = 100 and v = 0.03 m/s was run using the turbulent model.

Instead of the Laminar model, Turbulence and k-epsilon models were selected along with the auto-selected models.

                                                     Pressure Contour

                                                        Velocity Contour

                                                      Residuals Plot

Observations:

1. Contour Plots: As can be seen, the Pressure and Velocity contours for both the laminar and turbulent simulations are nearly identical with very minute distinctions. This can be attributed to the fact that a turbulent model won't affect an originally laminar flow as the turbulent viscosity computed by the turbulence model will be much lower as compared to the molecular viscosity (since turbulent viscosity is directly proportional to square of KE and KE for a laminar flow is very low in magnitude).

2. Residual Plots: Residuals for both the laminar and turbulent flows converge in the same range of iterations (60-70). However, the residuals for the turbulent model become unstable and vary sharply in a small value range. One reason to explain this unstable nature could be the round-off errors in calculating k and epsilon. Some amount of energy from the mean flow can be taken off by the KE and dissipation to cause this unstability. However, this small unstability doesn't cause any errors in the velocity and contour plots.

--------------------------------------------------------------------------------------------------------

Before proceeding with the further cases, let us take a look at a concept called y+.

Need for y+:

The above image shows the velocity variation from the wall boundary to the free flow. 

• As we can observe, there are steep gradients near the wall.
• In a second order FVM solver, variation across the cells is linear. hence, to capture this velocity profile accurately, we will need to place a large number of thin cells near the wall (middle figure) to resolve these gradients.
• Thin cells have high aspect ratio which increases the chances of cell skewness resulting in poor cell quality and ultimately divergence of the solution.  
• Instead of using large number of thin cells, we can use a single large cell and have non-linear variation (figure b) to capture the steep gradients.
• Advantages of such a cell would be less number of cells, bigger cells with no possibility of cell skewness and increased stability of the solution.

What is this non-linear function?

• Experimental measurements and DNS have shown the variation close to the wall, which can be used as a reference model. The function is as below.

               Non-linear function from a flow parallel to plane parallel channels

This is a function of dimensionless velocity U+, tangential velocity U(y) and distance normal to the wall y+.

Wall Functions:

Wall functions are empirical functions that are fitted to the observed behaviour close to the wall.

The standard wall functions for the viscous sublayer and the log-law region are clearly defined:

The variation is linear for y+ < 5 and logarithmic for y+ > 30.

[Note that the X-axis of the plot is a logarithmic scale. Hence the curvature in the viscous sublayer is actually linear and the linear variation in th log-law region is actually logarithmic]

Types of Wall Treatments:

1. Low y+ (y+ < 5):

The low y+ wall treatment resolves the viscous sublayer and needs little or no modeling to predict the flow across the wall boundary. The transport equations are solved all the way to the wall cell. The wall shear stress is computed as in laminar flows. To resolve the viscous sublayer, these models require a sufficiently fine mesh with near-wall cells located at y+ of around unity. The computational expense that is associated with this approach can be significant, particularly for large Reynolds number flows where the viscous sublayer can be very thin. Therefore this wall treatment is suitable only for low Reynolds number flows.

2. High y+ (y+ > 30):

The high-y+ wall treatment does not resolve the viscous sublayer. Instead wall functions are used to obtain the boundary conditions for the continuum equations. Wall shear stress, turbulent production, and turbulent dissipation are derived from equilibrium turbulent boundary layer theory.

This approach assumes that the near-wall cell lies within the log-law layer of the boundary layer at  y+> 30. The main advantage of the high-y+ wall treatment is therefore the significant savings in the number of near-wall cells.

3. All y+ (5 < y+ < 30):

The all-y+ wall treatment is a hybrid treatment that emulates the low-y+ wall treatment for fine meshes, and the high-y+ wall treatment for coarse meshes.

It is also formulated with the desirable characteristic of producing reasonable answers for meshes of intermediate resolution, that is, when the wall-cell centroid falls within the buffer region of the boundary layer. A blending function is then used to calculate turbulence quantities such as dissipation, production, and stress tensor.

 

In this project, we have used all y+ for all cases except the last one. The reason and result for that will be explained in that section.

--------------------------------------------------------------------------------------------------------

Case 2:

For Re = 1e5, v = 312.3 m/s. These simulations were run for 1000 iterations to attain the fully developed flow.

Here, the case was run using both the turbulence models - the realizable k-epsilon and the SST k-omega.

Results for realizable k-epsilon model-

                                                      Pressure Contour

                                                         Velocity Contour

                                                      Residuals Plot

Results for SST k-omega model-

                                                     Pressure Contour

                                                         Velocity Contour

                                                            Residuals Plot

Observations:

1. Contour Plots: The pressure contour plots are again almost identical for both the turbulence models. The velocity plots are also similar, however the k-omega model has better captured the shear layer and the recirculation region.

2. Residuals Plot: The residuals of both the models are observed to be stable. The residuals for k-epsilon models have stabilized close to 0, meaning the error is very low. For the k-omega model, however, the fluctuations are a bit too high. One reason for this could be the high sensitivity of the k-omega model to the free stream values. This can be resolved by using a finer mesh. Also, the residuals do not drop to zero. A reason for this could be the 'backstep disturbance' which means that the residuals do not drop until the disturbances have been transported out of the domain.

Case 3:

For the case Re = 1e7 v = 5205 m/s, the solution blew up and did not converge. One reason could be the very high velocity magnitude.

                                           The solution did not converge

So, the simulations were run for Re = 1e6 and v=520.5 m/s.

Results for realizable k-epsilon model-

                                                    Pressure Contour

                                                         Velocity Contour

                                                              Residuals Plot

                                             Surface Average Report for Pressure

                                            Surface Average Report for Velocity

Results for SST k-omega model (all y+ wall treatment)-

                                                       Pressure Contour

                                                       Velocity Contour

                                                            Residuals Plot

                                              Surface Average Report for Pressure

                                            Surface Average Report for Velocity

Observations:

1. Contour Plots:

The pressure and velocity contours for the k-epsilon and k-omega models are very different. The pressure plot for k-epsilon provides very high values of pressure. The pressure decreases after the step when compared to the inlet region. For the k-omega model, the pressure values are low and even negative in some parts of the domain. These negative pressure regions will actually restrict the flow in the positive direction. Thus the k-epsilon plot seems to be giving a close to actual consdition plot.

For the velocity contours, the SST k-omega all y+ wall treatment gives very poor resolution at the shear layer and the recirculation zone. However, the k-omega model resolves the flow in the free stream more accurately than the k-epsilon model. Thus, neither model gives the accurate results. As a compensation between accurate resolution of the near wall as well as the free stream, a simulation for SST k-omega (low y+ wall treatment) was done which gave fairly accurate results for both conditions. It will be discussed subsequently.

2. Residual Plots:

The residuals for the k-epsilon model rose almost linearly wheareas for the k-omega model, the residuals were almost constant with small spikes indicating sudden variation of variables. For the k-epsilon model, the residuals increase with progress of time as the values of variables keep on increasing leading to higher errors. For the k-omega, the fluctuations can be attributed to the constant switch between the low y+ and high y+ wall treatment. This can be concluded based on the fact that the spikes occur majorly in the Tke (Total KE) and Sdr (Specific Dissipation rate). This can also be verified with the low y+ case.

3. Surface Average Plots:

Surface average plots give us the plot of average value of a parameter over the entire selected face with iterations (or time). For both the models, the average velocities at the outlet and along the cut plane section stabilize after about 250 iterations. The pressure variation for both the models is different. This can be due to the different approaches in calculating the turbulent viscosity values.

Results for SST k-omega model (low y+ wall treatment)-

                                                     Pressure Contour

                                                  Velocity Contour

                                                   Residuals Plot

                                       Surface Average Report for Pressure

                                         Surface Average Report for Velocity

Observations:

1. Contour Plots:

The contour plots with low y+ wall treatment give the most accurate results for both pressure and temperature. The velocity contour perfectly captures the shear layer and recirculation zone as opposed to the all y+ wall treatment k-omega model.

2. Residuals Plot:

The residuals are quite similar to the all y+ wall treatment k-omega model. The only difference being the smooth curves for Tke and Sdr. As stated in the previous case, these are smoothened out due to the use of a single wall treatment (low y+) instead of fluctuating between low y+ and high y+.

3. Surface Average Plots:

Both the pressure and velocity plots have given  values which are between those of the k-epsilon and all y+ wall treatment k-omega models. This can be due to the fact that the k-epsilon model is quite unreliable when it comes to calculating the damping functions and probably under-estimates the damping. The all y+ SST k-omega on the other hand switches between the k-epsilon model at the boundary and the k-omega at the free stream. The fluctuations going from low y+ to high y+. In the low y+ condition, this constant transition is avoided giving the dissipation and KE values in a specific range.

CONCLUSION:

• The turbulent model works fine when applied for a laminar case (Case 1) with only the wavy nature of the residuals plot. The reason that the the turbulent model works well even for a laminar case is that the turbulent model solves the same basic flow equations as that of the laminar flow model and solves two more equations to find k and epsilon. Due to the laminar flow characteristics, the values of turbulent KE (k) and dissipation rate (epsilon) are very low and do not have a significant effect on the solution.
• When the flow is subsonic (Case 2), both the k-epsilon and k-omega turbulence models gave very close results. Thus, either of the two models can be used to simulate subsonic flows.
• For supersonic flows (Case 3), the SST k-omega model (low y+ wall treatment) gave the most accurate results. The SST k-omega model (all y+ wall treatment) gives accurate results for the freestream values but lacks accuracy at the boundary layers. On the contrary, the k-epsilon model resolves the boundary layer accurately but performs poorly for the freestream resolution. Thus, the SST k-omega low y+ wall treatment was the most accurate of the three and provides good balance between the resolution of boundary layer and freestream.

Google Drive link for simulation files:

https://drive.google.com/drive/folders/1RJn2Qe_ACouqKPqWG7tKkpR7ho3dWKr-?usp=sharing