A Low Dissipative and Stable Cell-Centered Finite Volume Method with the Simultaneous Approximation Term for Compressible Turbulent Flows

A simultaneous-approximation term is a non-reflecting boundary condition that is usually accompanied by summation-by-parts schemes for provable time stability. While a high-order convective flux based on reconstruction is often employed in a finite-volume method for compressible turbulent flow, finite-volume methods with the summation-by-parts property involve either equally weighted averaging or the second-order central flux for convective fluxes. In the present study, a cellcentered finite-volume method for compressible Naiver–Stokes equations was developed by combining a simultaneous-approximation term based on extrapolation and a low-dissipative discretization method without the summation-by-parts property. Direct numerical simulations and a large eddy simulation show that the resultant combination leads to comparable non-reflecting performance to that of the summation-by-parts scheme combined with the simultaneous-approximation term reported in the literature. Furthermore, a characteristic boundary condition was implemented for the present method, and its performance was compared with that of the simultaneous-approximation term for a direct numerical simulation and a large eddy simulation to show that the simultaneousapproximation term better maintained the average target pressure at the compressible flow outlet, which is useful for turbomachinery and aerodynamic applications, while the characteristic boundary condition better preserved the flow field near the outlet.


Introduction
Compressible turbulent flows are encountered in many engineering applications such as airplanes, high-speed vehicles, and turbomachinery. High-fidelity simulations of turbulence such as direct numerical simulations (DNS) or large eddy simulations (LES) require minimal numerical dissipation as upwind-biased schemes artificially dampen small-scale turbulent structures [1]. The adoption of low-dissipative schemes, however, leads to requiring non-reflecting boundary conditions at the subsonic compressible outlet, as the lack of numerical dissipation leads to the accumulation of incoming acoustic waves affecting the internal flow field unless exact outlet boundary conditions are imposed [2].
A simultaneous-approximation term (SAT) is a non-reflecting boundary condition that is often used with summation-by-parts (SBP) schemes, as the combination of an SBP scheme and an SAT (SBP-SAT approach) leads to provable time stability, meaning that the resultant numerical method is stable for long-time numerical integration [3][4][5][6].
The SBP-SAT approach was extensively studied for high-order finite-difference methods (FDMs) [5,6], but only a few studies regarding SBP-SAT approaches for finite-volume methods (FVMs) are available, while most research for FVMs concerns node-based FVMs. Nordström et al. [7] showed that the first derivative for unstructured node-centered FVM is on an SBP form given that equally weighted averaging is used for convective flux

Unfiltered Navier-Stokes Equations
Navier-Stokes equations for compressible flows are written as follows [14]: where t and x i are independent variables representing time and spatial coordinates. The summation convention over repeated indices is used. ρ, u i , p, E, σ ij , and q j denote density, velocity vector, pressure, total energy per unit mass, shear-stress tensor, and heat flux, respectively. E, σ ij , and q j are given as follows: where the heat-capacity ratio is given as γ = C p /C v . C p and C v are specific heat at constant pressure and constant volume, respectively. The rate of strain tensor S ij is given as Assuming the ideal gas model for air, the equation of state is given as where the gas constant is R = C p − C v . Variation of µ with T is approximated by a power law µ(T) = µ(T 0 )(T/T 0 ) 0.76 , which is valid between 150 and 500 K. κ is given as κ = µC p /Pr, where Pr is the Prandtl number.
Navier-Stokes equations are solved for DNS, which does not need a turbulence model given that flow structures are resolved by the Kolmogorov scale [15].

