[dlahaye]

Abstract
The objective of this humble endeavor is to describe how the celebrated Pitz-Daily incompressible tutorial (or any its variants) for simpleFoam can be extended with a Reynold stress tensor model for turbulence. This endeavor has been pursued by many (including the giants in this community) in the past and is by no means new. Our endeavor can in fact be stated as reproducing the tutorial [1,2] (among possible various others). Our motivation for stirring in this old pot once again is contribute to the documentation on how to set-up the case. We wish to complement existing tutorials with our insights. The Reynolds stress model solves for the six (in 2D four) components of the symmetric stress tensor R and the turbulent dissipation epsilon. The magnitude of the tensor R is equal to epsilon. Our biggest take-away lessons is that the over-relaxation of the equation (should this field?) for epsilon is essential to obtain convergence of the model.

1/ Introduction
We aim at extending the Pitz-Daily incompressible tutorial for simpleFoam with a Reynold Stress Tensor Model (RSM) for turbulent flow closure. Motivation for this endeavor is to be able to increase resolution in modeling flow with large recirculation. The Reynolds stress model is a six turbulence equation model for the six components of the symmetric Reynolds stress tensor. The Launders, Reece and Rodi variant of this model [3,4] is implemented as the LRR model in OpenFoam. Previous discussions on the use of the LRR model in this forum include [5,6]. Below we outline features of the the Pitz-Daily test case in Section 2, our understanding on how the model should be set up starting from the tutorial model in Section 3 and on how the model converges with the LRR RASModel activated.

2/ Pitz-Daily out of the box in OpenFoam-v1906
The Pitz-Daily test case as defined in OpenFoam-v1906 converges in 283 iterations (run simpleFoam) and meet the criterium for y+ on the mesh (run simpleFoam -postProcess -func yPlus -latestTime). The 283 iterations differs significantly from the tutorial [2]. This difference is likely due to different settings used.

3/ Case Set-Up
In this section we explain how the Pitz-Daily tutorial can be modified to run with the LRR RASModel in three steps.

Step 1/3: Setting up the system-folder
In the first step we set up the system-folder by editting the system-folder/fvSolution-file and the system-folder/fvSchemes-file. We edit the system-folder/fvSolution-file and add settings to solve for R including relaxation as described in e.g. [5]. We edit the system-folder/fvSchemes-file and add setting to discretize div(R) and div(phi,R) and settings for the wallDist dictionary. Sample settings for div(R) and div(phi,R) can be copied from other tutorials. To find such setting, use e.g. for i in `find /opt/OpenFOAM/OpenFOAM-v1906/tutorials/ -name fvSchemes`; do echo $i; grep $i "div(R)" ; done .

Step 2/3: Setting up the 0-folder
In the second step we edit the 0-folder. We do *not* define the boundary conditions for the components of the tensor R for inlet, outlet and wall patches directly. Instead, we resort to a two-substep procedure in which in the first substep, a simulation using a two-equation turbulence model runs for a number of iterations. In the second substep, the stress tensor R is computed in post-processing stage for the latest iteration. This computed stress tensor is used as boundary value setting for the simulation using the LRR model.

More specifically, in the first substep, we run simpleFoam on the out-of-the-box tutorial. At the last iteration (iteration 283 in our case) we compute the stress tensor R by post-processing computed results and running simpleFoam -postProcess -func R -latestTime. This command produce a file called turbulenceProperties: R in the latest time folder (folder 283). This file can be renamed to R and copied to the 0-folder. Observe that the settings for the turbulence equation used in the first run is used as input to compute R. One can thus *not* modify the constant-folder/turbulenceProperties-file prior to computing R. One can visualize components and magnitude of the tensor R in paraview.

Step 3/3: Running the simulations
In the third third we (prepare to) run simpleFoam with the LRR model. We edit the constant-folder/turbulenceProperties-file to change RASModel from current value into LRR. We edit the system-folder/fvSolution file and choose sufficiently low relaxation parameters for R and epsilon. We can now run simpleFoam with LRR. Alternatively, we can avoid copying the 283-folder/R-file to the 0-folder and run simpleFoam with LRR and startTime equal to the latestTime of the previous simulation.

4/ Discussion and Conclusions
We discussed how the incompressible tutorial Pitz-Daily for simpleFoam can be modified to include the Reynolds stress tensor wonder. We wonder whether the preprocessing using a two-equation turbulence model can be avoided.

