Code Verification for the SENSEI CFD Code

Xue, Weicheng; Wang, Hongyu; Roy, Christopher J.

doi:10.1115/1.4062609

Abstract

This work performs systematic studies for code verification for turbulence modeling in our research CFD code SENSEI. Turbulence modeling verification cases including cross term sinusoidal manufactured solutions and a few exact solutions are used to justify the proper Spalart–Allmaras and Menter’s SST turbulence modeling implementation of the SENSEI CFD code. The observed order of accuracy matches fairly well with the formal order for both the 2D/3D steady-state and 2D unsteady flows when using the cross term sinusoidal manufactured solutions. This work concludes that it is important to keep the spatial discretization error in a similar order of magnitude as the temporal error in order to avoid erroneous analysis when performing combined spatial and temporal order analysis. Since explicit time marching scheme typically requires smaller time-step size compared to implicit time marching schemes due to stability constraints, multiple implicit schemes such as the singly diagonally implicit Runge–Kutta multistage scheme and three point backward scheme are used in our work to mitigate the stability constraints.

Introduction

Computational fluid dynamics (CFD) is a method to mathematically model and numerically solve physical problems involving fluid flows. There are various types of flows, including laminar flows and turbulent flows. Laminar flow transitions to turbulent flow when the Reynolds number becomes large enough. Miscellaneous methods have already been developed including Reynolds-averaged Navier–Stokes (RANS) models, large Eddy simulation and direct numerical simulation to model turbulent flows. This work will only focus on the RANS method. RANS models can be categorized as zero equation models (Prandtl’s mixing length model for an example), one equation models (e.g., Spalart–Allmaras [1]), two equation models (⁠ $k - ϵ$ [2], $k - ω$ [3], $k - ω$ SST [4], etc.), etc. The Boussinesq approximation is usually introduced to close the RANS equations by ensemble averaging over the Navier–Stokes equations since there is no analytic expression for the Reynolds stresses. The Boussinesq approximation relates the Reynolds stresses with mean flow velocities through the use of eddy viscosity. The eddy viscosity is solved differently using different RANS models. For example, the Spalart–Allmaras (SA) model solves the eddy viscosity through a working variable, while the $k - ω$ SST model calculates the eddy viscosity by solving the kinetic energy and turbulent dissipation rate equations. The SA model is relatively cheap in terms of computational cost, and it is often used in flows with just slight separation such as transonic flow past airfoils, etc. The $k - ϵ$ model which is a two equation model performs well for free shear flows, and the $k - ω$ model is a good option for wall-bounded flows. The $k - ω$ SST model, also known as the Menter’s SST, is devised to combine the advantages of $k - ϵ$ (suitable for the free shear flow) and $k - ω$ model (suitable for the inner region of the boundary flow).

Code verification [5,6] is a set of procedures to ensure that the implementation of a code is numerically correct. Order of accuracy test is a more rigorous criterion than expert judgment, error quantification, and consistency/convergence [7]. The method of manufactured solutions or numerical benchmark solutions or exact solutions combined with order of accuracy verification is usually utilized to rigorously verify whether a code can achieve the expected order of accuracy or not. For steady-state flows, the spatial residuals are driven to machine zero and then the iteratively converged result is compared to the exact steady-state solution to obtain the discretization error. For unsteady flows, the iterative procedure is similar but the temporal discretization error also exists and should be accounted for, that is, the overall discretization error is comprised of the spatial and temporal discretization error. There are multiple sources of error affecting the code verification analysis [5], such as poor mesh quality, improper time-step size, large subiterative errors, round-off error accumulation, which is especially important for ill-conditioned systems and unsteady flows, and even implementation errors which can be easily overlooked.

Since not many exact solutions exist for the compressible Navier–Stokes equations, the method of manufactured solution (MMS) becomes a rigorous method for order of accuracy analysis [5]. The MMS can be applied to any codes using any numerical scheme as the verification process only deals with the correctness of the numerical solution in those codes, ignoring whether the solution itself is physical or not. However, MMS may require the manufactured solutions to be smooth so the analysis of nonsmooth solutions including shocks and discontinuities may not be well supported.

In fact, there are plenty of studies on 2D and 3D exact solutions to the unsteady incompressible Euler/Navier–Stokes equations. Al-Saif and Harfash [8] proposed a reduced differential transform method which is an iterative procedure to obtain some Taylor series solutions to the kinetically reduced local incompressible Navier–Stokes equations. Different isentropic Euler vortex problems were tested in Ref. [9], with characteristic and periodic boundaries being compared to obtain improved observed order of accuracy (OOA) for high order flux reconstructions. In Ref. [10], 2D steady and unsteady incompressible viscous flows in which the vorticity is proportional to the stream function transported by a uniform stream were studied. Shah et al. [11] presented a time-accurate numerical method using high order accurate compact finite difference to solve some incompressible Navier–Stokes problems. Dirichlet boundary condition was applied for general cases and periodic boundary condition was applied for periodic cases. Tavelli and Dumbser [12] proposed an arbitrarily high order accurate discontinuous Galerkin method for solving the three-dimensional incompressible Navier–Stokes equations on a staggered unstructured grid. Tsoutsanis et al. [13] compared structured and unstructured finite volume and Lagrange-remap methods and investigated the accuracy and efficiency of second- to ninth-order using the vortex pairing problem.

More analytical solutions to the incompressible Euler/Navier–Stokes equations can be also found in Refs. [14–19]. Plugging any of these solutions into the unsteady compressible Euler/Navier–Stokes equations only requires a source term to be derived for the energy equation [20,21], that is, the exact solution still satisfies the continuity and momentum equations. If the flow Mach number is small (0.1 or so), the magnitude of the source term is also very small, but it cannot be neglected when an OOA test is performed. For an incompressible flow for which only the exact initial solution exists, such as the 3D Taylor–Green vortex flow [22], a usually used approach is using an incompressible flow solver [12,23,24] or using a high order approximate process [25] to obtain a reference solution, and use the compressible solver solution to compare to the incompressible reference solution. Kokkinakis and Drikakis [26] did a comprehensive accuracy and computational cost comparison between different schemes (OOA ranging from 2 to 9) in a weakly compressible turbulent channel flow based on the direct numerical simulation result directly. Nonphysical MMS are also widely used for code verification of steady [27–29] and unsteady [30–32] CFD problems. In Ref. [27], a thorough code verification study using nonphysical MMS was performed on an unstructured finite volume CFD code. Different options including the governing equations, different boundary conditions, and different solvers were verified. Minion and Saye [31] applied various high order temporal time marching schemes for the incompressible Navier–Stokes equations and found that order reduction for the temporal accuracy can occur if applying time-dependent boundary conditions. Yu et al. [32] used MMS to verify a fluid–structure-interaction solver for both 2D and 3D cases. Instead of using a simultaneous refinement for the spatial and temporal spacing, a recursive Richardson extrapolation approach was used to subtract the spatial error from the total error, allowing the separation of the temporal error for unsteady problems.

Following up the work in Ref. [33] which compares the SENSEI results with NASA published numerical benchmark solutions by CFL3D [34,35] and FUN3D [36], this work is focused on the steady/unsteady Euler/Navier–Stokes solver code verification in the SENSEI CFD code [33,37,38]. SENSEI has the functionality to estimate the truncation error, which can be used to obtain high order discretization error estimates using error transport equations [39]. In this work, both exact solutions and manufactured solutions will be applied for code verification of the Euler/Navier–Stokes solvers in SENSEI.

The Computational Fluid Dynamics Code Base: Structured, Euler/Navier–Stokes Explicit–Implicit Solver

Overview of Structured, Euler/Navier–Stokes Explicit– Implicit Solver.

SENSEI is an acronym for structured, Euler/Navier–Stokes explicit–implicit solver, which is our in-house flow solver initially developed by Derlaga et al. [37]. The initial code was written in a very structured programing style which contained a lot of similar subroutines, with some modern Fortran features such as derived type data and pointers being used. Jackson et al. [38] refactored the SENSEI code completely using an object-oriented programing style and implemented the negative Spalart–Allmaras turbulence model into SENSEI for some 2D problems [40]. Turning SENSEI into an object-oriented code enables users to add new capabilities to SENSEI more easily due to the improvement of SENSEI’s modularity. Xue et al. [33] followed up the work of Jackson et al. [38] by implementing the 3D turbulence models including the negative Spalart–Allmaras [1] and Menter’s shear stress transport [4,41] and parallelized SENSEI using domain decomposition and message passing interface [42]. Now SENSEI is capable of solving the Euler, laminar, and Favre-averaged Navier–Stokes equations with Spalart–Allmaras and Menter’s shear stress transport models. Multiple turbulence cases have been compared by using SENSEI, CFL3D [35], and FUN3D [36].

Structured, Euler/Navier–Stokes explicit–implicit solver is a multiblock structured finite volume code and is embedded with several inviscid flux options including Roe’s flux difference splitting [43], Steger–Warming flux vector splitting [44], and Van Leer’s flux vector splitting [45]. Monotonic upstream-centered scheme for conservation laws (MUSCL) extrapolation and k-exact reconstruction are offered to achieve second- and high-order inviscid flux, respectively. The viscous and turbulent flux are second-order accurate using a central flux scheme after applying the Green–Gauss theorem. Various time marching schemes including explicit/implicit Runge–Kutta [46–48] and explicit/implicit Euler are offered. SENSEI has various boundary conditions such as slip/nonslip wall, supersonic/subsonic inflow/outflow, interblock boundaries, etc. Ghost cells as well as boundary face flux options are offered in SENSEI to enforce boundaries accurately, depending on the boundary type.