Spatially Filtered Navier-Stokes Equations
While DNS is prohibitively expensive for high Reynolds number flows encountered in practical applications, LES is an alternative high-fidelity simulation technique where large-scale turbulent eddies are directly solved, and small-scale eddies are modeled by a subgrid-scale (SGS) model.
A low-pass filter operation is defined as where x and ξ are spatial coordinate vectors in the flow domain of Ω. Filter kernel G satisfies the normalization condition and depends on the filter width defined as ∆. A Favre filter is also defined as φ ≡ ρφ/ρ. For the definition of filter width, the definition of an element size in three-dimensional space [16] is as follows: where ∆ i denotes the filter width in the i direction. Low-pass-filtered Navier-Stokes equations are given as follows [17]: where ρu i u i u j − ρu i u i u j ≈ 2τ ij u i is assumed following the work by Okong'o and Bellan [17]. E, σ ij , q j , S ij , τ ij , and ζ j are given as follows: Last, it is assumed that σ ij = σ ij , σ ij u i = σ ij u i , and q j = q j following the work by Okong'o and Bellan [17].

Subgrid-Scale Modeling
For LES, τ ij and ζ j are modeled by using eddy viscosity closures as follows: where τ kk and Π are modeled using Yoshizawa parameterization [18] and the Vreman kernel [19], given as follows: where B β , β ij , and α ij are defined as follows:

Numerical Methods, Verification, and Validation
The present numerical implementation uses Equations (1)-(3) for DNS, and Equations (11)-(13) for LES. As identical numerical methods were employed for both DNS and LES, descriptions of the numerical schemes are given for DNS.
Equations (1)-(3) are solved in cell-centered finite-volume form: where U = [ρ, ρu 1 , ρu 2 , ρu 3 , ρE] is the vector of conserved variables averaged in the control volume. F e and F v are the face-normal Euler and viscous flux vectors, respectively, corresponding to the right-hand sides of Equations (1)- (3). V cv is the volume of the control volume, ∑ f is the summation over the cell faces per cell, and A f is the area of the cell face. The volume and the surface integrals are approximated with the second-order midpoint quadrature rule. Following Khalighi et al. [10], biased polynomial reconstruction, which is maximal third-order accurate to the face centroid, is performed for primitive variables (ρ, u i , and p) on both sides of each cell face. Using the face-reconstructed values on the left-and the right-hand sides of a cell face, φ l and φ r , a central flux is calculated in skew-symmetric form [8], while dissipative essentially non-oscillatory (ENO) flux [20] is calculated by using the HLLC Riemman solver [13]. For more details regarding face reconstruction for fluxes, refer to Khalighi et al. [10].
A blended flux is expressed as follows: where 0 ≤ α ≤ 1 is a blending parameter. Unlike the implicit LES where a global constant of α = 0.25 [21] or α = 0.1 [22,23] is used, Khalighi et al. [10] computed α using the row norm of the symmetric part of the differencing operator for Euler fluxes D, as follows: where c = 2 was chosen from numerical tests [10]. The motivation to use Equation (29) was that a skew-symmetric operator (D = −D T ) is stable and non-dissipative (energyconserving) in the inviscid limit on a periodic domain [10]. The blending scheme by Khalighi et al. [10] is cost-effective, as α is precomputed using grid information prior to the main calculation while the use of high-order spatial filter requires filtering at each time step in addition to the calculation of fluxes. A relative solution switch [10] is employed to activate dissipative ENO fluxes in regions with shocks by setting α = 1. A third-order low-storage Runge-Kutta scheme [24] is employed to advance the equations in time.

Two-Dimensional Propagating Inviscid Vortex
The accuracy of the present numerical scheme was verified by using an exact vortex solution to the Euler equations, which is given as follows [8]: where c ∞ and ρ ∞ are the reference speed of sound and reference density, respectively. A vortex of radius one and Γ = 5.65 is initially located at (x 0 = −1, y 0 = −1) and is diagonally convected to (x f = 1, y f = 1) with the speed of M ∞ = 0.5. A 10 × 10 box is discretized on 50 × 50, 100 × 100, 200 × 200, and 400 × 400 cell uniform Cartesian grids. The exact solution is used for the initial condition, and periodic boundary conditions are imposed. The constant Courant-Friedrichs-Lewy (CFL) number of 0.3 is set for all calculations, such that both spatial and temporal convergence rates are measured [8].
L 2 and L ∞ norm error norms at (x f = 1, y f = 1) are shown against the grid width ∆ in Figure 1. L 2 and L ∞ norm errors both exhibit second-order convergence rates.

