CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Programming & Development (
-   -   Transport equation in chtMultiRegionSimpleFoam solver (

Aurelien Thinat July 13, 2010 05:55

Transport equation in chtMultiRegionSimpleFoam solver
Hello everybody,

First of all, I am a new user of OpenFoam.

I need to develop a scalar transport equation (which would be mass fraction) in the chtMultiRegionSimpleFoam solver in order to study a multispecies case. I'm using the reactingFoam solver as an example. Here are my problems :

1) I added the YEqn.H file in the "fluids" folder and I adapted the equation for a steady state solver and non-reactive species (no source term).
I have a problem with the constructor :
fluidRegions[i], // initally "mesh",

What is "fields" related to ? Should I modify "mesh.divScheme" into fluidRegions[i].divScheme ?

2) How can I initialize my scalars ? In the createFluidFields.H ? I didn't find the reader of the dictionnary "Ydefaut" in the reactingFoam source code.

Thank you very much for your help.


Aurelien Thinat July 15, 2010 11:03

I have added 3 transport equations in the chtMultiRegionSimpleFoam solver.

I defined scalar[i]Fluid, i = 1, 2 or 3 as a ptrList<volScalarField>

Then I define scalar[i]Fluid[j] , j= number of FluidRegion and I call


But then, when I launch the computation I got this message :

"Time = 0

Solving for fluid region fluid
DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 0.00068726968, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 0.00057863562, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 1, Final residual = 0.00097519241, No Iterations 1
#0 Foam::error::printStack(Foam::Ostream&) in "/home/trap/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/"
#1 Foam::sigFpe::sigFpeHandler(int) in "/home/trap/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/"
#2 ?? in "/lib64/"
#3 Foam::DILUPreconditioner::calcReciprocalD(Foam::Fi eld<double>&, Foam::lduMatrix const&) in "/home/trap/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/"
#4 Foam::DILUPreconditioner::DILUPreconditioner(Foam: :lduMatrix::solver const&, Foam::dictionary const&) in "/home/trap/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/"
#5 Foam::lduMatrix::preconditioner::addasymMatrixCons tructorToTable<Foam::DILUPreconditioner>::New(Foam ::lduMatrix::solver const&, Foam::dictionary const&) in "/home/trap/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/"
#6 Foam::lduMatrix::preconditioner::New(Foam::lduMatr ix::solver const&, Foam::dictionary const&) in "/home/trap/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/"
#7 Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/home/trap/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/"
#8 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/home/trap/OpenFOAM/OpenFOAM-1.7.0/lib/linux64GccDPOpt/"
#9 Foam::lduMatrix::solverPerformance Foam::solve<double>(Foam::tmp<Foam::fvMatrix<doubl e> > const&) in "/home/trap/OpenFOAM/OpenFOAM-1.7.0/applications/bin/linux64GccDPOpt/mychtMultiRegionSimpleFoam"
#10 main in "/home/trap/OpenFOAM/OpenFOAM-1.7.0/applications/bin/linux64GccDPOpt/mychtMultiRegionSimpleFoam"
#11 __libc_start_main in "/lib64/"
#12 _start at /usr/src/packages/BUILD/glibc-2.10.1/csu/../sysdeps/x86_64/elf/start.S:116
Floating point exception"

Do you have any idea ?

EDIT : I solve UEqn first, then YEqn1, YEqn2, YEqn3, hEqn...

Aurelien Thinat July 16, 2010 11:09

I'm still on it. I have a problem to declare the diffusivity of the scalar.

In a multiregion solver, the createFluidFields.H file is a loop on each fluid region and the variables are defined as ptrList with a component for each region.

I managed to declare correctly the scalar and to read the BCs' files. But when I'm trying to do quite the same for the diffusivity scalar, ie :
- defining a ptrList<dimensionedScalar> Diffusivity(fluidRegions.size());
- reading the transportProperties dictionnary,

there is a problem of "redeclaration" of the diffusivity.

Does anyone have an idea or any example ?

