[farzadmech]

Dear Foamers
Does anyone know how to generate Forced synthetic homogeneous isotropic turbulence(HIT) in OpenFOAM ?
I have found one example for decaying turbulence inOF2112;

OF2112/OpenFOAM-v2112/tutorials/incompressible/pimpleFoam/LES/decayIsoTurb

But I want a Force HIT. I will be thankful if anyone introduce me such an implementation.


Thanks,
Farzad

[farzadmech]

I was thinking for a simple method to resolve this complicated problem, and I am thinking about this approach;

    MRF.correctBoundaryVelocity(U);

    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(rho, U)
     ==
        fvOptions(rho, U) + (1-c)*(initial turbulent field) // where c = (current TKE)/(Initial TKE) 
    );

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                    mixture.surfaceTensionForce()
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                ) * mesh.magSf()
            )
        );

        fvOptions.correct(U);
    }

where initial filed is calculated using decayIsoTurb.

I will test this approach and let post it here.


Thanks,
Farzad

[farzadmech]

I think, I am the first person who wants to produce Forced HIT inside OpenFOAM(I have not seen any thread like that before), so I explain my problems here one by one;

1- what is the difference between blockMesh and createBlockMesh?
2- what are the parameters inside createBoxTurb especially the domain dimensions? I changed the Suggested box size L but it does not have any sensible change.

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

N           (64 64 64);
//N           (128 128 128);
//N           (256 256 256);

// Suggested box size of 9*2*pi [cm]
L           (0.56548667765 0.56548667765 0.56548667765); 


nModes      5000;
// Energy as a function of wave number
// Here using Comte-Bellot and Corrsin data at t.U_0/M = 42 (see Ref. table 3)
Ek          table
(
    (15 0)
    (20 0.000129)
    (25 0.00023)
    (30 0.000322)
    (40 0.000435)
    (50 0.000457)
    (70 0.00038)
    (100 0.00027)
    (150 0.000168)
    (200 0.00012)
    (250 8.9e-05)
    (300 7.03e-05)
    (400 4.7e-05)
    (600 2.47e-05)
    (800 1.26e-05)
    (1000 7.42e-06)
    (1250 3.96e-06)
    (1500 2.33e-06)
    (1750 1.34e-06)
    (2000 8e-07)
);


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

3- Can I use createBoxTurb for generating HIT for rectangular box instead of square box?


Thanks,
Farzad

After digging in the code, I am going to reply to my question myself;

1- what is the difference between blockMesh and createBlockMesh?
Answer 1: One of these tools must be used either blockMesh or createBlockMesh. createBlockMesh uses N(Nx,Ny,Nz) and L(Lx,Ly,Lz) to generate mesh.

2- what are the parameters inside createBoxTurb especially the domain dimensions? I changed the Suggested box size L but it does not have any sensible change.
Answer2: N(Nx,Ny,Nz) and L(Lx,Ly,Lz) are geometry and mesh parameters.

N           (64 256 64); //(64 64 64);
//N           (128 128 128);
//N           (256 256 256);

// Suggested box size of 9*2*pi [cm]
L           (0.1 0.4 0.1); //(0.56548667765 0.56548667765 0.56548667765);

3- Can I use createBoxTurb for generating HIT for rectangular box instead of square box?

Answer3: I believe so. First parameters are not restricted to Nx=Ny=Nz and Lx=Ly=Lz. Second, I do one initialization, and this is the quantitative result;

Generating kinetic energy field
min/max/average k = 1.03066e-06, 0.574491, 0.0546993
Generating div(U) field
min/max/average div(U) = -424.15, 363.329, -3.59639e-15

which shows that with rectangular box, it is still divergence free.
Also, I have integrate U components in whole domain, and no mean flow has been generated and it is just pure turbulence;

Ux=-0.00208645
Ux=+0.000625369
Ux=+0.000887475

Only thing remained is input and output powerSpectrum to analyze rectangular box turbulence. For input, I put the original decayIsoTurb tutorial for now;