LES of Supersonic Turbulent Channel Flow
To assess the predictive capability of the present method for compressible near-wall turbulence, LES is performed for supersonic turbulent channel flow at M b = 1.5 and . c w is the speed of sound at the wall and · denotes an average over time, and streamwise and spanwise directions. Periodic boundary conditions are imposed in the streamwise and the spanwise directions, while no-slip isothermal wall boundary conditions are imposed at y = δ and y = −δ, where δ is the channel half width. The spatially uniform body force is imposed in the streamwise direction to drive channel flow, and is adjusted each time-step to conserve the global mass rate.
To close the SGS models given in Equations (20) and (21), a C of 0.07 and Pr t of 0.5 obtained from an a priori study of Moin et al. [26] are used. C I in Equation (22) is assumed to be negligible due to the DNS study of supersonic channel flow [27], where the contribution of the dilatation to the turbulent kinetic energy dissipation rate was found to be insignificant.
Profiles of the mean and the root mean square (RMS) of velocity and temperature obtained from the present LES are shown against DNS data by Coleman et al. [27] in Figures 2 and 3, respectively. While the overprediction of the mean velocity in the core region and the peak of the RMS near the wall are well-known deficiencies of the eddyviscosity SGS models [28], both the mean and the RMS of temperature showed favorable agreement with those of the DNS by Coleman et al. [27]. All the results show more accurate results on the fine grid (80 × 64 × 64) compared with those on the coarse grid (40 × 48 × 32).
Next, consider a uniform grid with N nodes indexed from 1 to N for 0 ≤ x ≤ 1. For P is a diagonal positive definite matrix that defines a norm such that ||φ|| 2 P ≡ φ T Pφ. Q is a skew-symmetric matrix except at boundaries such that Q + Q T = diag(−1, 0, ..., 0, 1). Then, D x is called an SBP operator.

NSBP-SAT Approach for a Cell-Centered Finite-Volume Method
Usage of the SAT with non-SBP schemes (NSBP-SAT approach) was not reported in the literature, to the best of the authors' knowledge, presumably because SAT is coupled with the SBP for provable time stability, and partially due to successful applications of the CBC to compressible Navier-Stokes equations.
While both FDMs and nodal FVMs have data nodes at boundaries where the SAT can be imposed, the cell-centered FVM does not have data nodes where variables are advanced in time. Regarding this issue, Nordström and Björck [9] introduced data nodes at boundaries to prove time stability for the cell-centered FVM.
In the present work, the SAT is imposed on the reconstructed variables at boundaries using extrapolation. It is subsequently shown for numerical tests that the NSBP-SAT approach leads to comparable non-reflecting performance at outlets to that of the SBP-SAT approach by Shoeybi et al. [8].
While cell-volume-averaged variables are stored at cell centers for the cell-centered FVM, variable distribution within a cell can be reconstructed using gradients. For example, for the scalar variables φ L , φ P , and φ B shown in Figure 4, distributions of φ within cells can be reconstructed piecewise constantly (Figure 4a) or piecewise polynomially (Figure 4b). If an expression of the SAT at the boundary face b is given as 2σ , it is equivalent to assume a piecewise constant profile of φ for the cell B and impose an SAT at the boundary face b. For a more accurate representation of variables at the boundary, the piecewise polynomial reconstruction used for Equation (28) was employed. As variables outside boundaries are not available, linear extrapolations from the variables interior to the boundaries are used to calculate gradients required for the reconstruction. As a result, an SAT given as 2σ ] is used to drive the face-reconstructed variable φ b at boundary b to boundary data given as g(t).

