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/)
-   -   rhoCentralFoam transport equation (https://www.cfd-online.com/Forums/openfoam-programming-development/110274-rhocentralfoam-transport-equation.html)

JoaoDMiranda December 6, 2012 14:50

rhoCentralFoam transport equation
 
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 December 7, 2012 05:39

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 December 7, 2012 11:01

Thanks a lot! I will look into it!

Bryan S December 5, 2013 20:28

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

Henning86 December 6, 2013 08:03

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 December 6, 2013 08:36

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 December 6, 2013 10:59

2 Attachment(s)
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 December 6, 2013 14:42

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 December 6, 2013 18:53

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 December 6, 2013 19:02

1 Attachment(s)
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 December 6, 2013 19:09

Quote:

Originally Posted by Bryan S (Post 465170)
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:

Code:

   
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 December 6, 2013 20:59

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 December 6, 2013 22:16

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 December 7, 2013 02:02

Quote:

Originally Posted by Bryan S (Post 465194)
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 December 9, 2013 13:16

Quote:

Originally Posted by ChrisA (Post 465202)
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:

Quote:

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 December 9, 2013 13:29

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:

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 December 9, 2013 14:09

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:
Quote:

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 December 9, 2013 14:16

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:

Code:

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 December 9, 2013 14:49

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
Quote:

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 December 9, 2013 14:58

The first few lines are as follows:

Code:

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.

Bryan S December 9, 2013 15:46

Found it, I was missing #include "psiCombustionModel.H" in the main .C file. wmake worked successfully! Would you mind telling me what your chemistryProperties and combustionProperties would look like if you didn't want to model any reactions (i.e. passive scalar only)? I tried copying mine from the reactingFoam tutorial and changing "chemistry" to "off" in chemistryProperties and setting "active" to "false" in combustionProperties but the solver crashed immediately with
Quote:

terminate called after throwing an instance of 'std::bad_cast'
what(): std::bad_cast
Aborted (core dumped)
Any idea what's going on? Thanks!
-B

ChrisA December 9, 2013 15:50

1 Attachment(s)
Here's a copy of my chemistryProperties and combustionProperties files, they should get you started.

Bryan S December 9, 2013 19:48

Thanks Chris, that helped. I was able to work with my constant/ and system/ files and I successfully took 1 timestep forward before getting this:

Code:

...
Starting time loop

Mean and max Courant Numbers = 0.281002 0.96
deltaT = 0.000414938
Time = 0.000414938

diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal:  Solving for rhoUx, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal:  Solving for rhoUy, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal:  Solving for rhoE, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal:  Solving for rhoYi, Initial residual = 0, Final residual = 0, No Iterations 0
smoothSolver:  Solving for Gas1, Initial residual = 1.62351e-14, Final residual = 1.62351e-14, No Iterations 0
diagonal:  Solving for rhoYi, Initial residual = 0, Final residual = 0, No Iterations 0
smoothSolver:  Solving for Gas2, Initial residual = 1.62707e-14, Final residual = 1.62707e-14, No Iterations 0
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigSegv::sigHandler(int) at ??:?
#2  in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::operator=(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
#4 
 at ??:?
#5  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#6 
 at ??:?
Segmentation fault (core dumped)


This error looks familiar. I've gotten sigHandler errors before and I've never been able to fix them. Is there an easy fix or is it time to give up and try something else? Thanks.
-B

ChrisA December 9, 2013 19:52

Hard to say what's giving you a segmentation fault, looks like its solved all of your equations (your new Y equations included), try finding out what line is causing the fault (if you can) check your BC's etc. etc. Try reducing the courant number as well, things get shaky at Co = 1 with rhocentral from my experience. My cases run at 0.3-0.5.

Bryan S December 9, 2013 20:10

Yeah... It's time for me to hang up my hat on this one. I'm going to try to post-process the passive scalar equation in MATLAB. Thanks for all your help Chris. If you ever get a working version of your solver to the point that you'd feel like sharing it, I'd love to get a copy of it. Best,
-B

chriss85 July 18, 2014 08:52

I have been looking into related problems right now, and the YEqn posted by Chris has been quite helpful to me so far.

In my case I have one default gas and I'm producing two others. Do I also need a transport equation for my default gas then?
I have added the mass source terms to the total density as well as to the single transport equations, so the total mass should be conserved as far as I can see. Then I'm basically saying Y0 = 1 - Y1 - Y2. Is this approach sufficient or should I solve for the default gas as well?

On another topic, I'm currently in the process of switching from a sonicFoam-based solver to a rhoCentralFoam-based one, and it's unclear to me how I need to treat additional source terms in impulse and energy equations, i.e. forces and energy sources/sinks (from radiation and electrical current).

In sonicFoam I could just add them as source terms, but here I don't know if they should be added in predictor and/or corrector steps, and if some kind of central scheme is required for them.

I have also posted this in the following thread:
http://www.cfd-online.com/Forums/ope...tml#post502094

If anyone can shed some light on this I would be very grateful.

chriss85 July 24, 2014 10:21

No ideas? :(

reza_65 October 7, 2016 10:36

Quote:

Originally Posted by Henning86 (Post 465140)
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

Hi Henning,

I am wondering why you did not add any source term for combustion in the energy equation? if you solve just a cold flow it can be reasonable but for reacting flow this equations are missing combustion source term! So, it should not work for reacting flow!

Cheers,
Reza

TommyM January 15, 2020 09:12

Looking at the files provided by Henning, is the YEqn properly placed in the .C file?
As far as I know, the species concentration equation has to be solved between momentum and energy equation, while here it is solved last.
This is my opinion, but maybe I am missing something.


Tommy


All times are GMT -4. The time now is 03:44.