[TomKopp]

Is it possible to use

fvm::div(phi, U)
- fvc::div(tau)

instead of

fvm::div(phi, U)
- fvm::laplacian(fluid.nu(), U)
- (fvc::grad(U) & fvc::grad(fluid.nu()))

which is the standard formulation in non-Newtonian, steady-state solvers?

I would like to modify the stress tensor before the momentum equation is solved (i.e. include the elongational viscosity in normal direction).

For the energy equation this is already working:

// calculate stress tensor
    
volTensorField gradU = fvc::grad(U);
volTensorField tau = fluid.nu() * (gradU + gradU.T());
tau.component(tensor::XX) = tau.component(tensor::XX) / max(dimensionedScalar ("VSMALL", dimensionSet(0,2,-1,0,0,0,0), VSMALL), fluid.nu()) * elongFluid.nuElong();
tau.component(tensor::YY) = tau.component(tensor::YY) / max(dimensionedScalar ("VSMALL", dimensionSet(0,2,-1,0,0,0,0), VSMALL), fluid.nu()) * elongFluid.nuElong();
tau.component(tensor::ZZ) = tau.component(tensor::ZZ) / max(dimensionedScalar ("VSMALL", dimensionSet(0,2,-1,0,0,0,0), VSMALL), fluid.nu()) * elongFluid.nuElong();

// energy equation

tmp<fvScalarMatrix> TEqn
(
    fvm::div(phi, T)
  - fvm::laplacian(alpha, T)
  - 1.0/Cp * (tau && gradU)
);
      
TEqn().relax();

solve(TEqn());

However, when I use Gauss linear scheme for div((nu*(grad(U)+grad(U).T()))) the time step continuity error gets very high after only one step.

Does anyone have an idea how to solve the problem? Thank you in advance!

[TomKopp]

Even if I use

volTensorField gradU = fvc::grad(U);
volTensorField tau = fluid.nu() * (gradU + gradU.T());

without the elongational viscosity and then for the momentum equation

fvm::div(phi, U)
- fvc::div(tau)

it is not working (time step continuity error very high).

How can I tell OpenFOAM to solve

fvc::div(tau)

correctly?

[TomKopp]

Can it be due to the use of an explicit (fvc) instead of an implicit derivative (fvm)? Do I have to use an transient solver or should it also work for a steady-state solver like simpleFoam?