Two-Dimensional Acoustic Pulse
To test the non-reflecting capability of the NSBP-SAT approach, a two-dimensional acoustic pulse [29] was simulated using the present cell-centered FVM. A 10 × 10-sized two-dimensional box was uniformly discretized using 64 cells in each direction. The initial field is given as follows: where r = x 2 + y 2 , and origin (x, y) = (0, 0) is located at the center of the box for Cartesian coordinates. The amplitude of a Gaussian isentropic pulse A was set to 0.05, and α = 1.023 was used. While the SAT is applied to all boundaries, reference values are used for the boundary data required for the SAT (ρ b = ρ ∞ , u b = v b = 0, and p b = (ρ ∞ c 2 ∞ )/γ). The simulation was run until t = 10/c ∞ , when the acoustic wave entirely left the computational domain. Figure 5 shows instantaneous contours of p /p ∞ at t = 5/c ∞ and t = 10/c ∞ , while p /p ∞ = (p − p ∞ )/p ∞ . At t = 5/c ∞ , the acoustic pulse exited the domain without a noticeable sign of reflection. At t = 5/c ∞ , the acoustic pulse exited the domain and the remaining perturbations were observed. While Shoeybi et al. [8] reported the maximal value of p /p ∞ to be about 1% of the initial amplitude of the acoustic pulse at t = 10/c ∞ , the NSBP-SAT approach resulted in 1.6% of the initial amplitude of the acoustic pulse. Therefore, the non-reflecting performance of the NSBP-SAT approach was observed to be comparable to that of the SBP-SAT approach by Shoeybi et al. [8] for a two-dimensional acoustic pulse.

Low Reynolds Number Flow around a Two-Dimensional Circular Cylinder
To further test the performance of the NSBP-SAT approach, unsteady compressible laminar flow around a two-dimensional circular cylinder was simulated. A Reynolds number Re D = U ∞ D/ν of 100 was employed, where U ∞ is freestream velocity and D is the diameter of the cylinder.
Two grids of 260 × 150 and 400 × 200 cells in the circumferential and the radial directions were employed for the coarse-and fine-grid cases, respectively, and the coarse grid is shown in Figure 6. The computational domain was a circle of radius 30D, and the first grid spacing in the radial direction was 3 × 10 −3 D [8]. The no-slip isothermal wall boundary condition was imposed around the cylinder, and the SAT with the freestream boundary data was imposed at the far field.    Table 1 shows the time-averaged drag coefficient C D , maximal and minimal lift coefficients C L , and Strouhal number St for the present simulations with the coarse and the fine grids; simulations by Svärd and Nordström [30], and Shoeybi et al. [8]; and experimental data obtained by Fey et al. [31]. Svärd and Nordström [30] employed the 5th-order accurate scheme with a fine grid of 401 × 201, and Shoeybi et al. [8] employed a 2nd-order accurate scheme with 260 grids in the circumferential direction. Both Shoeybi et al. [8] and Svärd and Nordström [30] employed the SBP-SAT approach. As results from the coarse grid were as accurate as those reported by Shoeybi et al. [8], and results from the fine grid converged to results by Svärd and Nordström [30], the NSBP-SAT approach performed as well as the SBP-SAT approaches did in terms of non-reflecting performance at the outlet.  [31] 0.165 ± 0.001

