Benchmark for scale-resolving simulation with curved walls: the Taylor Couette flow

Flow between rotating concentric cylinders, or the Taylor Couette flow, has been studied extensively because of its rich physics, ranging from axisymmetric steady laminar flow, to fully developed turbulent flow. In the present study, we advocate the use of this problem as a benchmark case for scale-resolving simulation, such as large eddy simulation (LES) and direct numerical simulation (DNS). The problem is attractive because of its simple geometry, simple boundary conditions, and complex physics involving wall-shear induced and centrifugal instability. Unlike the well-known fully developed channel flow, this problem has a curved wall boundary, and it is unnecessary to add a source term to the governing equations to sustain the fully developed turbulent flow. A p-refinement study for Re = 4000 is performed first to establish DNS data, including the time history of enstrophy, which can be used as an accuracy and resolution indicator to evaluate numerical methods, and is orders of magnitude faster than using the mean flow quantities and Reynolds stresses to evaluate solution quality. Finally, an hp-refinement study is performed to establish the relative accuracy and efficiency of high-order schemes of various accuracy.


Introduction
Scale-resolving simulations (SRS) such as large eddy simulation (LES) [1] and direct numerical simulation [2] have been used extensively to study fundamental turbulent flow physics at relatively low Reynolds numbers. Because of the disparate length and time scales in a turbulent flow, the spectral method has been the method of choice for earlier DNS studies involving relatively simple geometries, such as fully developed channel flow [3], and the Taylor Couette flow [4].
In order to make LES feasible for real-world applications, it is necessary to develop high-order numerical methods capable of handling complex geometries. The last two decades saw significant development efforts in the CFD community for such methods [5][6][7][8][9]. Multiple international workshops on high-order CFD methods [10] have been organized to evaluate high-order methods on a wide variety of benchmark problems, from simple 2D flow problems to fully turbulent 3D problems involving complex geometries.
The research effort developing LES tools got a huge boost from NASA's CFD Vision 2030 Study [11], which states "The use of CFD in the aerospace design process is severely limited by the inability to accurately and reliably predict turbulent flows with significant regions of separation. Advances in Reynolds-averaged Navier-Stokes (RANS) modeling alone are unlikely to overcome this deficiency, while the use of Large-eddy simulation (LES) methods will remain impractical for various important applications for the foreseeable future, barring any radical advances in algorithmic technology. Hybrid RANS-LES and wall-modeled LES offer the best prospects …" With the recent development in high-order numerical methods and computer hardware, especially GPU architectures, LES has been increasingly used in industrial applications, especially in turbomachinery [12,13], for mission-critical problems. Flow problems with Reynolds numbers of several million have been tackled successfully.
With more LES applications in industry, the need to assess the solution error, mesh resolution and numerical methods becomes more critical. However, it is hard to evaluate numerical methods for LES because one needs to first establish a statistically steady turbulent flow. After that, it is still necessary to continue the simulation for a significant amount of time so that the mean flow and Reynolds stresses are sufficiently converged in time. If h or p-refinement studies are desired, it is necessary to repeat this process many times, making such evaluations very expensive.
Out of the many benchmark problems from the international high-order CFD workshops, the Green Taylor vortex problem [14] has become the most popular in evaluating numerical methods for SRS. There are several reasons for that popularity. First, the problem has the simplest geometry of a Cartesian box. Any numerical methods, based on any grids, structured or unstructured, can be easily evaluated with this case. Second, the problem has a well-defined initial condition, which triggers flow transition to a fully turbulent state through a well-known physical path of vortex stretching and eventual breakdown, but not through roundoff errors in the computation. Third, the use of the enstrophy history appears to distinguish the resolution capability of the numerical methods very well. One can even compute the effective numerical viscosity in the simulation.
The Green-Taylor vortex does have its deficiency. First and foremost, it does not have any wall-boundaries. As a result, it cannot assess performance for shear-induced turbulence. One could of course reuse the well-known fully developed turbulent channel flow problem. We did not select this case because of two reasons: the need to add a source term in the momentum equation to sustain the fully developed turbulent flow, and the wall is not curved. Curved wall boundaries can also be used to assess highorder boundary representations.
Our search has led to the Taylor-Couette flow, which has a very simple geometry, and curved wall boundaries, but very rich physics [15]. In addition, one can easily specify a computational mesh using one parameter, the mesh growth ratio in the radial direction, since uniform meshes are used in the spanwise and azimuthal directions. This is an important consideration because meshes in the same refinement family should be used when different meshes are employed to evaluate numerical methods.
The paper is organized as follows. In Section 2, we briefly review the numerical method used in the present LES and DNS of the Taylor-Couette flow. After that, the Taylor-Couette flow problem is described in Section 3. Section 4 presents the numerical results from the LES and DNS, and a hp-refinement study to evaluate the relative accuracy and efficiency of high-order schemes. The conclusions are summarized in Section 5.

