**2. Materials and Methods**

The present study made use of the high-fidelity implicit large-eddy simulation (iLES) approach based on spectral/hp element methods (SEM) [23]. Broadly speaking, this approach relies on numerical dissipation in lieu of a traditional turbulence model, no subgrid-scale model is employed, and has become very popular in recent decades. In this context, iLES has also been called under-resolved direct numerical simulation (uDNS) [24]. For a detailed study on why and how to use SEM-based iLES/uDNS, the reader is referred to references [25–36].

The SEM adopted was the high-order continuous Galerkin (CG) method coupled with a novel stabilization scheme particularly suited for the simulation of transitional and turbulent flows [37,38]. This technique is called gradient-jump penalty (GJP) and has even been proven superior to modern versions of the well-known spectral vanishing viscosity (SVV) operator [39]. Simulations were conducted with the Nektar++ open-source code [40,41], a mature platform that has been successfully employed multiple times for iLES/uDNS [42–45].

The adopted unsteady CG solver [46] is based on a spectral/hp element method that provides a continuous Galerkin projection as a linear combination of a set of a particular polynomial basis functions Φ*n*, such that:

$$\mu^{\delta} = \sum\_{n \in N} \Phi\_n \hat{\mathfrak{a}}\_n \tag{1}$$

where *u*ˆ*<sup>n</sup>* are the coefficients of our solution vector **uˆ** to the approximate solution *u<sup>δ</sup>* in the reduced finite-dimensional function space *U<sup>δ</sup>* ⊂ *U* for a linear differential equation type:

$$\mathbf{L}(\mathbf{u}) = f \,, \tag{2}$$

which can be replaced by our incompressible Navier–Stokes equations to be solved, i.e.,

$$\frac{\partial \mathbf{u}}{\partial t} = \mathbf{N}(\mathbf{u}) - \nabla p + \nu \mathbf{L}(\mathbf{u}) \, \text{ \,} \tag{3}$$

$$\nabla \cdot \mathbf{u} = 0 \,,\tag{4}$$

where:

$$\mathbf{N}(\mathbf{u}) = -(\mathbf{u} \cdot \nabla)\mathbf{u}\,,\tag{5}$$

$$\nu \mathbf{L}(\mathbf{u}) = \nu \nabla^2 \mathbf{u} \ . \tag{6}$$

They are solved along with the appropriate boundary conditions described later on. The pressure is solved by the following velocity correction scheme:

$$\int\_{\Omega} \nabla p^{n+1} \cdot \nabla \phi d\Omega = \int\_{\Omega} \nabla \cdot \left(\frac{-\hat{\mathbf{u}}}{\Delta t}\right) d\Omega + \int\_{\Gamma} \phi \left[\frac{\hat{\mathbf{u}} - \gamma\_0 \hat{\mathbf{u}}^{n+1}}{\Delta t} - \nu (\nabla \times \nabla \times \mathbf{u})^{\*}\right] \cdot \mathbf{n} \, dS \,, \tag{7}$$

where Ω is the domain and Γ is the boundary domain, with a backward approximation of the time derivative, such that:

$$
\hat{\mathbf{u}} = \mathbf{u}^{\ominus} + \Delta t \,\mathbf{N}^\* \,\mathrm{.}\tag{8}
$$

where (∗) denotes the time extrapolation, and (⊕) the backward differencing. It uses a second-order implicit–explicit (IMEX) time-integration scheme. The velocity is the solution of the Helmholtz problem given by:

$$\frac{\gamma\_0 \mathbf{u}^{n+1} - \hat{\mathbf{u}}}{\Delta t} = -\nabla p^{n+1} + \nu \mathbf{L}(\mathbf{u}^{n+1}) \,. \tag{9}$$

### *2.1. Computational Domain and Baseline Mesh*