Ek          table
(
    (15 0)
    (20 0.000129)
    (25 0.00023)
    (30 0.000322)
    (40 0.000435)
    (50 0.000457)
    (70 0.00038)
    (100 0.00027)
    (150 0.000168)
    (200 0.00012)
    (250 8.9e-05)
    (300 7.03e-05)
    (400 4.7e-05)
    (600 2.47e-05)
    (800 1.26e-05)
    (1000 7.42e-06)
    (1250 3.96e-06)
    (1500 2.33e-06)
    (1750 1.34e-06)
    (2000 8e-07)
);

Interestingly, energySpectrum function does not work properly for rectangular box and did not generate output where for the original square box decayIsoTurb energySpectrum works properly.

functions
{
    energySpectrum1
    {
        type            energySpectrum;
        libs            (randomProcessesFunctionObjects);
        writeControl    writeTime;
    }
}

Warning is as below;

--> FOAM Warning : 
energySpectrum1 function object is only appropriate for isotropic structured IJK meshes. Mesh extents: (0.1 0.4 0.1), computed IJK mesh extents: (0.1 0.1 0.1)

    From virtual bool Foam::functionObjects::energySpectrum::read(const Foam::dictionary&)
    in file energySpectrum/energySpectrum.C at line 179.
--> loading function object 'energySpectrum1'

I did several test and results are this:
1- For powerspectrum function, Lx/Nx=Ly/Ny=Lz/Nz is not the sufficient criteria and it need Nx=Ny=Nz(and basically Lx/Nx and Ly/Ny and Lz/Nz can be different), and I tested with Nx=Ny=Nz=64 and it was fine.
2- with the same Lx=Ly=Lz but different Nx=Ny=Nz, I got warning again;

N           (64 128 64); //(64 64 64);
// Suggested box size of 9*2*pi [cm]
L           (0.1 0.1 0.1); //(0.56548667765 0.56548667765 0.56548667765);

--> FOAM Warning : 
energySpectrum1 function object is only appropriate for isotropic structured IJK meshes. Mesh extents: (0.1 0.1 0.1), computed IJK mesh extents: (0.1 0.05 0.1)

    From virtual bool Foam::functionObjects::energySpectrum::read(const Foam::dictionary&)
    in file energySpectrum/energySpectrum.C at line 179.
--> loading function object 'energySpectrum1'

For powerSpectrum, below code is an easy code to implement in the solver;

https://github.com/xuhanwustl/WrayAg...ilities/calcEk

[farzadmech]

createBlockMesh.H

const cellModel& hex = cellModel::ref(cellModel::HEX);

cellShapeList cellShapes;
faceListList boundary;
pointField points;
{
    Info<< "Creating block" << endl;

    block b
    (
        cellShape(hex, identity(8)),
        pointField
        (
            {
                point(0, 0, 0),                   //1
                point(L.x(), 0, 0),               //2
                point(L.x(), L.y(), 0),          //3
                point(0, L.y(), 0),               //4
                point(0, 0, L.z()),               //5
                point(L.x(), 0, L.z()),          //6
                point(L.x(), L.y(), L.z()),     //7
                point(0, L.y(), L.z())           //8
            }
        ),
        blockEdgeList(),
        blockFaceList(),
        N
    );

    Info<< "Creating cells" << endl;
    cellShapes = b.shapes();

    Info<< "Creating boundary faces" << endl;

    boundary.setSize(b.boundaryPatches().size());
    forAll(boundary, patchi)
    {
        faceList faces(b.boundaryPatches()[patchi].size());
        forAll(faces, facei)
        {
            faces[facei] = face(b.boundaryPatches()[patchi][facei]);
        }
        boundary[patchi].transfer(faces);
    }

    points.transfer(const_cast<pointField&>(b.points()));
}

Info<< "Creating patch dictionaries" << endl;
wordList patchNames(boundary.size());
forAll(patchNames, patchi)
{
    patchNames[patchi] = polyPatch::defaultName(patchi);
}

