CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Solving N (Non-Coupled) Scalar Transport Equations (http://www.cfd-online.com/Forums/openfoam-programming-development/90374-solving-n-non-coupled-scalar-transport-equations.html)

joel.lehikoinen July 8, 2011 03:20

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

l_r_mcglashan July 8, 2011 04:00

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

joel.lehikoinen July 12, 2011 02:41

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

dinksy February 19, 2013 03:56

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

marc.immer October 1, 2013 05:47

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

Marc

olivierG October 1, 2013 08:36

hello,

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

regards,
olivier

marc.immer October 1, 2013 08:44

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

olivierG October 1, 2013 08:50

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

marc.immer October 1, 2013 10:16

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

marc.immer October 3, 2013 08:00

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]);
        }



All times are GMT -4. The time now is 17:58.