Turbulence Models Into the Favre-Averaged Navier–Stokes Equations Equations.

The governing equations used in the SENSEI CFD code can be generalized into a weak form as

\frac{\partial}{\partial t} \int_{Ω} U d Ω + \oint_{\partial Ω} (F_{i, n} - F_{ν, n}) d s = \int_{Ω} S d Ω

(1)

where

U

is the vector of conversed variables,

F_{i, n}

and

F_{ν, n}

are the inviscid and viscous flux normal components, respectively, (can be understood as the dot product of the second-order flux tensor and the normal vector of the surface) and

S

is the source term from either body forces, chemistry source terms, or the method of manufactured solutions. The Favre-averaged Navier–Stokes equations can be given in the framework shown in Eq. as

U = [\begin{matrix} \bar{ρ} \\ \bar{ρ} \tilde{u} \\ \bar{ρ} \tilde{v} \\ \bar{ρ} \tilde{w} \\ \bar{ρ} e_{t} \end{matrix}], F_{i, n} = [\begin{matrix} \bar{ρ} {\tilde{V}}_{n} \\ \bar{ρ} \tilde{u} {\tilde{V}}_{n} + {\hat{n}}_{x} p \\ \bar{ρ} \tilde{v} {\tilde{V}}_{n} + {\hat{n}}_{y} p \\ \bar{ρ} \tilde{w} {\tilde{V}}_{n} + {\hat{n}}_{z} p \\ \bar{ρ} h_{t} {\tilde{V}}_{n} \end{matrix}], F_{ν, n} = [\begin{matrix} 0 \\ {\hat{n}}_{x} {\tilde{τ}}_{x x} + {\hat{n}}_{y} {\tilde{τ}}_{x y} + {\hat{n}}_{z} {\tilde{τ}}_{x z} \\ {\hat{n}}_{x} {\tilde{τ}}_{y x} + {\hat{n}}_{y} {\tilde{τ}}_{y y} + {\hat{n}}_{z} {\tilde{τ}}_{y z} \\ {\hat{n}}_{x} {\tilde{τ}}_{z x} + {\hat{n}}_{y} {\tilde{τ}}_{z y} + {\hat{n}}_{z} {\tilde{τ}}_{z z} \\ {\hat{n}}_{x} {\tilde{Θ}}_{x} + {\hat{n}}_{y} {\tilde{Θ}}_{y} + {\hat{n}}_{z} {\tilde{Θ}}_{z} \end{matrix}], S = 0

(2)

where

\bar{ρ}

is density,

\tilde{u}, \tilde{v}, \tilde{w}

are the Cartesian velocity components, e_t is the total energy, h_t is the total enthalpy,

{\tilde{V}}_{n} = {\hat{n}}_{x} \tilde{u} + {\hat{n}}_{y} \tilde{v} + {\hat{n}}_{z} \tilde{w}

and the

{\hat{n}}_{i}

terms are the components of the outward-facing normal unit vector. The laminar and turbulent effects are combined in the definition for the viscous stresses which gives

\begin{matrix} {\tilde{τ}}_{x x} = 2 μ_{eff} (S_{x x} - \frac{1}{3} (\frac{\partial \tilde{u}}{\partial x} + \frac{\partial \tilde{v}}{\partial y} + \frac{\partial \tilde{w}}{\partial z})) - \frac{2}{3} \bar{ρ} k \\ {\tilde{τ}}_{y y} = 2 μ_{eff} (S_{y y} - \frac{1}{3} (\frac{\partial \tilde{u}}{\partial x} + \frac{\partial \tilde{v}}{\partial y} + \frac{\partial \tilde{w}}{\partial z})) - \frac{2}{3} \bar{ρ} k \\ {\tilde{τ}}_{z z} = 2 μ_{eff} (S_{z z} - \frac{1}{3} (\frac{\partial \tilde{u}}{\partial x} + \frac{\partial \tilde{v}}{\partial y} + \frac{\partial \tilde{w}}{\partial z})) - \frac{2}{3} \bar{ρ} k \\ {\tilde{τ}}_{x y} = 2 μ_{eff} S_{x y} \\ {\tilde{τ}}_{y z} = 2 μ_{eff} S_{y z} \\ {\tilde{τ}}_{z x} = 2 μ_{eff} S_{z x} \end{matrix}

(3)

where S_xx, S_xy, …, etc., are the mean strain rates, k is the turbulent kinetic energy, and the effective viscosity μ_eff is defined as

μ_{eff} = μ + μ_{T}

(4)

where μ is the dynamic viscosity and μ_T is the turbulent eddy viscosity.

\tilde{Θ}

in Eq. represents the heat conduction, work from the viscous stresses, and contribution from turbulent effects. Its components are given as

\begin{matrix} {\tilde{Θ}}_{x} = \tilde{u} {\tilde{τ}}_{x x} + \tilde{v} {\tilde{τ}}_{x y} + \tilde{w} {\tilde{τ}}_{x z} + k_{eff} \frac{\partial \tilde{T}}{\partial x} + {MDTT}_{x} \\ {\tilde{Θ}}_{y} = \tilde{u} {\tilde{τ}}_{x y} + \tilde{v} {\tilde{τ}}_{y y} + \tilde{w} {\tilde{τ}}_{y z} + k_{eff} \frac{\partial \tilde{T}}{\partial y} + {MDTT}_{y} \\ {\tilde{Θ}}_{z} = \tilde{u} {\tilde{τ}}_{x z} + \tilde{v} {\tilde{τ}}_{y z} + \tilde{w} {\tilde{τ}}_{z z} + k_{eff} \frac{\partial \tilde{T}}{\partial z} + {MDTT}_{z} \end{matrix}

(5)

where MDTT is the lumped term associated with molecular diffusion and turbulent transport (MDTT) in the energy equation of a given turbulence model and k_eff is effective thermal conductivity which is given by

k_{eff} = C_{p} (\frac{μ}{{Pr}_{L}} + \frac{μ_{T}}{{Pr}_{T}})

(6)

In Eq. (6), Pr_L is the (laminar) Prandtl number and Pr_T is the turbulent Prandtl number. In order to close the system, a turbulence model is needed to model the turbulent eddy viscosity μ_T and the turbulent kinetic energy k.

The Negative Spalart–Allmaras Turbulent Model.

The implementation of the original SA turbulence model can be given in the form shown in Eq. as

U = [\bar{ρ} \tilde{ν}], F_{i, n} = [\bar{ρ} \tilde{ν} {\tilde{V}}_{n}], F_{ν, n} = [{\hat{n}}_{x} {\tilde{Υ}}_{x} + {\hat{n}}_{y} {\tilde{Υ}}_{y} + {\hat{n}}_{z} {\tilde{Υ}}_{z}]

(7)

where

\tilde{ν}

is the turbulence working variable which is related to the kinematic eddy viscosity ν_T by

ν_{T} = \tilde{ν} f_{ν 1}, f_{ν 1} = \frac{χ^{3}}{χ^{3} + c_{ν 1}^{3}}, χ = \frac{\tilde{ν}}{ν}

(8)

and

{\tilde{Υ}}_{x} = \frac{\bar{ρ}}{σ} (ν + \tilde{ν}) \frac{\partial \tilde{ν}}{\partial x}, {\tilde{Υ}}_{y} = \frac{\bar{ρ}}{σ} (ν + \tilde{ν}) \frac{\partial \tilde{ν}}{\partial y}, {\tilde{Υ}}_{z} = \frac{\bar{ρ}}{σ} (ν + \tilde{ν}) \frac{\partial \tilde{ν}}{\partial z}

(9)

where ν is the kinematic viscosity. The source term for the SA model is given as

S = [\bar{ρ} (P - D) + \frac{c_{b 2}}{σ} \bar{ρ} [{(\frac{\partial \tilde{ν}}{\partial x})}^{2} + {(\frac{\partial \tilde{ν}}{\partial y})}^{2} + {(\frac{\partial \tilde{ν}}{\partial z})}^{2}] - \frac{1}{σ} (ν + \tilde{ν}) [\frac{\partial \bar{ρ}}{\partial x} \frac{\partial \tilde{ν}}{\partial x} + \frac{\partial \bar{ρ}}{\partial y} \frac{\partial \tilde{ν}}{\partial y} + \frac{\partial \bar{ρ}}{\partial z} \frac{\partial \tilde{ν}}{\partial z}]]

(10)

where P is a production term and D is a destruction term given by

P = c_{b 1} (1 - f_{t 2}) \tilde{S} \tilde{ν}, D = (c_{w 1} f_{w} - \frac{c_{b 1}}{κ^{2}} f_{t 2}) {[\frac{\tilde{ν}}{d}]}^{2}

(11)

The term d, defined as the distance to the nearest wall, and the modified vorticity

\tilde{S}

⁠, are given by

\tilde{S} = Ω + \frac{\tilde{ν}}{κ^{2} d^{2}} f_{ν 2}, f_{ν 2} = 1 - \frac{χ}{1 + χ f_{ν 1}}

(12)

where Ω is the magnitude of vorticity. The remaining terms are given as