PtrList<dictionary> boundaryDicts(boundary.size());
forAll(boundaryDicts, patchi)
{
    boundaryDicts.set(patchi, new dictionary());
    dictionary& patchDict = boundaryDicts[patchi];
    word nbrPatchName;
    if (patchi % 2 == 0)
    {
        nbrPatchName = polyPatch::defaultName(patchi + 1);
    }
    else
    {
        nbrPatchName = polyPatch::defaultName(patchi - 1);
    }

    patchDict.add("type", cyclicPolyPatch::typeName);
    patchDict.add("neighbourPatch", nbrPatchName);
}

Info<< "Creating polyMesh" << endl;
polyMesh mesh
(
    IOobject
    (
        polyMesh::defaultRegion,
        runTime.constant(),
        runTime,
        IOobject::NO_READ
    ),
    std::move(points),
    cellShapes,
    boundary,
    patchNames,
    boundaryDicts,
    "defaultFaces",
    cyclicPolyPatch::typeName,
    false
);

Info<< "Writing polyMesh" << endl;
mesh.write();

[farzadmech]

createBoxTurb.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2018 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Application
    createBoxTurb

Description
    Create a box of isotropic turbulence based on a user-specified
    energy spectrum.

    Based on the reference
    \verbatim
    Saad, T., Cline, D., Stoll, R., Sutherland, J.C.
    "Scalable Tools for Generating Synthetic Isotropic Turbulence with
    Arbitrary Spectra"
    AIAA Journal, Vol. 55, No. 1 (2017), pp. 327-331.
    \endverbatim

    The \c -createBlockMesh option creates a block mesh and exits, which
    can then be decomposed and the utility run in parallel.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "block.H"
#include "mathematicalConstants.H"

using namespace Foam::constant;

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

Foam::vector randomUnitVector(Random& rndGen)
{
    // Sample point on a sphere
    scalar t = rndGen.globalPosition<scalar>(-1, 1);
    scalar phim = rndGen.globalSample01<scalar>()*mathematical::twoPi;
    scalar thetam = Foam::acos(t);

    return vector
    (
        Foam::sin(thetam*Foam::cos(phim)),
        Foam::sin(thetam*Foam::sin(phim)),
        Foam::cos(thetam)
    );
}


int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Create a box of isotropic turbulence based on a user-specified"
        " energy spectrum."
    );

    argList::addBoolOption
    (
        "createBlockMesh",
        "create the block mesh and exit"
    );

    #include "setRootCase.H"

    #include "createTime.H"
    #include "createFields.H"

    if (args.found("createBlockMesh"))
    {
        // Create a box block mesh with cyclic patches
        #include "createBlockMesh.H"
        Info<< "\nEnd\n" << endl;
        return 0;
    }

    #include "createMesh.H"

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

    // Minimum wave number
    scalar kappa0 = mathematical::twoPi/cmptMin(L);

    // Maximum wave number
    scalar kappaMax = mathematical::pi/cmptMin(delta);

    Info<< "Wave number min/max = " << kappa0 << ", " << kappaMax << endl;

    Info<< "Generating velocity field" << endl;

    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionedVector(dimVelocity, Zero)
    );

    vectorField& Uc = U.primitiveFieldRef();
    const scalar deltaKappa = (kappaMax - kappa0)/scalar(nModes - 1);
    const vectorField& C(mesh.C());
    for (label modei = 1; modei <= nModes; ++modei)
    {
        // Equidistant wave mode
        scalar kappaM = kappa0 + deltaKappa*(modei-1);

        Info<< "Processing mode:" << modei << " kappaM:" << kappaM << endl;

        // Energy
        scalar E = Ek->value(kappaM);

        // Wave amplitude
        scalar qm = Foam::sqrt(E*deltaKappa);

        // Wave number unit vector
        const vector kappaHatm(randomUnitVector(rndGen));

        vector kappaTildem(0.5*kappaM*cmptMultiply(kappaHatm, delta));
        for (direction i = 0; i < 3; ++i)
        {
            kappaTildem[i] = 2/delta[i]*Foam::sin(kappaTildem[i]);
        }

        // Intermediate unit vector zeta
        const vector zetaHatm(randomUnitVector(rndGen));

        // Unit vector sigma
        vector sigmaHatm(zetaHatm^kappaTildem);
        sigmaHatm /= mag(kappaTildem);

        // Phase angle
        scalar psim = 0.5*rndGen.position(-mathematical::pi, mathematical::pi);

        // Add the velocity contribution per mode
        Uc += 2*qm*cos(kappaM*(kappaHatm & C) + psim)*sigmaHatm;
    }

    U.write();

    {
        Info<< "Generating kinetic energy field" << endl;
        volScalarField k("k", 0.5*magSqr(U));
        k.write();
        Info<< "min/max/average k = "
            << gMin(k) << ", " << gMax(k) << ", " << gAverage(k)
            << endl;
    }

    {
       Info<< "Generating div(U) field" << endl;
        volScalarField divU(fvc::div(U));
        divU.write();
        Info<< "min/max/average div(U) = "
            << gMin(divU) << ", " << gMax(divU) << ", " << gAverage(divU)
            << endl;
    }

    Info<< nl;
    runTime.printExecutionTime(Info);

    Info<< "End\n" << endl;
    return 0;
}


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

