A New “λ2” Term for the Spalart–Allmaras Turbulence Model, Active in Axisymmetric Flows

The new term belongs in the “basic,” free-shear-flow part of the Spalart–Allmaras (SA) model, and extends an idea of the Secundov team, incorporated in the νt-92 model. It detects transverse curvature in the distribution of the eddy viscosity ν~\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \tilde{\nu } $$\end{document}, so that it is passive in two-dimensional thin shear flows but potent especially in round jets. It eliminates the large over-prediction of the growth rate of such jets by the SA model, first detected by Birch in 1993. The originality is that the term is proportional to the middle eigenvalue λ2 of the Hessian operator of ν~\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \tilde{\nu } $$\end{document}. This is of course an empirical concept, but it is discriminating and rises when the distance r from the cylindrical axis becomes comparable with the length scale δ of the variations in the r direction. The inverted parabola is a prime example of such a distribution, and not unlike the ν~\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \tilde{\nu } $$\end{document} distribution in the round jet. The quantity λ2 is not infinitely differentiable, but it is free of singularities, and unlike the νt-92 version, is not dependent on two large quantities cancelling. The core term added to the Lagrangian derivative Dν~/Dt\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ D\tilde{\nu }/Dt $$\end{document} is simply cb3 λ2ν~\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \tilde{\nu } $$\end{document}, where cb3 is a new constant. The computing cost of calculating and ordering the eigenvalues is moderate. We have no proof of well-posedness for the new equation set, but the evidence so far is favorable, both in structured and unstructured grids. The λ2 term is calibrated on a fully-developed round jet, and tested in nine other cases, either 2D flows or flows in which r » δ, finding that in the latter it is negligible as expected. This is although the cb3 constant is rather large, namely 6. The λ2 term is not strong enough to make a mature vortex fully relaminarize as would be desirable, but the eddy viscosity drops by 74%. The raw λ2 term reduces the eddy viscosity in pipe flow, where that is detrimental; therefore, in the final model, it is multiplied by a function of the rSA parameter of the SA model, which is a measure of wall proximity. The λ2 term appears to be a safe addition to the SA model, and its application in different codes and to a variety of flows to be desirable.


Introduction
The difficulty in achieving accuracy for both planar and round flows has a long history in Reynolds-Averaged Navier-Stokes (RANS) modeling, and the SA model initially of 1992 is no exception (Spalart and Allmaras 1994). In fact, Birch first showed particularly poor results in the round jet for SA already in 1993 (Birch 1993). A recent review is given by Pollard (Pollard 2020). A well-regarded experiment is by Hussein, Capp, and George (Hussein et al. 1994) and produced a higher growth rate than earlier tests, but still far lower than the SA rate. This flow was not part of the calibration base of the SA model, and the fullydeveloped static jet, as opposed to the initial development from the nozzle, is not of high importance in the industry. Some adjustments to similar eddy-viscosity models, notably that of Wilcox (2008), rest on an empirical hypothesis of Pope (1978) based on mean vortex stretching and therefore on the velocity field. Secundov had a very different approach in the ν t -92 model, the correction being based on the eddy-viscosity field only (Shur et al. 1995). In such one-equation models, which have no exact terms, there seems to be no rigorous or even intuitive reason in favor of one or the other approach, or one that combines the two. Here, we first present Secundov's insight, and explain why we use different mathematics to inject it into the SA model, or any other.
The ν t -92 model includes a C 3 term, added to those which are active in 2D flows and are symbolized by the dots …: It combines the Laplacian of ̃ (using the variable name of the SA model for brevity), and an N 2 term which is the L 2 norm of the gradient of the L 2 norm of the gradient of ̃ . In a 2D distribution with a convex profile ̃(y), the two terms cancel because N 2 = − 2̃∕ y 2 . In contrast, in an axisymmetric distribution with only the dependence ̃(r), we have the Laplacian and, if it is convex, N 2 is the opposite of the first term in the right-hand side of (2), so that the C 3 term brings out the first derivative, divided by r: The 1/r factor is the transverse curvature of the isosurfaces, and is the key to the new term.
A useful example is an inverted parabola: In this field, the two terms in the Laplacian (2) are equal: As a result, the term is usable. The ν t -92 C 3 term however has some substantial disadvantages. It rests on the cancellation between two terms in a 2D flow, but numerically, this is (2) Δ̃= 2r 2 + 1 rr 1 3 not always achieved. One reason is that the most natural differentiation stencils for the two terms are not the same. In a simple structured grid, the second derivative at point i is easily calculated using the points (i − 1, i, i + 1), but the usual stencil for N 1 is (i − 1, i, i + 1) so that the stencil for N 2 is (i − 2, i − 1, i, i + 1, i + 2). This can cause significant numerical differences especially near a maximum. Another consideration is the absolute values, which prevent a fully smooth dependence of N 2 on the original distribution, again near an extremum. Essentially, N 2 is defined and has equal limits from the left and the right of the extremum, but is undefined at the extremum itself. The other issue with the ν t -92 model's C 3 term, namely that it is not 0 in 1D concave distributions, is not as clear-cut, considering that removing local minima is not obviously undesirable physically. Nevertheless, these considerations motivated us not to use exactly Secundov's approach.
Turning now to the Pope vortex-stretching term, conceived to explain the behavior of the round jet, we consider its empirical justification to be fairly weak (Pope 1978). In addition, the use of its absolute value as in the k − ω model (Wilcox 2008) is certainly not justified, because logically the effect of vortex shortening should be the opposite of that of vortex stretching. An additional argument against the Pope term is our keen interest in the rectilinear vortex (Spalart and Garbaruk 2019). In that flow, there is a clear need to reduce the eddy viscosity of the SA and most models, yet vortex stretching is absent. In contrast, a term which detects transverse variations as in (2) has much potential. In both the round jet and the vortex, the eddy-viscosity distributions are similar to inverted parabolas (although in the jet, the second derivative at r = 0 is 0). This leads us to the objective of approximately extracting the two terms in (2), in a robust manner, in general eddy-viscosity fields.