f_{t 2} = c_{t 3} exp (- c_{t 4} χ^{2}), f_{w} = g {[\frac{1 + c_{w 3}^{6}}{g^{6} + c_{w 3}^{6}}]}^{1 / 6}, g = r + c_{w 2} (r^{6} - r), r = min (\frac{\tilde{ν}}{\tilde{S} κ^{2} d^{2}}, r_{lim})

(13)

The coefficients used in the SA model are given as

\begin{matrix} c_{b 1} = 0.1355 σ = 2 / 3 c_{b 2} = 0.622 κ = 0.41 \\ c_{w 1} = \frac{c_{b 1}}{κ^{2}} + \frac{1 + c_{b 2}}{σ} c_{w 2} = 0.3 c_{w 3} = 2 c_{ν 1} = 7.1 \\ c_{t 1} = 1 c_{t 2} = 2 c_{t 3} = 1.2 c_{t 4} = 0.5 \\ r_{lim} = 10 \end{matrix}

The MDTT term for the SA model is zero. Since the model does not model the turbulent kinetic energy, k, all terms related to k in Eq. (2) are dropped for the SA model.

In the original SA model, the working variable,

\tilde{ν}

⁠, is not allowed to be negative. The initial transients in the solution can drive

\tilde{ν}

to negative especially near the edge of a wake. The SA model ran into many such issues in our initial testing. The negative Spalart– Allmaras model [1] is a modified version of the SA model which allows

\tilde{ν}

to be negative. The negative Spalart–Allmaras model is the same as the SA model when

\tilde{ν}

is greater than or equal to zero. The modification is only applied when

\tilde{ν}

is negative. μ_T is set to zero for negative

\tilde{ν}

and the definitions for

Υ_{x}, Υ_{y}

and

Υ_{z}

are modified to be

{\tilde{Υ}}_{x} = \frac{\bar{ρ}}{σ} (ν + f_{n} \tilde{ν}) \frac{\partial \tilde{ν}}{\partial x}, {\tilde{Υ}}_{y} = \frac{\bar{ρ}}{σ} (ν + f_{n} \tilde{ν}) \frac{\partial \tilde{ν}}{\partial y}, {\tilde{Υ}}_{z} = \frac{\bar{ρ}}{σ} (ν + f_{n} \tilde{ν}) \frac{\partial \tilde{ν}}{\partial z}

(14)

The source term is also modified to be

S = [\bar{ρ} (P_{n} - D_{n}) + \frac{c_{b 2}}{σ} \bar{ρ} [{(\frac{\partial \tilde{ν}}{\partial x})}^{2} + {(\frac{\partial \tilde{ν}}{\partial y})}^{2} + {(\frac{\partial \tilde{ν}}{\partial z})}^{2}] - \frac{1}{σ} (ν + f_{n} \tilde{ν}) [\frac{\partial \bar{ρ}}{\partial x} \frac{\partial \tilde{ν}}{\partial x} + \frac{\partial \bar{ρ}}{\partial y} \frac{\partial \tilde{ν}}{\partial y} + \frac{\partial \bar{ρ}}{\partial z} \frac{\partial \tilde{ν}}{\partial z}]]

(15)

where the modified production and destruction term P_n and D_n are given as

P_{n} = c_{b 1} (1 - c_{t 3}) S \tilde{ν}, D_{n} = - c_{w 1} {(\frac{\tilde{ν}}{d})}^{2}

(16)

The coefficient f_n that is introduced in the modification is given by

f_{n} = \frac{c_{n 1} + χ^{3}}{c_{n 1} - χ^{3}}, c_{n 1} = 16

(17)

The Menter’s Shear Stress Transport Turbulence Model.

The governing equations for the Menter’s shear stress transport turbulence (⁠

k - ω

SST) model can also be cast into the form given in Eq. as

U = [\begin{matrix} \bar{ρ} k \\ \bar{ρ} ω \end{matrix}], F_{i, n} = [\begin{matrix} \bar{ρ} k {\tilde{V}}_{n} \\ \bar{ρ} ω {\tilde{V}}_{n} \end{matrix}], F_{ν, n} = {\hat{n}}_{x} Ξ_{x} + {\hat{n}}_{y} Ξ_{y} + {\hat{n}}_{z} Ξ_{z}

(18)

where ω is the specific rate of dissipation and

Ξ_{x} = [\begin{matrix} (μ + σ_{k} μ_{T}) \frac{\partial k}{\partial x} \\ (μ + σ_{ω} μ_{T}) \frac{\partial ω}{\partial x} \end{matrix}], Ξ_{y} = [\begin{matrix} (μ + σ_{k} μ_{T}) \frac{\partial k}{\partial y} \\ (μ + σ_{ω} μ_{T}) \frac{\partial ω}{\partial y} \end{matrix}], Ξ_{z} = [\begin{matrix} (μ + σ_{k} μ_{T}) \frac{\partial k}{\partial z} \\ (μ + σ_{ω} μ_{T}) \frac{\partial ω}{\partial z} \end{matrix}]

(19)

For the Menter’s SST model, the source term is given by

S = [\begin{matrix} P_{s} - β^{*} \bar{ρ} ω k \\ \frac{γ}{ν_{T}} P_{s} - β \bar{ρ} ω^{2} + 2 (1 - F_{1}) \frac{\bar{ρ} σ_{ω 2}}{ω} (\frac{\partial k}{\partial x} \frac{\partial ω}{\partial x} + \frac{\partial k}{\partial y} \frac{\partial ω}{\partial y} + \frac{\partial k}{\partial z} \frac{\partial ω}{\partial z}) \end{matrix}]

(20)

where

ν_{T} = μ_{T} / \bar{ρ}

and the term

P_{s}

is a production term which is given by

P_{s} = {\tilde{τ}}_{x x} \frac{\partial \tilde{u}}{\partial x} + {\tilde{τ}}_{x y} \frac{\partial \tilde{u}}{\partial y} + {\tilde{τ}}_{x z} \frac{\partial \tilde{u}}{\partial z} + {\tilde{τ}}_{x y} \frac{\partial \tilde{v}}{\partial x} + {\tilde{τ}}_{y y} \frac{\partial \tilde{v}}{\partial y} + {\tilde{τ}}_{y z} \frac{\partial \tilde{v}}{\partial z} + {\tilde{τ}}_{x z} \frac{\partial \tilde{w}}{\partial x} + {\tilde{τ}}_{y z} \frac{\partial \tilde{w}}{\partial y} + {\tilde{τ}}_{z z} \frac{\partial \tilde{w}}{\partial z}

(21)

Here the

\tilde{τ}

terms do not take the dynamic viscosity into consideration, i.e., they only include the effect of the turbulent eddy viscosity, which is calculated as

μ_{T} = \frac{\bar{ρ} a_{1} k}{max (a_{1} ω, Ω F_{2})}

(22)

where Ω is the vorticity magnitude. The terms F₁ and F₂ are used to blend an inner (expressed with subscript 1) and outer (expressed with subscript 2) coefficients, which are given as

F_{1} = tanh ({arg}_{1}^{4}), {arg}_{1} = min [max (\frac{\sqrt{k}}{β^{*} ω d}, \frac{500 ν}{d^{2} ω}), \frac{4 \bar{ρ} σ_{ω 2} k}{C D_{k ω} d^{2}}]

(23)

where

C D_{k ω} = max [2 \bar{ρ} σ_{ω 2} \frac{1}{ω} (\frac{\partial k}{\partial x} \frac{\partial ω}{\partial x} + \frac{\partial k}{\partial y} \frac{\partial ω}{\partial y} + \frac{\partial k}{\partial z} \frac{\partial ω}{\partial z}), 10^{- 20}]

(24)

and

F_{2} = tanh ({arg}_{2}^{2}), {arg}_{2} = max (2 \frac{\sqrt{k}}{β^{*} ω d}, \frac{500 ν}{d^{2} ω})

(25)

The blending process is performed by

ϕ = F_{1} ϕ_{1} + (1 - F_{1}) ϕ_{2}

(26)

where

ϕ

can be any of the coefficients. A production limiter [49] is applied, which replaces the

P_{s}

term in the k-equation by

min (P_{s}, 20 β^{*} \bar{ρ} ω k)

⁠. The coefficients in the model are given as

\begin{matrix} γ_{1} = \frac{β_{1}}{β^{*}} - \frac{σ_{ω_{1}} κ^{2}}{\sqrt{β^{*}}} γ_{2} = \frac{β_{2}}{β^{*}} - \frac{σ_{ω_{2}} κ^{2}}{\sqrt{β^{*}}} \\ σ_{k 1} = 0.85 σ_{ω 1} = 0.5 β_{1} = 0.075 \\ σ_{k 2} = 1.0 σ_{ω 2} = 0.856 β_{2} = 0.0828 \\ β^{*} = 0.09 κ = 0.41 a_{1} = 0.31 \end{matrix}

The MDTT terms for the Menter’s SST model is computed as

{MDTT}_{x} = (μ + \frac{μ_{T}}{σ_{k}}) \frac{\partial k}{\partial x}, {MDTT}_{y} = (μ + \frac{μ_{T}}{σ_{k}}) \frac{\partial k}{\partial y}, {MDTT}_{z} = (μ + \frac{μ_{T}}{σ_{k}}) \frac{\partial k}{\partial z}

(27)

The Menter’s Shear Stress Transport Turbulence Model With Vorticity Source Term.

There are many variants of the Menter’s SST model, such as the Menter’s SST model with vorticity source term (SST-V). The SST-V model is slightly different from the standard SST model shown earlier in that the SST-V model uses the vorticity magnitude to compute the production term instead of using the shear stresses. The production term given in Eq. is replaced with