Numerical method
The Navier-Stokes solver is based on the flux reconstruction method originally developed by Huynh [16] for hyperbolic partial differential equations. This method was extended to mixed unstructured meshes in [17,18]. Many further developments were reviewed in [19]. This method belongs to discontinuous finite element methods, similar to the discontinuous Galerkin (DG) [20] and spectral difference [21] methods. However, this method also has some unique advantages. Here we present a brief introduction of the method starting from a hyperbolic conservation law governing inviscid flow with proper initial and boundary conditions, where the vector U consists of conservative variables, and F is the flux. By discretizing the computational domain with nonoverlapping elements, and introducing an arbitrary test function φ in each element, the weighted residual formulation of Eq. (1) on element V i can be expressed as Z The conservative variables inside one element are assumed to be polynomials, and expressed by nodal values at certain points called solution points (SPs). After applying integration by parts to the divergence of flux, replacing the normal flux term with a common Riemann flux F n com and integrating back by parts, we obtain Z Here, the common Riemann flux is computed with a Riemann solver where U i+ stands for the solution outside the current element, and n denotes the outward normal direction of the interface. The normal flux at the interface is: Note that if the face integral in Eq.
(3) can be transformed into an element integral then the test function will be eliminated. In order to do so, a "correction field" δ i is defined in each element as Z where ½ F n ¼ ½F n com −F n ðU i Þ is the normal flux jump. Eqs.
(3) and (6) result in The final formulation for each solution point j is where Π j denotes a projection to the polynomial space, and the subscript j denotes the j-th solution point in a certain element. For the viscous fluxes involving the gradient of conservative variables, we use the Bassi-Rebay 2 scheme [22]. No sub-scale stress models are used. Therefore, the present simulations are called implicit LES (ILES) or DNS when the simulation captures all relevant scales. For time integration, we employ the three-stage SSP Runge-Kutta scheme [23].