The above description is likely to be incomplete and fragile as time evolves. We theferore wonder whether a tutorial using the LRR model can be made available as tutorial in the OpenFoam distribution.

References
[1]: Turbulence Steady State, Tutorial Six of TU Vienna: https://www.cfd.at/sites/default/fil...xampleSix.pdf: shows change divSchemes in fvSchemes to take R into account: no settings explained for fvSolution: decrease in turbulent viscosity is shown;
[2]: Introduction to Stationary Turbulence Modeling, Jozsef Nagy, https://www.youtube.com/watch?v=-46pgYQYER8. The LRR model is discussed from minute 18 of 26 onwards
[3]: NASA Langley Turbulence Modeling Resource on Launders, Reece and Rodi model: https://turbmodels.larc.nasa.gov/rsm-ssglrr.html ;
[4]: Ansys Fluent Manual online
[5] Jasak on cfd.online.com on how to use the LRR Model: details of how to run the LRR model and showing convergence of the residual for the stress tensor components: Problems with the RSM in simpleFoam , 2006.
[6]: discussion regarding the LRR RASModel on cfdonline.com: settings for volSymmTensorField R in fvSolution are shown. No settings for fvSchemes are given unfortunately: Usage of R and LRR , 2010.

[HPE]

Hi,

- Finally a nice-looking thread unlike some of the questions here. Thanks.
- Could you please clarify if there is a question from your side, though? I didn't quite get the purpose of the post?

Thank you.

[dlahaye]

Hi,

Thank you for your words of appreciation. These words come as an encouragement that I much appreciate.

I indeed failed to phrase questions explicitly. I am wondering though:

1/ It took me a while to realize that the stress tensor model was indeed available. Can documentation be improved, e.g., by subdividing the documentation for RAS models into one-equation, two-equation and stress-tensor model?

2/ It took me a while to realize how to set up the model. Can documentation be improved, e.g., by providing in the OpenFoam release an Allrun script that performs all of the above via scripting?

3/ Can running the two-equation model to in the set-up phase be avoided? Does this require implementing wall functions for the stress tensor?

Thanks again and kind wishes. Domenico Lahaye

[HPE]

My apologies for my English skills, yet I still couldn't fully grasp what you mean by your questions, e.g.:

>> Can documentation be improved, e.g., by subdividing the documentation for RAS models into one-equation, two-equation and stress-tensor model?

Do you suggest/wish to improve the header file or extended code guide documentation on these models?

[vitor.geraldes@ist.utl.pt]

Can you explain why "The magnitude of the tensor R is equal to epsilon." . The units of R (in fact R/rho) are m2/s2 and the units of epsilon are m2/s3.

[klitz]

Hi, I have tried to use your steps to try and use LRR. However, I am facing the below issue:

--> FOAM FATAL ERROR:
[2] cannot be called for a calculatedFvPatchField
    on patch airfoil of field R in file "airfoil_aoa0/processor2/200/R"
    You are probably trying to solve for a field with a default boundary condition.

I know that my mesh and setup is ok because I have succesfully run simulations with SST and SA for 20 different angle of attack. I can see that when using -func R the airfoil and inlet boundary conditions changes to calculated. However, not the outlet it stays as the same. According to your steps I shouldnt need to change boundary conditions to start the simulation with LRR.

Thank you.

[dlahaye]

Originally Posted by vitor.geraldes@ist.utl.pt

Can you explain why "The magnitude of the tensor R is equal to epsilon." . The units of R (in fact R/rho) are m2/s2 and the units of epsilon are m2/s3.

Apologies for the late reply and the glitch.

Sincere thanks for pointing out an ambiguity.

Should be: "The trace of the tensor R divided by 2 is equal to k." See e.g. book Versteeg 2007, Section 3.7.3.

[dlahaye]

Originally Posted by klitz

Hi, I have tried to use your steps to try and use LRR. However, I am facing the below issue:

--> FOAM FATAL ERROR:
[2] cannot be called for a calculatedFvPatchField
    on patch airfoil of field R in file "airfoil_aoa0/processor2/200/R"
    You are probably trying to solve for a field with a default boundary condition.