LES of Flow over a Circular Cylinder at M = 0.4
To assess the non-reflecting performance of the NSBP-SAT approach for compressible turbulent flow, LES of flow over a circular cylinder at M ∞ = 0.4 and Re D = 3900 was performed. A circle of radius 35D with spanwise length of πD was used for the computational domain where the cylinder of a diameter of 1D was located at the center of the radial and circumferential plane. A sponge layer [32] with thickness of 15D was applied at the outer boundary to damp out vortices [33]. A grid of 288 × 400 × 64 cells in the radial, circumferential, and spanwise directions, respectively, was used. Grid resolution was similar to that of Mani et al. [33] except for the finer grid resolution in the spanwise direction, considering that Mani et al. [33] used a 6th-order compact finite-difference code [34].
For SGS models, C of 0.07 and Pr t of 0.9 were used for Equations (20) and (21), while C I in Equation (22) was assumed to be zero due to the low Mach number [35]. Figure 8 shows the contour of the instantaneous spanwise vorticity to indicate the separating laminar shear layers. While coarse spanwise resolution is known to result in transition to turbulence in the separating shear layers closer to the cylinder [36], sufficient grid resolution in the spanwise direction for the present simulation could be confirmed from the elongated separating laminar shear layers.  Table 2 shows a comparison of C D and St calculated from the present LES with LES results by Mani et al. [33]. Although the present simulation used a similar number of grids, except for the finer spanwise resolution, to those of Mani et al. [33], C D and St from the present simulation were quite close to those of Mani et al. [33] despite the lower order of accuracy for the present numerical method.

Average C D St
Present results 1.18 0.199 Mani et al. [33] 1.17 ± 0.01 0.200 ± 0.002 Figure 9 shows density spectra in the wake for both LES results. While good agreement of low-frequency components, including the location of the peak, explains the accurate prediction of the Stouhal number in Table 2, slight overpredictions at high frequencies are attributed to the low order of accuracy of the present solver. Note that the present solver predicted high-frequency components that are qualitatively similar to those predicted by the high-order low-dissipative solver [33]. As Mittal and Moin [1] showed that even 5th-order accurate upwind-biased schemes lead to significant underprediction for high-wavenumber components of the energy spectra compared with the second-order central difference, effects of the SAT are found to be limited near the boundary; therefore, the present numerical method remained low-dissipative for the main flow feature of interest. For reference, temporal and spatial scales for small turbulent eddies in large eddy simulation are implicitly associated [37].

Description of the Characteristic Boundary Condition
The characteristic boundary condition (CBC) is arguably the most well-known nonreflecting boundary condition for compressible Navier-Stokes equations. By applying one-dimensional characteristic analysis by Thompson [38] to the inviscid parts of the compressible Navier-Stokes equations, Poinsot and Lele [2] obtained ∂ρ ∂t ∂ ∂t where L i are the amplitudes of characteristic waves that are defined as follows: where λ i are the characteristic velocities associated with L i that are defined as follows: Because the characteristic velocity for Equation (49) is negative for subsonic flows, while upwind differencing requires data outside the computational domain, numerical calculations of spatial derivatives in Equation (49) involve downwind differencing, which is unstable [2]. To avoid this problem, Poinsot and Lele [2] modeled L 1 as follows: where σ is a tunable parameter, c is the speed of sound, M is the reference Mach number, and p t is the target pressure at the outlet. While the CBC by Poinsot and Lele [2] was limited to one-dimensional outward flows, Yoo et al. [39] proposed adding a transverse term to Equation (57) to account for the multidimensionality of outward flows (i.e., vortices): where β is another tunable parameter, and T 1 is defined as follows:

Characteristic Boundary Condition for a Cell-Centered FVM
Regarding the CBC for a cell-centered FVM, Gross and Fasel [40], and Motheau et al. [41] used spatial gradients modified by characteristic analysis to impose boundary conditions using ghost cells, which is referred to as a gradient-based CBC (GBCBC). For instance, an expression for a density gradient was given by Poinsot and Lele [2]: where L 1 is modified using Equation (58). Using such modified gradients of variables required at outlets, one can calculate all variables required at ghost cells. Alternatively, one can extrapolate variables to the boundary, and directly integrate Equations (44)-(48) at the boundary, which is referred to as the extrapolated original CBC (EOCBC) [42]. Kang and You [42] showed that the EOCBC performed better than the GBCBC for a convecting two-dimensional vortex and LES of a subsonic turbine cascade when the EOCBC and the GBCBC are implemented in a cell-centered FVM. Therefore, the EOCBC was implemented in the present study.
For the rest of the Section 5, the SAT refers to the SAT in the NSBP-SAT approach.

