[cyln]

Hello OF users,

I am doing a URANS simulation and trying to understand the difference between div(phi,U) linearUpwind default; and div(phi,U) linearUpwind grad(U);.

linearUpwind is a bounded (with limiters) second order upwind scheme. When I study the OF documentation, I see that linearUpwind is usually used with grad(U) term at the end for the convection term. However, the turbulence terms such as k and omega, which have the same format as the convection term, are used with the default option. So, I did a laminar benchmark study of flat-plate boundary layer against the Blasius solution.

For the convection term, I used the followings, respectively.
1: div(phi,U) Gauss linearUpwind default;
2: div(phi,U) Gauss linearUpwind grad(U);
Both generated almost the same results. I also calculated the shape factor for the two and they are as follows: H_1=2.69, H_2=2.69.

my fvSchemes file looks like this:

 gradSchemes 
 { 
     default         Gauss linear; 
  } 
  divSchemes 
  { 
      default         none; 
      div(phi,U)      Gauss linearUpwind default; 
      //div(phi,U)      Gauss linearUpwind grad(U); 
      div((nuEff*dev2(T(grad(U))))) Gauss linear; 
  }

This suggests that linearUpwind does not necessarily require grad(U) term at the end for the convection term. Remember, for k and omega, Gauss linearUpwind default; is commonly used in openfoam documentations.

My main question, can I use linearUpwind with default option for the convection term? More specifically, what is the difference between using the default and grad(U)? My gradient terms are already discretized using second order centered differencing.

[tomf]

Hi,

  gradSchemes 
  { 
     default         Gauss linear; 
  } 
  divSchemes 
  { 
      default         none; 
      div(phi,U)      Gauss linearUpwind default; 
      //div(phi,U)      Gauss linearUpwind grad(U); 
      div((nuEff*dev2(T(grad(U))))) Gauss linear; 
  }

As you have only a default gradScheme, grad(U) equals default. You may see some different results if you use gradSchemes and divSchemes as:

  gradSchemes 
  { 
     default         Gauss linear; 
     grad(U)        leastSquares;
  } 
  divSchemes 
  { 
      default         none; 
      //div(phi,U)      Gauss linearUpwind default; 
      div(phi,U)      Gauss linearUpwind grad(U); 
      div((nuEff*dev2(T(grad(U))))) Gauss linear; 
  }

Regards,
Tom

[cyln]

Hi Tom,

Thanks for your reply. But what is the default (none)?

Cyln

[tomf]

Hi Cyln,

The default none just means that there is no default divScheme. That means that you would have to provide one for every term that you are solving for. If you comment both

div(phi,U)

lines you should get an error from OpenFOAM that

div(phi,U)

is not defined. If you then change

default none

to

default Gauss linearUpwind default

all terms that are not explicitly specified will use

Gauss linearUpwind default

.

Regards,
Tom

[cyln]

Thanks for the clear explanation. It makes a lot more sense now. What I am trying to really understand is the meaning of "default" added at the end of div(phi,U) discretization scheme.

for example:

div(phi,U) Gauss linearUpwind default;

What does default mean here? Does it mean I am using the second-order upwind as unbounded?

Please correct me if I am wrong, but to my understanding from my readings, when we replace default with grad(U), then linearUpwind becomes bounded scheme through the limiters defined in gradSchemes.

In the case below, grad(U) term is used rather than default. Moreover, there is no limiter defined in gradSchemes as they are second-order centered. Yet, I am not getting exactly the same results for div(phi,U) Gauss linearUpwind default; and div(phi,U) Gauss linearUpwind grad(U);

  gradSchemes 
  { 
     default         Gauss linear; 
  } 
  divSchemes 
  { 
      default         none; 
      //div(phi,U)      Gauss linearUpwind default; 
      div(phi,U)      Gauss linearUpwind grad(U); 
      div((nuEff*dev2(T(grad(U))))) Gauss linear; 
  }