[farzadmech]

I checked the energySpectrum.C code, and I believe it is able to calculate powerspectrum in a "structured box". Only limitation in the code is that it assumes "// Assume all cells are the same size..." which means no expansion rate can be used. For Safety, it is better to use Lx/Nx=Ly/Ny=Lz/Nz which means mesh cell are cubic and have the same size.

https://www.openfoam.com/documentati...8C_source.html

/*---------------------------------------------------------------------------*\
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
  
     OpenFOAM is free software: you can redistribute it and/or modify it
     under the terms of the GNU General Public License as published by
     the Free Software Foundation, either version 3 of the License, or
     (at your option) any later version.
  
     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
     for more details.
  
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  
 \*---------------------------------------------------------------------------*/
  
 #include "energySpectrum.H"
 #include "fft.H"
 #include "fvMesh.H"
 #include "boundBox.H"
 #include "OFstream.H"
 #include "mathematicalConstants.H"
 #include "volFields.H"
 #include "addToRunTimeSelectionTable.H"
  
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  
 namespace Foam
 {
 namespace functionObjects
 {
     defineTypeNameAndDebug(energySpectrum, 0);
     addToRunTimeSelectionTable(functionObject, energySpectrum, dictionary);
 }
 }
  
  
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
  
 void Foam::functionObjects::energySpectrum::writeFileHeader(Ostream& os)
 {
     writeHeader(os, "Turbulence energy spectra");
  
     writeCommented(os, "kappa E(kappa)");
  
     os  << endl;
 }
  
  
 void Foam::functionObjects::energySpectrum::calcAndWriteSpectrum
 (
     const vectorField& U,
     const vectorField& C,
     const vector& c0,
     const vector& deltaC,
     const Vector<int>& N,
     const scalar kappaNorm
 )
 {
     Log << "Computing FFT" << endl;
     const complexVectorField Uf
     (
         fft::forwardTransform
         (
             ReComplexField(U),
             List<int>({N.x(), N.y(), N.z()})
         )
        /scalar(cmptProduct(N))
     );
  
  
     Log << "Computing wave numbers" << endl;
     const label nMax = cmptMax(N);
     scalarField kappa(nMax);
     forAll(kappa, kappai)
     {
         kappa[kappai] = kappai*kappaNorm;
     }
  
     Log << "Computing energy spectrum" << endl;
     scalarField E(nMax, Zero);
     const scalarField Ec(0.5*magSqr(Uf));
     forAll(C, celli)
     {
         point Cc(C[celli] - c0);
         label i = round((Cc.x() - 0.5*deltaC.x())/(deltaC.x())*(N.x() - 1));
         label j = round((Cc.y() - 0.5*deltaC.y())/(deltaC.y())*(N.y() - 1));
         label k = round((Cc.z() - 0.5*deltaC.z())/(deltaC.z())*(N.z() - 1));
         label kappai = round(Foam::sqrt(scalar(i*i + j*j + k*k)));
  
         E[kappai] += Ec[celli];
     }
  
     E /= kappaNorm;
  
     Log << "Writing spectrum" << endl;
     autoPtr<OFstream> osPtr = createFile(name(), time_.value());
     OFstream& os = osPtr.ref();
     writeFileHeader(os);
  
     forAll(kappa, kappai)
     {
         os  << kappa[kappai] << tab << E[kappai] << nl;
     }
 }
  
  
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  
 Foam::functionObjects::energySpectrum::energySpectrum
 (
     const word& name,
     const Time& runTime,
     const dictionary& dict
 )
 :
     fvMeshFunctionObject(name, runTime, dict),
     writeFile(mesh_, name),
     cellAddr_(mesh_.nCells()),
     UName_("U"),
     N_(Zero),
     c0_(Zero),
     deltaC_(Zero),
     kappaNorm_(0)
 {
     read(dict);
 }
  
  
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  
 bool Foam::functionObjects::energySpectrum::read(const dictionary& dict)
 {
     fvMeshFunctionObject::read(dict);
     writeFile::read(dict);
  
     dict.readIfPresent("U", UName_);
  
     const boundBox meshBb(mesh_.bounds());
  
     // Assume all cells are the same size...
     const cell& c = mesh_.cells()[0];
     boundBox cellBb(boundBox::invertedBox);
     forAll(c, facei)
     {
         const face& f = mesh_.faces()[c[facei]];
         cellBb.add(mesh_.points(), f);
     }
  
     const vector L(meshBb.max() - meshBb.min());
     const vector nCellXYZ(cmptDivide(L, cellBb.max() - cellBb.min()));
  
     N_ = Vector<int>
     (
         round(nCellXYZ.x()),
         round(nCellXYZ.z()),
         round(nCellXYZ.z())
     );
  
     // Check that the mesh is a structured box
     vector cellDx(cellBb.max() - cellBb.min());
     vector expectedMax(N_.x()*cellDx.x(), N_.y()*cellDx.y(), N_.z()*cellDx.z());
     vector relativeSize(cmptDivide(L, expectedMax));
     for (direction i = 0; i < 3; ++i)
     {
         if (mag(relativeSize[i] - 1) > 1e-3)
         {
             FatalErrorInFunction
                 << name() << " function object is only appropriate for "
                 << "isotropic structured IJK meshes. Mesh extents: " << L
                 << ", computed IJK mesh extents: " << expectedMax
                 << exit(FatalError);
         }
     }
     Log << "Mesh extents (deltax,deltay,deltaz): " << L << endl;
     Log << "Number of cells (Nx,Ny,Nz): " << N_ << endl;
  
     // Map into i-j-k co-ordinates
     const vectorField& C = mesh_.C();
     c0_ = returnReduce(min(C), minOp<vector>());
     const vector cMax = returnReduce(max(C), maxOp<vector>());
     deltaC_ = cMax - c0_;
  
     forAll(C, celli)
     {
         label i = round((C[celli].x() - c0_.x())/(deltaC_.x())*(N_.x() - 1));
         label j = round((C[celli].y() - c0_.y())/(deltaC_.y())*(N_.y() - 1));
         label k = round((C[celli].z() - c0_.z())/(deltaC_.z())*(N_.z() - 1));
  
         cellAddr_[celli] = k + j*N_.y() + i*N_.y()*N_.z();
     }
  
     kappaNorm_ = constant::mathematical::twoPi/cmptMax(L);
  
     return true;
 }
  
  
 bool Foam::functionObjects::energySpectrum::execute()
 {
     return true;
 }
  
  
 bool Foam::functionObjects::energySpectrum::write()
 {
     const auto& U = mesh_.lookupObject<volVectorField>(UName_);
     const vectorField& Uc = U.primitiveField();
     const vectorField& C = mesh_.C();
  
     if (Pstream::parRun())
     {
         PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
  
         UOPstream toProc(0, pBufs);
  
         toProc << Uc << C << cellAddr_;
  
         pBufs.finishedSends();
  
         if (Pstream::master())
         {
             const label nGlobalCells(cmptProduct(N_));
             vectorField Uijk(nGlobalCells);
             vectorField Cijk(nGlobalCells);
  
             for (const int proci : Pstream::allProcs())
             {
                 UIPstream fromProc(proci, pBufs);
                 vectorField Up;
                 vectorField Cp;
                 labelList cellAddrp;
                 fromProc >> Up >> Cp >> cellAddrp;
                 UIndirectList<vector>(Uijk, cellAddrp) = Up;
                 UIndirectList<vector>(Cijk, cellAddrp) = Cp;
             }
  
             calcAndWriteSpectrum(Uijk, Cijk, c0_, deltaC_, N_, kappaNorm_);
         }
  
     }
     else
     {
         vectorField Uijk(Uc, cellAddr_);
         vectorField Cijk(C, cellAddr_);
  
         calcAndWriteSpectrum(Uijk, Cijk, c0_, deltaC_, N_, kappaNorm_);
     }
  
     return true;
 }
  
  
 // ************************************************************************* //

