tushargh29/transonic_shock_buffet

This repository is dedicated to the studienarbeit project of "Numerical investigation of 2D transonic shock-buffet around a NACA 0012-34 airfoil using OpenFOAM and flowTorch" undertaken by Mr. Tushar Gholap under the guidance of Dr. Andre Weiner.

transonic_shock_buffet's Issues

Literature review

Goal:
Create a literature review concerned with the numerical simulation of buffeting at NACA0012(-34) and similar airfoils. Information of interest might be:

• investigated range in terms of Reynolds and Mach number as well as the angle of attack; conditions for buffeting
• type of turbulence modeling, e.g, RAS, LES, hybrid
• turbulence model
• application of wall modeling
• mesh type (block-structured, unstructured) and mesh resolution
• 2D or 3D simulation

There are some sources in the shared cloud folder to start with. A good strategy might be to investigate the cited articles and to look at newer articles citing the article in question.

Mesh improvements

Hi Tushar,

your next task will be to make the mesh of the base case more suitable for buffet simulations. I suggest the following changes for a start:

• create a workflow without the execution of restore0Dir, the second decomposePar, and rhoPimpleFoam (only meshing)
• change the workflow to use 4 instead of 2 processes
• in the snappyHexMeshDict (reference)
  • increase the surface refinement to reduce the y-plus value
  • add more surface layers to reach a y-plus value close to one (you computed a theoretical estimate)
  • add distance-based refinement to increase the region refinement close to the airfoil
  • add a refinement box to increase the wake resolution
  • add a refinement box to increase the resolution around the shock

Suggestion: initially, be moderate the changes to see how meshing time and memory consumption are affected be the different meshing parameters.

Best, Andre

Perform test simulation on Phoenix

Hi Tushar,
please, try to follow these instructions to perform a first test simulation on Phoenix. Let me know if you need more detailed instructions. You can also get in touch with Darshan since he is working currently on a similar step.
Best, Andre

Obtaining buffet onset angle of attack for constant Ma and Re numbers for NACA 0012 airfoil

Hello sir,

As per our discussion in the last meeting, I created a new simulation case setup as follows:
Reduced Y Grading of 5000, using wall functions and 450 cells chordwise.

Then, I tried to obtained the onset buffet angle of attack for constant Re and Ma numbers specifically for NACA 0012 airfoil using the bisection method.

Thanks.

Best,
Tushar

Comparison with experimental data for NACA0012 profile

Hi Tushar,

I've created a pull request with a new setup. Changes are:

• improved blockMeshDict; have a look at the new mesh
• standard NACA 0012 airfoil
• files and options for Spalart-Allmaras turbulence model
• two function objects to compute and write the surface pressure coefficient

I suggest taking a step back. The article

Reynolds-Averaged Navier–Stokes Study of the Shock-Buffet Instability Mechanism

contains comparisons of the surface pressure coefficient for pre-buffet and buffet conditions (standard NACA 0012 airfoil). To see how the new setup performs, I suggest the following steps:

• create the setups 1, 2, 4, 6 from table 1 in the article; consider the nominal angle of attack and Mach number (subscript N); keep the Reynolds number at 10^7 (current base setup)
• from this article, digitize the corresponding experimental data (1: page 23 (26 in the PDF); 2: page 29 (32 in the PDF); 4: page 43 (46 in the PDF); 6: page 55 (58 in the PDF)); tables lower and upper airfoil surface; columns X/C and C_P
• run the simulations until steady-state; you can take the lift or drag coefficient as a reference
• compare the pressure coefficient between experiment and simulation
  • C_P in the simulation is written to postProcessing/surfaces/
  • load both the experimental and numerical data in a jupyter notebook and plot them comparatively
  • typically, -C_P (negative) is plotted over the length coordinate instead of C_P
  • the x-coordinate should be normalized with the chord length such that x in [0, 1]
  • the plots should look similar to the ones in figure 1 of the article mentioned above

