[ChrisA]

Firstly, thanks everyone who contributed to this thread, it's been very helpful with my development.

The issue I'm having is that I would like to handle the species conservation equation at the solver level (not use mstm().correct() ) for a variety of reasons.

I'd like to use diffusion coefficients from the Maxwell Stefan mass transport model and Chapman Enskogg binary diffusion coefficients.

From what I've managed to figure out I'm going to need to be able to assign a variable to my diffusion coefficients and be able to call updatecoefficients() from the mass transport model (maxwell steffan in this case). It's easy enough to say but I've had little luck, I think I have pulling the diffusion coefficients up to my solver from the mass transport model done properly but I can't seem to figure out how to call updatecoefficients().

I feel like this is a really simple thing to program with only minor changes to the code but I'm stuck as to what to do.

What I've tried is adding the function declaration for updateMolarFractions() (void updateMolarFractions(); ) to the multispeciestransportmodel header file and then calling updateMolarFractions using mstm().updateMolarFractions() (the mstm() implementation is the same as in the modified reacting foam). This compiles with /wmake libso but when i go to make the application with /wmake I get the following error: undefined reference to `Foam::multiSpeciesTransportModel::updateCoefficie nts()' Like I said, I feel this is simple and I'm just missing something basic.

Edit: managed to get it to solve for the maxwell coefficients at the solver level using the mstm pointer for the Dij information. I would still be interested to see if there's a more eloquent solution out there though.

[olivierG]

Hello mhsn,

Originally Posted by mhsn

Hey guys, thanks to your help, I could make modifiedReactingFoam work with openFoam2.1.1.
you may have to change some parts to be adapted to your openFoam especially in modifiedReactingFoam/Make/options and modifiedReactingFoam/Make/files

let me know if you needed that. I couldn't upload them regarding the maximum limit for the file size.

I would use the multi diffusion with 2.1.x, so can you please show the modification you have done ?

regards,
olivier

[mhsn]

Olivier, I couldn't attach the files here. If you give me your email, I can email you the modified files.

[olivierG]

hello,

Thanks mhsn , Valerio, and Co .....
I put the 2.1.x version here if other people want it.

regards,
olivier

[ldeng]

hello everyone,

i am now using also this multispeciesLibary (champan-Enskog+ Fick) in OpenFoam-2.1.0 . i simulated a One Dimensional premixed laminar flame with reduced mechanisms (methane oxidation, 15 species, 21 reactions) in OpenFoam, and got two problems when I compared the results with cantera and original reactingFoam.

1. temperature does not match
2. production (H2O, CO2) does not match

The operating conditions and boundary conditions in cantera and openfoam are same.

i already compared all coefficient in OpenFoam and cantera, such as mixture diffusion coefficient, thermoal conductivity, viscosity, ... They are all same. I also directly replace the term -fvm::lapacian(D_[i], yi) with the unity lewis Nr. term -fvm::laplacian(turbulence.muEff(), yi) in the subclass Fick.C, and remove mass diffusion term in hs equation, the it performs just like reactingFoam. It proved that, the connection between this multispeciesLibary and themoPhysical packet in Openfoam is fine.

I am wondering the mixture mass diffusion coefficient may break the species conservation equation. I really appreciate that anyone who are using this library can share your experiences with me.

I checked the boundary condition at the inlet used in cantera for burner stabilized flame, there is a mass-flux balance boundary condition (tau_i = y_i +j_i/phi) . I think i need to apply such boundary condition in OpenFoam. i will try mix of fixed value + fixed Gradient.

[olivierG]

hello,

I think you are right.
There is a common issue with inlet diffusion: the total flux is not only convective , but also diffusive, so if you set a fixedValue for a concentration, whit a know velocity, you will be wrong with your species flux calculation.
Take a look at this thread:
http://www.cfd-online.com/Forums/ope...chemistry.html

So try to set diffusivity to 0 at inlet, or set a mixed condition.
P.S: if you could set up a mixed condition for this, fell free to share !

regards,
olivier

[ChrisA]

Not sure if this will help, but I've managed to get around the "extra" mass due to diffusion at an inlet by correcting the inflow conditions for the mass diffusion.