That's why I am seeking the meaning of the default written in div(phi,U) Gauss linearUpwind default; .

[tomf]

I guess you would have to specify what you mean with not exactly the same results, can you quantify?

The default in div(phi,U) linearUpwind default; points to the default gradScheme.

So you should not see any difference. Otherwise the implementation is different from what I said and I would like to know how this differs as well.

Some information can be found on the online cpp-guide.

[cyln]

Regarding the flat-plate boundary layer benchmark, even though I obtain the same shape factor for both div(phi,U) linearUpwind default; and div(phi,U) linearUpwind default;, the shape factors are not exactly equal:

i.e. H_1=2.652, H_2=2.654. Moreover, the velocity values should be extactly the same if default in div(phi,U) linearUpwind default; refers to the default gradScheme which is Gauss linear. However, I am not getting exactly the same results. More specifically, we expect that the following two fvSchemes should give us exactly the same results, but they dont.

  gradSchemes 
  { 
     default         Gauss linear; 
  } 
  divSchemes 
  { 
      default         none; 
      div(phi,U)      Gauss linearUpwind default; 
      div((nuEff*dev2(T(grad(U))))) Gauss linear; 
  }

  gradSchemes 
  { 
     default         Gauss linear; 
  } 
  divSchemes 
  { 
      default         none; 
      div(phi,U)      Gauss linearUpwind grad(U); 
      div((nuEff*dev2(T(grad(U))))) Gauss linear; 
  }

When I try these two in another bluff-body simulation, I get higher numerical diffusion for div(phi,U) Gauss linearUpwind grad(U);, suggesting second-order upwind scheme is bounded in this case.

[tomf]

If you are absolutely sure that this scheme is the only difference between the two cases (setup, software, hardware), than I think there may be something wrong indeed. You could consider filing a bug report. If there is a reply on the bug report, one may know if the implementation is wrong, or if this behavior is as expected and the way I believe it works is wrong.

Of course, if someone sees this thread on the forum and has another explanation, I would be more than interested to know.

[pete20r2]

I just tried the simpleCar tutorial case from 1712 except I changed fvSchemes to:

/*--------------------------------*- C++ -*----------------------------------*\
ddtSchemes
{
    //default         steadyState;
    ddt(epsilon)    steadyState;
    ddt(k)          steadyState;
}

gradSchemes
{
    //default        none;
    grad(U) Gauss linear;
    grad(p) Gauss linear;
    grad(epsilon) Gauss linear;
    grad(k) Gauss linear;
}

divSchemes
{
    //default         none;
    div(phi,U)     Gauss linearUpwind grad(U);
    //div(phi,U)     Gauss linearUpwind grad(p); //can also use this
    div(phi,k)      bounded Gauss upwind;
    div(phi,epsilon) bounded Gauss upwind;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    //default         Gauss linear corrected;
    laplacian(nuEff,U)      Gauss linear corrected;
    laplacian((1|A(U)),p)   Gauss linear corrected;
    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
    laplacian(DkEff,k)      Gauss linear corrected;
}

interpolationSchemes
{
    //default         linear;
    flux(U) linear;
    flux(HbyA)  linear;
}

snGradSchemes
{
    //default         corrected;
}


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

I tested the fact that the grad(U) part is simply to define a gradient scheme by sending it to an irrelevant reference (grad(p)). I get exactly the same results with both, same number of iterations, same residuals, same sensor readings.

Changing div(phi,U) to:

div(phi,U)     Gauss linearUpwind default;

in the above file results in:

--> FOAM FATAL IO ERROR: 
keyword default is undefined in dictionary "/home/xxxxx/OpenFOAM/xxxxx-v1712/run/simpleCar/system/fvSchemes.gradSchemes"

Supplying a default grad scheme:

default   Gauss linear;

allows the case to run again, with exactly the same results.

Unless I've missed something, there is something else you are doing to cause this issue.

When do you stop the jobs?