Best, Andre

Change fvSchemes

Hi Tushar,

I may have found the issue with the setup. Please, try the following change in fvSchemes:

// current setting
    div(phiv,p)     Gauss upwind;
    div(phid,p)     Gauss upwind;
// new setting
    div(phiv,p)     Gauss limitedLinear 1.0;
    div(phid,p)     Gauss limitedLinear 1.0;

You do not have to stop your currently running simulations. OpenFOAM dictionaries are modifiable at runtime. You can:

1. login on the cluster
2. go to the run directory
3. for each running case
   3.1 open system/fvSchemes
   3.2 modify the two lines mentioned above and save the file
   3.3 move on to the next case

Let me know if this works for you.

Best, Andre

Test of different inlet values for k and omega

Hi Tushar,
I got feedback from a colleague that the inlet values for k and omega might be influencing the solution. In particular, high values for k might stabilize the solution too much. Therefore I suggest testing the following settings:

• take one case with the base mesh and fixed Re, Ma, and angle of attack
• vary k from a very small number, say 1.0E-5, up to the current one, which is about 320; the results for k=320 you already have; I suggest adding 3-4 cases with values in between
• the formula to compute omega you find here in the section Initialization; for the mixing length L, take 5% of the chord length
• visualize the fields nut, k, omega and observe the difference

Best, Andre

Compute estimate for cell-width at airfoil surface

Hi Tushar,

