Looking at Figure 1 of Park et al. (2017), or even opening the mesh file from the supplemental material of the same paper with GMSH, it is evident that the mesh is built starting with a particular spacing near the wall (which, following the reference, is chosen in order to obtain a y+ close to 1) which grows with a given factor up to a certain height from the airfoil. However, from this length forward, there seem to be some other layers growing with other factors, greater than the previously mentioned. Farther away from the airfoil (e.g. at the boundaries), the length of the sides of the elements in the stream- and crosswise seems to be very big (definitely more than one chord - I suppose that the growth of elements far away from the boundary should have big growth factor).

My question is: what is the best practice for defining what the distances between the different layers should be when building meshes for PyFR? From what I understood, the direction of the refinement behind the airfoil should depend on the angle of attack with respect to the incoming flow (this is also one intuitive concept). Moreover, it is clear that the actual grid length should be defined in accordance with the polynomial order to be used in the code (in fact each element should accommodate possibly more than one node). However, I am more interested in understanding how much space one should care about in the refinement. I imagine that the number of layers and their â€śheightsâ€ť from the airfoil should be defined by the physics of the problem (Ma and rho at free-stream, viscosity, desired y+, and the Re) but I need to understand how to build a good mesh and its refinement. Are there any empirical relations that I should care about or it only depends on trial-and-error?

For the boundary layer I would start with having the first interior point inside the first element at a y+ = 1. With Gauss-Legendre solution points at p = 4 this works out at the first element having y+ = 20. Then, proceed to add layers with an expansion ratio of 1.1 to 1.2. Usually 10 layers is a good first guess. Once youâ€™ve got your simulation up and running you can compute the numerical y+ using the real skin friction rather than an estimate. This can be used to see if additional refinement is required.

The size and number of elements in the wake depends very much on the physics you are looking to capture. If you are mostly interested in quantities on the surface then you can often get away with less resolution than if you want to take statistics in the wake itself. Since the elements in the wake tend to be much larger than those in the boundary layer they are comparatively cheap and so there is usually no harm in adding extra resolution here.

To see how you are doing in the wake, as in how close you are to DNS levels of resolution, you can compute the numerical Kolmogorov length scale. This will let you know what kind of physics can be resolved in various parts of the wake.

The first solution point is at a y+ of 1. This gives an overall element size of approximately y+ = 20. This is a trivial consequence of the fact that at p = 4 the first Gaussâ€“Legendre point (over [-1,1]) is at 0.906 or ~5% of the overall domain.

Thanks for the answer, but i still dont fully understand this. If i understood correctly to get the first point of the gauss-legendre quadrature to be at 0.906 in an interval of [-1,1] one would need 5 points inbetween this interval. This would result in a hexahedral element to have 125 solution points which seems to be correct. But i cannot seem to find a source that says that the gauss-legendere quadrature at p=4 is with n = 5 (points). Im only able to find this:

"Antialiasing via an L2 projection on the P3 and P4 meshes, with third-order and fourth-order solution polynomials, respectively, refers to use of 9th-order and 11th-order Gaussâ€“Legendre quadrature rules,
"
I am using antialiasing, this would mean that an 11th order quadrature is used, but are those points solution points ? If so, the first solution point would be at ~ Â± 0.978. Thanks for you help.

A degree four Lagrange interpolating polynomial (p = 4) has five points. Anti-aliasing does not impact the ultimate ability of the scheme to resolve the solution as this is set by the number of solution points rather than the degree of anti-aliasing.

Alright, so the default behaviour of the solver is to place solution points with the Gaussâ€“Legendre quadrature rule depending on the solution polynomial order. Thanks for youre help, very much appreciated.