marupio July 17, 2010 15:30

So your problem is now with defining a ptrList < dimensionedScalar > ?

Not sure why it's complaining about redeclaration... are you using the ptrList object correctly? I wrote a solver that did something similar. My scalars are not volume fractions, though... and I called it "gamma" instead of "diffusivity". Heres the code:


    // scalars is a subDictionary containing admScalar data
    dictionary admScalarsDict(admDict.subDict("scalars"));
    wordList admScalarKeywords(admScalarsDict.toc());
    label nScalars(admScalarKeywords.size());

    if (admDebug > 0)
        Info << "Found " << nScalars << " scalars listed under keywords:\n" << endl;
        for (label i=0; i<(nScalars-1); i++)
            Info << "[" << admScalarKeywords[i] << "], ";
        Info << "[" << admScalarKeywords[nScalars-1] << "]\n" << endl;
    PtrList<volScalarField> admVar(nScalars);
    PtrList<dimensionedScalar> gamma(nScalars);
    admScalar tempScalarID;
    forAll(admScalarKeywords, i)
        tempScalarID = admScalarsDict.lookup(admScalarKeywords[i]);
        Info << "Reading field named " <<;
        Info << " from keyword " << admScalarKeywords[i] << " in admDict. \n" << endl;

            new dimensionedScalar
                "gamma(" + + ")",
                dimensionSet(0, 2, -1, 0, 0),

            new volScalarField

Bufacchi July 17, 2010 19:55


Are you working on diffusion of species from one region to another? I'm not sure I understood your model correctly.


marupio July 17, 2010 21:12

This section of code is only in the initial setup (ie createFields). Since you were running into duplicate definition problems, I suspected you were using the ptrList wrong. The code is an example for that. As for the model, these are scalars that are used in a parallel biochemical model... so they are being transported by convection and diffusion, but they are not volume fractions because the biochemical model isn't mass preserving.

In other words, as far as OpenFOAM is concerned, these are just scalars.

Aurelien Thinat July 19, 2010 02:40

Thank you Marupio. It was a problem with the prtList. It seems it's working now, I just have a problem of dimension conflict.

Edit :
I might need your help again.

The equation I'm solving is based on the scalarTransport solver ie :

- fvm::laplacian(Dscalar1Fluid,scalar1Fluid)

but the equation I need to solve is :
div(U*c) - laplacian(Dc * c) = 0
with Dc = Dscalar1Fluid, c = scalar1Fluid, U = velocity.

I saw in the user guide that phi = rho * U,

The error I get is :
incompatible dimensions for operation
[scalar1Fluid[1 -3 -1 0 0 0 0] ] - [scalar1Fluid[0 0 -1 0 0 0 0] ]

It confirms my problem with phi. So I tried to solve "fvm::div(U,scalar1Fluid)" which is not working of course :

"no matching function for call to 'div(Foam::GeometricField<Foam::Vector<double>,Foa m::fvPatchField, Foam::volMesh>&,...)"


marupio July 19, 2010 09:19

phi is actually velocity flux interpolated to the volume faces, and appears frequently in the discretization process.

All the incompressible solvers I've looked at use transport equations that are divided by density... this changes the units you'd expect to see. What dimensions are you giving Dscalar1Fluid? In my solver, it turns out to be [0 2 -1 0 0 0 0].

Aurelien Thinat July 19, 2010 09:28

I had a look to the createFluidFields and phi is defined as :

new surfaceScalarField
& fluidRegions[i].Sf()

Well, I'm adding rho in my transport equation.

marupio July 19, 2010 10:07

Ah, so phi actually matches convention for once. I'm used to working with incompressible. You don't want to define another phi. I think the problem is simply a units issue with your diffusion coefficient. The code I gave you (above) was for incompressible... so maybe your diffusion coefficient needs to be multiplied by density.

It's something simple like that, I'm sure. Take a moment to work out the units on paper - make sure you are defining the diffusion coefficient properly.