The simulation domain is initially based on a two-dimensional mesh of quadrilaterals, as shown in Figure 3, which is then extruded along the spanwise direction in a way that allows for the desired waviness to be incorporated. The dimensions of the domain, i.e., the distances upstream and downstream of the body, the height and the width were based on dimensions previously established as sufficient to avoid blockage, with a blockage ratio (diameter/height) less than 3.4%, and allow the proper development of the wake [47,48].

**Figure 3.** Two-dimensional mesh. It is worth noting that this mesh has not yet received p-type refinement (*Np* = 2), that is, in this case, in each direction the elements still will be "divided" in half, then increasing its resolution.

Note the mesh is of a structured pattern around the cylinder and along the wake, whereas elsewhere it is of an unstructured pattern. In SEM, the solution inside each element is represented through a polynomial expansion. Here, the element-wise polynomial order has been chosen as *P* = 2 (parabolic solution per element). Note that CG's nominal order of accuracy is *P* + 1, hence all the solutions presented here are nominally third-order accurate.

The circular cylinder is placed at the center of the domain (origin of Cartesian coordinates) and has a diameter *D* of unit size. The outermost boundary is made up of a semi-circle (inflow boundary of radius 15*D*) and a rectangle whose downstream closure coincides with the outflow boundary (at *x* ≈ 21*D*). The intermediate part of the mesh, as also shown in Figure 3, has another semi-circular section (of radius 4*D*) that connects to two symmetric triangular sections. These exist as means for the mesh to transition into the wake region, which extends downstream until the outflow boundary.

Around the cylinder, there is an outer circle (of radius 1.5*D*) and an inner circle (of radius *D*/2 + *δBL*) placed very close to the cylinder's surface, as shown in Figure 4. The latter

is devoted to capturing the (laminar) boundary layer of the simulation, whereby its thickness is adjusted based on the cylinder's Reynolds number *ReD* according to the estimate

$$
\delta\_{\rm BL} \approx D \left( \text{Re}\_D \right)^{-1/2} \text{ }, \tag{10}
$$

which was found to approximately match the boundary layer thickness next to its separation point. This innermost circular layer is subdivided into 5 equispaced layers that are deemed sufficient to resolve the cylinder's laminar boundary layer, given the high-order nature of the elements within the equispaced layers. For example, when using *P* = 2 elementwise polynomials in the solution, the innermost circular region will have 2 × 5 = 10 DOF (degrees of freedom) across the radial direction.

It is also worth mentioning that, in terms of resolution power, a SEM-type DOF is somewhat superior to the DOF of classical CFD schemes, as can be seen in dispersion/diffusion analyses [25,27]. In fact, these considerations apply to all the elements in the domain, which is why the mesh in Figures 3 and 4 may seem too coarse for those unfamiliar with SEM. If the mesh was shown with "subelements" representing the element-wise polynomial DOFs, it would look much finer.

**Figure 4.** Close-up view of the mesh next to the cylinder. It is worth noting that this mesh has not yet received p-type refinement (*Np* = 2), that is, in this case, in each direction the elements still will be "divided" in half, increasing their DOF.

From the edge of the boundary layer region, the radial mesh size grows geometrically toward the outer circle around the cylinder, when the radial mesh size achieves the same mesh spacing employed horizontally along the wake (≈0.25*D*). Along the cylinder's wall, in the azimuthal direction, 54 equispaced elements (or angular subdivisions) are used to cover the whole circle. This also defines the azimuthal size of the mesh along the outer circle, which in turn defines the vertical mesh size along the wake (≈0.35*D*).

As mentioned previously, the regions within the cylinder's outer circle, along with the two symmetric triangular regions plus the wake region follow a structured mesh pattern. The other regions of the domain have a mesh with an unstructured pattern. All meshes have been generated with GMSH [49], an open-source finite element mesh generator. A thirdorder mesh representation was used in all the cases, i.e., one order above the polynomial order adopted for the numerical solution near the surface, as per typical SEM guidelines for simulating flows around curved geometries.
