[JoaoDMiranda]

Hello everybody,

I am trying to put a new transport equation to see the mixture of flow (non-reactive) in rhoCentralFoam. I found a very simple model, http://openfoamwiki.net/index.php/Ho...sport_equation, but this does not work for rhoCentralFoam.

Has anyone implemented a transport equation for the purpouse of mixture in rhoCentralFoam?
Can you give me some advice on how to do it?

Thanks in advance for your help.

João

[treima]

Hello,

I didn´t add a equation for mixture in rhoCentralFoam, but I´ve worked with this solver for some time.

In rhoCentralFoam you need a special discretisation of the equation, not the "simple one" which is described in the wiki. You can have a look at the paper form Greenshield, Weller, Gasparini and Reese,
"Implementation of semi-discrete, non-staggered central schemes in a colocated, polyhedral, finite volume framework, for high-speed viscous flows".
Perhaps this helps you. The discretisation is based on an upwind sheme form Kurganov, Noelle and Petrova.


regards

treima

[JoaoDMiranda]

Thanks a lot! I will look into it!

[Bryan S]

Did you figure out a discretization of the transport equation that works with rhoCentralFoam? I'm trying to do the same thing. Thanks.

[Henning86]

i solved the problem for version 2.1.1 but it wont work on the new version.


i justed link the psireactionthermo libary in rhoCentralFoam. (see rhoreactionfoam for details)

this should be pretty easy in v2.2.x

[immortality]

Hi Henning,
I before added a transport equation to rhoCentralFoam that worked well but not a mixture one,could you tell us what equation you used and how?could you put the equation in this thread? that could be very helpful to OpenFOAMers.

[Henning86]

This Implementation should work for v2.2.2 it was easier than i tought

in version 2.1.1 you have to create an new class

[Bryan S]

Thank you Henning, those are very helpful! I notice that the mass fraction of species 2 in the forward step case is greater than 1 in some places in the domain. Do you know the cause of this issue and if it can be fixed? Thanks!
-B

[Bryan S]

Actually never mind. I'm an experimentalist so I didn't realize that this can happen with linear equations like the scalar transport equation in certain numerical schemes. I had a colleague explain it to me and it all makes sense now. Thanks for your help!
-B

[ChrisA]

Henning,