I know that my mesh and setup is ok because I have succesfully run simulations with SST and SA for 20 different angle of attack. I can see that when using -func R the airfoil and inlet boundary conditions changes to calculated. However, not the outlet it stays as the same. According to your steps I shouldnt need to change boundary conditions to start the simulation with LRR.

Thank you.

Pls. elaborate. What solver and solver settings are you using to obtain this error message?

[klitz]

Originally Posted by dlahaye

Pls. elaborate. What solver and solver settings are you using to obtain this error message?

I am using simpleFoam.
fvSchemes:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
    grad(U)         cellLimited Gauss linear 1;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linearUpwind grad(U);
    div(phi,k)      Gauss upwind;
    div(phi,omega)  Gauss upwind;
    div(phi,nuTilda) Gauss upwind;
    div(phi,epsilon) Gauss upwind;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
    div((nu*dev2(T(grad(U)))))  Gauss linear;
    div(R)          Gauss linear;
    div(phi,R)      Gauss upwind;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method meshWave;
}


// ************************************************************************* //

fvSolution

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-7;
        relTol          0.1;
        smoother        GaussSeidel;
    }

    "(U|omega|nuTilda|epsilon|R)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-8;
        relTol          0.01;
    }

    k
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-9;
        relTol          0.01;
    }


   /* pFinal
    {
        solver          GAMG;
        tolerance       1e-04;
        relTol          0.1;
        smoother        GaussSeidel;
    }
    pcorrFinal
    {
        solver          GAMG;
        tolerance       1e-04;
        relTol          0.1;
        smoother        GaussSeidel;
    }
    "(U|k|omega)Final"
    {
        $U;
        tolerance       1e-04;
        relTol          0;
    }
*/

}

SIMPLE
{
    nNonOrthogonalCorrectors 2;
    consistent yes;

    residualControl
    {
       /* p               1e-6;
        U               1e-6;
        "(omega|f|v2)" 1e-6;
        k               1e-6;
        nuTilda			1e-6;
        epsilon 		1e-6;
        */
    }
}

/*{
    nOuterCorrectors 2;
    nCorrectors     3;
    nNonOrthogonalCorrectors 3;
    pRefCell        0;
    pRefValue       0;

}

PISO
{
    nCorrectors     2;
    nOuterCorrectors 3;
    nNonOrthogonalCorrectors 2;
    pRefCell        0;
    pRefValue       0;

    residualControl
    {
        p               1e-5;
        U               1e-5;
        "(k|epsilon|omega|f|v2)" 1e-5;
    }
}

potentialFlow
{
    nNonOrthogonalCorrectors 0;
}
*/
relaxationFactors
{
    fields
    {
        p               0.9;
        R               0.3;
    }
    equations
    {
        U               0.9;
        k               0.8;
        omega           0.65;
        nuTilda         0.7;
        epsilon         0.7;
        
    }
}

// ************************************************************************* //

[dlahaye]

How many iterations does simpleFoam perform before giving the error message?

[klitz]

Originally Posted by dlahaye

How many iterations does simpleFoam perform before giving the error message?

First iteration U,p,k and epsilon are done. But reynolds stress iteration blows up.

[dlahaye]

Not sure. What about lowering relaxation factor?

[klitz]

Originally Posted by dlahaye

Not sure. What about lowering relaxation factor?

They didnt work unfortunately. Anymore ideas? Because I cannot find more.

[dlahaye]

Pls. elaborate. What happens in case of lowering relaxation factors?

Also, pls. explain the starting guess for the LRR model. Did you run the tow-equation model to full convergence as a starting guess?

Furthermore, do the results for the Reynolds stress model obtained from post-processing the two-equation model look reasonable?

[agustinvo]

Originally Posted by klitz

They didnt work unfortunately. Anymore ideas? Because I cannot find more.

The calculated BC is applied on a wall? Have you tried to use another BC there? The R utility calculates values on the mesh, but at the walls it uses calculated (from nut, k and U fields, I think). A fixed value equal to zero should be a good approach. What are the values you get on this patch?

[crubio.abujas]

Thank you Domenico for the initiative. I also think it is extrange that no example with LRR is included in the tutorials. Indeed, I've found your post looking for a precise way to understand this model and its implementation in OpenFOAM, so I think the post is quite valuable.

However, I've not been able to replicate the results satisfactorily and I found the tutorial quite vague in some specific details that may be critical to successfully complete the simulation.