P_{s} = μ_{T} Ω^{2} - \frac{2}{3} \bar{ρ} k (\frac{\partial \tilde{u}}{\partial x} + \frac{\partial \tilde{v}}{\partial y} + \frac{\partial \tilde{w}}{\partial z})

(28)

The production limiter used in the standard Menter’s SST model is still employed in the k-equation.

Finite Volume Discretization.

For finite volume discretization, the domain of interest, Ω, is partitioned into a sequence of nonoverlapping control volumes, Ω_i, such that

Ω = \cup_{i = 1}^{N_{ν}} Ω_{i}

⁠. The weak form in Eq. can be rewritten for each control volume as

\frac{\partial}{\partial t} \int_{Ω_{i}} U d Ω + \oint_{\partial Ω_{i}} (F_{i, n} - F_{ν, n}) d s = \int_{Ω_{i}} S d Ω

(29)

Denote the discrete solution of the finite volume method as

U_{h}

which is assumed to be constant in each control volume and approximates the control volume average of the exact solution. With the discrete steady-state residual given as

R_{h}

⁠, the discrete version of Eq. can be given in a semidiscrete form as

| Ω_{i} | \frac{\partial}{\partial t} U_{h} + R_{h} = 0

(30)

where

| Ω_{i} |

is the volume of Ω_i,

U_{h}

is the cell averaged solution vector, and

R_{h}

is the spatial residual vector, which is given in the following equation:

R_{h} = \sum_{1}^{f} (F_{i, n} - F_{ν, n}) Δ s - | Ω_{j} | S_{h}

(31)

where f is the cell face number (there are six faces for each hexahedral cell), $Δ s$ is the face area, and S_h is the cell averaged source term vector. It should be noted that we have not yet discretized the temporal term, which will be covered next.

Temporal Discretization.

In this work, explicit/implicit Runge–Kutta [46–48] and three point backward [50] temporal schemes are implemented. The temporal discretization deals with the numerical approximation of the

(\partial / \partial t) U_{h}

term in Eq. . Runge–Kutta is a family of time marching schemes, and SENSEI has explicit Runge–Kutta 2/4 stage schemes and singly diagonal implicit Runge–Kutta (SDIRK) 1/2/3 stage schemes implemented. A general form of the Runge–Kutta schemes can be seen in the following equation:

\begin{matrix} U_{h}^{n + 1} = U_{h}^{n} - Δ t \sum_{i = 1}^{s} b_{i} R_{h}^{i} \\ U_{h}^{i} = U_{h}^{n} - Δ t \sum_{j = 1}^{s} a_{i j} R_{h}^{j} \end{matrix}

(32)

where s is the number of stages, b_i and a_ij are constants, which can be found in Ref. [51]. For explicit Runge–Kutta schemes, all a_ij = 0 when $i \leq j$ ⁠, and for SDIRK, all a_ij = 0 when i < j. Many Runge–Kutta schemes can achieve high order accuracy for the temporal terms, but this work mainly uses second-order unless otherwise specified. Newton’s iteration is used for implicit Runge–Kutta schemes in this work.

A general form of the three point backward scheme used in this paper is given in the following equation:

\frac{3 U_{h}^{n + 1} - 4 U_{h}^{n} + U_{h}^{n - 1}}{2} = R_{h}^{n + 1}

(33)

As can be seen, three point backward requires the solutions from the previous two steps (no intermediate substep required) and it can achieve second-order accuracy for the temporal terms. SDIRK is used for the first step in SENSEI so that the three point backward can be initiated.

Nondimensionalization.

In SENSEI, the governing equations are nondimensionalized so that the convergence of the implicit solver performs better. Reference quantities for the density, temperature, length, speed of sound, and time are applied. Turbulence variables can also be nondimensionalized based on these reference quantities. The reference variables are given in the following equation:

\begin{matrix} length : & l_{ref} \\ time : & \frac{l_{ref}}{a_{ref}} \end{matrix} \begin{matrix} density : ρ_{ref} \\ velocity : a_{ref} \\ pressure : ρ a_{ref}^{2} \\ temperature : t_{ref} \end{matrix} \begin{matrix} dynamic viscosity : μ_{ref} \\ turbulence working variable : \frac{μ_{ref}}{ρ_{ref}} \\ turbulence kinetic energy : {(s f_{k} a_{ref})}^{2} \\ turbulence eddy dissipation ratio : \frac{μ_{ref}}{ρ_{ref} {(s f_{ω} a_{ref})}^{2}} \end{matrix}

(34)

where sf_k and

s f_{ω}

are scaling factors of the reference velocity a_ref for the turbulence kinetic energy k and the turbulence eddy dissipation ratio ω, respectively. sf_k and

s f_{ω}

are user-setting parameters to obtain a good nondimensionalization for the linear system of equations, especially for the turbulence equations.

Observed Order of Accuracy Test

For steady flows, the observed order of accuracy can be obtained from performing a systematic refinement over a series of spatial meshes. The process of performing code verification for unsteady flow and steady flow is generally similar but different in the following aspects. First, the steady flow does not have any temporal discretization error, while the unsteady flow does and the temporal discretization error should be quantified. Second, the order of the spatial discretization error can be different from the order of the temporal discretization error; therefore, attention should be paid to ensure these two errors are in the same order if possible. For some explicit schemes, it is difficult to use a large time-step size due to the numerical stability constraints. Due to this reason, implicit time marching schemes may be required. Also, usually different primitive variables for the Euler/Navier–Stokes equations may have different frequencies in time, and a common time-step size may not guarantee the spatial discretization errors be close to the temporal discretization errors for all variables. The subiterative error and round-off error need to be driven down small enough compared to the discretization error, in order to obtain convincing results.

Spatial Discretization Error.

Systematic mesh refinement should be performed on all spatial dimensions. Using a 1D problem as an example for simplicity, for a general pth order accurate numerical scheme, the spatial discretization error (spatial DE) can be written in the following equation:

ϵ_{h} = g_{x} Δ x^{p} + O (Δ x^{p + 1})

(35)

where h stands for the grid size in computational coordinates. When the spatial mesh is refined enough to be in the asymptotic range, the difference of the error magnitude on two consecutively refined meshes (refined by a factor of 2) is

2^{p}

times. In fact, the “consecutive” condition is not required to compute the OOA. A more general way of computing the observed OOA is shown in the following equation:

\hat{p} = \frac{ln (\frac{‖ ϵ_{r h} ‖}{‖ ϵ_{h} ‖})}{ln (r)}

(36)

where r represents the mesh refinement factor of one mesh over the other. Note that this refinement factor r is assumed to be the same in all spatial directions (as well as in time for unsteady flows).

Temporal Discretization Error.

The verification procedure for a temporal scheme with no spatial discretization is almost the same as that for the spatial scheme. The only difference is that the time-step is refined instead of the mesh spacing. The temporal discretization is given in Eq. . For explicit time marching schemes, the time-step size can be much smaller than using implicit time marching schemes

ϵ_{h} = g_{t} Δ t^{q} + O (Δ t^{q + 1})

(37)

where q is the formal order of accuracy for the temporal discretization error (temporal DE).

Spatial and Temporal Discretization Error.

Based on the spatial discretization error and temporal discretization error aforementioned, a simple form of the overall discretization error should combine both [52] and can be written in the following equation:

ϵ_{h} = g_{x} Δ x^{p} + g_{t} Δ t^{q} + O (Δ x^{p + 1}) + O (Δ t^{q + 1})

(38)

A more complicated form of the overall discretization error has spatial-temporal mixed terms [53], which is not considered in this paper, as the feature only exists for some specific schemes.

It should be mentioned that the temporal discretization error may be much smaller compared to the spatial discretization error for some unsteady problems, so attention should be paid to ensure that it does not adversely affect the final observed order of accuracy [5]. The temporal domain can be regarded as one extra dimension to the spatial domain, requiring a reasonable time-step size so that the temporal discretization error is of a similar order to the spatial discretization error.

Results

We tested several cases in this work, including 2D Euler convecting vortex flow, 2D Taylor–Green vortex flow, and cross-term sinusoidal (CTS) manufactured solutions. Both the steady and unsteady solutions for CTS are considered for the purpose of verification of the steady/unsteady solvers in SENSEI. All the test cases were run on a Dell Precision 5820 with an Intel Xeon Processor W2265 12 C 3.5 GHz 4.8 GHz Turbo-HT 19.25 MB (165 W) DDR4-2933, with 128 GB 8 × 16 GB DDR4 ECC (Error correction code). It should be noted that mesh refinement of a 3D unsteady case is easily stuck in an out of memory issue since we need to refine the mesh in all three dimensions, and unsteady cases need more memory space to store intermediate step results, only 2D cases are used for unsteady flow verification. SENSEI can output different norms including L₁, L₂, and $L_{\inf}$ norms, but their behaviors are quite similar for the test cases in this work so we will only show L norms.

Two-Dimensional Euler Convecting Vortex Flow.

The nondimensional solution of the 2D Euler convecting vortex flow [54] is given in the following equation:

\begin{matrix} u = & u_{\infty} - \frac{β}{2 π} exp \frac{1 - r^{2}}{2} (y - y_{0}) \\ v = & v_{\infty} + \frac{β}{2 π} exp \frac{1 - r^{2}}{2} (x - x_{0}) \\ T = & 1 - \frac{(γ - 1) β^{2}}{8 γ π^{2}} exp (1 - r^{2}) \\ ρ = & T^{\frac{1}{γ - 1}} \\ p = & ρ^{γ} \end{matrix}

