CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM (https://www.cfd-online.com/Forums/openfoam/)
-   -   Steady state multi species (https://www.cfd-online.com/Forums/openfoam/80130-steady-state-multi-species.html)

 Hrushi September 16, 2010 05:41

Steady state multi species

Hi, I am trying to develop liquid phase multi species solver with/without reaction (no combustion).
Can anybody explain me how to define species, species properties in OpenFOAM?
How to create lookup table for species so that each species is assigned some number?
So that I can solve the transport equation for each species like that done in Yeqn.H code in reactingFoam solver.
Mass fraction given by the above YEqn can then be used to calculate mass weighted average mixture properties like:
rhomix = rho1*Y[1] + rho2*Y[2] + rho3*Y[3]

This average mixture properties then can be used to solve the momentum equation.

Reaction effect can be introduced by adding the source term in Yeqn.H later on.

Any suggestions , Please...:)

 marcbest September 16, 2010 08:27

as a beginnig you can use this report http://projekter.aau.dk/projekter/fi...784/Report.pdf

 Hrushi September 27, 2010 05:06

individual species properties in for loop to calculate rho_mix and nu_mix

hi,

I would like to use mixture viscosity in the solving momentumm equation. For that I need to calculate mixture density (mass weighted). From that I can calculate mixture kinematic viscosity (nu_mix) . To do this, I need access to species mass fraction, species density, species viscosity.
rho_mix = rho1*Y[1]+rho2*Y[2]+rho3*Y[3]
nu_mix= (rho1*nu1*Y[1]+rho2*nu2*Y[2]+rho3*nu3*Y[3])/rho_mix

My question is how can I call individual species properties in for loop to calculate rho_mix and nu_mix?

for (lable i=0, i<Y.size(),i++)
(
rho_mix=rho_mix+rho1:(*Y[i]+rho2:(*Y[i]+rho3:(*Y[i]
)

 Cyp September 27, 2010 09:09

hi Hrushi!

I think you should defined your rho1,rho2,rho3 as a list of scalar : rho=(rho1,rho2,rho3)

hence, rho[0]=rho1 ; rho[1|=rho2 ;rho[2]=rho3

and you can use the following loop :

Code:

```rho_mix=0.0 for (int i=0; i<Y.size(); i++)     {     rho_mix=rho_mix+rho[i]*Y[i];     }```

else, you can define your mixture density as :

Code:

`rho_mix = rho1*Y[0] + rho2*Y[1] + rho3*Y[2]`
(keep in mind that C++ starts the numerotation from 0)

@++

 Hrushi October 24, 2010 08:03

error: no match for 'operator='

Hi Cyprien,

I have tried your suggestion. But when I compile the code, I get the following error message.
SOURCE=simpleMyFoam.C ; g++ -m64 -Dlinux64 -DWM_DP -Wall -Wno-strict-aliasing -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3 -DNoRepository -ftemplate-depth-40 -I/export/apps/OpenFOAM/OpenFOAM-1.6/src/turbulenceModels -I/export/apps/OpenFOAM/OpenFOAM-1.6/src/turbulenceModels/incompressible/RAS/RASModel -I/export/apps/OpenFOAM/OpenFOAM-1.6/src/transportModels -I/export/apps/OpenFOAM/OpenFOAM-1.6/src/transportModels/incompressible/singlePhaseTransportModel -I/export/apps/OpenFOAM/OpenFOAM-1.6/src/finiteVolume/lnInclude -IlnInclude -I. -I/export/apps/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude -I/export/apps/OpenFOAM/OpenFOAM-1.6/src/OSspecific/POSIX/lnInclude -fPIC -c \$SOURCE -o Make/linux64GccDPOpt/simpleMyFoam.o
In file included from simpleMyFoam.C:66:
YEqn.H:21: warning: left-hand operand of comma has no effect
YEqn.H:21: warning: right-hand operand of comma has no effect
YEqn.H:21: error: no match for 'operator=' in 'rho = (0, rho3)'
/export/apps/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/List.C:447: note: candidates are: void Foam::List<T>:operator=(const Foam::UList<T>&) [with T = double]
/export/apps/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/List.C:479: note: void Foam::List<T>::operator=(const Foam::List<T>&) [with T = double]
/export/apps/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/List.C:494: note: void Foam::List<T>::operator=(const Foam::SLList<T>&) [with T = double]
/export/apps/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/List.C:522: note: void Foam::List<T>::operator=(const Foam::IndirectList<T>&) [with T = double]
/export/apps/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/List.C:541: note: void Foam::List<T>::operator=(const Foam::UIndirectList<T>&) [with T = double]
/export/apps/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/List.C:560: note: void Foam::List<T>::operator=(const Foam::BiIndirectList<T>&) [with T = double]
/export/apps/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/ListI.H:134: note: void Foam::List<T>:operator=(const T&) [with T = double]
/export/apps/OpenFOAM/OpenFOAM-1.6/src/finiteVolume/lnInclude/readSIMPLEControls.H:6: warning: unused variable 'momentumPredictor'
/export/apps/OpenFOAM/OpenFOAM-1.6/src/finiteVolume/lnInclude/readSIMPLEControls.H:12: warning: unused variable 'transonic'
make: *** [Make/linux64GccDPOpt/simpleMyFoam.o] Error 1

The error repeats for line 41 and 61 also.

Here is my YEqn.H:
/*tmp<fv::convectionScheme<scalar> > mvConvection
(
fv::convectionScheme<scalar>::New
(
mesh,
fields,
phi,
mesh.divScheme("div(phi,Yi_h)")
)
);*/

{
label inertIndex = -1;
volScalarField Yt = 0.0*Y[0];

const dimensionedScalar& rho1 = speciesProperties.lookup("rho1");
const dimensionedScalar& rho2 = speciesProperties.lookup("rho2");
const dimensionedScalar& rho3 = speciesProperties.lookup("rho3");

List <scalar> rho(3);
rho = (rho1,rho2,rho3);

volScalarField rho_mix
(
IOobject
(
"rho_mix",
runTime.timeName(),
mesh,
IOobject::AUTO_WRITE
),
mesh
);

const dimensionedScalar& mu1 = speciesProperties.lookup("mu1");
const dimensionedScalar& mu2 = speciesProperties.lookup("mu2");
const dimensionedScalar& mu3 = speciesProperties.lookup("mu3");

List <scalar> mu(3);
mu = (mu1,mu2,mu3);

volScalarField mu_mix
(
IOobject
(
"mu_mix",
runTime.timeName(),
mesh,
IOobject::AUTO_WRITE
),
mesh
);

const dimensionedScalar& nu1 = speciesProperties.lookup("nu1");
const dimensionedScalar& nu2 = speciesProperties.lookup("nu2");
const dimensionedScalar& nu3 = speciesProperties.lookup("nu3");

List <scalar> nu(3);
nu = (nu1,nu2,nu3);

for (label i=0; i<Y.size(); i++)
{
if (Y[i].name() != inertSpecie)
{
volScalarField& Yi = Y[i];

solve
(
fvm::div(phi,Yi)
// mvConvection->fvmDiv(phi, Yi)
- fvm::laplacian(turbulence->nuEff(), Yi)
==
source_r
);

Yi.max(0.0);
Yt += Yi;
rho_mix += rho [i] * Yi;
mu_mix += mu [i] * Yi;

}
else
{
inertIndex = i;
}
}

Y[inertIndex] = scalar(1) - Yt;
Y[inertIndex].max(0.0);
rho_mix += rho[inertIndex] * Y[inertIndex];
mu_mix += mu[inertIndex] * Y[inertIndex];

volScalarField nu_mix
(
IOobject
(
"nu_mix",
runTime.timeName(),
mesh,
IOobject::AUTO_WRITE
),
mu_mix/rho_mix
);

}

I am not able to understand the error message. I tried to look into list.H and list .C but could not get the clue.
Can you (anyone) please help me out in how to define the list of scalars and their assignment?

Thanks and Regards,

Hrushikesh

 Cyp October 24, 2010 14:54

did you try ptrList instead of List ?

 Hrushi October 27, 2010 06:24

Error while execution

Hi,

I have modified the simpleFoam solver. I have included specie equation. The code compiles giving no compliation error. But when I run my code, it gives following error message:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting RAS turbulence model kEpsilon
kEpsilonCoeffs
{
Cmu 0.09;
C1 1.44;
C2 1.92;
alphaEps 0.76923;
sigmaEps 1.3;
}

Starting time loop

Time = 1

DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 0.000312273, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 0.000245102, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.779043, Final residual = 0.000142658, No Iterations 1
DICPCG: Solving for p, Initial residual = 1, Final residual = 0.00865508, No Iterations 87
time step continuity errors : sum local = 0.00383834, global = 3.55093e-06, cumulative = 3.55093e-06
DILUPBiCG: Solving for epsilon, Initial residual = 0.868617, Final residual = 0.0142968, No Iterations 1
DILUPBiCG: Solving for k, Initial residual = 1, Final residual = 0.0286599, No Iterations 1

hanging pointer, cannot dereference#0 Foam::error::printStack(Foam::Ostream&) in "/export/apps/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#1 Foam::error::abort() in "/export/apps/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#2 main in "/home/na16116/run/run/simpleMyFoam/bin"
#3 __libc_start_main in "/lib64/libc.so.6"
#4 Foam::regIOobject::writeObject(Foam::IOstream::str eamFormat, Foam::IOstream::versionNumber, Foam::IOstream::compressionType) const in "/home/na16116/run/run/simpleMyFoam/bin"

From function PtrList::operator[]
in file /export/apps/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/PtrListI.H at line 123.

FOAM aborting

Abort

Can anyone explain what is the reason for this error message?

Thanks
Hrushikesh

 Hrushi October 28, 2010 06:02

Hanging pointer error

Here are the details of my case:
createFields.H

Info<< nl << "Reading thermophysicalProperties" << endl;
PtrList<volScalarField> Y(5);

IOdictionary speciesProperties
(
IOobject
(
"speciesProperties",
runTime.constant(),
mesh,
IOobject::NO_WRITE
)
);

wordList speciesNames
(
speciesProperties.lookup("speciesNames")
);

//speciesTable (speciesNames);

dimensionedScalar source_r
(
speciesProperties.lookup("source_r")
);

word inertSpecie(speciesProperties.lookup("inertSpecie" ));

Info << "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::AUTO_WRITE
),
mesh
);

Info << "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::AUTO_WRITE
),
mesh
);

# include "createPhi.H"

label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);