The setup of the Taylor Couette flow
Taylor-Couette (TC) flow, the name given to the flow between two coaxial co-or counter-rotating cylinders, is one of the classical problems used to study the physics of fluids. The problem carries the name of G.I. Taylor because of his landmark analysis of the instability and the discovery of the counter-rotating Taylor vortices when the Reynolds number is sufficiently large [15]. Theoretical, experimental and computational studies have been performed for the TC flow for many decades and interested readers can refer to a recent review article [24] for details.
In the present study, we assume that the inner cylinder rotates with an angular velocity of ω, while the outer cylinder is stationary, as shown in Fig. 1. The geometry of the flow is characterized by the radius ratio, η = r i /r o , with r i and r o being the inner and outer cylinder radiuses respectively. The wall velocity of the inner cylinder is U = ωr i . Let's first assume that the flow is incompressible with density ρ, dynamic viscosity of μ, and kinematic viscosity ν = μ/ρ. The width of the annulus is d = r or i . The Reynolds number is defined as Instead of using the Reynolds number, one can alternatively characterize the driving of TC flow by the Taylor number Another important geometric parameter is the aspect ratio, which is defined as where L is the length of the cylinder. In the present study, Γ is selected to be 2π based on an earlier DNS study [4] and our own tests. Since d is chosen as the primary length scale of the TC flow, other parameters such as the radius ratio and the aspect ratio determine the computational cost once the Reynolds number is specified. The domain size in the azimuthal direction is roughly π(r i + r o ), which can be expressed as πdð1þηÞ 1−η . It is easy to see that when η goes to 1, the domain size (with respect to d) approaches infinity. In other words, a narrower gap width incurs higher computational cost assuming the same Reynolds number. In the present study, we consider the wide gap case with η = ½. There have been many existing experimental and computational studies of this geometry [4].
When the Reynolds number is sufficiently low, an axisymmetric steady laminar solution exists. In fact, for an incompressible flow, the analytical velocity field in a polar coordinate system can be written as: The momentum equation reduces to A simple integration yields the following formula where C is an arbitrary constant. In the case of compressible flow, such simple analytical solutions have not been found. However, we could assume a low Mach number to approximate an incompressible flow. In addition, we need to fix at least the temperature or pressure of one wall to make the problem well-posed. Let's assume the inner wall temperature is T (with a Prandtl number of 0.71), with a local speed of where γ is the ratio of specific heats, and R the gas constant. The inner wall Mach number is The outer wall is assumed adiabatic. To approximate an incompressible flow, we assume M = 0.1. It is well known that the flow is fully turbulent when Re > 3000. In the present study, we select a Reynolds number of 4000 to make the computational cost reasonable for benchmark studies. In the axial direction (or spanwise direction), we employ a periodic boundary condition to approximate infinite cylinder length. The TC flow has no inlet or exit. Therefore, the mean density remains a constant forever. A fixed wall temperature will also fix the pressure at the inner wall, given the density.
In summary, we compute the laminar low Mach compressible TC flow in the following manner. Given a constant density, the wall temperature and Mach number, we can obtain the inner wall velocity U, and the angular velocity. In addition, the wall pressure is computed from the ideal gas law. From the wall pressure, we can determine C in (15). After that the pressure in the entire domain is known. Given density and pressure, we can compute the temperature. The velocity field is computed from (12) and (13). Finally, the viscosity is computed from (9).
The inner wall velocity can be used to define a time scale which is used to non-dimensionalize the physical time, i.e.,t ¼ t t Ã . We can of course start a LES/DNS from this laminar initial condition with no variation in the axial direction, and no radial flow velocity. The interactions between the round-off error, the truncation error and the physical instabilities at Re = 4000 will trigger the eventual transition from a laminar flow to a fully turbulent one. Once the flow becomes fully developed, we can extract mean forces, Reynolds stresses, and power spectral densities of certain flow variables to evaluate solution accuracy. Unfortunately, this approach is expensive, and often has uncertainty associated with the convergence of time-averaging.
Here is what we can learn from the popular Green-Taylor vortex problem. A welldefined initial condition appears to drive the transition process following a physical path of vortex stretching, tilting and breakdown. This process appears to severely test the resolution of the numerical methods. After some trial and error, we employ the following initial condition by adding variations to the steady laminar flow in all three directions, i.e., with ε = 0.1. We have assumed that z varies between 0 and L. The time histories of several global integral variables can be used for benchmarking purpose. The most used are the kinetic energy and enstrophy, defined for incompressible flow as where v ! is the velocity vector, and ω If the domain is periodic, one can derive the following well-known identity [25][26][27].
From the kinetic energy dissipation rate, we can estimate the Kolmogorov scale where ϵ is the dissipation rate of the turbulent kinetic energy, which can be computed as In (27), the overbar denotes time-averaging. The friction torque generated by the inner cylinder, T i , is computed according to with τ i the skin friction at the inner wall. The torque coefficient is defined as For a fully developed TC flow, the total kinetic energy approaches a constant because of the work done by the inner cylinder, ωT i . The effective dissipation rate becomes For the present benchmark problem, non-dimensional quantities are preferred. In addition, we include the density to take into account the density variation. The quantities with a hat denote non-dimensional quantities, i.e., where ρ is the average density. Then identities (25) become If we express the Kolmogorov scale in a non-dimensional form, we obtain 4 Numerical results and discussions