(Boussinesq doesn't use density)

Aurelien Thinat July 19, 2010 10:21

I have added the density inside all my term of the equation. It still computes and the scalar is still unbounded.

The Y1Eqn.H files :

for (label i=0; i>scalar1Fluid.size(),i++)
- fvm::laplacian(rhoFluid[i]*Dscalar1Fluid[i], scalar1Fluid[i])

Do I have to add a test to detect the max between scalar1Fluid and 0/1 ?

marupio July 19, 2010 10:31

You could try it. I'm not entirely sure. If you want to bound the solution, try using:


for (label i=0; i>scalar1Fluid.size(),i++)
        - fvm::laplacian(rhoFluid[i]*Dscalar1Fluid[i], scalar1Fluid[i])



Aurelien Thinat July 20, 2010 05:57


My scalars transport equations are now working well. I didn't add the max and min function. Now I will have a look at the rho/Cp/lambda/mu state equations.

Thank you again.


Aurelien Thinat July 21, 2010 05:53

Hello everybody,

I need some help to clear my mind. My purpose is to adapt the chtMultiRegionSimpleFoam solver for multispecies case.

- So I started by adding 3 scalar transport equations, which will be use for the mass fraction of each species.

The next steps are :

- Adapt the EoS for the density and thermal diffusivity,
- Add the condition : Sum(Yi) = 1.

Am I right ?

siavash_abadeh December 4, 2010 13:09

dear Aurélien
i have a problem and i think you can help me.
my problem is mass transfer from solid to liquid with electrochemical method and i want simulate this problem with openfoam and i cant find solver for it.
how resolvent your problem?
use a solver or write a code for it.
if you write a code : can i have your code?:o
very very thank in advance

mm.abdollahzadeh April 10, 2012 13:51

Hi all,

I am also want to add a scaler equation to Chtmultiregionfoam.
I have two question:

1- How is it possible to add a dimenedScaler in creatsfeild which be readed from a dictionary. ( for example the diffissity )

2-is the boundary condition for conjucate heat transfer is general enough to be used for other scalar equations with out any modification?


nimasam April 11, 2012 01:17


Originally Posted by mm.abdollahzadeh (Post 354029)

1- How is it possible to add a dimenedScaler in creatsfeild which be readed from a dictionary. ( for example the diffissity )

to learn how to program in OpenFOAM, you should dig one solver for example icoFoam or .....
how ever i guess you want to make a variable diffusivity! thats why you want to use setField to assign it values! so first
create a volScalarField variable with the name diffusivity in creatFields.H in your solver, then you should compile your modified solver with wmake command (read user guide i think chapter 4 how to compile a solver)

then make a diffusivity file in folder 0 , here in this file you can assign uniform or non uniform value ( with setField) for diffusivity


2-is the boundary condition for conjucate heat transfer is general enough to be used for other scalar equations with out any modification?

whats your mean??

mm.abdollahzadeh April 11, 2012 03:51

Dear Nima

Thanks for your kind reply. I am a little bit familer with icoFoam. in icoFoam I can define a dimeionscaler and read it from a dictionary.
the thing is that I want to solve another equations like electric field in both solid and fluid zones. so i want to difine electric diffusitivity in a dictionary and read it seprately for both solid and fluid zone.
my second question is related to this fact that my electricpotential feild should be solved in both solid and fluid zones. I was thinking that can i use the turbulentTemperatureCoupledBaffleMixed boundary condition for another electricfield???
there is two option in that boundary ? K and KName which i dont know if can i change them to make this boundary applicable for electricfield?


nimasam April 11, 2012 05:36

so you should create your own boundary condition, im not familiar with electricpotential feild, but for first step to implement your bounday you may want to look
compressible::turbulentTemperatureCoupledBaffleMix ed
to see how it is implemented

mm.abdollahzadeh April 11, 2012 06:35

yes, I an on it, but a guy from chalmers where doing the same with privious versions and claim that its enough general!

All times are GMT -4. The time now is 04:43.