I simply found the difference between my expected mass flow rate and the mass flow rate that openfoam gave me (done by integrating rho*U* inlet species concentration over a slice). We know that the difference is due to diffusion so if we assume the mass flow due to diffusion is a constant correct the inflow for that difference we arrive at a better mass flow rate. It won't be exact as the diffusion will change slightly as you change the concentration by changing the inflow velocity.. but with a few iterations you could conceivably approach the correct inflow mass flow rate.

This of course only works for specific conditions where you can assume and calculate the necessary values... and isn't ideal because it potentially requires iterative simulations... but its working for my research and might be applicable to a few other people's cases.

[olivierG]

hello Chris,

Yes your method should work, but is case specific.
You should also take care at outlet, if you use inletOutlet condition: I've found that we can have local backflow over some face, which then lead to diffusive flux. The solution is to use zeroGradient, not inletOutlet.

regards,
olivier

[hg2lf]

In the file diffusivity.C:

void Foam::diffusivity::update()
{
for(label i=0;i=species_.size();i++)
{
for(label j=i;j<species_.size();j++)
{
label k=species_.size()*i+j-0.5*i*(i+1);
Dij_[k]=eps_/tau_*DijModels_[k].D();
}

}
}
"j=i" means that the diffusivity is calculated for species themselves? Maybe there are something that I haven't consider?

[wyldckat]

Greetings hg2lf,

Just to confirm: do you mean the code that is in the file "diffusivityModels/diffusivityModel/diffusivityModel.C"? (I've opened the one for 2.1.x.)

More specifically this piece of code?

void Foam::diffusivityModel::update()
{    
    for(label i=0; i < species_.size(); i++)
    {
        for(label j=i; j < species_.size(); j++)
        {
            label k = species_.size()*i+j-0.5*i*(i+1);
            Dij_[k] = eps_/tau_*DijModels_[k].D();
        }
    }
}

If so, here's what I can figure out:

1. This kind of dual for loop is for going through the upper triangle of the matrix, including the diagonal.
2. The "label k" line is for a special storage mechanism they've used, given the size of the "DijModels_" array, which is defined in the constructor of "diffusivityModel":
   

   DijModels_.setSize(0.5*species.size()*(species.size()+1));

3. If we also look at the constructor, we can see how the "DijModels_" array is initialized, namely with the following code:
   

               DijModels_.set
               (
                   k,
                   binaryDiffusivityModel::New
                   (
                       species[i],
                       species[j],
                       dic_,
                       p,
                       T
                   )
               );

4. So in conclusion, when "j=i", it should be related to the "binaryDiffusivityModel" of the "i" item of the species... which means that it's "binaryDiffusivityModel" that should tell the remaining story...
5. Looking at "diffusivityModels/binaryDiffusivityModel/binaryDiffusivityModel.C", the following can be found:
   

   Foam::binaryDiffusivityModel::binaryDiffusivityModel
   (
       const word& name1,
       const word& name2,
       const dictionary& dic,
       const volScalarField& p,
       const volScalarField& T
   )

   Which means that the "species" array is only carrying names of the species in question (I was thinking it was another model itself...). Therefore, it makes some sense that a model for the diffusion of the same species can come in handy, specially if it's the only species around

That's as far as I've looked. I'm not familiar with this library, so I don't know any specific details about it. But my estimate is that the "binaryDiffusivityModel" could refer to at least two extremes: two particles of the same "species" can repel each other or can attract each other.

Best regards,
Bruno

[hg2lf]

Thanks,Bruno,
But I still can't understand. I need to modify this library to calculate the diffusivity for gas mixtures. So maybe the physical problem is not the same with the exist one. I'll do more work for this library, so that I could solve this problem.
Best regards.
hg2lf

[wyldckat]

Hi hg2lf,

From your description, I can only guess that you need either one of the following possible solutions:

• Perhaps you need to expand from binary to ternary? I.e. to take into account how three particles of different species interact with each other at the same time?
• Or you need to change how the binary reactions are being used?
• Or change completely the model for diffusion?

Either way, I'm only able to suggest that you have a look into the folders nearby the one for the code I referred, because I vaguely remember seeing that there are already some other diffusion models that can be used instead of this one... but I haven't studied them properly.

Good luck! Best regards,
Bruno

[ChrisA]

Hi hg2lf,

I can't comment on it much because I stopped using this library for my work (I'm using rhoCentralFoam so I needed a predictor step in my species equation, which this library does not allow for) but I can say that the maxwell-steffan model is a full multi-component model (handles many gas species). Ansys FLUENT uses this model, although I'm not sure if its the same implementation as is in the library.