[joshwilliams]

Originally Posted by farzadmech

I was thinking for a simple method to resolve this complicated problem, and I am thinking about this approach;

    MRF.correctBoundaryVelocity(U);

    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(rho, U)
     ==
        fvOptions(rho, U) + (1-c)*(initial turbulent field) // where c = (current TKE)/(Initial TKE) 
    );

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                    mixture.surfaceTensionForce()
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                ) * mesh.magSf()
            )
        );

        fvOptions.correct(U);
    }

where initial filed is calculated using decayIsoTurb.

I will test this approach and let post it here.


Thanks,
Farzad

There is a solver called dnsFoam which is for "". Source code can be found here. If I recall though, there is a bug which restricts it to be run in serial (one core only). On OpenFOAM 6 I think the bug was due to the function to calculate FFT not being suitable for parallelisation. Unsure if it has been fixed on later versions.

[farzadmech]

Thanks Josh
I can't use dnsFOAM since I need to simulate multiphase flow using VOF(interFOAM), and basically I am interested in synthetic turbulence instead of DNS approach. I searched a lot, but I did not find anyone who previously used synthetic turbulence inside OpenFOAM framework. What do you think about the approach I have suggested?

fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(rho, U)
     ==
        fvOptions(rho, U) + (1-c)*(initial turbulent field) // where c = (current TKE)/(Initial TKE) 
    );