singlePhaseTransportModel laminarTransport(U, phi);

autoPtr<incompressible::RASModel> turbulence
(
incompressible::RASModel::New(U, phi, laminarTransport)
);

YEqn.H

/*tmp<fv::convectionScheme<scalar> > mvConvection
(
fv::convectionScheme<scalar>::New
(
mesh,
fields,
phi,
mesh.divScheme("div(phi,Yi_h)")
)
);*/

label inertIndex = -1;
volScalarField Yt = 0.0*Y[0];

const dimensionedScalar& rho1 = speciesProperties.lookup("rho1");
const dimensionedScalar& rho2 = speciesProperties.lookup("rho2");
const dimensionedScalar& rho3 = speciesProperties.lookup("rho3");

const dimensionedScalar rlist[3] = {rho1,rho2,rho3};
/*FixedList<dimensionedScalar,3> FixedRList(rlist);

List <dimensionedScalar> rho(FixedRList);*/

volScalarField rho_mix
(
IOobject
(
"rho_mix",
runTime.timeName(),
mesh,
IOobject::AUTO_WRITE
),
mesh
);

const dimensionedScalar& mu1 = speciesProperties.lookup("mu1");
const dimensionedScalar& mu2 = speciesProperties.lookup("mu2");
const dimensionedScalar& mu3 = speciesProperties.lookup("mu3");