(39)

where (x₀, y₀) is the vortex center, $r = \sqrt{({(x - x_{0})}^{2} + {(y - y_{0})}^{2})}$ ⁠, β is the vortex strength, and γ is the ratio of specific heats. For the Euler case in this work, β = 5, $γ = 1.4, u_{\infty} = 3, v_{\infty} = 0$ ⁠, with $l_{ref} = 1 m, a_{ref} = 1 m / s, t_{ref} = 1 K, ρ_{ref} = 1 {kg / m}^{3}$ ⁠. Therefore, this problem has a Mach number of 2.54. The reference values are unrealistic but it does not matter as the OOA test does not require realistic solutions. The domain of interest is $[0, 10 m] \times [- 5 m, 5 m]$ ⁠. The initial vortex center is at (5,0). Periodic boundaries are applied for the four faces.

As discussed earlier, proper time-step sizes should be determined so that the temporal discretization errors are of similar orders as the spatial discretization errors. Although SENSEI applies nondimensionalization to keep different equations on a similar order of magnitude, errors for different primitive variables are still likely to be not similar. A global consideration for the time-step size should be seriously taken. For this 2D Euler convecting vortex flow, a systematic separate OOA study is performed to determine a proper time-step size before a combined OOA study. The nondimensional discretization errors of the density and u velocity component applying the SDIRK 2 scheme are given in Fig. 1. It can be seen that a good nondimensional time-step size on the 128 × 128 mesh is about 0.025 (since $l_{ref} / a_{ref} = 1$ ⁠, the dimensional time-step size is 0.025 s), since similar orders of magnitude reductions can be achieved as the time-step is coarsened (moving to the right on the 128 × 128 curve) and as the mesh is coarsened (moving from the 128 × 128 curve to the 64 × 64 curve at the same dt = 0.025). However, due to the fact that Newton’s iteration requires the time-step size to be small to be stable, a nondimensional time-step size of 0.0125 is used, which will make the temporal discretization error smaller than the spatial discretization error, but they still have the same order of magnitude. When refining systematically in the combined order analysis, time-step sizes on other levels of meshes can be easily derived. Figure 2 shows the nondimensional error contours at $t = 2 s$ with $Δ t = 0.0125 s$ on the 128 × 128 mesh. The largest errors occur near the convecting vortex center (⁠ $x = 1 m$ ⁠), which is reasonable. The error contour plots can help detect potential implementation issues if there are some strange behaviors in the OOA test.

Fig. 1

View large Download slide

Two-dimensional Euler vortex flow separate order analysis (SDIRK 2): (a) discretization analysis (ρ) and (b) discretization analysis (u)

Fig. 2

Error contours of 2D Euler vortex flow (M∞=2.54, SDIRK 2): (a) density and (b) u velocity

View large Download slide

Error contours of 2D Euler vortex flow (⁠ $M_{\infty} = 2.54$ ⁠, SDIRK 2): (a) density and (b) u velocity

Figures 3 and 4 show the combined OOA results of applying the SDIRK 2 scheme and the three point backward scheme(⁠ $Δ t = 0.0125 s$ on the 128 × 128 mesh), respectively. For both schemes, an observed order of around 2 is achieved for all primitive variables if the refinement is enough, indicating the implementation should be correct. However, there are some discrepancies. First, the SDIRK 2 scheme has lower discretization errors compared to the three point backward scheme, which is expected as SDIRK 2 has more intermediate stages than three point backward. However, SDIRK 2 requires more memory space. Second, their behaviors for the OOA plots are different on coarse meshes, which suggests that SDIRK 2 enables a larger time-step size to balance the spatial and temporal discretization error, since the OOA is higher than 2 using SDIRK 2 on coarse level meshes. This 2D Euler convecting vortex flow can also be used to justify the correct implementation of the Euler solver and periodic boundary condition in SENSEI.

Fig. 3

Two-dimensional Euler vortex flow (M∞=2.54): SDIRK 2: (a) DE and (b) OOA

View large Download slide

Two-dimensional Euler vortex flow (⁠ $M_{\infty} = 2.54)$ ⁠: SDIRK 2: (a) DE and (b) OOA

Fig. 4

Two-dimensional Euler vortex flow (M∞=2.54): three point backward: (a) DE and (b) OOA

View large Download slide

Two-dimensional Euler vortex flow (⁠ $M_{\infty} = 2.54)$ ⁠: three point backward: (a) DE and (b) OOA

Two-Dimensional Taylor–Green Vortex.

The Taylor–Green vortex is a classic incompressible solution to the laminar Navier–Stokes equations. It has already been simulated using subgrid scale models and high resolution numerical methods [55]. The Taylor–Green vortex case is a good test case to verify whether the laminar NS solver is implemented correctly in SENSEI. A nondimensional 2D solution can be found in the following equation:

\begin{matrix} ρ = & ρ_{0} \\ u = & M cos (x) sin (y) exp - 2 t M / Re \\ v = & - M sin (x) cos (y) exp - 2 t M / Re \\ p = & p_{0} - \frac{ρ}{4} M^{2} (cos (2 x) + cos (2 y)) exp - 4 t M / Re \\ T = & \frac{p}{ρ R} \end{matrix}

(40)

where $ρ_{0} = 1, p_{0} = 0.713$ ⁠, Re is the Reynolds number, M is the Mach number, and R is the gas constant. For this laminar flow case, $l_{ref} = 1 m, a_{ref} = 374.17 m / s, ρ_{ref} = 1 {kg / m}^{3}$ ⁠. A case using Reynolds number of 1 and Mach number of 0.027 is tested. In fact, other configurations are also tried and show similar behavior. The domain of interest is $[- π m, π m] \times [- π m, π m]$ ⁠, which is one period spatially. Dirichlet boundaries are applied for the four boundaries. Note that the pressure term decays faster compared to the velocity terms, so attention may be needed when choosing the time-step size.

The first thing to verify an unsteady case is to determine a proper time-step size. Figure 5 shows the separate order analysis results for the pressure and density (Re = 1). It can be seen that a good nondimensional time-step size on the 256 × 256 is around 0.7. However, still being subject to the stability constraints of Newton’s iterative method, a nondimensional time-step size of 0.375 is used. This implies that the temporal discretization error is smaller than the spatial discretization error but their difference is not large (still in the same order). Similarly, when refining or coarsening systematically, time-step sizes on other levels of meshes can be easily calculated.

Fig. 5

View large Download slide

Two-dimensional Taylor–Green vortex flow separate order analysis (Re = 1, SDIRK 2): (a) discretization analysis (p) and(b)discretization analysis (u)

Figure 6 shows the dimensional solution contours obtained at $t^{*} = 12$ ⁠. The units for the pressure and u velocity contours are Pa and m/s, respectively. The decaying factors for the pressure and u velocity component at $t^{*} = 12$ are 0.28 and 0.53, respectively, where $t^{*}$ is the nondimensional time (dimensional time t over $l_{ref} / a_{ref}$ ⁠). Figure 7 shows the nondimensional error contours for the pressure and u velocity component at $t^{*} = 12$ ⁠. For the pressure, as it damps out 72% of its initial amplitude, larger errors only exist near the boundaries. In contrast, the u velocity has larger regions with large errors. The errors look reasonable with no strange behaviors shown. One interesting thing we found about this case is that on the coarsest mesh (8 × 8), only one time-step is used to iterate to the final time-step, since there is enough viscosity (Re = 1) to guarantee a stable iterating process.

Fig. 6

View large Download slide

Solution contours of 2D Taylor–Green decaying vortex (⁠ $t^{*} = 12$ ⁠, SDIRK 2): (a) pressure and (b) u velocity

Fig. 7

View large Download slide

Error contours of 2D Taylor–Green decaying vortex (SDIRK 2): (a) pressure and (b) u velocity

Figure 8 shows the OOA results of applying the SDIRK 2 scheme. For this case, the nondimensional discretization errors for the u and v velocity are quite close so their lines in Fig. 8 are highly overlapped. It is surprising to see that the discretization errors drop faster than second-order on some intermediate levels of mesh, and it causes the OOA to reduce on fine meshes. The overall observed OOA for the pressure and density variables are around 1.5, which is not perfect. One reason of this order reduction may be that the temporal discretization errors of SDIRK 2 are still smaller than the spatial discretization errors on the meshes of 64 × 64 and 128 × 128, and then return to comparable order of magnitude when refining further. Another potential reason would be that some components of the discretization errors of SDIRK 2 for this case have opposite signs on some intermediate meshes and causes the order increase on these meshes. It can be seen from Fig. 8 in which the line connecting the 32 × 32 and 256 × 256 points directly are close to second-order slope, especially taking into account that we use consecutive data points to calculate the OOA, instead of using a base data point as the reference. The 2D Taylor–Green vortex case is used to verify the NS implementation in SENSEI.

Fig. 8

View large Download slide

Two-dimensional Taylor–Green vortex flow (Re = 1, SDIRK 2): (a) DE and (b) OOA

Cross-Term Sinusoidal Manufactured Solutions.