I think it is the only way I can change decaying turbulence to forced turbulence. Do you know any other approach? or another implemented method inside OpenFOAM?

Thanks,
Farzad

[aghsin]

Hi Farzad,

your question is a little unclear for me. do you need continuously synthetic turbulent content at the inlet of your computational domain? or you need a volumetric turbulence filed?

Regards

[farzadmech]

Dear Aghsin
I have no inlet or outlet. It is a water tank, and I want to generate synthetic isotropic turbulence inside it. I want to use the method I explained in the second post. Do you know any other way to generate forced HIT inside a tank?


Thanks,
Farzad

[farzadmech]

I have searched for original paper used in createBoxTurbDict for powerSpectrum input named "G. Comte-Bellot and S. Corrsin. Simple Eulerian time correlation of full- and narrow-band velocity signals in grid-generated, `isotropic' turbulence. Journal of Fluid Mechanics, 48:273–337, 1971". I have attached table 3 to this post. All input in createBoxTurbDict must be in SI units (see figure1). For small box (0.09*2Pi)*(0.09*2Pi)*(0.09*2Pi) = 0.5654*0.5654*0.5654 = 0.1808 m3, prescribed powerSpectrum and output powerspectrum is very similar as shown in figure2(note: it is decaying turbulence).

[farzadmech]

However, when I use 2Pi*2Pi*2Pi instead of (0.09*2Pi)*(0.09*2Pi)*(0.09*2Pi), my prescribed powerSpectrum is different from output powerspectrum which has been shown in below figure. Reason for this lies in the calculating the kappaNorm_ where smaller box size with the same number of mesh leads to higher wavenumbers.

     kappaNorm_ = constant::mathematical::twoPi/cmptMax(L);

So, in original decayIsoturb where we have used "G. Comte-Bellot and S. Corrsin. Simple Eulerian time correlation of full- and narrow-band velocity signals in grid-generated, `isotropic' turbulence. Journal of Fluid Mechanics, 48:273–337, 1971", createBoxTurb cannot generate the same powerSpectrum when I use 2Pi*2Pi*2Pi instead of (0.09*2Pi)*(0.09*2Pi)*(0.09*2Pi) box. Also, remember Comte-Bellot(1971) data is experimental.