First of all I want to say I've tried the model in OF7, so many details may be different on other versions, but I think the essentials are preserved. Please correct me if I',m wrong in that assumption.

Stage 1/3:

system/fvSolution:
I've followed the link to the e.g [5], but it is not specifically clear which example to use. The blend case link is dead so no information can be extracted for the original question. Some years later some users post settings for another case "blayer.tar.gz" which is indeed accessible, but as the user received no response to his problem I'm not sure if the values are safe to use.
As you said in the beginning of the post Under relaxation factors are a critical point, so I think a more solid recommendation should be present for the tutorial.
Finally I've used the following values:

fields
{
    p           0.3;
    R           0.3;
}
equations
{
    U           0.6;
    epsilon     0.1;
}

system/fvSchemes:

You mention that settings for wallDist dictionaries are required, but you did not mention any specific setting or how to deal with it. Moreover this dictionary is already defined in the pitzDaily settings so no further modification shall be required, right?
In the tutorials the only method employed is meshWave, although two other models are available (Poisson and advectionDiffusion).

I've used the following command to explore the schemes used in other tutorials.

grep -rh "div(R)" --include fvSchemes $FOAM_TUTORIALS | sort -u
grep -rh "div(phi, R)" --include fvSchemes $FOAM_TUTORIALS | sort -u

For div(R) there is just one option and is Gauss linear.

For div(phi,R) there are 5 options:

• div(phi,R) Gauss limitedLinear 1;
• div(phi,R) Gauss upwind;
• div(phi,R) bounded Gauss limitedLinear 1;
• div(phi,R) bounded Gauss linear;
• div(phi,R) bounded Gauss upwind;

I've tried bounded Gauss linear but was too inestable, so I've decided to use upwind instead.

When I tried to solve the case it complained about lacking a scheme for div((nu*dev2(T(grad(U))))). I don't know if you forgot to mention in the tutorial or was not required for your case. Anyway I used Gauss linear as well for this field.

Stage 2/3:
Here the follow up is quite straight. I was glad with the comment on not changing the RASModel prematurely to avoid problems with the postProcessing comment

State 3/3:
This is also a straightforward step.

The simulation runs 2000 iterations, but the residuals don't look so good. I've tried to extend it to 5000 iterations but I seem to be stuck.

I've checked the velocity profile and it seems like a transient behaviour is occurring. I don't think this is the expected behavior. It may be something about the URF that is limiting it for fully converging or other point that I'm missing, but I really think that a tutorial covering a difficult converging model should have more details concerning its numerics configuration.

It would be great if you can share the specific configuration of the case and some of the results so it can be checked that the model is ok and I'm not missing something essential.

For those who are interested in replicating my exact steps I've written an Allrun file to do so. Just remind you that I'm using OF7, so it may not work in other OF versions.

PD: The forum don't let me attach the Allrun file, so I share it here as a code extract.

Allrun

#!/bin/bash
# Using OpenFOAM v7