const dimensionedScalar mulist[3] = {mu1,mu2,mu3};
/*FixedList<dimensionedScalar,3> FixedMuList(mulist[3]);

List <dimensionedScalar> mu(FixedMuList);*/

volScalarField mu_mix
(
IOobject
(
"mu_mix",
runTime.timeName(),
mesh,
IOobject::AUTO_WRITE
),
mesh
);

/*const dimensionedScalar& nu1 = speciesProperties.lookup("nu1");
const dimensionedScalar& nu2 = speciesProperties.lookup("nu2");
const dimensionedScalar& nu3 = speciesProperties.lookup("nu3");

const dimensionedScalar nulist[3] = {nu1,nu2,nu3};
FixedList<dimensioned<double>,3> FixedNuList(nulist);

List <dimensionedScalar> nu(FixedNuList);*/

for (label i=0; i<Y.size(); i++)
{
if (Y[i].name() != inertSpecie)
{
volScalarField& Yi = Y[i];

solve
(
fvm::div(phi,Yi)
// mvConvection->fvmDiv(phi, Yi)
- fvm::laplacian(turbulence->nuEff(), Yi)
==
source_r
);

Yi.max(0.0);
Yt += Yi;
rho_mix += rlist [i] * Yi;
mu_mix += mulist [i] * Yi;

}
else
{
inertIndex = i;
}
}