Convecting Two-Dimensional Inviscid Vortex
A convecting two-dimensional inviscid vortex is a well-known benchmark test for non-reflecting boundary conditions [2,[39][40][41][43][44][45]. A single vortex is superimposed on a uniform flow aligned with the x 1 direction. The initial velocity field is given as follows [45]: where Γ is the strength of the vortex, and (x 1c , x 2c ) is the coordinate for the center of a computational domain. The computational domain is a square of dimension L of 0.013 that is discretized by 80 × 80 uniform cells. Initially, uniform pressure and temperature of 101,300 Pa and 300 K, respectively, were assumed. While the SAT was imposed at boundaries in the streamwise direction, periodic boundary conditions were imposed for the rest of the boundaries. While a vortex case at a low Mach number (close to 0) presents a greater challenge than that of a case at a high Mach number (close to 1) for the non-reflecting outlet [43,44], a case with M ∞ = u 1 /c ∞ = 0.028 (U 0 = 10 (m/s)) and Γ = 0.011 [45] was chosen for the present test. Figure 10 shows the exact solution, and obtained results using the CBC and the SAT from the top to the bottom row, respectively. For each case, streamwise velocity isocontours with normalized pressure field p * at t * = 0.45, 1.12, and 1.43 are shown, where t * = t/τ, τ = L/(2U 0 ), and p * is given as follows: Compared with the reference solution, the CBC allowed for the vortex to pass the outlet with small distortion, which was confirmed by the similar minimal and maximal pressure values at t * = 1.12 and 1.43, respectively, along with vortex structures shown by the streamwise velocity isocontours.
The SAT, on the other hand, seemed to destroy the vortex at the outlet, along with significant deviations of the streamwise velocity and the pressure fields from those of the reference solution. However, at t * = 1.12 and 1.43, spatially averaged pressure p * was maintained at 0.0, which is the target pressure imposed by the SAT. Results obtained using the CBC showed slightly larger deviations of p * = 0.03 from the reference value.
These observations can be understood from the difference between the CBC and the SAT. While the CBC utilizes linear characteristic analysis to model the local incoming characteristic wave, the SAT utilizes a linear characteristic analysis to impose a given boundary dataset, which is usually given to be either reference values or averaged values at the outlet [4]. In other words, while the CBC attempts to preserve flow structures at outlets by modeling the incoming wave, the SAT ignores such details and enforces variables at outlets to become closer to what is prescribed as boundary values. The observation that the SAT slightly better keeps the spatially averaged pressure field close to the target pressure is encouraging as the accurate imposition of the mean pressure at the outlet is a particular interest for turbomachinery and aerodynamic applications. For example, while the performance of a compressor is often investigated through a thermodynamic-efficiency versus mass-flow rate graph for a multiple of operation conditions [46] obtained from a Reynolds-averaged Navier-Stokes technique [47], the mean mass-flow rate is determined by the mean pressure ratio at the inlet and the outlet [46]. Because LES is often performed at a specific operation condition due to computational cost [47], keeping the mean pressure at the outlet as close to a prescribed value as possible is desirable to simulate a certain operation condition. Time-averaged drag and lift coefficients are also often the primary interests for aerodynamic applications such as an airfoil [48] where the accurate imposition of the mean pressure is closely related to the accurate calculation of the drag and the lift.
The CBC, on the other hand, produces a more accurate flow field near the boundary along with more accurate pressure fluctuations from Figure 10. This is beneficial for flow simulations in which an accurate description of unsteady flow physics is more important than mean quantities such as the drag and the lift.
Another interesting aspect for the NSBP-SAT approach is that, unlike the SBP-SAT approach, which guarantees linear time stability, the NSBP-SAT approach does not allow for such provable stability, but it was found to be stable for the vortex test.
Last, the dissipative characteristic of the SAT observed in Figure 10 does not necessarily mean that it renders the present numerical method dissipative for the entire flow field if the outlet is sufficiently far from the flow feature of interest. For instance, the density spectrum shown in Figure 9 shows that the internal field was not noticeably affected by the SAT, as high-frequency components are qualitatively well-captured.

LES of Flow through the VKI Subsonic Turbine Cascade
To further test the NSBP-SAT approach, subsonic flow over a VKI turbine cascade at a Reynolds number of 2.8 × 10 6 based on chord length was simulated by employing the SAT at the outlet. The LES result with the SAT was compared with the LES result obtained using the CBC and the experimental data by [49]. The VKI turbine cascade was also used by Granet et al. [45] to assess the non-reflecting capabilities of several variants of CBCs.
A blade with a chord length C of 140 mm was located from the inlet and the outlet by 0.5 and 2C, respectively. The three-dimensional domain was discretized by 8.16 × 10 6 hexahedral cells in total, and the mesh is shown in Figure 11. The thickness of the domain in the spanwise direction was 0.1C, for which 64 cells were used. While synthetic turbulence [50] is added to the uniform incoming flow at the Mach number of 0.19, uniform thermodynamic variables were assumed at the inlet. Either the CBC or the SAT was imposed at the outlet. Periodic conditions were applied to the rest of the boundaries. Statistics were collected for 10C/U b , where U b is bulk streamwise velocity. For further details regarding the geometry of the cascade, refer to Sieverding et al. [49]. The isentropic Mach number is defined as follows: where p 0 is the total pressure given at the inlet and p w denotes pressure along the wall of the cascade. For SGS models, C of 0.07 and Pr t of 0.9 were used for Equations (20) and (21), while C I in Equation (22) was assumed to be zero due to a low Mach number [35].
The isentropic Mach number distribution along the cascade, which is averaged in time, and the spanwise direction for simulations using the CBC and the SAT are shown in Figure 12. While both the SAT and the CBC produced reasonably accurate results, the SAT performed as well as the CBC did as a non-reflecting boundary condition. At the peak of the isentropic Mach number, the SAT result even showed slightly improved performance, which may have been due to the capability of the SAT to keep the time and spatially averaged pressure field closer to the target pressure, which was also observed in the two-dimensional inviscid vortex case.  [49]; red squares, obtained results using CBC; green circles, obtained results using SAT.

Concluding Remarks
A cell-centered finite-volume method for compressible turbulent flow was developed by combining a non-SBP low-dissipative discretization method and an SAT based on extrapolation, leading to an NSBP-SAT approach. Although the NSBP-SAT approach does not lead to provable time stability, it showed comparable non-reflecting performance to that of the SBP-SAT approach reported in the literature for a two-dimensional acoustic pulse, flow around a two-dimensional circular cylinder at a low Reynolds number, and LES of compressible flow over a circular cylinder at a Mach number of 0.4. The CBC was also implemented in the present finite-volume method and compared with the SAT in the NSBP-SAT approach against a convecting two-dimensional inviscid vortex and a subsonic turbine cascade in terms of non-reflecting performance. The SAT maintained the mean pressure field near the outlet slightly closer to the prescribed value than the CBC did, rendering the SAT slightly more beneficial for turbomachinery where accurate mean pressure imposition at the outlet is important for simulating a certain operation condition and airfoils where time-averaged drag and lift are the primary interests. The CBC, on the other hand, better preserved flow structures, leading to more accurate predictions of pressure fluctuations near the outlet. This is attractive for flow simulations where accurate descriptions of unsteady flow features are more important than accurate predictions of mean quantities are, such as the drag and the lift.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflicts of interest.

Abbreviations
The