# Run from this directory
cd ${0%/*} || exit 1        

echo "Creating case pitzDailyKE"
cp -r $FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily pitzDailyKE

echo "Meshing pitzDailyKE"
blockMesh -case pitzDailyKE > pitzDailyKE/log.blockMesh

echo "Running simpleFoam on pitzDailyKE"
simpleFoam -case pitzDailyKE > pitzDailyKE/log.simpleFoam

echo "Postprocessing pitzDailyKE"
simpleFoam -case pitzDailyKE \
    -postProcess -latestTime -func yPlus > pitzDailyKE/log.post-yPlus
simpleFoam -case pitzDailyKE \
    -postProcess -latestTime -func R > pitzDailyKE/log.post-R

echo "Creating case pitzDailyLRR"
cp -r $FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily pitzDailyLRR

echo "Setting pitzDailyLRR"

echo " * Setting fvSolution"
foamDictionary \
    -entry "solvers.R" -set "{
        solver PBiCG;
        preconditioner DILU;
        tolerante      1e-6;
        relRol         0.1;
    }" pitzDailyLRR/system/fvSolution;

foamDictionary \
    -entry "relaxationFactors" -set "
    {
fields
{
    p           0.3;
    R           0.3;
}
equations
{
    U           0.6;
    epsilon     0.1;
}
    }
    " \
    pitzDailyLRR/system/fvSolution

echo " * Setting fvSchemes"
foamDictionary  \
    -entry "divSchemes.div(R)" -set "Gauss linear" \
    pitzDailyLRR/system/fvSchemes

echo " * Setting fvSchemes"
foamDictionary  \
    -entry "divSchemes.div(R)" -set "Gauss linear" \
    pitzDailyLRR/system/fvSchemes
foamDictionary  \
    -entry "divSchemes.div(phi,R)" -set "bounded Gauss upwind" \
    pitzDailyLRR/system/fvSchemes
foamDictionary  \
    -entry "divSchemes.div((nu*dev2(T(grad(U)))))" -set "Gauss linear" \
    pitzDailyLRR/system/fvSchemes
foamDictionary  \
    -entry "wallDist.method" -set "meshWave" \
    pitzDailyLRR/system/fvSchemes

echo " * Setting initial conditions"
latestTime=$(foamListTimes -case pitzDailyKE -latestTime)
cp -r pitzDailyKE/$latestTime/{U,p,R,epsilon,nut} pitzDailyLRR/0

echo " * Adding more simulation time"
foamDictionary  \
    -entry "endTime" -set "2000" \
    pitzDailyLRR/system/controlDict

echo " * Setting turbulence Model"
foamDictionary \
    -entry "RAS.RASModel" -set "LRR" \
    pitzDailyLRR/constant/turbulenceProperties

echo "Meshing pitzDailyLRR"
blockMesh -case pitzDailyLRR > pitzDailyLRR/log.blockMesh

echo "Solving case pitzDailyLRR"
simpleFoam -case pitzDailyLRR > pitzDailyLRR/log.simpleFoam

[dlahaye]

Dear Carlos,

Thank you much for your wonderful input.

1/ I have tried using OpenFoam-v1906 only. I agree that particular details might be version dependent;

2/ For the underrelaxation, I suggest to use values closer to 0 than to 1 in initial runs. Once the model runs fine, you can increase the underrelaxation factors to increase the speed of convergence. Possibly it is beneficial to lower the overrelazation factor for the velocity from 0.6 to even lower values;

3/ Possibly it is valuable to check whether the model does converge with the k-epsilon model;

Kind wishes, Domenico.

[lpz456]

Thank you very much!

[lpz456]

Why the results are unsteady?

[lpz456]

Originally Posted by crubio.abujas

Thank you Domenico for the initiative. I also think it is extrange that no example with LRR is included in the tutorials. Indeed, I've found your post looking for a precise way to understand this model and its implementation in OpenFOAM, so I think the post is quite valuable.

However, I've not been able to replicate the results satisfactorily and I found the tutorial quite vague in some specific details that may be critical to successfully complete the simulation.

First of all I want to say I've tried the model in OF7, so many details may be different on other versions, but I think the essentials are preserved. Please correct me if I',m wrong in that assumption.

Stage 1/3:

system/fvSolution:
I've followed the link to the e.g [5], but it is not specifically clear which example to use. The blend case link is dead so no information can be extracted for the original question. Some years later some users post settings for another case "blayer.tar.gz" which is indeed accessible, but as the user received no response to his problem I'm not sure if the values are safe to use.
As you said in the beginning of the post Under relaxation factors are a critical point, so I think a more solid recommendation should be present for the tutorial.
Finally I've used the following values:

fields
{
    p           0.3;
    R           0.3;
}
equations
{
    U           0.6;
    epsilon     0.1;
}

system/fvSchemes:

You mention that settings for wallDist dictionaries are required, but you did not mention any specific setting or how to deal with it. Moreover this dictionary is already defined in the pitzDaily settings so no further modification shall be required, right?
In the tutorials the only method employed is meshWave, although two other models are available (Poisson and advectionDiffusion).

I've used the following command to explore the schemes used in other tutorials.

grep -rh "div(R)" --include fvSchemes $FOAM_TUTORIALS | sort -u
grep -rh "div(phi, R)" --include fvSchemes $FOAM_TUTORIALS | sort -u

For div(R) there is just one option and is Gauss linear.

For div(phi,R) there are 5 options:

• div(phi,R) Gauss limitedLinear 1;
• div(phi,R) Gauss upwind;
• div(phi,R) bounded Gauss limitedLinear 1;
• div(phi,R) bounded Gauss linear;
• div(phi,R) bounded Gauss upwind;

I've tried bounded Gauss linear but was too inestable, so I've decided to use upwind instead.

When I tried to solve the case it complained about lacking a scheme for div((nu*dev2(T(grad(U))))). I don't know if you forgot to mention in the tutorial or was not required for your case. Anyway I used Gauss linear as well for this field.

Stage 2/3:
Here the follow up is quite straight. I was glad with the comment on not changing the RASModel prematurely to avoid problems with the postProcessing comment

State 3/3:
This is also a straightforward step.

The simulation runs 2000 iterations, but the residuals don't look so good. I've tried to extend it to 5000 iterations but I seem to be stuck.

I've checked the velocity profile and it seems like a transient behaviour is occurring. I don't think this is the expected behavior. It may be something about the URF that is limiting it for fully converging or other point that I'm missing, but I really think that a tutorial covering a difficult converging model should have more details concerning its numerics configuration.

It would be great if you can share the specific configuration of the case and some of the results so it can be checked that the model is ok and I'm not missing something essential.

For those who are interested in replicating my exact steps I've written an Allrun file to do so. Just remind you that I'm using OF7, so it may not work in other OF versions.

PD: The forum don't let me attach the Allrun file, so I share it here as a code extract.

Allrun

#!/bin/bash
# Using OpenFOAM v7

# Run from this directory
cd ${0%/*} || exit 1        

echo "Creating case pitzDailyKE"
cp -r $FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily pitzDailyKE

echo "Meshing pitzDailyKE"
blockMesh -case pitzDailyKE > pitzDailyKE/log.blockMesh

echo "Running simpleFoam on pitzDailyKE"
simpleFoam -case pitzDailyKE > pitzDailyKE/log.simpleFoam

echo "Postprocessing pitzDailyKE"
simpleFoam -case pitzDailyKE \
    -postProcess -latestTime -func yPlus > pitzDailyKE/log.post-yPlus
simpleFoam -case pitzDailyKE \
    -postProcess -latestTime -func R > pitzDailyKE/log.post-R

echo "Creating case pitzDailyLRR"
cp -r $FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily pitzDailyLRR

echo "Setting pitzDailyLRR"

echo " * Setting fvSolution"
foamDictionary \
    -entry "solvers.R" -set "{
        solver PBiCG;
        preconditioner DILU;
        tolerante      1e-6;
        relRol         0.1;
    }" pitzDailyLRR/system/fvSolution;

foamDictionary \
    -entry "relaxationFactors" -set "
    {
fields
{
    p           0.3;
    R           0.3;
}
equations
{
    U           0.6;
    epsilon     0.1;
}
    }
    " \
    pitzDailyLRR/system/fvSolution

echo " * Setting fvSchemes"
foamDictionary  \
    -entry "divSchemes.div(R)" -set "Gauss linear" \
    pitzDailyLRR/system/fvSchemes

echo " * Setting fvSchemes"
foamDictionary  \
    -entry "divSchemes.div(R)" -set "Gauss linear" \
    pitzDailyLRR/system/fvSchemes
foamDictionary  \
    -entry "divSchemes.div(phi,R)" -set "bounded Gauss upwind" \
    pitzDailyLRR/system/fvSchemes
foamDictionary  \
    -entry "divSchemes.div((nu*dev2(T(grad(U)))))" -set "Gauss linear" \
    pitzDailyLRR/system/fvSchemes
foamDictionary  \
    -entry "wallDist.method" -set "meshWave" \
    pitzDailyLRR/system/fvSchemes

echo " * Setting initial conditions"
latestTime=$(foamListTimes -case pitzDailyKE -latestTime)
cp -r pitzDailyKE/$latestTime/{U,p,R,epsilon,nut} pitzDailyLRR/0

echo " * Adding more simulation time"
foamDictionary  \
    -entry "endTime" -set "2000" \
    pitzDailyLRR/system/controlDict

echo " * Setting turbulence Model"
foamDictionary \
    -entry "RAS.RASModel" -set "LRR" \
    pitzDailyLRR/constant/turbulenceProperties

echo "Meshing pitzDailyLRR"
blockMesh -case pitzDailyLRR > pitzDailyLRR/log.blockMesh

echo "Solving case pitzDailyLRR"
simpleFoam -case pitzDailyLRR > pitzDailyLRR/log.simpleFoam

Hello, I meet the same problem in openfoamv2106, have you solve the problem?
Thank you in advance!