We're also using the maxwell-steffan model for our work, although directly implemented in the solver rather than through a library.

The way it works in this library (and the model in general), to my limited understanding, is that the maxwell-steffan model uses binary diffusion coefficients to build an interaction matrix and then gets a multi component diffusion coefficient from that matrix for each species by summing the relevant terms. You can specify which binary diffusion coefficient to method to use (eg. chapman enskogg) and the maxwell-steffan model (should) reduce to the binary diffusion coefficient when only 2 gases are present.

Hope that helps,
Chris

[tmu]

Hi Friends
I use this library for simulation laminar combustion in tube. but I have some problem:
maximum temperature in simulation differ about 100k with maximum temp in Reference. for this reason mass fraction of species have difference.
I read library an manual and compare equations in this library with Ref (for example: Theoretical and Numerical Simulation : Poinsot ). I found there are some difference.
1. In YEqn: we have the correction velocity that is added to convection velocity
to ensure conservation of species.although this term is disappeared in library.
2 . In hsEqn: we need to add one term in hs equation. although in library two terms are added.
please help me if there are any reference for this solver that speak about equation.

[hg2lf]

Hello,
everyone, in the file Fick.C, in the function updateCoefficient(), the code like this:

forAll(species(), j)
{
if (j != i)
{
tmpGamma += x_[j] / Dij(i,j);
}
}

I think I could understand the question above.

[hg2lf]

Hello,
everyone,
In the file multiSpeciesTransportModel.C, in the function updateMolarFractions(), there is a piece of code, like this:

const scalarField& psiCells = thermo_.psi().internalField();
const scalarField& TCells = thermo_.T().internalField();
const scalarField& yCells =
thermo_.composition().Y(i).internalField();
scalarField& xCells = x_[i].internalField();
forAll(xCells, cellI)
{
xCells[cellI] = yCells[cellI] * psiCells[cellI] * TCells[cellI]
* R / W(i);
}

The code:"xells[cellI] = yCells[cellI] * psiCells[cellI] * TCells[cellI]* R / W(i);" express the idea:

x=Y*(1/RT)*T*Rcom/w(i);

the bug is w(i), which should be the molecular weight of the mixtures.

So, I think we should modified this piece of code.

[ldeng]

hi hg2lf,
x(i) = p(i)/p
p(i) = rho(i) *R/w(i) *T
p = rho/psi
then x(i) = rho(i) *R/w(i) *T *psi/rho = y(i)* psi *R *T/w(i), w(i) is for pure specie
so code is correct

[ChrisA]

Originally Posted by tmu

Hi Friends

1. In YEqn: we have the correction velocity that is added to convection velocity
to ensure conservation of species.although this term is disappeared in library.

I believe that the method being used in this code is the Y(inert) = 1 - Y(everything else). You'll have to verify this though, I haven't looked at the library in a long time. The method is mentioned briefly in "Theoretical and Numerical Combustion" by Poinsot (Which I think is the book you meant to refer to in your post?) Anyway, as is mentioned there, there is an error associated with the method but if the inert species concentration is dominant in the simulation this error becomes negligible. It also disappears if the diffusion coefficients are all equal.

For reference, reactingFoam also uses this method in it's species equations.

I cannot comment too much on the heat transfer terms, in my own energy equation I'm only including the enthalpy transport term (sum of h_i*J_i, still ironing out the formulation for that myself.)

Regarding references for the library and it's equations, there should be a PDF in the earlier posts describing the equations used, there's also a paper by Valerio to look at.

[ldeng]

Hi everyone,

as a conclusion to my previous diffusion flux problem:

The diffusion flux at inlet is a common issue. To avoid that we have two options:

1. specify the total flux( diffusion + convection) equals to a fixed value. it is a kind of mixture boundary condition, but difficult(my opinion) to implement into OpenFOAM.

2. close the diffusion flux(set neumann bc to zeroGradient at inlet), this is not elegant but the easiest.

because i also simulated a 1-D burner stabilized flame in Cantera. In order to obtain a better results specially near the inlet, I specified a fixed value to all species as an initial condition based on the simulated results. After the combustion proceeded, i change the BC of intermediate species back to zero. Anyway, by such approach, i got a nearly perfect results (see attachment) compared to Cantera.

regards,

lei

[idjibb]

hello All,

Did anybody manage to get the install script to work for OF211 ?

regards,