Spherical part of a tensor is wrongly computed in less than three dimensions

Summary

Spherical part of a tensor is wrongly computed in less than three dimensions

Steps to reproduce

Use turbulence models in 2D and 1D cases and other applications computing the deviatoric part of a tensor. The bug comes from using the definition for sph(T) in three dimensions for 2D and 1D cases.

D = T - 1/n*tr(T)*I

where n is the number of solved dimensions in the problem. We can also write,

D = T - sph(T)

with sph(T) = 1/n*tr(T)*I

Example case

Let's see the code in TensorI.H

//- Return the spherical part of a Tensor
template<class Cmpt>
inline SphericalTensor<Cmpt> sph(const Tensor<Cmpt>& t)
{
    return SphericalTensor<Cmpt>
    (
        (1.0/3.0)*tr(t)
    );
}

as we can see n is hardcoded as 3, where it should be linked to the actual number

of dimensions being resolved.

What is the current bug behaviour?

Spherical part of tensors are formally wrongly computed, the practical importance

of this in 2D and 1D has to be assessed.

What is the expected correct behavior?

The correction should be introduced in the spherical tensor method

sph(T) = 1/n*tr(T)*I

where n is the number of solved dimensions in the problem.

Environment information

• OpenFOAM version : any
• Operating system : any

Possible fixes

I was tracing the way how sph(T) is computed and since for tensor fields it is done through macros

I didn't find an easy way to make the changes. The macros are thought to use one-parameter methods

and here we need to pass the number of solved dimensions. Computing the dimensions in each calling

of sph(T) is unaffordable since it requires running solutionD() method for each element of the field. Additionally

in the scope of sph method you don't have access to the mesh since you are inside a plain Field.