P-refinement study
We first conduct a p-refinement study to establish the resolution for DNS. The computational mesh is shown in Fig. 2, which has 64, 80 and 20 elements in the axial, azimuthal and radial directions respectively, resulting in a total of 102,400 elements. The mesh is uniform in the axial and azimuthal directions, and clustered in the radial direction with a minimum element size of 0.007d, and a growth ratio of 1.5 near the wall. The histories of the non-dimensional kinetic energy and enstrophy are displayed in Fig. 3. Except p2, the kinetic energy predicted by all other schemes agrees very well until t = 40. After that, even the p4 and p5 predictions have visible differences, indicating that the flow has become chaotic and turbulent by that time. On the other hand, the enstrophy histories show significant differences after t = 20, especially for p2 and p3. It is obvious that the enstrophy histories are approaching p-independence, and the difference between the p4 and p5 results is very small before t = 25. The p4 and p5 enstrophy histories have two distinctive peaks and one valley between t = 20 and 30. The first peak is at t 1 = 22.058, and the first valley is at t 2 = 24.016, and the second peak is at t 3 = 27.200.
Since there is a small difference between the p4 and p5 enstrophy histories in the neighborhood of t 3 , a finer mesh with 227,136 elements is used for a p5 simulation. The enstrophy histories on the base and the finer meshes are displayed in Fig. 4. Note that there is an excellent agreement between the results until t 4 = 33.5. After t 4 , the difference in the enstrophy histories is significant, indicating the flow has become irregular and turbulent.
Based on the enstrophy histories in the hp-refinement study, we can conclude that the flow appears to follow a definitive path from a smooth initial condition to a fully turbulent state. From t = 0 to t 1 , the enstrophy grows by more than an order of magnitude, indicating a strong vortex stretching/tilting process, which is a severe test for the resolution of numerical methods. To visualize the process, the distributions of the vorticity magnitude at t 1 from p2 to p5 simulations on the r-z plane (x = 0 cutting plane) are displayed in Fig. 5. Clearly thin shear layers are generated near the cylinder walls, and inside the flow field, indicating strong vortex stretching/tilting. Note that the flow field displays strong symmetric patterns, which indicate a laminar flow field at t 1 . It is obvious that the p2 and p3 simulations are not capable of fully resolving the thin vortex sheets, while the p4 and p5 simulations produce p-independent results. Between t 1 and t 2 , the enstrophy decreases, signaling a breakdown of vortex sheets into smaller eddies. Because of the close agreement between the p4 and p5 (on both the base and finer meshes) enstrophy histories until t 2 , we believe that the flow is still largely laminar at t 2 . To verify this claim, the distributions of the vorticity magnitude at t 2 from p2 to p5 simulations are displayed in Fig. 6. The strong symmetry in the distribution of vorticity magnitude suggests indeed a laminar flow at t 2 . One can clearly see the breakdown of the vortex sheets into smaller vortices. After the vortex breakdown, we expect the flow to transition into a turbulent one soon.
The increase in enstrophy after t 2 is a surprise because the phenomenon did not appear in the Green-Taylor vortex flow. We do have the following explanation. Soon after t 2 , the flow starts to transition into a nonsymmetric and chaotic state. The boundary layers near the inner and outer cylinders also start to transition into a turbulent one, creating strong shear layers which briefly drive up the enstrophy until t 3 . The distributions of the vorticity magnitude at t 3 from p2 to p5 simulations are presented in Fig. 7. We indeed see more clearly the stronger shear layer generated near the outer cylinder wall at t 3 than at t 2 . Note that nonsymmetric flow patterns start to appear, indicating transition to a chaotic and turbulent flow. It is obvious that vortex breakdown continues between t 2 and t 3 inside the flow domain, and visible differences appear between the p4 and p5 results.
After the second enstrophy peak at t 3 , the flow becomes more and more chaotic and turbulent. However, the enstrophy histories between the two p5 simulations on the Fig. 6 The vorticity distributions at t 2 computed with different schemes base and finer meshes remain quite close until t 4 . After that, we believe the flow becomes very irregular, and it is meaningless to compare instantaneous states. The distributions of the vorticity magnitude at t 4 from p2 to p5 simulations are displayed in Fig. 8. Note clearly that the flow has become non-symmetric, irregular, and the flow patterns are very different for all schemes.
We now suggest the following parameters to compare for this scale-resolving benchmark problem: Set the initial condition according to (18)- (21) Run the simulation until t = t 4 = 33.5. Compare the enstrophy history with the DNS result at p5. Compare the distribution of the vorticity magnitude with the DNS p5 result at t = t 2 .
After the establishment of the p5 simulation on the base mesh as a DNS, we continued the simulation for another 1400 time units. The histories of kinetic energy, enstrophy, and inner wall torque coefficient are displayed in Fig. 9. Although the enstrophy and torque appear to reach a statistically steady turbulent state around t = 100, it takes another 500 time units for the kinetic energy to reach a fully developed state. One can also estimate the Kolmogorov scale to be λ = 0.011d. The minimum element size in the wall normal direction is 0.007d. At this resolution for p = 5, there are roughly 10 DOFs per Kolmogorov scale, clearly enough to capture the smallest turbulent scale. To illustrate the flow vortical structures, the instantaneous and mean velocity vectors are plotted in Fig. 10. From the mean velocity vector plot, we clearly see six vortices (or three pairs), with each vortex having an aspect ratio close to 1. A similar vortical structure was also observed for Reynolds number of 8000 in [4]. This indicates that the shape and size of the vortices are not strongly dependent on the Reynolds number. Note that many smaller vortices exist in the instantaneous velocity plot, highlighting the chaotic and turbulent nature of the flow field.

