Abstract
We present FLEW, an in-house high-fidelity solver for direct numerical simulation (DNS) of turbulent compressible flows over arbitrary shaped geometries. FLEW solves the Navier–Stokes equations written in a generalized curvilinear coordinate system, in which the surface coordinates are non-orthogonal, whereas the third axis is normal to the surface. Spatial discretization relies on high-order finite-difference schemes. The convective terms are discretized using an hybrid approach, combining the near-zero numerical dissipation provided by central approximations with the robustness of weighted essentially non-oscillatory (WENO) schemes, required to capture shock waves. Central schemes are stabilized using a skew-symmetric-like splitting of convective derivatives, endowing the solver with the energy-preserving property in the inviscid limit. The maximum order of accuracy is eighth for central schemes (also used for viscous terms discretization) and seventh for WENO. The code is oriented to modern high-performance computing (HPC) platforms thanks to message passing interface (MPI) parallelization and the ability to run on graphics processing unit (GPU) architectures. Reliability, accuracy and robustness of the code are assessed in the low-subsonic, transonic and supersonic regimes. We present the results of several benchmarks, namely the inviscid Taylor–Green vortex, turbulent curved channel flow, transonic laminar flow over a NACA 0012 airfoil and turbulent supersonic ramp flow. The results for all configurations proved to be in excellent agreement with previous studies.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
1 Introduction
Compressible flows over curved surfaces are ubiquitous in the aerospace field, both in external (e.g., aircraft wings) and internal (e.g., turbomachinery) configurations. The dynamics of these flows involves complex phenomena, such as shock wave/boundary layer interactions, non-uniform pressure gradients and centrifugal instabilities. Common feature of these phenomena is the generation of a highly unsteady flow field. Besides requiring further fundamental research, these phenomena make it challenging to predict turbulence intensities, skin friction, wall pressure fluctuations and heat transfer. Concerning technology applications, mispredicted flow instabilities can be detrimental to structural integrity and aerodynamic performance. A key tool in this context is the direct numerical simulation (DNS), whose predictive capabilities can aid both fundamental research and technological development. However, numerical algorithms for compressible flows simulations are less standardized then those for incompressible flows [1]. The hyperbolic nature of compressible Navier–Stokes equations allows for the presence of propagating disturbances and discontinuities, which interact with the underlying turbulence. Accurate simulations of such problems are challenging due to the contradictory requirements of numerical methods suited for turbulence, which must minimize any numerical dissipation that would otherwise overwhelm the small scales, and shock-capturing schemes, which must introduce numerical dissipation to stabilize the solution across discontinuites [2].
Further complications arise from geometric complexities. These are typically addressed using local mesh refinement in Cartesian frameworks or unstructured meshes, which, in turn, causes high computational cost or decreased accuracy, respectively [3]. In recent years, the immersed boundary method (IBM) has gained popularity and it has been successfully used in conjunction with turbulence models to simulate high-Reynolds-number flows [4]. Within DNS, though, the IBM struggles to provide a sufficiently accurate description of the boundary layer, especially in presence of separated flow [5]. A different approach to deal with complex geometries relies on the use of body-fitted meshes in a generalized curvilinear coordinate framework. This technique ensures most accurate simulation of the fluid dynamics near the wall and straightforward implementation of high-order finite-difference (FD) schemes.
Among FD approximations, central schemes are preferable as they do not introduce any artificial dissipation, but require special care for their susceptibility to numerical instability. This non-linear instability is due to the accumulation of aliasing errors resulting from discrete evaluation of the convective terms [6]. To address this problem, Pirozzoli [7] proposed an energy-preserving formulation of convective terms, based on fully splitting the derivatives of triple products through the product rule. This technique allows to capture the whole range of turbulence scales without needing any additional stabilization, such as artificial diffusivity [8], explicit filtering of physical variables [9] or relying solely on WENO schemes for convective terms discretization [10]. For years these stabilization strategies have set the standards for DNS of compressible flows, especially in the case of body-fitted grids, as most of the existing energy-preserving schemes are only formulated in Cartesian coordinates [11]. Another appealing property of energy-preserving formulations is that they can be efficiently combined with shock-capturing methods, yielding hybrid schemes that currently represent an optimal strategy for the computation of shocked flows [12].
Within this framework, we developed the high-fidelity solver FLEW, tailored for DNS of turbulent compressible flows over complex geometries. The solver relies on a third-order Runge–Kutta scheme for time integration and high-order FD methods for spatial discretization. Using generalized curvilinear coordinates, a body-fitted grid is transformed from the physical space, where its cells are skewed, to the computational space, where it appears as a regular hexahedron with cubic cells. This process allows to approximate spatial derivatives with FD methods originally developed for Cartesian frameworks. The switch between physical and computational space is carried out through the metric terms, which are computed by numerically deriving the physical mesh coordinates. FLEW employs message passing interface (MPI) parallelization and supports graphics processing unit (GPU) architectures, essential requirements to access the latest generation of high-performance computing (HPC) clusters. This paper is organized as follows. In Sect. 2 we describe the key points of the implemented algorithms, providing a thorough explanation of how energy-preserving schemes are tailored to a curvilinear framework in a computationally efficient manner. In Sect. 3 we present an extensive validation of the code, demonstrating its capabilities in a wide range of flow regimes and configurations. Results of the inviscid Taylor–Green vortex, turbulent curved channel flow, transonic laminar airfoil flow and supersonic turbulent ramp flow are discussed.
2 Methodology
2.1 Balance Equations
FLEW solves the compressible Navier–Stokes equations for a perfect gas written in generalized curvilinear coordinates. The Cartesian coordinates \((x_i)\) of a body-fitted grid are transformed into boundary-conforming curvilinear coordinates \((\xi _j)\) through the mapping \(x_i(\xi _j)\), with \(i,j=1,2,3\). Thus, the boundary surface is composed of segments of coordinate surfaces \(\xi _j=\text {const}\). In such a framework, assuming a stationary grid, the Navier–Stokes equations read
Here and throughout the paper the repeated index convention is used to indicate summation. The vectors of conservative variables, convective and viscous fluxes are
where \(E=c_vT+u_iu_i/2\) is the total energy, \(H=E+p/\rho\) the total enthalpy, \(u_i\) the velocity component in the \(i^{th}\) Cartesian direction, \(\hat{u}_j=u_i\hat{J}_{ji}\) the velocity component (normalized by J) in the \(j^{th}\) curvilinear direction, with \(\hat{J}_{ji}=J_{ji}/J\), being \(J_{ji}=\partial \xi _j / \partial x_i\) the Jacobian matrix of the transformation \(\xi _j(x_i)\), and J its determinant. The constitutive equations defining the viscous stress tensor and heat flux vector are
where \(\mu\) is the dynamic viscosity and \(\kappa =c_p\mu /Pr\) the thermal conductivity. Due to the perfect gas assumption, the Prandtl number is constant and fixed to \(Pr=0.72\). The system is completed by the equation of state \(p=\rho RT\). For the viscosity–temperature relationship, both the power law and Sutherland’s law are available. The Navier–Stokes equations are solved in non-dimensional form. From the minimal set of four independent reference variables, \(\rho _0, T_0, l_0, R\), the other reference quantities are derived as \(p_0=\rho _0 R T_0,~u_0=\sqrt{RT_0},~\mu _0=\rho _0 u_0 l_0,~t_0=l_0/u_0\), where \(l_0\) and \(t_0\) are the reference length and time.
2.2 Spatial Discretization
FLEW relies on high-order FD schemes designed for uniform stencils. By recasting the governing equations from Cartesian into curvilinear form (1), skewed input cells are stretched into cubic uniform cells, allowing the application of standard FD schemes. In the following, we will refer to the Cartesian system \((x_i)\) as ‘physical space’ and to the curvilinear one \((\xi _j)\) as ‘computational space’.
2.2.1 Convective Terms
The convective derivative in the \(j^{th}\) direction (temporarily omitting pressure forces) can be expressed as
where \(\varphi\) is the transported scalar quantity, i.e., 1 for the continuity equation, \(u_i\) for the momentum equation, H for the energy equation. Convective derivatives are manipulated according to the work of Pirozzoli [13] who, after assessing several energy-preserving types of splitting, proposed the following as the most robust:
The numerical discretization relies on identifying a numerical flux at the intermediate nodes, such that
An approximation of convective derivatives that is locally conservative (and computationally efficient) is based on the following expression for the numerical flux:
where
is the two-point, three-variable discrete averaging operator, \(a_{\ell }\) are the standard coefficients for central approximations of the first derivative, and 2L is the order of accuracy, which can be set by the user from second to eighth. The scaling of the order of accuracy was verified by Pirozzoli [13], who performed a grid sensitivity study using the same numerical scheme employed here.
2.2.2 Metric Terms
The metric terms are the elements of the Jacobian matrix, \(J_{ji} = \partial \xi _j / \partial x_i\), of the coordinate transformation \(\xi _j(x_i)\). Being generally unknown, \(J_{ji}\) is obtained by inverting the Jacobian matrix of the inverse transformation \(x_i(\xi _j)\). That is
For the present purposes, we assume that the coordinate transformation is locally invertible and does not depend on time. The inverse matrix \(\partial x_i / \partial \xi _j\) is computed by numerically deriving the Cartesian coordinates of the mesh along each \(\xi _j\) direction. A crucial aspect in evaluating metrics is to preserve the freestream. To address this problem, we use the same approximations as employed for the convective terms. Thus, metric terms are evaluated at each node \(\mathcal {N}\) as
where \((x_i)_{\mathcal {N}}\) are the Cartesian coordinates of the grid node \(\mathcal {N}\), and \(\hat{{f}}_{j;1/2}\) is the numerical flux taken at the intermediate node. The subscript j; 1/2 denotes a shift by 1/2 about \(\mathcal {N}\) in the positive \(j^{th}\) direction. Spacings in the computational space are set to be all unity, thus grid nodes correspond to integer values of \(\xi _j\). The numerical flux is computed as
where the accuracy order must match that of convective derivatives. Using these scheme, the following metric identities are numerically satisfied:
which is necessary and sufficient condition to ensure freestream preservation [9].
2.2.3 Shock Capturing
The locally conservative formulation (7) allows straightforward hybridization of central approximations with shock-capturing methods, based in FLEW on Lax–Friedrichs flux vector splitting and weighted essentially non-oscillatory (WENO) reconstruction. To implement such schemes, the characteristic form of the Euler equations must be expressed in the curvilinear system \(\xi _j\). This is possible by introducing a metric normalization such that the physical velocity components projected along \(\xi _j\) appear in the vectors of convective fluxes. That is
where \(\tilde{J}_{ji}\) are the metric terms normalized by \(m_j\) (e.g., \(\tilde{J}_{1i}=J_{1i}/m_1\)), defined as
With such normalization, \(\tilde{u}_j=u_i \tilde{J}_{j i}\) is the \(j^{th}\) component of the contravariant velocity, i.e., the projection of the velocity vector along \(\xi _j\). The Jacobian matrix of the normalized flux vector can be diagonalized as \({\partial \tilde{\textbf{F}}_j}/{\partial \textbf{Q}}= \textbf{R}_j\mathbf {\Lambda }_j\textbf{L}_j\), where \(\textbf{R}_j\) and \(\textbf{L}_j\) are the right and left eigenvectors matrices, and \(\mathbf {\Lambda }_j\) is the diagonal matrix of the eigenvalues. For each \(j^{th}\) direction, left and right eigenvectors are evaluated at the interface (denoted by the subscript j; 1/2), using the Roe’s average for flow variables and the arithmetic mean for metric terms. Positive and negative flux components are projected along the characteristic directions using a local Lax–Friedrichs flux splitting:
where \(\left| \lambda _{\max }\right|\) is the maximum eigenvalue over the entire stencil. First, the characteristic fluxes are consistently normalized as
Then, the \(\hat{\textbf{g}}_{j}^{\pm }\) are reconstructed by non-linearly combining the interpolation polynomials on WENO sub-stencils. Finally, the WENO reconstructions \(\hat{\textbf{g}}_{j}^{\textrm{WENO}\pm }\) are projected back into physical space, obtaining the numerical flux at the interface:
The accuracy of WENO schemes reaches up to seventh order. The transition between energy-preserving and shock-capturing discretization is controlled by an improved version of the Ducros shock sensor [14]:
where \(\omega _i\) is the vorticity and \(u_0\) is the reference velocity. The sensor is designed to be \(\theta \approx 0\) in smooth flow regions and \(\theta \approx 1\) in the presence of shock waves. Fixed a customized threshold \(\theta _0\), the shock-capturing scheme is locally activated if \(\theta >\theta _0\).
2.2.4 Viscous Terms
Spatial derivatives of viscous terms are approximated using central schemes, which are applied in the computational space. Our goal is to find the analytic relation to reconstruct derivatives from the computational into physical space. First, let us consider the viscous terms in Cartesian coordinates and expand them to Laplacian form,
to avoid odd–even decoupling phenomena. As for the energy equation
The relationship between the first derivative of a generic variable \(\varphi\) expressed in physical (\(x_i\)) and computational (\(\xi _j\)) space is obtained applying the chain-rule:
while the Laplacian, assuming that the \(\zeta\) axis is normal to the surface, is written as
where
Equations (21) and (22) allow us to reconstruct in physical space derivatives that are approximated in computational space, where central schemes can be applied. First and second derivatives in the \(j^{th}\) direction are computed at the node \(\mathcal {N}\) as
where \(b_{\ell }\) are the coefficients for second derivatives. The accuracy order ranges from second to eighth.
2.3 Characteristic Boundary Conditions
In addition to physical boundary conditions, FLEW features dynamic conditions based on characteristic wave relations [15], namely non-reflecting for far-field boundaries and reflecting for no-slip walls. Considering the boundaries at the coordinate surface \(\xi =\textrm{const}\), the convective derivative in the \(\xi\) direction reads
where we split the contribution of the convective flux in Cartesian coordinates, \({\textbf{F}}_i\), from that of the metric terms, \(\hat{J}_{1i}\), associated to the \(\xi\) direction. Expanding the product derivative, we get
The first term of the right-hand side can be recast as
where the Jacobian matrix of the flux, \({\partial \tilde{\textbf{F}}_1}/{\partial \textbf{Q}}\), has been diagonalized as outlined in paragraph 2.2.3. Exploiting the local one-dimensional approximation, the conservative variables are updated at the boundary according to the following equation:
where \(\mathcal {L}_{1}=\mathbf {\Lambda }_{1}\textbf{L}_{1} {\partial \textbf{Q}}/{\partial \xi }\) has the physical interpretation of wave amplitudes associated with the characteristic field. Analogous considerations can be made for boundaries at the surface \(\eta =\textrm{const}.\)
2.4 Time Integration
Spatial discretization leads to a semi-discrete system of ordinary differential equations:
where \(\textbf{R}\) is the vector of the residuals. The system is advanced in time using a three-stage, third-order Runge–Kutta scheme [16]:
where \(\textbf{Q}^{(0)}=\textbf{Q}^n\) and \(\textbf{Q}^{(3)}=\textbf{Q}^{n+1}\) are the solution vectors at the time-step n and \(n+1\), respectively. The integration coefficients are \(\alpha _{\ell }=(0,17 / 60,-5 / 12)\) and \(\beta _{\ell }=(8 / 15,5 / 12,3 / 4)\). The time-step \(\Delta t\), which is updated at each iteration by default, is chosen as the minimum value among the convective, viscous and thermal limitation:
which are computed, in turn, as the minimum value among all cells. Respectively
where \(j=1,2,3\) is the curvilinear direction, c is the speed of sound, \(\nu =\mu /\rho\) is the kinematic viscosity, \(\alpha = \kappa /\rho\) is the molecular conductivity. The Courant–Friedrichs–Lewy \((\textrm{CFL})\) number can be set by the user (a value close to unity is usually optimal).
3 Validation
FLEW is designed to perform DNS for various flow configurations, including the curved channel, compression ramp and wing airfoil. In this section, we validate the solver for these configurations by comparing the results to experimental and numerical data available in the literature. In addition, we use the inviscid Taylor–Green vortex on a wavy mesh as a benchmark.
3.1 Inviscid Taylor–Green Vortex
The energy-preserving properties of the solver are tested using the inviscid Taylor–Green flow, proposed by Duponcheel et al. [17] as time-reversibility test. The solution is computed in a triply periodic domain consisting of a wavy two-dimensional mesh, depicted in Fig. 1, which is extruded in the third dimension. The mesh is generated according to the following mapping [8]:
where \(L_x=L_y=L_z=2\pi\) and \(N_x=N_y=N_z=32\). The flow is initialized as follows:
where \(u_0=M_0c_0\) is the reference velocity (the Mach number is \(M_0=0.01\)) and \(c_0,~p_0,~\rho _0\) are the reference speed of sound, pressure and density.
Due to the lack of viscosity, from the ordered initial conditions small-scale turbulence structures are rapidly formed with incurred growth of vorticity. The solution is advanced up to time \(t/t_0 = 8\), at which the velocity field is reversed, and then further advanced up to \(t/t_0 = 16\) (the reference time and length are \(t_0=l_0/u_0\) and \(l_0=1\)). Based on the time-reversibility property of the Euler equations, the initial conditions should be exactly recovered.
In Fig. 2, we display a sequence of snapshots of the flow field, which qualitatively show that the initial coherent structures are faithfully reproduced at the end of the simulation. Simulation results are depicted in Fig. 3, where we report the time evolution of total kinetic energy and enstrophy. The kinetic energy is perfectly conserved over time, showing that no numerical dissipation is introduced. The total enstrophy grows markedly up to the reversal time (\(t/t_0 = 8\)), at which it begins to decrease following a symmetrical trend, with exact recovery of the initial condition at the final time (\(t/t_0 = 16\)).
3.2 Turbulent Curved Channel Flow
We focused on validating FLEW in the presence of pressure gradients due to surface curvature. For this purpose, we considered the classical study by Moser & Moin [18], where pressure-driven flow in a mildly-curved channel was simulated. The same flow configuration was recently used by Brethouwer [19] as benchmark.
The computational domain is bounded by sectors of concentric cylinders (a scale drawing appears in Fig 4). The curvature parameter is \(\delta /r_c=1/79\), where \(\delta\) is the channel half-width and \(r_c\) is the radius of curvature measured at the centreline. The domain has a length of \({4}\pi \delta /3\) in the spanwise (z) direction and subtends an angle of \(0.16~\textrm{rad}\) in the azimuthal (\(\theta\)) direction, which yields a length of \(12.64\delta\) along the centreline. An imposed mean-pressure gradient in the azimuthal direction drives the flow. No-slip isothermal conditions are imposed at the walls, while periodic boundary conditions are enforced in the azimuthal and spanwise directions, so that the simulated flow is the fully evolved turbulent flow in a curved channel. The bulk Reynolds number is \(\textrm{Re}_b=\rho _b u_b \delta / \mu _0 = 2600\) and the bulk Mach number is \(M_b=u_b/c_0=0.1\), where \(\mu _0\) and \(c_0\) are the viscosity and speed of sound at the wall temperature. The bulk density and velocity are defined as
where \(\mathcal {V}\) is the domain volume. The computation is carried out with \(N_{\theta }\times N_{r}\times N_{z}=216\times 72\times 144\), leading to grid spacings of \(\Delta z^+=5\), \(r_{\textrm{c}}^{+} \Delta \theta =10\), \(\Delta r^+=0.2\) at the wall. The mesh in the redial direction is staggered such that the wall coincides with an intermediate node, where the convective fluxes are identically zero. This approach guarantees correct telescoping of the numerical fluxes and exact conservation of the total mass, with the further benefit of doubling the time-step [20].
The surface curvature generates centrifugal instability, which breaks the symmetry in wall-normal direction. Specifically, the convex curvature of the inner wall has a stabilizing effect, leading to turbulence attenuation and smaller wall shear stress. Conversely, the concave curvature of the outer wall has a destabilizing effect, leading to turbulence amplification and larger wall shear stress. This trend is confirmed by the local values of friction Reynolds number, which resulted \(\textrm{Re}_{\tau ,i}=155\) at the convex (inner) wall and \(\textrm{Re}_{\tau ,{o}}=179\) at the concave (outer) wall. The local friction velocities are computed as
where \(\overline{u}_\theta\) is the mean streamwise velocity. A global friction velocity can be defined as follows. The mean pressure gradient in the curved channel is \((\textrm{d} P / \textrm{d} \theta ) / r\), which must be equal to \(-\rho u_{\tau }^2 / \delta\). From the streamwise momentum balance we obtain \(\textrm{d} P / \textrm{d} \theta = -\rho \left( r_i^2 u_{\tau ,i}^2+r_o^2 u_{\tau ,o}^2\right) /\left( 2 r_c \delta \right)\), hence
A comparison with the results of Moser & Moin and Brethouwer’s is shown in Table 1, where we report the local and global values of friction Reynolds number. The flow can be visualized in Fig. 5. The instantaneous field of velocity magnitude (5a) is depicted alongside two wall-parallel slices (5b), revealing a predominance of high-speed streaks near the concave wall and of low-speed streaks near the convex one.
One-point statistics are shown in Fig. 6, where we report the mean velocity profile in local wall units and turbulence intensities. The mean velocity in the log-law region is higher at the convex wall (black line in Fig. 6a) due to normalization with the local friction velocity, which is lower at the convex wall (i.e., lower wall shear stress). Velocity fluctuations, especially in the azimuthal direction (red line in Fig. 6b), highlight amplification of turbulence near the concave wall due to centrifugal instability. In general, we can appreciate the excellent agreement of the results.
One-point statistics for the turbulent cuved channel flow. Mean velocity profile in local wall units (a): black line refers to the convex wall, blue line to the concave one. Root mean square of velocity fluctuations (Instantaneous field of the velocity magnitude normalized by the bulk velocity ): red line refers to \(u_{\theta ,\text {rms}}\), blue line to \(u_{r,\text {rms}}\), black line to \(u_{z,\text {rms}}\)
3.3 Transonic Laminar Airfoil Flow
Laminar flows over an airfoil are commonly employed to assess accuracy, stability and convergence of solution algorithms for the Navier–Stokes equations. In recent years, these flows have also served as test cases for high-order numerical schemes, as in the study by Swanson [21], which we reference. Our investigation involves a two-dimensional laminar flow over a NACA 0012 airfoil at free-stream Mach number \(M_{\infty }=0.8\), Reynolds number \(\textrm{Re}=500\) (based on the chord length, c) and angle of attack \(\alpha =10^{\circ }\). The fore part of the domain has a radius of curvature 20c, and the outlet is positioned 20c downstream of the trailing edge. Four different grid densities are considered to analyze grid sensitivity. In all cases, structured meshes with C-type topology are generated using the open-source program Construct2D [22]. The finest mesh consists of \(N_x \times N_y = 1280\times 512\) grid points, with 768 cells on the airfoil clustered near the leading and trailing edges. The number of grid points was successively halved in both coordinate directions, to the coarsest mesh consisting of \(N_x \times N_y = 160\times 64\). Figure 7 shows the near-field region of the \(320\times 128\) grid.
No-slip adiabatic conditions are applied at the airfoil surface, while non-reflective boundary conditions are employed at the far field. The C-mesh closure generates a computational boundary in the wake region. Wake boundary conditions are handled by acting only on the ghost nodes, in which the solution of the corresponding internal nodes is set. This approach ensures that the boundary cutting through the domain in the wake region remains perfectly transparent to the flow.
Inflow conditions are such that the flow accelerates sufficiently to create a small supersonic region, without forming a shock wave. The supersonic zone of the flow is revealed in the Mach contours displayed in Fig. 8, where streamline patterns around the airfoil are also presented. We observe two recirculation regions with large extent both on the suction side of the airfoil and in the wake.
The comparison of the aerodynamic coefficients (Table 2), as well as the distribution of surface pressure and skin-friction (Fig. 9), reveals a very close agreement between our results and Swanson’s. These results also indicate that reasonably good estimates of the \(C_L\) and \(C_D\) can be achieved on fairly coarse grids.
3.4 Supersonic Turbulent Ramp Flow
A canonical configuration to study shock wave/boundary layer interactions is the compression ramp. Supersonic turbulent flows over a ramp has been extensively investigated through wind-tunnel experiments, whereas only a few DNS are available in the literature for this configuration.
One of the most relevant computational study is by Wu & Martin [23], who performed a DNS of supersonic turbulent flow over a \(24^{\circ }\) compression ramp. Our results are compared to their study, as well as to the experimental work of Bookey et al. [24] for the same configuration. A scale drawing of the computational domain is presented in Fig. 10. Dimensions are \(29\delta\) upstream and \(18\delta\) downstream of the corner in the streamwise direction, \(2.2\delta\) in the spanwise direction and \(6\delta\) in the wall-normal direction, where \(\delta\) is the \(99\%\) thickness of the incoming boundary layer. The mesh is generated analitically following the procedure used by Martin et al. [25], who provides full details of the mapping in the relevant work. The grid is clustered near the wall in the wall-normal direction, and a corner singularity is avoided by prescribing a finite curvature at the ramp corner. The DNS is performed with \(N_{x}\times N_y\times N_{z}=2432\times 256\times 160\) points, leading to grid spacings of \(\Delta x^+ \simeq 5\), \(\Delta z^+ \simeq 3\) and \(\Delta y^+ = 0.5\) at the wall. The incoming flow conditions are \(M_{\infty }=2.9\) and \(\textrm{Re}_{\delta }=31,500\). A recycling–rescaling method is applied to feed the inflow turbulence [26], with the recycling station located \(15\delta\) downstream of the inlet. Non-reflecting boundary conditions are used at the outlet and top boundary, while at the wall is enforced a no-slip isothermal condition with reflecting treatment based on characteristic wave relations.
The approaching turbulent boundary layer is deflected at the corner, where the shock system originates. Turbulence intensity amplifies as the flow crosses the shock, resulting in large density and pressure fluctuations. This interaction is depicted in Fig. 11, where we present an instantaneous field of numerical Schlieren. The coherent flow structures, visible in Fig. 11, indicate the presence of an extended region of separated flow.
Flow separation and reattachment can be precisely located as the points at which the mean friction coefficient changes sign. In our DNS, they are respectively predicted at \(x_s=-3.2\delta\) and \(x_r=1.3\delta\) (the ramp corner is located at \(x=0\)). A comparison with the results of Wu & Martin and Bookey et al. is shown in Table 3. To asses the accuracy of our simulation, we compare the wall pressure distribution and mean velocity profile across the interaction. Mean wall-pressure distribution (Fig. 12a) is predicted within the experimental uncertainty of \(5\%\) [24]. Two velocity profiles are shown in Fig. 12b, one in the incoming boundary layer and the other \(4\delta\) downstream of the corner. In both cases, the agreement is within \(5\%\).
4 Conclusions
We presented FLEW, a high-fidelity solver for DNS of compressible flows over complex geometries. The code is HPC-oriented thanks to MPI parallelization and the ability to run on multi-GPU architectures. Various aspects of extending the energy-preserving schemes proposed by Pirozzoli [13] to simulate flows over curvilinear meshes are presented. In Sect. 2, we describe a locally conservative approximation of convective derivatives suitable for generalized curvilinear coordinates. We also discuss how shock waves, viscous terms and characteristic boundary conditions are treated. Special attention is devoted to the metric terms, which are computed with the same scheme used for convective derivatives, ensuring the freestream preservation.
Demonstrative examples of the efficacy of the implemented algorithms are provided in Sect. 3. We started our assessment from the time-reversibility benchmark of the inviscid Taylor–Green vortex, aiming to to test the energy-preserving capability of our solver on curvilinear meshes. The results demonstrate that the solver can perfectly recover the initial conditions, indicating the absence of any artificial dissipation. Various viscous test cases, including the curved channel flow, airfoil flow and ramp flow, further substantiate the accuracy and robustness of our solver across a wide range of flow regimes, spanning from low-subsonic, through transonic to supersonic conditions. Comparison of the results with experimental and numerical data available in the literature revealed excellent agreement for all tested configurations.
While employing the generalized curvilinear coordinates in two directions already allows us to address numerous significant problems, our next goal is to remove the constraint that the third axis be normal to the surface. This extension will enable the code to handle fully arbitrary geometries.
Data Availability
All data that support the findings of this study are available from the corresponding author upon request.
References
Bernardini, M., Modesti, D., Salvadore, F., Pirozzoli, S.: Streams: a high-fidelity accelerated solver for direct numerical simulation of compressible turbulent flows. Comput. Phys. Commun. 263, 107906 (2021)
Johnsen, E., Larsson, J., Bhagatwala, A.V., Cabot, W.H., Moin, P., Olson, B.J., Rawat, P.S., Shankar, S.K., Sjögreen, B., Yee, H.C.: Assessment of high-resolution methods for numerical simulations of compressible turbulence with shock waves. J. Comput. Phys. 229(4), 1213–1237 (2010)
Piquet, A., Zebiri, B., Hadjadj, A., Safdari Shadloo, M.: A parallel high-order compressible flows solver with domain decomposition method in the generalized curvilinear coordinates system. Int. J. Numer. Methods Heat Fluid Flow 30(1), 2–38 (2020)
Bernardini, M., Modesti, D., Pirozzoli, S.: On the suitability of the immersed boundary method for the simulation of high-Reynolds-number separated turbulent flows. Comput. Fluids 130, 84–93 (2016)
Johnson, J.P., Iaccarino, G., Chen, K.-H., Khalighi, B.: Simulations of high Reynolds number air flow over the NACA-0012 airfoil using the immersed boundary method. J. Fluids Eng. 136(4), 040901 (2014)
Coppola, G., Capuano, F., Pirozzoli, S., Luca, L.: Numerically stable formulations of convective terms for turbulent compressible flows. J. Comput. Phys. 382, 86–104 (2019)
Pirozzoli, S.: Generalized conservative approximations of split convective derivative operators. J. Comput. Phys. 229(19), 7180–7190 (2010)
Kawai, S., Lele, S.K.: Localized artificial diffusivity scheme for discontinuity capturing on curvilinear meshes. J. Comput. Phys. 227(22), 9498–9526 (2008)
Visbal, M.R., Gaitonde, D.V.: On the use of higher-order finite-difference schemes on curvilinear and deforming meshes. J. Comput. Phys. 181(1), 155–185 (2002)
Sun, Z.-S., Ren, Y.-X., Zhang, S.-Y., Yang, Y.-C.: High-resolution finite difference schemes using curvilinear coordinate grids for DNS of compressible turbulent flow over wavy walls. Comput. Fluids 45(1), 84–91 (2011)
Kuya, Y., Kawai, S.: High-order accurate kinetic-energy and entropy preserving (keep) schemes on curvilinear grids. J. Comput. Phys. 442, 110482 (2021)
Pirozzoli, S.: Numerical methods for high-speed flows. Annu. Rev. Fluid Mech. 43, 163–194 (2011)
Pirozzoli, S.: Stabilized non-dissipative approximations of Euler equations in generalized curvilinear coordinates. J. Comput. Phys. 230(8), 2997–3014 (2011)
Ducros, F., Ferrand, V., Nicoud, F., Weber, C., Darracq, D., Gacherieu, C., Poinsot, T.: Large-eddy simulation of the shock/turbulence interaction. J. Comput. Phys. 152(2), 517–549 (1999)
Poinsot, T.J., Lele, S.: Boundary conditions for direct simulations of compressible viscous flows. J. Comput. Phys. 101(1), 104–129 (1992)
Spalart, P.R., Moser, R.D., Rogers, M.M.: Spectral methods for the Navier–Stokes equations with one infinite and two periodic directions. J. Comput. Phys. 96(2), 297–324 (1991)
Duponcheel, M., Orlandi, P., Winckelmans, G.: Time-reversibility of the Euler equations as a benchmark for energy conserving schemes. J. Comput. Phys. 227(19), 8736–8752 (2008)
Moser, R.D., Moin, P.: The effects of curvature in wall-bounded turbulent flows. J. Fluid Mech. 175, 479–510 (1987)
Brethouwer, G.: Turbulent flow in curved channels. J. Fluid Mech. 931, 21 (2022)
Modesti, D., Pirozzoli, S.: Reynolds and Mach number effects in compressible turbulent channel flow. Int. J. Heat Fluid Flow 59, 33–49 (2016)
Swanson, R., Langer, S.: Steady-state laminar flow solutions for NACA 0012 airfoil. Comput. Fluids 126, 102–128 (2016)
Prosser, D.: Construct2D. https://sourceforge.net/projects/construct2d/
Wu, M., Martin, M.P.: Direct numerical simulation of supersonic turbulent boundary layer over a compression ramp. AIAA J. 45(4), 879–889 (2007)
Bookey, P., Wyckham, C., Smits, A., Martin, P.: New experimental data of STBLI at DNS/LES accessible Reynolds numbers. In: 43rd AIAA Aerospace Sciences Meeting and Exhibit, p. 309 (2005)
Martin, P., Xu, S., Wu, M.: Preliminary work on DNS and LES of STBLI. In: 33 rd AIAA Fluid Dynamics Conference and Exhibit (2003)
Ceci, A., Palumbo, A., Larsson, J., Pirozzoli, S.: Numerical tripping of high-speed turbulent boundary layers. Theor. Comput. Fluid Dyn. 36(6), 865–886 (2022)
Acknowledgements
Computational resources kindly made available by CINECA are thankfully acknowledged. G.S. acknowledges the financial support from ICSC—Centro Nazionale di Ricerca in ‘High Performance Computing, Big Data and Quantum Computing’, funded by European Union—NextGenerationEU. A.C. acknowledges the financial support from TEAMAero Horizon 2020 research and innovation programme under grant agreement no. 860909. S.P. acknowledges the financial support from the European High Performance Computing Joint Undertaking and Germany, Italy, Slovenia, Spain, Sweden and France under grant agreement no. 101092621.
Funding
Open access funding provided by Università degli Studi di Roma La Sapienza within the CRUI-CARE Agreement.
Author information
Authors and Affiliations
Contributions
GS contributed to the code development, wrote the main manuscript text, and supervised the simulations of the curved channel flow and airfoil flow. AC contributed to the code development and supervised the simulations of the Taylor–Green vortex and compression ramp flow. SP conceived the code, wrote its core, and supervised its development and validation throughout all simulations. All authors reviewed the manuscript.
Corresponding author
Ethics declarations
Conflict of interest
The authors declare that they have no conflict of interest.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Soldati, G., Ceci, A. & Pirozzoli, S. FLEW: A DNS Solver for Compressible Flows in Generalized Curvilinear Coordinates. Aerotec. Missili Spaz. 103, 413–425 (2024). https://doi.org/10.1007/s42496-024-00199-4
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s42496-024-00199-4