Repository navigation
add support for axisymmetric jets with SA turbulence model #2923
Description
Activity
I had a look at this. Reporting what I found and what I have so far, rather than opening a PR blind the issue has no description, so please correct me if this is not what you meant.
What is already there
SA with
AXISYMMETRIC= YESis not blocked: the guard "Axisymmetry is currently only supported for KIND_TURB_MODEL chosen as SST" was removed in d90e8e1 (Jan 2024), which addedResidualAxisymmetricDiffusiontoCSourceBase_TurbSA. The axisymmetric nozzle case runs withKIND_TURB_MODEL= SAwithout error.What looks incomplete
SA adds only the diffusion part of the axisymmetric contribution. The two comparable transported scalars add convection as well:
scalar axisym. convection axisym. diffusion Jacobian SST k,omega(turb_sources.hpp)-yinv*V*rho*v*phiyes yes species Y(species_sources.cpp)-yinv*V*rho*v*phiyes yes SA nu_tildemissing yes none SU2 discretizes the SA convective term in divergence form (
Flux = a0*nue_i + a1*nue_j, no density). In cylindrical coordinates that divergence isd(u*nue)/dx + (1/r) d(r*v*nue)/dr, so the 2D planar fluxes miss(1/r)*v*nue. Equivalently: for a uniformnuein an incompressible axisymmetric flow the convective residual should vanish, becausedu/dx + dv/dr + v/r = 0; without the term it evaluates to-nue*v/rper unit volume. That error is proportional to the radial velocity, so it is small in a boundary layer and not small in a jet or a diverging nozzle.Change I have on a branch
In
CSourceBase_TurbSA, mirroring the SST naming:ResidualAxisymmetricDiffusion->ResidualAxisymmetricConvectionDiffusion, addingResidual[0] -= yinv*Volume*v*nueandJacobian_i[0][0] -= yinv*Volume*v(velocity only, no density, to match SA's density-free convective flux).- The call is moved to after
Jacobian_i[0][0] *= Volume;, otherwise the new Jacobian term would be scaled by the volume twice. Residual behaviour of the existing diffusion term is unchanged by the move. Coord_i[1]andV_i[idx.Velocity()+1]are registered as preaccumulation inputs when axisymmetric.Coord_i[1]was already used by the existing diffusion term without being registered, so the discrete adjoint of axisymmetric SA was missing that dependency.
Verification so far
- Non-axisymmetric is untouched: everything new sits behind
if (axisymmetric), andturb_naca0012_saat iteration 5 still gives-12.037372, -16.384159, 1.080346, 0.018385, 20, -3.456391, 20, -4.641246, 0against the committed-12.037323, -16.384159, 1.080346, 0.018385, 20, -3.457510, 20, -4.641262, 0(serial GCC 15.2 build, differences only in the linear-solver columns). - The term is active in the axisymmetric nozzle case (
axisymmetric_rans/air_nozzleswitched to SA, cold start, 1000 iterations): rms[nu_tilde] -7.327 before, -7.155 after.
That second point only shows the term does something, not that it is right. Before I open a PR I would like to validate it against a 3D wedge of the same geometry, which is the only check that really settles the sign and the factor, and add an axisymmetric SA regression case so this path is covered in CI. Two questions:
- Is this the gap you had in mind, or were you after a validation/tutorial jet case?
- For the regression entry, would you prefer an SA variant of the existing
air_nozzlecase (mesh already in the TestCases repo), or a dedicated jet case (new mesh)?
SA cannot simulate axisymmetric jets, it needs a separate submodel:
https://www.researchgate.net/publication/349801309_A_New_l2_Term_for_the_Spalart-Allmaras_Turbulence_Model_Active_in_Axisymmetric_FlowsThanks, that is the piece I was missing I had been looking at the discretization rather than the model. For anyone reading later, that paper is documented on the NASA TMR as SA-TC (transverse curvature free-shear correction), Spalart & Garbaruk, Flow, Turbulence and Combustion 107:245-256 (2021).
The term added to the RHS of the SA equation is
max(1 - r, 0)^3 * c_b3 * lambda_2 * nu_hat / sigma, c_b3 = 6, sigma = 2/3with
rthe same as infw, so it switches off near walls and is fully active in free shear, andlambda_2the middle eigenvalue of the Hessian ofnu_hat. The mechanism is that for a plane jet the out-of-plane Hessian component is zero and the term vanishes, while for a round jet the azimuthal component of the Hessian of an axisymmetric field isnu_hat_r / r, which is not - that is what fixes the round-jet/plane-jet anomaly.I would like to implement it. From reading the code, the pieces are:
- A new
SA_OPTIONS= TRANSVERSE_CURVATUREentry and the corresponding flag inSAParsedOptions. - The Hessian of
nu_tilde, which does not exist anywhere in SU2 today. The route that does not require new infrastructure is the existing auxiliary-variable gradient: setnAuxVar = nDimfor the SA variable, store the components ofgrad(nu_tilde)inAuxVar, and callSetAuxVar_Gradient_GG/LS, which already has an MPI communication type. The result is the Hessian at each point. lambda_2itself: in 2D planar, the eigenvalues of the 2x2 in-plane Hessian plus a zero third one; in 2D axisymmetric, the 2x2 block plusnu_tilde_r / ras the out-of-plane eigenvalue; in 3D the analytic eigenvalues of a symmetric 3x3. Middle = sum - min - max in all cases.- The source term and its Jacobian, which is straightforward since the term is linear in
nu_tilde.
One thing I would like your opinion on before I get too far. Neither the paper nor the TMR page prescribes how to discretize the second derivatives, and a Green-Gauss gradient of a Green-Gauss gradient is noisy on unstructured meshes, which is exactly where
lambda_2would be evaluated in a shear layer. If you know how this was done in the reference implementation, or have a preference (Green-Gauss twice, least-squares on the gradients, something else), that would save me guessing.For validation I would target the round jet spreading rate against the plane jet, which is what the paper is about. That needs a jet mesh, which means something in the TestCases repo; tell me if you would rather I propose a case there first.
Separately, and not part of this: while reading
CSourceBase_TurbSAI noticed the axisymmetric contribution only includes the diffusion part. SST and the species solver both also add the convection term-(1/y) * rho * v * phi, and SU2 writes the SA convective flux in divergence form, so the cylindrical divergence appears to be missing(1/y) * v * nu_tilde. Equivalently, for a uniformnu_tildein an incompressible axisymmetric flow the convective residual should vanish becausedu/dx + dv/dr + v/r = 0, and without that term it does not. I have a patch for it but have not validated it against a 3D wedge yet, so I am not claiming it is right. Should I open a separate issue for that, or is it deliberate?- A new
Related finding, independent of the convection term discussed above: the axisymmetric source of the flow equations (
CEulerSolver.cpp/flow_sources.cpp) took turbulence variable 0 askfor any RANS model. With SA that variable isnu_tilde, so the2/3 rho kterms usednu_tildeas if it were a kinetic energy (wrong quantity and units). This is fixed in #2329 (commit 771b26e):kis used only with SST, and the2/3 rho kterm only wherekis part of the stress tensor.Two practical consequences for the wedge validation: the SA axisymmetric nozzle (
axisymmetric_rans/air_nozzle) results change with that fix, independently of any SA-equation change, so compare against a 3D wedge with both in place, and expect the reference values of that case to move. I have not checked the convection term itself against a wedge.Thanks, that's useful. One correction for anyone following: #2329 was closed unmerged, and I think the fix is now in your draft #2946 (tkeNeeded = rans && SST in CEulerSolver, plus the SA air_nozzle variant). Develop still passes nu_tilde as k at CEulerSolver.cpp:2207. I'll validate the SA convection term against a 3D wedge on top of #2946 so the two effects are separated.
Add support for axisymmetric jets with SA turbulence model. SA needs an additional submodel:
https://www.researchgate.net/publication/349801309_A_New_l2_Term_for_the_Spalart-Allmaras_Turbulence_Model_Active_in_Axisymmetric_Flows