For one extra test I have used DNS data in "Dizaji, Farzad F., and Jeffrey S. Marshall. "An accelerated stochastic vortex structure method for particle collision and agglomeration in homogeneous turbulence." Physics of Fluids 28, no. 11 (2016): 113301." which resulted in good comparison between input and output powerSpectrum since both has used where in both simulation 2Pi*2Pi*2Pi box has been used which has been shown in figure 4.

Originally Posted by farzadmech

I have searched for original paper used in createBoxTurbDict for powerSpectrum input named "G. Comte-Bellot and S. Corrsin. Simple Eulerian time correlation of full- and narrow-band velocity signals in grid-generated, `isotropic' turbulence. Journal of Fluid Mechanics, 48:273–337, 1971". I have attached table 3 to this post. All input in createBoxTurbDict must be in SI units (see figure1). For small box (0.09*2Pi)*(0.09*2Pi)*(0.09*2Pi) = 0.5654*0.5654*0.5654 = 0.1808 m3, prescribed powerSpectrum and output powerspectrum is very similar as shown in figure2(note: it is decaying turbulence).

[joshwilliams]

Originally Posted by farzadmech

However, when I use 2Pi*2Pi*2Pi instead of (0.09*2Pi)*(0.09*2Pi)*(0.09*2Pi), my prescribed powerSpectrum is different from output powerspectrum which has been shown in below figure.

Hi Farzad,


I had this issue a few months ago. You can visualise various spectra using this online tool by University of Utah. For 6.28 x 6.28 x 6.28 domain, I used the Kang-Chester-Meneveau spectra. There is a limitation to the resolution that can be generated using the online tool, for finer resolution, you can use their python script. I would use 10k Fourier modes for accurate representation of the spectra.



J

[joshwilliams]

This is my createBoxTurb file to achieve Re_\lambda = 20 (nu = 0.005 on 128 cell mesh)

[farzadmech]

Thanks a lot Josh
I was in the middle of editing my post. 2Pi*2Pi*2Pi with 128^3 mesh cant generate desired output since its mesh is too large and kappaNorm_ do not let this happen;

   kappaNorm_ = constant::mathematical::twoPi/cmptMax(L);

but for another inputs code works fine. I just have a question, what is lambda in Re_\lambda = 20 ?
Also, do you know can I use createboxTurb to generate turbulence in rectangular box? I have used it to generate and I did not get any error, but powerSpectrum has not been calculated since it is defined for square box and not a rectangular box.

Thanks,
Farzad

Originally Posted by joshwilliams

This is my createBoxTurb file to achieve Re_\lambda = 20 (nu = 0.005 on 128 cell mesh)

[joshwilliams]

Originally Posted by farzadmech

Thanks a lot Josh
I was in the middle of editing my post. 2Pi*2Pi*2Pi with 128^3 mesh cant generate desired output since its mesh is too large and kappaNorm_ do not let this happen;

   kappaNorm_ = constant::mathematical::twoPi/cmptMax(L);

but for another inputs code works fine. I just have a question, what is lambda in

Re_\lambda = 20

?

Thanks,
Farzad

\lambda is Taylor microscale. See equation 5 in the following paper: Elghobashi, S. and Truesdell, G.C., 1992. Direct simulation of particle dispersion in a decaying isotropic turbulence. Journal of Fluid Mechanics, 242, pp.655-700.

[Smyrna]

Farzad,
Would you able to generate the forced HIT?

I was trying to do similar thing (forced HIT to use for multiphase flow) couple of years ago using Saad's tool as well.

I was able to initialize the turbulence but was not able to force it through the entire simulation.

[Alish1984]

Hi,
This is late but for anyone interested...
You can simply use linear forcing methods,
see:

https://pubs.aip.org/aip/pof/article...ren-s-physical

for forcing scalar you can use:
https://link.springer.com/article/10...94-021-00314-6

or the extention of Carrol for scalar.

[ishtaiaq9]

Hi Farzad, were you able to create a working code for the forced isotropic turbulence?