Since the most rigorous way of verifying a code is the order of accuracy test, this work is focused on applying manufactured solutions. It should be emphasized that physical solutions are not required as the verification process only deals with the mathematics of a problem and thus can be used to determine the correctness of the implementation. After applying the governing equations to the manufactured solutions, SENSEI calculates the manufactured source terms numerically. The most important reason why the numerical form for the source terms is used is that there is no need to derive the analytic expression for the source terms (usually the source term expressions are extremely long and implementing them into the code is not easy). In this way, new various manufactured solutions can be implemented easily and flexibly in SENSEI. The manufactured solutions are three sine functions in each direction and cross term ones as well. The equation takes the general form

\begin{matrix} f = a_{1} + a_{2} sin (a_{3} π \frac{x}{l_{ref}} + a_{4} π) + a_{5} sin (a_{6} π \frac{y}{l_{ref}} + a_{7} π) + a_{8} sin (a_{9} π \frac{z}{l_{ref}} + a_{10} π) \\ + a_{11} sin (a_{12} π \frac{x y}{l_{ref}^{2}} + a_{13} π) + a_{14} sin (a_{15} π \frac{y z}{l_{ref}^{2}} + a_{16} π) + a_{17} sin (a_{18} π \frac{x z}{l_{ref}^{2}} + a_{19} π) \\ + a_{20} sin (a_{21} π \frac{t a_{ref}}{l_{ref}} + a_{22} π) \end{matrix}

(41)

where a_i are coefficients (mean flow and turbulence equations use different coefficients). For the manufactured solution, $l_{ref} = 1 m, a_{ref} = 340 m / s$ ⁠, which has a Reynolds number of 500 and a Mach number of 0.15. We run all cases to a final time of $0.01 s$ so that we can have enough variations in time. In our tests, the coarsest mesh has eight cells in each dimension of a curvilinear domain of about $[0, 1 m] \times [0, 1 m]$ ⁠. For unsteady CTS, we choose the coefficients a₂₁ intentionally so that $a_{21} π (t a_{ref} / l_{ref})$ is the same order as $a_{3} π (x / l_{ref})$ or $a_{6} π (y / l_{ref})$ ⁠, or equivalently $a_{21} π (Δ t a_{ref} / l_{ref})$ is the same order as $a_{3} π (Δ x / l_{ref})$ ⁠. Our choice ensures that the variations in space and time are similar. We will apply a separate order analysis to further show the discretization errors are similar in space and time.

Steady-State Cross-Term Sinusoidal Method of Manufactured Solution.

For the turbulence model implementation verification in SENSEI, the first thing to do is verify steady-state flows. The 3D steady-state manufactured solutions for some of the primitive variables in the Menter’s SST model are shown in Fig. 9. For the steady-state cases, since temporal schemes do not affect the verification results as the final solution is steady-state, only results using Euler implicit are shown for the steady-state CTS MMS. The Reynolds number for this CTS MMS case is about 78 × 10⁶ (low Reynolds numbers range from 5000 to 7000 are also tried and show similar behavior). If a lower dimension problem is solved then the relevant coefficients are set to 0. The grid size in the paper ranges from $4^{d}$ to $256^{d}$ (d is the dimension, i.e., two for 2D and three for 3D).

Fig. 9

Three-dimensional cross term sinusoidal manufactured solutions (Re=78 × 106): (a) density, (b) pressure, (c) turbulent kinetic energy, and (d) Eddy dissipation rate

View large Download slide

Three-dimensional cross term sinusoidal manufactured solutions (⁠ $Re = 78$ × 10⁶): (a) density, (b) pressure, (c) turbulent kinetic energy, and (d) Eddy dissipation rate

Prior to verifying Menter’s SST models in SENSEI, we have already tested the SA model implementation by comparing to some published results and using CTS MMS. The relevant results can be found in Ref. [33]. In this work, the observed order of accuracy results for the 2D and 3D Menter’s SST models are given in Fig. 10. For both the 2D and 3D CTS MMS cases, the values for the turbulent kinetic energy and turbulent eddy dissipation rate are chosen carefully so that the turbulent eddy viscosity has the same order O(1) as the viscous kinematic viscosity. Although the chosen value of the turbulent eddy viscosity and the viscous kinematic viscosity is $0.01 m^{2} / s$ ⁠, which is much higher than the actual value (about $1.8 \times 10^{- 5} m^{2} / s$ ⁠), the purpose here is to make the influence of the mean flow and turbulent flow approximately equal to verify the two systems together. Note that it is only for the code verification purpose. Actually, it is easier to converge the MMS cases using a physical value (e.g., $1.8 \times 10^{- 5} m^{2} / s$ ⁠) for the viscosity than using a nonphysical value (e.g., $0.01 m^{2} / s$ ⁠). Also, the observed order of accuracy using the physical value for the viscosity is still near two on the 2D and 3D CTS MMS cases. For the turbulent kinetic energy and turbulent eddy dissipation rate, the observed OOA is slightly greater than 2 on fine meshes for the 2D case, which may imply that some components of the discretization errors are canceled out. The Reynolds number for the 2D case and the 3D case shown are about 5600 and 7000, respectively. For the 3D Menter’s SST, the grid size of 256³ is not applicable as it encounters an out of memory issue. The steady-state turbulent CTS MMS is used to verify the steady-state turbulent flow implementation in SENSEI.

Fig. 10

View large Download slide

CTS MMS: observed order of accuracy: (a) 2D curvilinear CTS (Re = 5600) and (b) 3D curvilinear CTS (Re = 7000)

Unsteady Cross-Term Sinusoidal Method of Manufactured Solution.

For unsteady solutions, because different primitive variables have different periodic frequencies (their largest difference may be five times, still less than an order of magnitude), different primitive variables have different magnitudes of temporal discretization errors. Similar to what we have done earlier in this work for the unsteady Euler and laminar NS test cases, we need to find a proper time-step size for all variables on a given mesh and then coarsen or refine to get time-step sizes on other meshes. Figure 11 shows a systematic separate order analysis for the pressure variable and u velocity component using the $k - ω$ SST turbulence model. It can be seen that on the 256 × 256 mesh, a dimensional time-step size of $3.125 \times 10^{- 5}$ s is reasonable to guarantee that the temporal discretization errors are of the same order of magnitude as the temporal discretization errors for both variables. Coarsening the mesh to the 8 × 8 level, we can obtain the coarsest time-step size to be $0.001 s$ ⁠.

Fig. 11

View large Download slide

Two-dimensional CTS MMS separate order analysis (SDIRK 2): (a) discretization analysis (p) and (b) discretization analysis (u)

Figure 12 shows the 2D CTS OOA results using the Menter’s SST model. After enough refinement, the OOA is close to second-order for all the primitive variables. It is not surprising to see that the ω variable performs poorly on some intermediate meshes, but eventually improves to be second-order accurate on finer meshes. In this OOA test, some components of the discretization errors may cancel out each other because having the opposite signs or boost each other because having the same sign, and this causes the OOA oscillates between the 1.5 and 2.0 range. This behavior is related to the specific values in the CTS MMS and can disappear with some other values. This CTS MMS is used to verify the Menter’s SST model implementation in SENSEI.

Fig. 12

View large Download slide

Two-dimensional CTS MMS combined order analysis (Menter’s SST, SDIRK 2): (a) DE and (b) OOA

For the SA model, since the 1024 × 1024 mesh cannot be fully converged in the Newton’s iteration when using the same time-step size as the Menter’s SST model, we reduced the time-step sizes by 37.5% (now using 0.625 of the time-step size for the Menter’s SST model), which will make the temporal discretization errors smaller but still comparable to the spatial discretization errors. Figure 13 shows the OOA results using the SA model. For all the variables, the OOA is close to second-order when sufficiently refined. This SA CTS MMS is used to verify the SA turbulence implementation in SENSEI.

Fig. 13

View large Download slide

Two-dimensional CTS MMS combined order analysis (SA, SDIRK 2): (a) DE and (b) OOA

Figure 14 shows the nondimensional error contours using the SA model at $t = 0.01 s$ with a prescribed $Δ t = 1.953 \times 10^{- 5}$ s on the 256 × 256 mesh. It can be seen that the errors are uniformly distributed in the whole domain for both the u velocity component and the turbulence working variable.

Fig. 14

View large Download slide

Error contours of 2D CTS MMS (SA, SDIRK 2): (a) u velocity and (b) turbulence working variable

Conclusions

In this work, the SENSEI compressible flow solver was verified using the cross-term sinusoidal manufactured solutions and a few exact solutions. For all the solution verification cases presented in this work, fairly good observed order of accuracy results were obtained for the discretization error. Based on the results for all the test cases in this work, we believe that both the temporal and spatial terms for the Navier–Stokes equations (laminar or turbulent) have been implemented correctly in SENSEI. Based on the results in this work, we want to remind that special attention should be paid so that the temporal discretization errors are of the same order as the spatial discretization errors. When performing a systematic OOA test for unsteady flow solutions. Otherwise, some implementation errors for the temporal terms may not be easily found because the temporal discretization errors are usually smaller for general cases.

Acknowledgment

This material is partially based on research sponsored by the U.S. Air Force.

Funding Data

U.S. Air Force (Agreement No. FA865019-2-2204; Funder ID: 10.13039/100006831).

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

References

1.

Allmaras

,

S. R.

, and

Johnson

,

F. T.

,

2012

, “

Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model

,” Proceedings of the Seventh International Conference on Computational Fluid Dynamics (

ICCFD7

),

Big Island, HI

, July 9–13, pp.

1

–

11

.https://www.iccfd.org/iccfd7/assets/pdf/papers/ICCFD7-1902_paper.pdf

2.

Chien