Your Yeqn.H implementation is missing the central upwind scheme, you're going to notice some strange behavior in your solutions as a result. I've been using the following yeqn which works quite well (you're going to have to change the diffusion coefficient bits back to Schmidt number as this is developed for laminar diffusion problems.)

It's a bit of a mess as I've just finished developing it to fit my purposes and cut a good chunk of that stuff out before posting it... I haven't had time to clean it up as well, as a result there's some spurious comments that should be deleted/ignored. Credit for the original upwind scheme implementation to the yeqn should go to tatu, who provided it to me.

[ChrisA]

Originally Posted by Bryan S

Thank you Henning, those are very helpful! I notice that the mass fraction of species 2 in the forward step case is greater than 1 in some places in the domain. Do you know the cause of this issue and if it can be fixed? Thanks!
-B

This is an issue with the method being used. Mass conservation isn't satisfied with the way the Y equation is written. There are a couple solutions to this which are available in many combustion textbooks but the approach reactingfoam takes is to use an inert species, N2 for example, and have any errors be absorbed by that. This is done in the following lines:

     
Yt += Yi; (inside the loop species loop)

Y[inertIndex] = scalar(1) - Yt;
Y[inertIndex].max(0.0);

Basically you force conservation by taking Yinert = 1 - all other species, the error is then "absorbed" (mostly) by the inert species concentration. This is applicable for certain situations, others you will need to implement a proper transport equation formulation that satisfies mass conservation without this trick.

[Bryan S]

Thanks Chris, I was wondering about that since the previous code submitted used the same discretization I was trying before with no success.

Could you please post a .zip file with all the files required for your solver? I am very interested in using it. I tried simply inserting your .H file into Henning's solver directory and recompiling, but your file includes a DijEqn.H file that I don't have. Thanks!
-B

[Bryan S]

Another issue that comes up with this is that I'd like to use a thermotype with hePsiThermo, multiComponentMixture, sutherlandTransport, hConstThermo, perfectGas, specie, and sensibleEnthalpy. This isn't one of the default supported options and there is a thread somewhere about how to create a new set of options and link that library to the solver, but I'm not clear on how to link libraries and compile a new thermo model etc. because as I mentioned I'm an experimentalist so I'm new to C++ and don't really know how to do this sort of thing.

Do any of you know of a nice set of step-by-step instructions for creating and compiling my own thermo model and then linking it to my solver? Thanks!
-B

[ChrisA]

Originally Posted by Bryan S

Thanks Chris, I was wondering about that since the previous code submitted used the same discretization I was trying before with no success.

Could you please post a .zip file with all the files required for your solver? I am very interested in using it. I tried simply inserting your .H file into Henning's solver directory and recompiling, but your file includes a DijEqn.H file that I don't have. Thanks!
-B

Ah, yes, Dijeqn.h was for calculating the laminar binary diffusion coefficient for a case I was running with binary laminar diffusion, not terribly interesting and not really ready for any sort of public viewing. Also, my full solver will likely be much more of a hindrance than a help in its current form, its very much coded for a niche situation currently.

I changed the original file I uploaded to use the standard turbulent Schmidt number formulation. That should work when combined with the other files hennings provided... I think, completely untested though (Sorry I haven't been able to go through it more thoroughly, busy busy)

[Bryan S]

Originally Posted by ChrisA

I changed the original file I uploaded to use the standard turbulent Schmidt number formulation. That should work when combined with the other files hennings provided... I think, completely untested though

It almost works with Henning's files, I think you must be calling some functions that aren't linked in the Make/options file or something. When I try to wmake I get the following errors:

In file included from rhoCentralReactionFoam.C:247:0:
YEqn.H: In function â:
YEqn.H:9:9: error: â was not declared in this scope
YEqn.H:24:29: error: â was not declared in this scope
YEqn.H:26:35: error: â was not declared in this scope
YEqn.H:62:16: error: â was not declared in this scope

I'm not really sure what those mean, I tried looking at that file using a couple different editors but I can't figure out what function or functions it's talking about and if I knew I don't know how to fix the problem. Any advice? Thanks,
-B

[ChrisA]

That's and interesting and slightly uninformative error message that i've never seen before... not sure why it's not giving you actual functions and variable names. Anyway, based on the lines it trips up at it's missing a few declarations related to the multispecies bits. I took a quick look at henning's createFields.H, it looks like you're missing the following lines of code:

multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;

forAll(Y, i)
{
    fields.add(Y[i]);
}

That should do it. Probably wouldn't have been an issue had I given you all my files but we can file this one under teach a man to fish and whatnot

Also, you can look at at reactingfoam and see what files it has linked, what initial headers are included at the top of the code and any general yeqn code differences. That's what the yeqn is based on with the addition of the central-upwind scheme so any additional headers and linked libraries that are required can be found there.

[Bryan S]

Thanks Chris, both tips were very helpful. The error message is frustratingly vague, and it turns out it's referring to several different undeclared functions but calling them all â. So that's fun. I've fixed a couple of them by adding lines to createFields.H based on what I saw in rhoReactingFoam.

I've gone down from 4 errors to 2, and they may be coming from the same thing. I'm pretty sure it involves the variable rhoY in YEqn.H, which I assume is a volScalarField that is just rho*Y. I tried declaring it like this:

volScalarField rhoY
(
IOobject
(
"rhoY",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
rho*Y
);

but I get an error that says "invalid initialization of reference of type â from expression of type â" which I assume means that rhoY is supposed to be a list with 2 elements and I'm treating it like a single scalar value. How do I declare it properly? Thanks!
-B

[ChrisA]

Yeah, somethings wrong with your compiler's error messages, I've never ever seen that happen; and I've seen my fair share of errors :P.

So, the error here comes in because Y is a list (or a list of pointers? this is where my knowledge gets vague) containing mass fraction information for each individual species Y[i]. So Y is sort of a 4D array if you want to think of it like that, its a list of all the 3D arrays of species mass fraction. Anyway, rhoY needs to be a list as well, containing rho*Y[i], the following is the proper declaration for createFields.H:

multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;

forAll (Y, i)
{
    fields.add(Y[i]);
}
PtrList<volScalarField> rhoY(Y.size());
forAll (Y, i)
{
         rhoY.set
         (
             i,
             new volScalarField
             (
                 IOobject
                 (
                     "rhoYi",
                     runTime.timeName(),
		     mesh,
                     IOobject::NO_READ,
                     IOobject::NO_WRITE
                 ),
                 rho*Y[i]
             )
         );
    	 fields.add(Y[i]);
}

Pretty sure the fields.add(Y[i]) is redundant but... yeah.. :P. Sorry for forgetting that part, that's what I get for not looking at my own code to fix your issue. But, now you know how to add lists of other things which comes in very handy when dealing with multispecies things!

[Bryan S]

Thanks Chris, that fixed it. Turns out the last error is independent from that one, it comes from using "reaction" in YEqn.H, which isn't declared in Henning's createFields.H file. I tried copying the lines from rhoReactingFoam and reactingFoam's createFields.H files, neither seems to work. I also made sure to include the necessary libraries and files in my Make/options file, but that's not helping. I realize I could just take the reactions bit out of the equation YEqn.H because I don't have a source term for the application I have in mind, but it might be nice to have that utility later on.

What do the first few lines of your createFields file look like? At the moment I have

autoPtr<combustionModels::rhoCombustionModel> reaction
(
combustionModels::rhoCombustionModel::New(mesh)
);

rhoReactionThermo& thermo = reaction->thermo();
thermo.validate(args.executable(), "h", "e");

My Make/options file includes both -I$(LIB_SRC)/combustionModels/lnInclude and -lcombustionModels

Thanks,
-B

[ChrisA]

The first few lines are as follows:

autoPtr<combustionModels::psiCombustionModel> reaction
(
    combustionModels::psiCombustionModel::New(mesh)
);

psiReactionThermo& thermo = reaction->thermo();
thermo.validate(args.executable(), "h", "e");
basicMultiComponentMixture& composition = thermo.composition();

Which is similar to yours, just with psi combustion model instead of rho, it should work, I'm not sure why it isn't picking up on reaction being declared, it looks to be done correctly.