currently, the basic setup uses wall-functions to correct mu, k, and omega at the airfoil surface. Ideally, we would like to avoid the usage of wall functions because they introduce uncertainty to our simulation (we don't know how well assumptions made in derivation apply to our case). To see if it is feasible to resolve the boundary layer down to yPlus=1, it is useful to make an estimate for the necessary cell-width at the airfoil surface:

1. Follow this procedure to compute an estimate for the cell-width; report also intermediate results; you can use this syntax to render Latex formulas
2. check the yPlus field using paraview; select only the airfoil patch; report the yPlus values around the thickest part of the airfoil; get an estimate for the width of the first cell layer (load internal mesh; pick one cell that is roughly horizontal; use the extractSelection filter; look at the bounds in the information tab)
3. Compare the current cell-width with the width estimated in 1.

Best, Andre

Incorrect inlet velocity vector

Hi Tushar,

I noticed in one of your setups that the inlet velocity vector is not set correctly. E.g. here, the inlet velocity has a non-zero z-component.

The compute_inlet_velocity.py script prints the x and y components of the velocity vector in square brackets. The value in the next line is the angle of attack and not the z-component (it's a consistency check). The z-component is always zero in our simulations.

Please check your simulations.

Best, Andre

Adoption of final IDDES simulation setup and obtaining buffet onset angle of attack for constant Ma and Re numbers for NACA 0012-34 airfoil

Hello sir,

I did adopt the final IDDES simulation case setup of
Reduced Y Grading of 5000, using wall functions and 450 cells chordwise for NACA 0012-34 airfoil.

Then, I tried to obtained the onset buffet angle of attack for constant Re and Ma numbers specifically for NACA 0012 airfoil using the bisection method.

Thanks.

Best,
Tushar

Implementing IDDES approach for NACA0012 airfoil

Hello Sir,

As per your instructions in the pull request, I started cases of 2D and 3D setup based on IDDES approach.
I will let you know the results of the simulation soon.
Thanks.

Best regards,
Tushar

Extract reference data for drag and lift

Hi Tushar,

in [1], you find reference data for drag and lift coefficients of the airfoil investigated in this project:

• Fig 175 on page 284 provides drag coefficients for varying angle of attack and Mach number
• Fig 161 on page 274 provides lift coefficients for varying angle of attack and Mach number

To work with the reference data,

• extract the measurements using digitizer or a similar tool; after this step, you should have two or more textfiles holding the numerical values for drag and lift
• re-create the original plots using Python (the plots don't have to look exactly the same; only the provided information should be the same); Matplotlib is a great package to create visualizations; Pandas or Numpy can help you to load the data

Get in touch if you need help :-)
Best, Andre

reference 1

Trial of OnSet-Buffet cases

Hello Sir,

As per our discussion, I created the onset cases for dataset1 and started running the simulations on cluster.
The cases are as follows:

1. Case1: OnSet_Set1_Ma-0.75_Re-10Mio_GaussLinear_alpha-4_OriginalMesh
2. Case2: OnSet_Set1_Ma-0.75_Re-10Mio_GaussLinear_alpha-4_FineMesh
3. Case3: OnSet_Set1_Ma-0.75_Re-10Mio_GaussLinear_alpha-6_OriginalMesh
4. Case4: OnSet_Set1_Ma-0.75_Re-10Mio_GaussLinear_alpha-6_FineMesh

I am changing the Courant number after some time interval as well to accelerate the simulations.
I will plot the Cp, Cl and yPlus distribution graphs for onset as well as prebuffet cases very soon.
Thanks.

Best,
Tushar

Getting started

Hi Tushar,

to get familiar with the workflow and the simulation setup, here are some tasks to get on track:

• merge the pull request
• follow the instructions in the README file to run your first simulation
• while the simulation is running, use tail -f log.rhoPimpleFoam to follow the solver output inside the simulation folder; you can interrupt the simulation or the output of tail using Ctrl+C
• open the airfoil geometry in ParaView and take a look: paraview test_cases/naca0012-34-base/constant/triSurface/naca0012-34.stl
• add a file called naca_simulation_setup.md to the repository and create a list of all files in the setup; create a short description of the files' purpose; if you're uncertain, try to find a description online or post open questions in this issue; for basic Markdown syntax, see here

Best, Andre

Variation of Reynolds number

Hi Tushar,
great job starting the new setup on the cluster. Unfortunately, there is still no buffeting. To make the computation easier, I suggest going to lower Reynolds numbers. After merging the latest pull request, I would like you to try the following:

• create a refined version of the base case; to do so, double the number of cells per block in the blockMeshDict, something like

    // Number of cells
    yCells  320; // airfoil to far field
    yGap     16;  // resolution of trailing edge
    xUCells 120;  // upstream
    xMCellsLower 200; // middle low
    xMCellsUpper 200; // middle up
    xDCells 420; // downstream

• for default and fine mesh, vary the Reynolds number while keeping the Mach number and the angle of attack fixed at 0.85 and 5°; for the Reynolds number, I suggest values of Re={0.375, 0.75, 1.5, 3} Mio (8 simulations in total)
• for the fine mesh, use 8 MPI tasks; see here how to do this
• set the final simulation times to {2.4, 1.2, 0.6, 0.3}s and the writeInterval to {0.008, 0.004, 0.002, 0.001} (controlDict)
• to scale the mesh, e.g., to adjust the Reynolds number, change the scale parameter in blockMeshDict, e.g. scale 0.157637;
• if you change the scale parameter or the inlet velocity, you also have to update the parameters magUInf, lRef, and Aref in FOforces; lRef is the chord length (the same value as you use for scaling); Aref is the projected area in lift-direction (the chord length times the depth of our profile (depth is always 0.02)); magUInf is the magnitude of the inlet velocity vector
• the angle of attack in this setup is implemented by changing the direction of the inlet velocity vector rather than rotating the geometry; the script compute_inlet_vector.py helps you to determine the velocity vector for new speeds
• set the time parameter in jobscript to 48:00:00

How to create a new setup:

• make a copy of the base case in test_cases; give the folder a sensible name, e.g., naca-Re-0.375-Ma-0.85-fine
• adjust the settings in the new case
• make a copy of the jobscript and adjust the folder names
• submit the new jobscript

Let's see if we can find a suitable setup. Let me know if the instructions are not clear enough.
Best, Andre