,

K.-Y.

,

1982

, “

Predictions of Channel and Boundary-Layer Flows With a Low-Reynolds-Number Turbulence Model

,”

AIAA J.

,

20

(

1

), pp.

33

–

38

.10.2514/3.51043

Google Scholar

Crossref

3.

Wilcox

,

D. C.

,

2008

, “

Formulation of the k-w Turbulence Model Revisited

,”

AIAA J.

,

46

(

11

), pp.

2823

–

2838

.10.2514/1.36541

Google Scholar

Crossref

4.

Menter

,

F. R.

,

1994

, “

Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications

,”

AIAA J.

,

32

(

8

), pp.

1598

–

1605

.10.2514/3.12149

Google Scholar

Crossref

5.

Oberkampf

,

W. L.

, and

Roy

,

C. J.

,

2010

,

Verification and Validation in Scientific Computing

,

Cambridge University Press

, Cambridge, UK.

Google Scholar

Crossref

6.

Oberkampf

,

W. L.

, and

Trucano

,

T. G.

,

2002

, “

Verification and Validation in Computational Fluid Dynamics

,”

Prog. Aerosp. Sci.

,

38

(

3

), pp.

209

–

272

.10.1016/S0376-0421(02)00005-2

Google Scholar

Crossref

7.

Roy

,

C. J.

,

2005

, “

Review of Code and Solution Verification Procedures for Computational Simulation

,”

J. Comput. Phys.

,

205

(

1

), pp.

131

–

156

.10.1016/j.jcp.2004.10.036

Google Scholar

Crossref

8.

Al-Saif

,

A.-S. J.

, and

Harfash

,

A. J.

,

2019

, “

A New Approximate Analytical Solutions for Two-and Three-Dimensional Unsteady Viscous Incompressible Flows by Using the Kinetically Reduced Local Navier–Stokes Equations

,”

J. Appl. Math.

,

2019

, pp.

1

–

19

.10.1155/2019/3084394

9.

Spiegel

,

S. C.

,

Huynh

,

H.

, and

DeBonis

,

J. R.

,

2015

, “

A Survey of the Isentropic Euler Vortex Problem Using High-Order Methods

,”

AIAA

Paper No. 2015-2444.10.2514/6.2015-2444

10.

Hui

,

W.

,

1987

, “

Exact Solutions of the Unsteady Two-Dimensional Navier-Stokes Equations

,”

Z. Angew. Math. Phys. ZAMP

,

38

(

5

), pp.

689

–

702

.10.1007/BF00948290

Google Scholar

Crossref

11.

Shah

,

A.

,

Yuan

,

L.

, and

Khan

,

A.

,

2010

, “

Upwind Compact Finite Difference Scheme for Time-Accurate Solution of the Incompressible Navier–Stokes Equations

,”

Appl. Math. Comput.

,

215

(

9

), pp.

3201

–

3213

.10.1016/j.amc.2009.10.001

12.

Tavelli

,

M.

, and

Dumbser

,

M.

,

2016

, “

A Staggered Space–Time Discontinuous Galerkin Method for the Three-Dimensional Incompressible Navier–Stokes Equations on Unstructured Tetrahedral Meshes

,”

J. Comput. Phys.

,

319

, pp.

294

–

323

.10.1016/j.jcp.2016.05.009

Google Scholar

Crossref

13.

Tsoutsanis

,

P.

,

Kokkinakis

,

I. W.

,

Könözsy

,

L.

,

Drikakis

,

D.

,

Williams

,

R. J.

, and

Youngs

,

D. L.

,

2015

, “

Comparison of Structured-and Unstructured-Grid, Compressible and Incompressible Methods Using the Vortex Pairing Problem

,”

Comput. Methods Appl. Mech. Eng.

,

293

, pp.

207

–

231

.10.1016/j.cma.2015.04.010

Google Scholar

Crossref

14.

Panton

,

R. L.

,

2013

,

Incompressible Flow

,

Wiley

, Hoboken, NJ.

Google Scholar

Crossref

15.

White

,

F. M.

, and

Corfield

,

I.

,

2006

,

Viscous Fluid Flow

, Vol.

3

,

McGraw-Hill

,

New York

.

16.

Tryggeson

,

H.

,

2007

, “

Analytical Vortex Solutions to Navier-Stokes Equation

,” Ph.D. thesis,

Växjö University Press

, Växjö, Sweden.

17.

Wang

,

C.

,

1989

, “

Exact Solutions of the Unsteady Navier–Stokes Equations

,”

ASME Appl. Mech. Rev.

,

42

(

11S

), pp.

S269

–

S282

.10.1115/1.3152400

Google Scholar

Crossref

18.

Deissler

,

R. G.

,

1965

, “

Unsteady Viscous Vortex With Flow Toward the Center

,” NASA Lewis Research Center, Cleveland, OH,

Report No. NASA-TN-D-3026.

19.

Sengupta

,

T. K.

,

Sharma

,

N.

, and

Sengupta

,

A.

,

2018

, “

Non-Linear Instability Analysis of the Two-Dimensional Navier-Stokes Equation: The Taylor–Green Vortex Problem

,”

Phys. Fluids

,

30

(

5

), p.

054105

.10.1063/1.5024765

Google Scholar

Crossref

20.

Wang

,

H.

,

Tyson

,

W. C.

, and

Roy

,

C. J.

,

2019

, “

Discretization Error Estimation for Discontinuous Galerkin Methods Using Error Transport Equations

,”

AIAA

Paper No. 2019-2173.10.2514/6.2019-2173

21.

Tyson

,

W. C.

,

Yan

,

G. K.

,

Roy

,

C. J.

, and

Ollivier-Gooch

,

C. F.

,

2019

, “

Relinearization of the Error Transport Equations for Arbitrarily High-Order Error Estimates

,”

J. Comput. Phys.

,

397

, p.

108867

.10.1016/j.jcp.2019.108867

Google Scholar

Crossref

22.

Taylor

,

G. I.

, and

Green

,

A. E.

,

1937

, “

Mechanism of the Production of Small Eddies From Large Ones

,”

Proc. R. Soc. London, Ser. A

,

158

(

895

), pp.

499

–

521

.10.1098/rspa.1937.0036

23.

DeBonis

,

J.

,

2013

, “

Solutions of the Taylor–Green Vortex Problem Using High-Resolution Explicit Finite Difference Methods

,”

AIAA

Paper No. 2013-0382.10.2514/6.2013-0382

24.

Wang

,

Z. J.

,

Fidkowski

,

K.

,

Abgrall

,

R.

,

Bassi

,

F.

,

Caraeni

,

D.

,

Cary

,

A.

,

Deconinck

,

H.

, et al.,

2013

, “

High-Order CFD Methods: Current Status and Perspective

,”

Int. J. Numer. Methods Fluids

,

72

(

8

), pp.

811

–

845

.10.1002/fld.3767

Google Scholar

Crossref

25.

Bo

,

Y.

,

Wang

,

P.

,

Guo

,

Z.

, and

Wang

,

L.-P.

,

2017

, “

DUGKS Simulations of Three-Dimensional Taylor–Green Vortex Flow and Turbulent Channel Flow

,”

Comput. Fluids

,

155

, pp.

9

–

21

.10.1016/j.compfluid.2017.03.007

Google Scholar

Crossref

26.

Kokkinakis

,

I.

, and

Drikakis

,

D.

,

2015

, “

Implicit Large Eddy Simulation of Weakly-Compressible Turbulent Channel Flow

,”

Comput. Methods Appl. Mech. Eng.

,

287

, pp.

229

–

261

.10.1016/j.cma.2015.01.016

Google Scholar

Crossref

27.

Veluri

,

S. P. K.

,

2010

, “

Code Verification and Numerical Accuracy Assessment for Finite Volume CFD Codes

,” Ph.D. thesis,

Virginia Tech

, Blacksburg, VA.

28.

Salari

,

K.

, and

Knupp

,

P.

,

2000

, “

Code Verification by the Method of Manufactured Solutions

,”

Sandia National Laboratories

,

Albuquerque, NM

, Technical Report No. SAND2000-1444.

29.

Roy

,

C. J.

,

Ober

,

C. C.

, and

Smith

,

T. M.

,

2002

, “

Verification of a Compressible CFD Code Using the Method of Manufactured Solutions

,”

AIAA

Paper No. 2002-3110.10.2514/6.2002-3110

Google Scholar

Crossref

30.

Etienne

,

S.

,

Garon

,

A.

, and

Pelletier

,

D.

,

2009

, “

Code Verification for Unsteady Flow Simulations With High Order Time-Stepping Schemes

,”

AIAA

Paper No. 2009-169.10.2514/6.2009-169

Google Scholar

Crossref

31.

Minion

,

M. L.

, and

Saye

,

R.

,

2018

, “

Higher-Order Temporal Integration for the Incompressible Navier–Stokes Equations in Bounded Domains

,”

J. Comput. Phys.

,

375

, pp.

797

–

822

.10.1016/j.jcp.2018.08.054

Google Scholar

Crossref

32.

Yu

,

K. R.

,

Étienne

,

S.

,

Hay

,

A.

, and

Pelletier

,

D.

,

2015

, “

Code Verification for Unsteady 3-D Fluid–Solid Interaction Problems

,”

Theor. Comput. Fluid Dyn.

,

29

(

5–6

), pp.

455

–

471

.10.1007/s00162-015-0367-4

33.

Xue

,

W.

,

Wang

,