Discussion of the Eigenvalues
In a simple axisymmetric distribution with only the dependence ̃(r), the two terms in the right-hand side of (2) give the two eigenvalues of the Hessian operator, which motivated us to involve these eigenvalues in the model (the third eigenvalue is 0, this field having no dependence on x). The reason the eigenvalues can be identified with particular directions is symmetry when ̃ is a function only of r: the eigenvectors must be in the r, θ and x directions. Thus, the axisymmetric case is well-understood. In a general 3D flow, the r direction is not defined, but the eigenvalues are still accessible and potentially helpful, assuming they are well-behaved. If the ̃ distribution versus r is convex, which is the common case, the two eigenvalues are similar, and therefore the middle eigenvalue is representative of the transverse-curvature effect. For a 2D thin shear flow, the middle eigenvalue is close to 0, and therefore the term essentially vanishes, as desired. The middle eigenvalue, by construction, is smaller in magnitude than the larger components of the Hessian, which are already present in the existing diffusion term of the SA model. This makes it look like a safe quantity to use. The eigenvalue, called λ 2 , has the dimension 1/T, which led instantly to the form c b3 λ 2 ̃∕ . The c b3 constant will be positive, to bring about a reduction of ̃ in typical cases, with r < 0.
The trigonometric solution of the cubic equation provides the real eigenvalues of a 3-by-3, symmetric, real matrix such as the Hessian of ̃ . The computing cost is not out of line with the exponential and trigonometric functions the SA model already involves.
The edge of the turbulent region is a critical aspect of all turbulence models. Recall that the eddy viscosity drops to 0 linearly (Spalart and Allmaras 1994). This is part of the inverted parabola (4), extended with ̃ = 0 for r > R. The second derivative then has a Dirac distribution at r = R: whereas (̃∕ r)∕r merely jumps from the finite value −2̃0∕R 2 to 0. Therefore, (6) gives the largest eigenvalue λ 3 , and λ 2 is finite so that the new term is of the same order as the existing diffusion terms of the SA model; in contrast with (6), it is not singular (which would be unacceptable). The edge behavior is altered only quantitatively, with the turbulent ramp propagating outwards more slowly than it does without the new term.
Fields with singular λ 2 on the axis would be those proportional to r or to R − r, in two or three dimensions. However, such conical apexes are immediately suppressed by the diffusion terms, which makes them irrelevant.

