[CFDanielGER]

Hello everybody,

for a new hybrid turbulence model I want to code a new subgrid length scale delta where delta is defined as follows:

delta = (1 / sqrt(3)) * max | (I_n - I_m) | where n,m = 1,8 and stand for the vertices of a cell.

I_n is defined as follows and is the cross product of the vorticity vector and the vector from the cell center to the cell vertex:

I_n = n x r_n

Now I am looking for an efficient way of defining the new subgrid length scale. Since the vectors of the cell center to the cell vertices are constant I thought about saving them in a list in the createFields part and then call them in the solver. But I don't know how to realize it yet why I am asking you for your help. How do you think could the code be implemented and be structured? I am looking forward for your help. If you need additional information please let me know.

Many thanks and best regards,
Daniel

[CFDanielGER]

Anyone with an idea how to create vectors from cell center to cell vertices?

[mAlletto]

For each cell you can assess the labels of the points building the cell. See e.g. Point labels associated with a cell

The rest is up to you

[CFDanielGER]

My current implementation looks like this but I get some errors when creating the vectors r_n and I_n and r_m / I_m respectively.

template<class BasicTurbulenceModel>
tmp<volScalarField> SAEDES<BasicTurbulenceModel>::deltaOmega
(

) const
{
  volVectorField vorticity = fvc::curl(this->U());
  volScalarField magVorticity = mag(vorticity);
  volVectorField vorticityUnit = vorticity / magVorticity;

  const polyMesh& pMesh = this->mesh();
  const pointField& p = pMesh.points();
  volVectorField cellCenter = this->mesh().C();

  forAll(cellCenter, celli)
  {
    const labelList& cellPts = this->mesh().cellPoints()[celli];

    forAll(cellPts, pointi)
    {
      vector rn[pointi](cellCenter[celli].x() - p[cellPts[pointi]].x(), cellCenter[celli].y() - p[cellPts[pointi]].y(), cellCenter[celli].z() - p[cellPts[pointi]].z());

      volVectorField In[pointi] = (vorticityUnit[celli]  ^ rn[pointi]);

      forAll(cellPts, pointj)
      {
        vector rm[pointj](cellCenter[celli].x() - p[cellPts[pointj]].x(), cellCenter[celli].y() - p[cellPts[pointj]].y(), cellCenter[celli].z() - p[cellPts[pointj]].z());

        volVectorField Im[pointj] = (vorticityUnit[celli]  ^ rn[pointj]);
      }
    }

    volScalarField maxI[celli] = max(mag(In - Im));
  }

  return 1.0 / sqrt(3.0) * maxI;
}

The first loop loops over all cells, the second loop loops over all vertices for each cell and creates the vectors r_n (vector from the cell center to the cell vertices) and I_n (cross product of r_n with the vorticity vector). The third loop loops again over all vertices of each cell and creates r_m and I_m respectively. At the end of the first loop then the maximum of the difference between I_n and I_m is calculated.

I get following errors when creating I_m:
- array must be initialized with a brace-enclosed initializer

And some more errors when creating r_m and r_n respectively:
- array must be initialized with a brace-enclosed initializer
- expression list treated as compound expression in initializer [-fpermissive]
- value computed is not used [-Wunused-value]

Do you have an idea what the problem could be?

Many thanks for your help.

Best regards,
Daniel

[mAlletto]

Try to Google for the constructors of a vector and vectorfield in openfoam.

[CFDanielGER]

Thanks for your help, it is working now!

[mAlletto]

Can you post your solution here to help others who have the same problem

[CFDanielGER]

Sure, maybe not the "best" implementation but it is working

template<class BasicTurbulenceModel>
tmp<volScalarField> kOSSTEDES<BasicTurbulenceModel>::deltaOmega
(

) const
{
  volVectorField vorticity = fvc::curl(this->U());
  volScalarField magVorticity = mag(vorticity);
  volVectorField vorticityUnit = vorticity / magVorticity;

  const polyMesh& pMesh = this->mesh();
  const pointField& p = pMesh.points();
  volVectorField cellCenter = this->mesh().C();

  List<vector> In(8, vector::zero);
  List<vector> Im(8, vector::zero);

  volScalarField maxI
  (
    IOobject
    (
      "magI",
      this->mesh().time().timeName(),
      this->mesh(),
      IOobject::NO_READ,
      IOobject::NO_WRITE
    ),
    this->mesh(),
    dimensionedScalar("magI", dimensionSet(0, 1, 0, 0, 0, 0, 0), 0)
  );

  scalar previousMaxI = 0.0;

  forAll(cellCenter, celli)
  {
    const labelList& cellPts = this->mesh().cellPoints()[celli];

    forAll(cellPts, pointi)
    {
      label cellP = cellPts[pointi];
      vector pnt = p[cellP];
      vector rn(cellCenter[celli] - pnt);

      forAll(cellPts, pointj)
      {
        label cellP = cellPts[pointj];
        vector pnt = p[cellP];
        vector rm(cellCenter[celli] - pnt);

        In[pointj] = vorticityUnit[celli]  ^ rn;
        Im[pointj] = vorticityUnit[celli]  ^ rm;
      }

      maxI[celli] = max(max(mag(In - Im)), previousMaxI);

      previousMaxI = maxI[celli];

    }

  }

  return 1.0 / sqrt(3.0) * maxI;