H.

, and

Roy

,

C. J.

,

2020

, “

Code Verification for 3D Turbulence Modeling in Parallel SENSEI Accelerated With MPI

,”

AIAA

Paper No. 2020-0347.10.2514/6.2020-0347

Google Scholar

Crossref

34.

Krist

,

S. L.

,

1998

, “

CFL3D User’s Manual (Version 5.0)

,”

National Aeronautics and Space Administration, Langley Research Center

, Hampton, VA.

35.

Bartels

,

R. E.

,

Rumsey

,

C. L.

, and

Biedron

,

R. T.

,

2006

, “

CFL3D Version 6.4-General Usage and Aeroelastic Analysis

,” NASA Langley Research Center, Hampton, VA,

Report No. NASA/TM-2006-214301.

36.

Biedron

,

R. T.

,

Carlson

,

J.-R.

,

Derlaga

,

J. M.

,

Gnoffo

,

P. A.

,

Hammond

,

D. P.

,

Jones

,

W. T.

,

Kleb

,

B.

, et al.,

2019

, “

FUN3D Manual: 13.5

,”

Langley Research Center

,

Hampton, VA

,

Report No. NASA/TM–2019–220271.

37.

Derlaga

,

J. M.

,

Phillips

,

T.

, and

Roy

,

C. J.

,

2013

, “

SENSEI Computational Fluid Dynamics Code: A Case Study in Modern Fortran Software Development

,”

AIAA

Paper No. 2013-2450.10.2514/6.2013-2450

38.

Jackson

,

C. W.

,

Tyson

,

W. C.

, and

Roy

,

C. J.

,

2019

, “

Turbulence Model Implementation and Verification in the SENSEI CFD Code

,”

AIAA

Paper No. 2019-2331.10.2514/6.2019-2331

39.

Wang

,

H.

,

Xue

,

W.

, and

Roy

,

C. J.

,

2020

, “

Error Transport Equation Implementation in the SENSEI CFD Code

,”

AIAA

Paper No. 2020-1047.10.2514/6.2020-1047

40.

Rumsey

,

C.

,

Smith

,

B.

, and

Huang

,

G.

,

2010

, “

Description of a Website Resource for Turbulence Modeling Verification and Validation

,”

AIAA

Paper No. 2010-4742.10.2514/6.2010-4742

41.

Menter

,

F. R.

,

Kuntz

,

M.

, and

Langtry

,

R.

,

2003

, “

Ten Years of Industrial Experience With the SST Turbulence Model

,”

Turbul., Heat Mass Transfer

,

4

(

1

), pp.

625

–

632

.https://cfd.spbstu.ru/agarbaruk/doc/2003_Menter,%20Kuntz,%20Langtry_Ten%20years%20of%20industrial%20experience%20with%20the%20SST%20turbulence%20model.pdf

42.

Gropp

,

W.

,

Lusk

,

E.

, and

Skjellum

,

A.

,

1999

,

Using MPI: Portable Parallel Programming With the Message-Passing Interface

, Vol.

1

,

MIT Press

, Cambridge, MA.

43.

Roe

,

P. L.

,

1981

, “

Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes

,”

J. Comput. Phys.

,

43

(

2

), pp.

357

–

372

.10.1016/0021-9991(81)90128-5

Google Scholar

Crossref

44.

Steger

,

J. L.

, and

Warming

,

R.

,

1981

, “

Flux Vector Splitting of the Inviscid Gasdynamic Equations With Application to Finite Difference Methods

,”

J. Comput. Phys.

,

40

(

2

), pp.

263

–

293

.10.1016/0021-9991(81)90210-2

Google Scholar

Crossref

45.

Van Leer

,

B.

,

1997

, “

Flux-Vector Splitting for the Euler Equation

,”

Upwind and High-Resolution Schemes

,

Springer

, Berlin, pp.

80

–

89

.

46.

Ascher

,

U. M.

,

Ruuth

,

S. J.

, and

Spiteri

,

R. J.

,

1997

, “

Implicit-Explicit Runge–Kutta Methods for Time-Dependent Partial Differential Equations

,”

Appl. Numer. Math.

,

25

(

2–3

), pp.

151

–

167

.10.1016/S0168-9274(97)00056-1

47.

Kennedy

,

C. A.

, and

Carpenter

,

M. H.

,

2016

, “

Diagonally Implicit Runge-Kutta Methods for Ordinary Differential Equations. A Review

,” NASA Langley Research Center, Hampton, VA,

Report No. NASA/TM-2016-219173.

48.

Jameson

,

A.

,

Schmidt

,

W.

, and

Turkel

,

E.

,

1981

, “

Numerical Solution of the Euler Equations by Finite Volume Methods Using Runge Kutta Time Stepping Schemes

,”

Proceedings of the 14th Fluid and Plasma Dynamics Conference

, Palo Alto, CA, June 23–25, p.

1259

.

Google Scholar

Crossref

49.

Menter

,

F.

,

1993

, “

Zonal Two Equation k-w Turbulence Models for Aerodynamic Flows

,”

AIAA

Paper No. 93-2906.10.2514/6.93-2906

50.

Wu

,

J.

,

Fan

,

L.

, and

Erickson

,

L.

,

1990

, “

Three-Point Backward Finite Difference Method for Solving a System of Mixed Hyperbolic-Parabolic Partial Differential Equations

,”

Comput. Chem. Eng.

,

14

(

6

), pp.

679

–

685

.10.1016/0098-1354(90)87036-O

Google Scholar

Crossref

51.

Ferracina

,

L.

, and

Spijker

,

M.

,

2008

, “

Strong Stability of Singly-Diagonally-Implicit Runge–Kutta Methods

,”

Appl. Numer. Math.

,

58

(

11

), pp.

1675

–

1686

.10.1016/j.apnum.2007.10.004

Google Scholar

Crossref

52.

Richards

,

S. A.

,

1997

, “

Completed Richardson Extrapolation in Space and Time

,”

Commun. Numer. Methods Eng.

,

13

(

7

), pp.

573

–

582

.10.1002/(SICI)1099-0887(199707)13:7<573::AID-CNM84>3.0.CO;2-6

Google Scholar

Crossref

53.

Kamm

,

J.

,

Rider

,

W.

, and

Brock

,

J.

,

2003

, “

Combined Space and Time Convergence Analysis of a Compressible Flow Algorithm

,”

AIAA

Paper No. 2003-4241.10.2514/6.2003-4241

54.

Yee

,

H. C.

,

Sandham

,

N. D.

, and

Djomehri

,

M. J.

,

1999

, “

Low-Dissipative High-Order Shock-Capturing Methods Using Characteristic-Based Filters

,”

J. Comput. Phys.

,

150

(

1

), pp.

199

–

238

.10.1006/jcph.1998.6177

Google Scholar

Crossref

55.

Drikakis

,

D.

,

Fureby

,

C.

,

Grinstein

,

F. F.

, and

Youngs

,

D.

,

2007

, “

Simulation of Transition and Turbulence Decay in the Taylor–Green Vortex

,”

J. Turbul.

,

8

(

8

), p.

N20

.10.1080/14685240701250289

Code Verification for the SENSEI CFD Code

Abstract

Introduction

The Computational Fluid Dynamics Code Base: Structured, Euler/Navier–Stokes Explicit–Implicit Solver

Overview of Structured, Euler/Navier–Stokes Explicit– Implicit Solver.

Turbulence Models Into the Favre-Averaged Navier–Stokes Equations Equations.

The Negative Spalart–Allmaras Turbulent Model.

The Menter’s Shear Stress Transport Turbulence Model.

The Menter’s Shear Stress Transport Turbulence Model With Vorticity Source Term.

Finite Volume Discretization.

Temporal Discretization.

Nondimensionalization.

Observed Order of Accuracy Test

Spatial Discretization Error.

Temporal Discretization Error.

Spatial and Temporal Discretization Error.

Results

Two-Dimensional Euler Convecting Vortex Flow.

Two-Dimensional Taylor–Green Vortex.

Cross-Term Sinusoidal Manufactured Solutions.

Steady-State Cross-Term Sinusoidal Method of Manufactured Solution.

Unsteady Cross-Term Sinusoidal Method of Manufactured Solution.

Conclusions

Acknowledgment

Funding Data

Data Availability Statement

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

Code Verification for the SENSEI CFD Code

Abstract

Introduction

The Computational Fluid Dynamics Code Base: Structured, Euler/Navier–Stokes Explicit–Implicit Solver

Overview of Structured, Euler/Navier–Stokes Explicit– Implicit Solver.

Turbulence Models Into the Favre-Averaged Navier–Stokes Equations Equations.

The Negative Spalart–Allmaras Turbulent Model.

The Menter’s Shear Stress Transport Turbulence Model.

The Menter’s Shear Stress Transport Turbulence Model With Vorticity Source Term.

Finite Volume Discretization.

Temporal Discretization.

Nondimensionalization.

Observed Order of Accuracy Test

Spatial Discretization Error.

Temporal Discretization Error.

Spatial and Temporal Discretization Error.

Results

Two-Dimensional Euler Convecting Vortex Flow.

Two-Dimensional Taylor–Green Vortex.

Cross-Term Sinusoidal Manufactured Solutions.

Steady-State Cross-Term Sinusoidal Method of Manufactured Solution.

Unsteady Cross-Term Sinusoidal Method of Manufactured Solution.

Conclusions

Acknowledgment

Funding Data

Data Availability Statement

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only