[K62]

Dear OpenFOAM-users,

My aim is to run a RANS simulation with a given Reynolds-stress tensor field . This will come from a previous RANS simulation of the same case but will be (slightly) modified, i.e. . In order to write a solver, which is able to do that, I used simpleFoam (OF2.3.0). The UEqn.C file contains the following lines

tmp<fvVectorMatrix> UEqn
    (
        fvm::div(phi, U)
      + turbulence->divDevReff(U)
      ==
        fvOptions(U)
    );

The pointer

turbulence->divDevReff(U)

is defined e.g. in LRR.C

tmp<fvVectorMatrix> LRR::divDevReff(volVectorField& U) const
{
...

return
        (
            fvc::div(R_)
          + fvc::laplacian(nut_, U, "laplacian(nuEff,U)")
          - fvm::laplacian(nuEff(), U)
        );
}

This leads to the following implementation of my new solver

tmp<fvVectorMatrix> UEqn
    (
        fvm::div(phi, U)
      + fvc::div(R)
      + fvc::laplacian(nut, U, "laplacian(nuEff,U)")
      - fvm::laplacian(nu+nut, U)
      ==
        fvOptions(U)
    );

Also the line

turbulence->correct();

is commented, since it is not necessary anymore.

My workflow is like that.

STEP 1: I run a simulation using the un-modified simpleFoam solver until the residuals go to machine accuracy or drop a value of 1e-10 or so. I then run the utility R, which computes the Reynolds stresses using the Boussinesq approximation

tmp<volSymmTensorField> kEpsilon::R() const
{
    return tmp<volSymmTensorField>
    (
        new volSymmTensorField
        (
            IOobject
            (
                "R",
                runTime_.timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
            k_.boundaryField().types()
        )
    );
}

STEP 2: The output field R is used as an input for the modified solver described above, so that no additional transport equation for is solved and it stays (so far) unmodified. Since, the field is exactly the one which is the output of the step 1, I naively (?) assumed that using it as an input for the modified solver will lead to the same U and p fields as in STEP 1.

For the cases I am looking at I observe that for a smooth geometry (2D hill), the modified solver converges but gives velocity and pressure fields which show differences at the wall of up to 5%. I guess this has to do with the fact how boundaries are treated by the different differential operators (laplacian, div and grad) in OpenFOAM, but I wasn't able to get to the core of this problem. Does anybody have an idea about that?

Best wishes,
Martin

[zym604]

Thanks you Martin, your thought is exactly the same as mine. Regarding the 5% difference, did you use the same system files in both case? For example, fvscheme.

[zym604]

Hi Martin, although I didn't work out your method, I found another way to deal with it. I want to share it with you:

1. Calculate nut according to R;
2. Add nut to nu as an internal field;
3. Solve the problem as laminar flow;
4. Done.