CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   How to address a certain region in solver? (https://www.cfd-online.com/Forums/openfoam-programming-development/185951-how-address-certain-region-solver.html)

mizo April 6, 2017 09:32

How to address a certain region in solver?
 
Hello,

say I create a simple mesh with blockMesh. Then I use topoSet and splitMeshRegions to create two regions Reg1 and Reg2. Now, I would like to calculate a certain parameter, but differently for each region. Is there a way to "tell" the solver that for region Reg1 he should use a certain equation and for region Reg2 he should use a different one?

I know that multiple regions are used in chtMultiRegionFoam tutorial. The solver handles it in a very sophisticated way:
Code:

forAll(fluidRegions, i)
        { ...

I have tried to understand this tutorial and to use something similar in my case but sadly, I have not been successful. I do not have C++ language knowledge, but it is my impression that this case is connected to the heat transfer case very strongly. I have not been able to understand, how exactly should I create fields. The solver does it via #include "createFields.H", which includes the creation of fluid and solid fields.

The createFluidFields.H file looks like this:

Code:

// Initialise fluid field pointer lists
PtrList<rhoThermo> thermoFluid(fluidRegions.size());
PtrList<volScalarField> rhoFluid(fluidRegions.size());
PtrList<volVectorField> UFluid(fluidRegions.size());
PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
PtrList<uniformDimensionedScalarField> hRefFluid(fluidRegions.size());
PtrList<volScalarField> ghFluid(fluidRegions.size());
PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
PtrList<volScalarField> p_rghFluid(fluidRegions.size());
PtrList<radiation::radiationModel> radiation(fluidRegions.size());

List<scalar> initialMassFluid(fluidRegions.size());
List<label> pRefCellFluid(fluidRegions.size(), 0);
List<scalar> pRefValueFluid(fluidRegions.size(), 0.0);
List<bool> frozenFlowFluid(fluidRegions.size(), false);

PtrList<dimensionedScalar> rhoMax(fluidRegions.size());
PtrList<dimensionedScalar> rhoMin(fluidRegions.size());

PtrList<IOMRFZoneList> MRFfluid(fluidRegions.size());
PtrList<fv::options> fluidFvOptions(fluidRegions.size());

// Populate fluid field pointer lists
forAll(fluidRegions, i)
{
    Info<< "*** Reading fluid mesh thermophysical properties for region "
        << fluidRegions[i].name() << nl << endl;

    Info<< "    Adding to thermoFluid\n" << endl;

    thermoFluid.set
    (
        i,
        rhoThermo::New(fluidRegions[i]).ptr()
    );

    Info<< "    Adding to rhoFluid\n" << endl;
    rhoFluid.set
    (
        i,
        new volScalarField
        (
            IOobject
            (
                "rho",
                runTime.timeName(),
                fluidRegions[i],
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            thermoFluid[i].rho()
        )
    );

    Info<< "    Adding to UFluid\n" << endl;
    UFluid.set
    (
        i,
        new volVectorField
        (
            IOobject
            (
                "U",
                runTime.timeName(),
                fluidRegions[i],
                IOobject::MUST_READ,
                IOobject::AUTO_WRITE
            ),
            fluidRegions[i]
        )
    );

    Info<< "    Adding to phiFluid\n" << endl;
    phiFluid.set
    (
        i,
        new surfaceScalarField
        (
            IOobject
            (
                "phi",
                runTime.timeName(),
                fluidRegions[i],
                IOobject::READ_IF_PRESENT,
                IOobject::AUTO_WRITE
            ),
            linearInterpolate(rhoFluid[i]*UFluid[i])
                & fluidRegions[i].Sf()
        )
    );

    Info<< "    Adding to gFluid\n" << endl;
    gFluid.set
    (
        i,
        new uniformDimensionedVectorField
        (
            IOobject
            (
                "g",
                runTime.constant(),
                fluidRegions[i],
                IOobject::MUST_READ,
                IOobject::NO_WRITE
            )
        )
    );

    Info<< "    Adding to hRefFluid\n" << endl;
    hRefFluid.set
    (
        i,
        new uniformDimensionedScalarField
        (
            IOobject
            (
                "hRef",
                runTime.constant(),
                fluidRegions[i],
                IOobject::READ_IF_PRESENT,
                IOobject::NO_WRITE
            ),
            dimensionedScalar("hRef", dimLength, 0)
        )
    );

    dimensionedScalar ghRef
    (
        mag(gFluid[i].value()) > SMALL
      ? gFluid[i]
          & (cmptMag(gFluid[i].value())/mag(gFluid[i].value()))*hRefFluid[i]
      : dimensionedScalar("ghRef", gFluid[i].dimensions()*dimLength, 0)
    );

    Info<< "    Adding to ghFluid\n" << endl;
    ghFluid.set
    (
        i,
        new volScalarField
        (
            "gh",
            (gFluid[i] & fluidRegions[i].C()) - ghRef
        )
    );

    Info<< "    Adding to ghfFluid\n" << endl;
    ghfFluid.set
    (
        i,
        new surfaceScalarField
        (
            "ghf",
            (gFluid[i] & fluidRegions[i].Cf()) - ghRef
        )
    );

    Info<< "    Adding to turbulence\n" << endl;
    turbulence.set
    (
        i,
        compressible::turbulenceModel::New
        (
            rhoFluid[i],
            UFluid[i],
            phiFluid[i],
            thermoFluid[i]
        ).ptr()
    );

    p_rghFluid.set
    (
        i,
        new volScalarField
        (
            IOobject
            (
                "p_rgh",
                runTime.timeName(),
                fluidRegions[i],
                IOobject::MUST_READ,
                IOobject::AUTO_WRITE
            ),
            fluidRegions[i]
        )
    );

    // Force p_rgh to be consistent with p
    p_rghFluid[i] = thermoFluid[i].p() - rhoFluid[i]*ghFluid[i];

    fluidRegions[i].setFluxRequired(p_rghFluid[i].name());

    radiation.set
    (
        i,
        radiation::radiationModel::New(thermoFluid[i].T())
    );

    initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();

    const dictionary& simpleDict =
        fluidRegions[i].solutionDict().subDict("SIMPLE");

    setRefCell
    (
        thermoFluid[i].p(),
        p_rghFluid[i],
        simpleDict,
        pRefCellFluid[i],
        pRefValueFluid[i]
    );

    simpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);

    rhoMax.set
    (
        i,
        new dimensionedScalar
        (
            dimensionedScalar::lookupOrDefault
            (
                "rhoMax",
                simpleDict,
                dimDensity,
                GREAT
            )
        )
    );

    rhoMin.set
    (
        i,
        new dimensionedScalar
        (
            dimensionedScalar::lookupOrDefault
            (
                "rhoMin",
                simpleDict,
                dimDensity,
                0
            )
        )
    );

    Info<< "    Adding MRF\n" << endl;
    MRFfluid.set
    (
        i,
        new IOMRFZoneList(fluidRegions[i])
    );

    Info<< "    Adding fvOptions\n" << endl;
    fluidFvOptions.set
    (
        i,
        new fv::options(fluidRegions[i])
    );

    turbulence[i].validate();
}

Right at the beginning, all the lists are created using functions pretty specific for the heat transfer case. I do not know how could I include something similar in my case.

My question is, is there a way around this? My case is much simpler, I just need two, maximum three, regions. Is there a way I could just assign equations to a specific region?

Thanks for reading!
mizo

p.b April 19, 2017 21:06

OpenFOAM uses a hierarchical system in terms of object registries and its registered (named) objects, where each object is associated with a case directory:
https://openfoamwiki.net/index.php/O...objectRegistry
http://foam.sourceforge.net/docs/Gui...df#section.4.1
https://github.com/OpenFOAM/OpenFOAM...jectRegistry.H
https://github.com/OpenFOAM/OpenFOAM.../regIOobject.H

In short and simplified words:
- The topmost registry is runTime:
https://github.com/OpenFOAM/OpenFOAM...e/createTime.H
Code:

Foam::Info<< "Create time\n" << Foam::endl;

Foam::Time runTime(Foam::Time::controlDictName, args);

- A new mesh is normally registered as object of runTime (third parameter of IOobject below):
https://github.com/OpenFOAM/OpenFOAM...e/createMesh.H
Code:

Foam::Info
    << "Create mesh for time = "
    << runTime.timeName() << Foam::nl << Foam::endl;

Foam::fvMesh mesh
(
    Foam::IOobject
    (
        Foam::fvMesh::defaultRegion,
        runTime.timeName(),
        runTime,
        Foam::IOobject::MUST_READ
    )
);

- if there is more than one mesh, every mesh simply gets its own name from a list of solidNames, fluidNames or whateverNames (first parameter of IOobject below):
https://github.com/OpenFOAM/OpenFOAM...eFluidMeshes.H
Code:

const wordList fluidNames(rp["fluid"]);

PtrList<fvMesh> fluidRegions(fluidNames.size());

forAll(fluidNames, i)
{
    Info<< "Create fluid mesh for region " << fluidNames[i]
        << " for time = " << runTime.timeName() << nl << endl;

    fluidRegions.set
    (
        i,
        new fvMesh
        (
            IOobject
            (
                fluidNames[i],
                runTime.timeName(),
                runTime,
                IOobject::MUST_READ
            )
        )
    );
}

https://github.com/OpenFOAM/OpenFOAM...eSolidMeshes.H
Code:

const wordList solidsNames(rp["solid"]);

PtrList<fvMesh> solidRegions(solidsNames.size());

forAll(solidsNames, i)
{
    Info<< "Create solid mesh for region " << solidsNames[i]
        << " for time = " << runTime.timeName() << nl << endl;

    solidRegions.set
    (
        i,
        new fvMesh
        (
            IOobject
            (
                solidsNames[i],
                runTime.timeName(),
                runTime,
                IOobject::MUST_READ
            )
        )
    );
}

- every field is registered as an object (third parameter of IOobject below) of ONE selected mesh, e.g.:
https://github.com/OpenFOAM/OpenFOAM...createFields.H
Code:

Info<< "Reading field p\n" << endl;
volScalarField p
(
    IOobject
    (
        "p",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);


Info<< "Reading field U\n" << endl;
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

- every equation "belongs" to mesh the first field which is used in the context of a finite volume method (e.g. fvm::ddt, fvm::laplacian, ...), e.g.:
https://github.com/OpenFOAM/OpenFOAM...Foam/icoFoam.C
Code:

fvVectorMatrix UEqn
(
    fvm::ddt(U)
  + fvm::div(phi, U)
  - fvm::laplacian(nu, U)
);

So you simply have to do something like:
Code:

Time runTime(Time::controlDictName, args);

fvMesh meshA
(
    Foam::IOobject
    (
        word("regionA"),
        runTime.timeName(),
        runTime,
        IOobject::MUST_READ
    )
);

fvMesh meshB
(
    Foam::IOobject
    (
        word("regionB"),
        runTime.timeName(),
        runTime,
        IOobject::MUST_READ
    )
);

volVectorField fieldA
(
    IOobject
    (
        "fieldA",
        runTime.timeName(),
        meshA,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    meshA
);

volVectorField fieldB
(
    IOobject
    (
        "fieldB",
        runTime.timeName(),
        meshB,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    meshB
);

fvVectorMatrix fieldAinRegionAEqn
(
    fvm::laplacian(fieldA)
);

fvVectorMatrix fieldBinRegionBEqn
(
    fvm::laplacian(fieldB)
);

- this is linked to a directory structure like
Code:

0
- regionA
  - fieldA
- regionB
  - fieldB
1
- regionA
  - fieldA
- regionB
  - fieldB
constant
- regionA
  - polyMesh
    - points
    - faces
    - owner
    - neighbour
    - boundary
- regionB
  - polyMesh
    - points
    - faces
    - owner
    - neighbour
    - boundary
system
- controlDict
- regionA
  - fvSolution
  - fvSchemes
- regionB
  - fvSolution
  - fvSchemes

- I suggest you make yourself familiar with C++ and refer to the following guides for further reading:
http://foam.sourceforge.net/docs/Gui...erGuide-A4.pdf
http://foam.sourceforge.net/docs/Gui...mmersGuide.pdf

Best, Pascal

mizo April 20, 2017 05:25

Thank you very much for the elaborate answer, it will surely help!

xiaoliang_CFD September 11, 2019 12:03

This is really helpful for a new OpenFoam user. Thanks a lot!


All times are GMT -4. The time now is 09:56.