HP-refinement studies
In this subsection, we perform h and p refinement studies to evaluate the relative accuracy and efficiency between different solution polynomial orders. We employ the enstrophy history until t = t 4 as the resolution/accuracy indicator. In order to guarantee that all the meshes are in the same refinement family, the finer meshes are generated through uniform refinements of the coarser meshes. The coarsest mesh has 12,800 elements (with 10 elements in the wall normal, 40 elements in the azimuthal and 32 in the axial direction respectively), as shown in Fig. 11. The mesh has a minimum element size of 0.032d in the wall normal direction with a growth ratio of 1.8. The coarsest mesh is called the level 1 mesh. Uniform refinements are then used to generate three extra levels of finer meshes. Therefore, the level 2 mesh has 102,400 elements, the level 3 mesh has 819,200 elements, while the level 4 mesh has 6,553,600 elements.
We then perform a p1 simulation on the level 4 mesh. The total number of DOFs per equation for this simulation is roughly 52.4 M. Based on earlier benchmark results [10], the resolution of this simulation is similar to a 2nd-order FV simulation with 52.4 M control volumes. After that, we perform higher p simulations on the coarser meshes to obtain enstrophy histories to estimate the solution resolution and accuracy. The  Fig. 12. Note that the p2 simulation on the level 2 mesh has higher resolution than the p1 simulation on the Level 4 mesh. The p2 simulation has a total of 2.76 M DOFs, which is roughly 1/19 of the DOFs in the p1 simulation. It is interesting to see that the p3 simulation on the level 1 mesh is comparable to the p2 simulation on the level 2 mesh. The number of DOFs in the p3 simulation is only 0.82 M, which is 1/64 of that in the p1 simulation.  At p > 2, refining the mesh once produces results more accurate than by increasing the order by 1. For example, the p3 (p4) simulation on the level 2 mesh is more accurate than the p4 (p5) simulation on the level 1 mesh. We do note that the p5 simulation on the level 1 mesh is more accurate than the p3 simulation on the level 2 mesh. In other words, the p3 simulation with 2.4 times the DOFs of the p5 simulation is less accurate than the p5 simulation. Summarizing all the results together, we can conclude that the p5 simulation with 1/150 the DOFs of the p1 simulation is more accurate than the p1 simulation. To obtain the DNS resolution using the p1 simulation, at least 3.3 billion DOFs are needed for this relatively simple benchmark problem.

Conclusions
In the present study, the Taylor-Couette flow at the Reynolds number of 4000 is suggested as a benchmark test case for scale resolving simulation with curved wall boundaries. Similar to the well-known Taylor Green vortex problem, the enstrophy history is used as the resolution and accuracy indicator. The use of enstrophy can dramatically reduce the cost of performing benchmark studies for LES and DNS. The usual approach for such benchmarks involves comparing the mean flow and Reynolds stresses, which are very time-consuming to obtain, and dependent on the duration of the simulations. With a well-defined smooth initial condition, it appears a definitive path can be identified for the flow to undergo vortex stretching/ tilting, breakdown, and finally transition to a fully turbulent flow. Through a prefinement study, we have identified one such smooth initial condition, and the resolution needed for DNS.
This benchmark is then used to perform an hp-refinement study to evaluate the relative accuracy and efficiency of FR/CPR schemes of various orders of accuracy. A family of meshes generated through uniform refinements are used in this hp-refinement study. It is shown that the p2 simulation on the level 2 mesh is more accurate than the p1 simulation (2nd order accurate) on the level 4 mesh. The total number of DOFs of the p1 simulation on the level 4 mesh is roughly 19 times of that of the p2 simulation on the level 2 mesh. It is also shown that the p3 simulation on the level 1 mesh has a similar accuracy as that of the p2 simulation on the level 2 mesh. In other words, to achieve a similar accuracy as the p3 simulation, the p1 simulation should have at least 64 times of the DOFs in a p3 simulation. Results at higher p orders indicate that the p1 simulation needs at least 150 times the number of DOFs to achieve a similar accuracy to a p5 simulation. In other words, a p5 FR/CPR simulation on a mesh of 47,000 elements should be more accurate than a 2nd-order FV simulation on a mesh of the same family with 1.5 billion control volumes.

Availability of data and materials
The benchmark data will be available for others to download.