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)

 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 23:46.