Y[inertIndex] = scalar(1) - Yt;
Y[inertIndex].max(0.0);
//rho_mix += rlist[inertIndex] * Y[inertIndex];
//mu_mix += mulist[inertIndex] * Y[inertIndex];

volScalarField nu_mix
(
IOobject
(
"nu_mix",
runTime.timeName(),
mesh,
IOobject::AUTO_WRITE
),
mu_mix/rho_mix
);

simplemyFoam.C
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
This file is part of OpenFOAM.

OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.

OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

Application
simpleFoam

Description
Steady-state solver for incompressible, turbulent flow

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "RASModel.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createFields.H"
#include "initContinuityErrs.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nStarting time loop\n" << endl;

while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;

#include "initConvergenceCheck.H"

p.storePrevIter();

// Pressure-velocity SIMPLE corrector
{
#include "UEqn.H"
#include "pEqn.H"
}

turbulence->correct();

#include "YEqn.H"

runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;

#include "convergenceCheck.H"
}

Info<< "End\n" << endl;

return 0;
}

// ************************************************** *********************** //

The code compiles. But it crashes when executing during 1st iteration with hanging pointer error described in previous reply.

Hrushikesh

 Cyp November 2, 2010 11:09

Code:

```createFields.H Info<< nl << "Reading thermophysicalProperties" << endl; PtrList<volScalarField> Y(5);```
Is it normal that you have Y(5) instead of Y(3) ??

 Hrushi November 3, 2010 00:13

Yes, it doesnot make any difference. The problem lies in how to fill up the ptrList? I looked at ptrList code. The set function is used for the same.
I am trying with that now. Do you know the correct way/alternate way to fill up the ptrList?

Hrushi

 nakul November 3, 2010 02:34

Hi Hrushi,

I am also trying to implement something similar. Can you tell me what thermodynamics package you are using for your mixture in your thermoPhysicalProperties dictionary?

I am using hsCombustionThermo with "multiComponentmixture" but I am not being able to get it right.

There is no example as to how multiple species are read by OpenFOAM while using multicomponent mixture.

Could you also tell how are you reading in the multiple species of the mixture?

Thanx

 Hrushi November 14, 2010 06:40

Multi species solver

Hi,

Sorry for late reply (was on Diwali vacation). The problem is solved now. Code is compiling and running. I will update you soon with next development.

Currently, I am trying to introduce heat transfer effects in above developed solver.

Thanks

Hrushikesh

 alfa_8C December 1, 2010 10:05

Hello Hrushi,

I followed your thread and it seems you've done exactly the same as I am starting to do right now. Is your solver anyway available as a contribution?

