CFD Online Discussion Forums - Solving N (Non-Coupled) Scalar Transport Equations

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Solving N (Non-Coupled) Scalar Transport Equations (https://www.cfd-online.com/Forums/openfoam-programming-development/90374-solving-n-non-coupled-scalar-transport-equations.html)

Solving N (Non-Coupled) Scalar Transport Equations

Hi,

I need to create a solver based on buoyantPimpleFoam which solves an arbitrary (around 10) scalar transport equations. I have written a solver which solves one passive scalar transport equation. I could just copy-paste that code enough many times, but I would like to know if a more elegant solution is possible: Namely, how easy it would be to implement a for-loop that solves N scalar transport equations, where N is specified by the user? The problem is, I don't know C++ and I'm not too familiar with the classes present in OpenFOAM.

As far as I know, the vectors in OpenFOAM are always 3-dimensional. But could I use a tensor of rank (1,N) or (N,1) to store the scalars? Or is using a tensor a bad idea, if I want to have different BCs for the scalars? Tensor would be nice because then I wouldn't have to worry about how to get OF to read/write files with indices in their names (scalar1, scalar2, etc.).

This is not a critical problem, as I said I can just copy-paste the code snippet that solves one scalar transport equation N times, but if I want to change the number of scalars later it becomes cumbersome to recompile the solver every time. I was just wondering if anyone with more knowledge of C++ and OF source code knows these things.

Regards,
Joel

Sure, you can use a PtrList. Here's how you initialise them:

Code:

//  Create a list of pointers to the mass fraction field of each species.

    PtrList<volScalarField> chemicalSpecies(numberOfSpecies);

    for (label i=0; i<numberOfSpecies; ++i)

    {

        Info << "Creating Species " << namesOfSpecies[i] << endl;

        chemicalSpecies.set

        (

            i,

            new volScalarField

            (

                IOobject

                (

                    namesOfSpecies[i],

                    runTime.timeName(),

                    mesh,

                    IOobject::NO_READ,

                    IOobject::AUTO_WRITE

                ),

                mesh,

                dimensionedScalar("zero", dimless, 0.0),

                mixtureFraction.boundaryField().types()

            )

        );

    }

Then use a for loop to solve the equations, as done in, for example, YEqn.H in solvers/combustion/chemFoam

Thank you very much for your help, I managed to create the solver I will need.

I am also trying to solve a similar problem of solving transport equation of n scalars, say c, and would like to define them as an array, such as c[n] instead of defining n separate scalarfields. The above thread was of some help but if you can tell me which solver it is a part of then I could look into it to clarify my doubts. The main problem I face is with the definition of c as an array in the initial time folder [0].
Regards

Hi,
I'm asking myself the same question (how to initialize the array in the 0 folder).
Did you make any progress?

Marc

hello,

Why you don't use reactingFoam or rhoReacting(Buoyant)Foam ?
At least take a look at this solver.

regards,
olivier

Thanks, reactingFoam is helpful. But I guess there is no easier way than duplicating/renaming the files in the 0 Folder for each species. But this at least can be automated with a script.

Cheers
Marc

yes and no ;)

If you take a look at reactingFoam tutorial, for product of reaction, you have a "Ydefault" file, in which you specify the defaut boundary conditions for the species. So you may try to use this.

regards,
olivier

I took a look at the reactingFoam solver. Since I'm only interested in a passive scalar transport, having a whole thermodynamic model might be a bit overkill... I'm not sure how I could only extract the classes I need, since it all seems to dependend on the classes basicMultiComponentMixture and chemistryReader.

Regards
Marc

Hi,

I managed to get it to work. I used parts of the reactingFoam for the loop structure,

Code:

tmp<fv::convectionScheme<scalar> > mvConvection

        (

            fv::convectionScheme<scalar>::New

            (

                mesh,

                fields,

                phi,

                mesh.divScheme("div(phi,si_h)")

            )

        );

        volScalarField kappaEff

        (

            "kappaEff",

            turbulence->nu()/Pr + turbulence->nut()/Prt

        );



        forAll(s, i)

            {

            volScalarField& si = s[i];



            tmp<fvScalarMatrix> siEqn

            (

                    fvm::ddt(si)

                  + mvConvection->fvmDiv(phi, si)

                  - fvm::laplacian(kappaEff, si)

            );



            sources.constrain(siEqn());

            solve(siEqn() == sources(si),mesh.solver("si"));

        }

and copy pasted parts of the chemistryReader to get the initial files.
To define the species, I added a list into the transportProperties dict. (makes it easier to read in the createFields.H

Code:

wordList names(transportProperties.lookup("scalars"));

    PtrList<volScalarField> s(names.size());



    forAll(s, i)

    {

        IOobject header

        (

            names[i],

            mesh.time().timeName(),

            mesh,

            IOobject::NO_READ

        );



        // check if field exists and can be read

        if (header.headerOk())

        {

            s.set

            (

                i,

                new volScalarField

                (

                    IOobject

                    (

                        names[i],

                        mesh.time().timeName(),

                        mesh,

                        IOobject::MUST_READ,

                        IOobject::AUTO_WRITE

                    ),

                    mesh

                )

            );

        }

        else

        {

            volScalarField sdefault

            (

                IOobject

                (

                    "sdefault",

                    mesh.time().timeName(),

                    mesh,

                    IOobject::MUST_READ,

                    IOobject::NO_WRITE

                ),

                mesh

            );



            s.set

            (

                i,

                new volScalarField

                (

                    IOobject

                    (

                        names[i],

                        mesh.time().timeName(),

                        mesh,

                        IOobject::NO_READ,

                        IOobject::AUTO_WRITE

                    ),

                    sdefault

                )

            );

        }

    }

    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;



    forAll(s, i)

        {

            fields.add(s[i]);

        }