Full Formulation of the New Term, and Adjustment of the Constant
As mentioned, the basic model is now D̃∕Dt = ⋯ + c b3 λ 2 ̃∕ . The division by σ gives the message through the notation that it belongs with the diffusion terms, but the distinction is not rigorous. As mentioned below, a reduction of production could, physically, have been an option. However, arguments of simplicity and robustness are important and motivated the choice we made over altering the production.
The fully-developed round jet is the driving case for calibration. Figures 1 and 2 show that the far-field is well-reproduced with c b3 = 6, and we employ this value from here on. The experimental datasets are by Albertson et al. (1948), while the range for spreading rate is taken from Wilcox (2006).
Incompressible round jet with uniform inlet (Albertson et al. 1948) and (Wilcox 2006) Fig. 2 Incompressible round jet with developed pipe inlet (Albertson et al. 1948;Wilcox 2006) Figure 3 details the Hessian components H(i, j) in (x, r, θ) axes, and the eigenvalues in the fully-developed jet. H(2, 2), in the r direction, is the first term of the right-hand side of (2), while H(3, 3) in the circumferential direction is the second term in (2), which is always negative. The components involving the streamwise direction x are small, as expected. The singularity in (6) is visible near r/x = 0.29, and becomes λ 3 as expected so that it has no impact on the model. The eigenvectors (not shown) include one in the circumferential direction, for the eigenvalue equal to H(3, 3); this is required by symmetry. In the bulk of the flow, the other two are close to the x and the r direction, respectively.
The budget of ̃ is shown in Fig. 4. The λ 2 term is the dominant negative one over much of the flow, offsetting most of the production; it then drops to small values near the edge of the turbulent region, where the basic SA diffusion term dominates as usual.

Tests in Other Free Shear Flows
The first tests involve the initial evolution of round jets, shown in Figs. 5 and 6. The standard SA (with Compressibility Correction, which is not strong at this Mach number) model is not inaccurate in the mixing-layer region, but it terminates the potential core too fast, which the λ 2 term remedies well. Thus, the calibration in the fully-developed jet is successful in the entry region also. This case shows that the λ 2 term is not problematic in regions where ̃∕ r > 0 , in the inside of the mixing layer, although it was aimed at the opposite sign. Its role is rather weak, compared with the basic diffusion terms, the mixing layer being thin compared with radius.
Next, Fig. 7 covers the Mature Vortex case, studied by the authors in Spalart and Garbaruk (2019). Recall the reasoning in that paper that the best possible behavior for a model at  vortex maturity is to fully suppress the eddy viscosity (or Reynolds stresses in other types of models). We find that the new term, without any adjustment, reduces it by almost a factor of 4; in other words, it is partly successful. This response could be expected, the eddyviscosity profile being quite close to an inverted parabola as seen in Fig. 7a. The vortex core radius, accordingly, is reduced by almost a factor of 2; in the developed state, lengths scale like √̃m ax t . Curiously, the circulation overshoot, precisely the flow feature we argued in Spalart and Garbaruk (2019) is unphysical, is not reduced at all as seen in Fig. 7b. It just happens at a smaller radius. This remaining flaw is fully corrected by the SARC model, which can be combined with the new term. The SARC model extinguishes the eddy viscosity in a mature vortex (Spalart and Garbaruk 2019), resulting a laminar vortex as seen in Fig. 7b.
We now turn to the vortex generated by a NACA0012 wing with a round tip, using experimental data of Chow et al. (1997). The robustness of the new model in this fully three-dimensional flow and the weak change at the trailing edge (X/C = 0) are positive signs, but the accuracy along the vortex is poor. The new term visibly improves the crossflow velocity, but nowhere as much as the SARC models does. The improvement for the axial velocity is also disappointing. The SARC version is considered to be among the current best for this and other vortex flows, and with it, the λ 2 term makes essentially no difference; this is consistent with the mature-vortex findings (Fig. 8).
The next case again shows a sizeable effect for the λ 2 term, but again a weaker one than from the RC correction. The slat area of the DLR-F15 high-lift aerofoil (LEISA configuration) is shown (Wild et al. 2006). The λ 2 term reduces the eddy viscosity by about half, which could be expected based on near-circular contours, and then thins out the slat wake. We have no quantitative flow and turbulence measurements, which would allow us to determine which model is more accurate; in fact, eddy viscosity is not a directly measurable quantity. The effect is much weaker for the flap-cove eddy viscosity in which the contours are oval rather than near-circular (not shown), and also for the pressure and skin-friction distributions (not shown). The SA and SA-λ 2 models are closer to experiment than the SARC model, which disturbs the separation on the flap, leading to an excess of circulation around the flap and the entire system Fig. 9.
The final case in this section is a supersonic axisymmetric base flow. Figure 10 compares the SA and SA-λ 2 models each with the experiment of Herrin and Dutton (1994) in real space. The new model provides a substantial improvement, at least once the wake is well-developed. Figures, not shown, indicate that the eddy viscosity is lower by a factor of about 2. The wake width is consequently smaller, and closer to the experiment. In contrast, the velocity near the base shows no improvement.
The improvement is confirmed in Fig. 11, showing that the return of the centerline velocity to freestream is noticeably slower, although not very accurate. In this figure, "CC" refers to the SA compressibility correction (http:// turbm odels. larc. nasa. gov/). Unfortunately, the combination of two corrections reduces the eddy viscosity more than is desirable. The poor accuracy of CFD in the nether-region, x/R < 1.5, is also confirmed.   Herrin and Dutton (1994) 1 3

Wall-Bounded Flows
The new term was applied to fully-developed pipe flow, in which the eddy-viscosity distribution is very similar to that in round free shear flows. Unfortunately, the consequent reduction of the eddy viscosity is not desirable at all, since the SA model is quite accurate in raw form; this is shown in Fig. 12. Therefore, in order to disable the correction in wall-bounded flows, it is multiplied by a function of the r SA parameter of the SA model, which is 1 near a wall and 0 in free shear flows. The final form of the correction is This correction was present, in particular, in Figs. 9 and 10.

Further Tests
A number of flows were simulated with and without the λ 2 correction, and found to have very little sensitivity to it. Most are listed at the TMR (http:// turbm odels. larc. nasa. gov/). First is the plane jet. In 2D flows, usually λ 2 is not the 0 that comes from the third direction, but rather the small eigenvalue associated with the slow streamwise evolution. Next is the boundary layer in zero pressure gradient, which is unresponsive for the same reasons as the plane jet, and is of course an essential flow module. The fully-developed flow in a square duct, which contains weak streamwise vortices and is a key test case for nonlinear eddy-viscosity models, is similarly unaffected. Two flows with axisymmetric geometries, but boundary-layer thicknesses much smaller than the radius of curvature, show very little effect. They are the CS0 low-speed diffuser and the Bachalo-Johnson transonic bump (http:// turbm odels. larc. nasa. gov/). The NACA 4412 aerofoil at high angle of attack 13.87°, with a small recirculation zone near the trailing edge, has a similar weak response. The backward-facing step is also immune, probably for the reason invoked about regarding the flap cove, namely that the eddy-viscosity contours are elongated ovals, so that the λ 2 correction is near-passive over most of the region. The effect is also weak both for pressure and skin friction over the Common Research Model, which is a clean wing in transonic conditions. Finally, a flow field with Vortex Generators in a supersonic boundary layer was re-calculated. In such a flow, the SARC version of the model again has far more control than the new term does.

Discussion of Stability, and of other Options
Numerical stability is an essential requirement for a PDE. We have no mathematical proof of Hadamard well-posedness, but we have performed numerical tests centred on oscillations of the eddy viscosity. In one dimension, any ̃(y) distribution makes λ 2 zero, so that there is no effect. In two dimensions, the effect on distributions of the type sin(k x x) sin(k y y) is stabilizing, as was verified by calculating the effect on the integral of ̃′ 2 . This is the third graph in Fig. 13: the contribution to the variance of ̃ is zero or negative. This would also apply with different wavenumbers in x and y, no matter how large; in other words, short waves are not amplified, which is strongly related to Hadamard's criterion. Incidentally, the λ 2 contours in the center graph illustrate the nonlinear character of this quantity, and the discontinuities in its gradient.
In terms of iterative convergence, a linearization of λ 2 appears difficult, so that true Newton methods may not be possible, but in the NTS code we used, in general the convergence was not significantly degraded.
The spatial smoothness of the eigenvalues is also of interest. If the two-dimensional Hessian is h ij , the eigenvalues are h 11 + h 22 ± √ h 11 − h 22 2 + h 2 12 ∕2 , so that they are smooth functions of the h ij entries unless the matrix is proportional to the identity matrix (that is, h ij ∝ ij ), which gives the square root an argument of zero, and therefore infinite slope. However, the argument of the square root being quadratic in h ij , the slope of the eigenvalue versus the h ij entries themselves is finite; this was visible in Fig. 13b.
Another source of non-smoothness is the situation when two eigenvalues that are each a smooth function switch roles, between being the middle and the largest-magnitude eigenvalue; this is seen in Fig. 3b. In summary, λ 2 is a continuous but not everywhere differentiable function of space. We note that the original SA model is built on differentiable functions, but its solutions are not differentiable everywhere, and particularly at the edge of the turbulent region where the diffusion coefficient falls to 0. Therefore, the non-smoothness of λ 2 is of the same degree.
The conception of the new term was fairly arbitrary, as was the case with the original SA terms even within the usual constraints such as dimensional analysis. A possibility that was entertained was to form a non-dimensional combination of the eigenvalues, such as the ratio of λ 2 to the largest eigenvalue, or to the trace of the Hessian matrix. Various ratios of the coefficients of the characteristic polynomial were also considered. A function of that quantity could then be used to modify any term in the model, whether production or diffusion. An empirical motivation would be sought. It seems that the only candidate is the coefficient of λ in the characteristic polynomial, which is dominated by the product of H 22 and Fig. 13 Fields for the distribution ̃= 10 + sin (x) sin (y). No z dependence was assumed. Left, eddy viscosity; centre, λ 2 ; Right, effect of the λ 2 term on the variance of ̃ 1 3 H 33 (in the (y, z) plane). However, the dimensional analysis is unfavorable, as the dimension of this quantity is 1/T 2 . An additional issue with non-dimensional quantities is that in a region of uniform values, such a quantity is undefined, essentially 0/0. Such regions are very common, especially outside the turbulent region. In contrast, in those regions, λ 2 is simply equal to 0.
As it is, the new term is formally grouped with the diffusion terms inside the 1/σ bracket, but this classification is not rigorous (any more than presenting the c b2 term as a diffusion term was). The argument that the term derives only from the eddy-viscosity distribution, and not the velocity field, has the merit of conceptual simplicity, and we believe using subtle features of the velocity field and its derivatives would carry more risk of unintended consequences in flow fields very different from those considered in the design and testing.
Another possibility which was considered was to argue that the clear need for the new term is to reduce the eddy viscosity. Thus, an option would be to take min(λ 2 , 0). However, λ 2 is positive only in small volumes, and that has not degraded the results. Thus, simplicity was preferred, in the absence of any strong argument, although the r SA -dependent factor in (7) does add a little complexity.
A consideration is the application to other empirical turbulence models. The λ 2 quantity can be a useful tool with any model, to isolate transverse-curvature effects. However, the Hessian of a quantity such as k or ε would not have the dimension 1/T, which is very convenient in the Hessian of ̃ . Therefore, further work would be needed even if the modeler accepted λ 2 as a legitimate measure.

Summary
A correction to the Spalart-Allmaras turbulence model, aimed at eliminating its wellknown failure in round jets, was created. We used the notation SA-λ 2 for it, but the name SA-TC for Transverse Curvature is more descriptive and will be used in the long run, including in the TMR (http:// turbm odels. larc. nasa. gov/). It has the novel feature of using an eigenvalue of the Hessian operator of the eddy viscosity, which has favorable properties in axisymmetric fields, is formally convenient, rather easy to introduce, and generally well-behaved including at the abrupt edge of the turbulent region. Its heuristic character is consistent with the rest of the model. The new term is fully successful in the round jet, and partly in the vortex [where the SARC model is fully successful (Spalart and Garbaruk 2019)]. There is little evidence so far that it severely degrades convergence, let alone cause instabilities, although any reduction in the eddy viscosity has the potential of precluding steady solutions (which we observed in the case of the LEISA high-lift aerofoil). The term is rather simple to calculate, and not likely to invite coding errors. It can be combined with any of the other corrections to SA, such as RC or QCR (http:// turbm odels. larc. nasa. gov/). An unexpected fact is that the value c b3 = 6 is significantly larger than 1, and also than the original constants such as c b2 = 0.622. This echoes values such as c r2 = 12 in the SARC model, but this latter constant renders streamwise curvature effects, which have long been known to be surprisingly strong. The transverse-curvature effect appears similarly powerful, and we recall that other RANS models have suffered from low accuracy in the static round jet.
We note that the new term was adjusted in the fully-developed jet, but also makes a favorable difference in the early development of the jet. It also reduces the eddy viscosity in recirculating regions, such as that under the slat of a modern wing, but we have not been in a position to determine whether this raises or lowers the model's accuracy. The effect is favorable in the wake of an axisymmetric base. A fair number of other flows were simulated, and showed very little effect.