If not, what did you change then from the last post on, to make the solver not only compiling but also running?

Best Regards form Zurich to Mumbai,
Tony

 Hrushi December 2, 2010 00:42

Quote:
 Originally Posted by alfa_8C (Post 285575) Hello Hrushi, I followed your thread and it seems you've done exactly the same as I am starting to do right now. Is your solver anyway available as a contribution? If not, what did you change then from the last post on, to make the solver not only compiling but also running? Best Regards form Zurich to Mumbai, Tony
Hi Tony,

No, it is not. The current solver does only the multi species part without reaction and heat transfer. It just transports mixture of species frpom one point to another. Also, the solver is not stable when it is running. The initial residual values do not come below 0.1. They fluctuate between 0.12 to 0.73. The solution is not converged. Do you know any reason for this?:(:)

Also now I am trying to include heat transfer in the above solver by introducing enthalpy equation. I am confused with thermophysical model framework available in OF.:confused: Can you help me out in this or give suggestions on it?:)

Hrushi

 alfa_8C December 2, 2010 06:23

Hello Hrushi,

what is the reason for you to create a solver with multicomponent mixture model on your own? You could easily use reactingFoam. I think this solver would fulfill all your needs, as it provides mixture, heat transfer, reactions on/off, but no radiation by the way. But this is not really a problem since it is very easy to implement a radiation model in OpenFoam. If you need help regarding this just let me know.

Tony

 the king 2 December 2, 2010 06:33

hi

sry...my reply nt regarding thread..i am new to the forum....iwant to learn gambit and fluent.do u knw some institute for training?

 Hrushi December 2, 2010 07:37

Quote:
 Originally Posted by alfa_8C (Post 285703) Hello Hrushi, what is the reason for you to create a solver with multicomponent mixture model on your own? You could easily use reactingFoam. I think this solver would fulfill all your needs, as it provides mixture, heat transfer, reactions on/off, but no radiation by the way. But this is not really a problem since it is very easy to implement a radiation model in OpenFoam. If you need help regarding this just let me know. Tony
Hi Tony,

Thanks for giving your time. The reactingFoam solver updates the species properties using gas law. I don't want that. I needed a solver to solve the reaction, heat transfer in liquid phase with steady state approach and which would update physical properties as mass weighted average. I think reactingFoam doesn't provide this option. If it can, please let me know the way for this. The reason behind developing this solver is that I first wanted to introduce YEqn and calculate mixture density and use this in UEqn followed by heat transfer addition and reaction stuff.

Thanks again for sharing your thought on this. I really appreciate that. Will contact you if I need any help.;):D

Hrushi

 Hrushi December 2, 2010 07:38

Quote:
 Originally Posted by alfa_8C (Post 285703) Hello Hrushi, what is the reason for you to create a solver with multicomponent mixture model on your own? You could easily use reactingFoam. I think this solver would fulfill all your needs, as it provides mixture, heat transfer, reactions on/off, but no radiation by the way. But this is not really a problem since it is very easy to implement a radiation model in OpenFoam. If you need help regarding this just let me know. Tony
Hi Tony,

Thanks for giving your time. The reactingFoam solver updates the physical and thermophysical properties using gas law. I don't want that. I needed a solver to solve the reaction, heat transfer in liquid phase with steady state approach and which would update physical properties as mass weighted average. I think reactingFoam doesn't provide this option. If it can, please let me know the way for this. The reason behind developing this solver is that I first wanted to introduce YEqn and calculate mixture density and use this in UEqn followed by heat transfer addition and reaction stuff.

Thanks again for sharing your thought on this. I really appreciate that. Will contact you if I need any help.;):D

Hrushi

 prasanth August 11, 2011 02:07

Multi Species in OpenFOAM

Hello Hrushi,

I want to solve the multSpecies problem using OpenFOAM. I have seen your posting in forum regarding this solver. Have you success on this solver? If So, Could you please tell me what are the changes did you made for this solver

Regards
Prasanth.

Quote: