CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Adding transport equation to twoPhaseEulerFoam (OF231)

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By vkoppejan

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   February 7, 2015, 09:04
Question Adding transport equation to twoPhaseEulerFoam (OF231)
  #1
Member
 
A. Bernath
Join Date: Jun 2011
Location: Karlsruhe, Germany
Posts: 39
Rep Power: 14
derkermit is on a distinguished road
Hello all,
I'm struggling with the implementation of a scalar transport equation in twoPhaseEulerFoam. I've done that in various solvers before but this one seems to be "special".
Also, the explanation in the wiki (http://openfoamwiki.net/index.php/Ho...sport_equation) isn't sufficient for this solver.

At the moment, the psiEqn.H looks like this:

Code:
{
    fvScalarMatrix psiEqn
    (
        fvm::ddt(alpha1, rho1, psi) + fvm::div(alphaRhoPhi1, psi)
      - fvm::Sp(contErr1, psi)

      - fvm::laplacian
        (
            fvc::interpolate(alpha1)
           *fvc::interpolate(thermo1.alphaEff(phase1.turbulence().mut())),
            psi
        )
      ==
    psiSource*rho1*alpha1
    );

    psiEqn.relax();
    psiEqn.solve();

    Info<< "min " <<psi.name()
        << " " << min(psi).value() << endl;
}
In createFields.H, I added the lines
Code:
    volScalarField psi
    (
        IOobject
        (
            "psi",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedScalar("zero", dimless, 0.0),
    "zeroGradient"
    );
    volScalarField psiSource
    (
        IOobject
        (
            "psiSource",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionedScalar("zero", dimless/dimTime, 0.01),
    "zeroGradient"
    );
and in twoPhaseEulerFoam.C Ii included the psiEqn.H right after the energy equations.

Also the divScheme for "div\(alphaRhoPhi.*,psi\)" was set to "Gauss limitedLinear 1;" (upwind doesn't help) and an entry in fvSolution was done according to the one of the energy equation (smoothSolver etc.). I also tried different solvers but that didn't make a change at all.

When I run the solver on the injection tutorial case, it crashes in the moment, when it should solve psiEqn:

Code:
Starting time loop

Courant Number mean: 0 max: 0
Max Ur Courant Number = 0
Time = 0.005

PIMPLE: iteration 1
MULES: Solving for alpha.air
MULES: Solving for alpha.air
alpha.air volume fraction = 0.293333  Min(alpha1) = 0  Max(alpha1) = 1
smoothSolver:  Solving for e.air, Initial residual = 1, Final residual = 7.75344e-14, No Iterations 1
smoothSolver:  Solving for e.water, Initial residual = 0.00402093, Final residual = 0.00367231, No Iterations 1000
min T.air 300
min T.water 300
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::symGaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int) at ??:?
#4  Foam::symGaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
#5  Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#6  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
#7  Foam::fvMatrix<double>::solve(Foam::dictionary const&) at ??:?
#8  Foam::fvMatrix<double>::solve() at ??:?
#9  
 at ??:?
#10  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11  
 at ??:?
Does anyone have an idea what is wrong here? From what I can read from the error message it complains about signs or possibly division by zero/square root of a negative number,... but I can't figure out the origin of the problem.

Would be very glad to get any help on this! Thanks in advance.

Regards, Alex
derkermit is offline   Reply With Quote

Old   February 9, 2015, 10:12
Default
  #2
New Member
 
Ramon
Join Date: Feb 2014
Location: Eindhoven
Posts: 25
Rep Power: 12
RjwV is on a distinguished road
Hello Alex,

Perhaps I can help you. I am not sure why your method is wrong, but I can push you in the right direction.

twoPhaseEulerFoam makes use of two phases, so your scalar has to be defined for one or both of the phases, i.e. in phaseModel.C and .H, you have to add your scalar similar to 'alpha'. So for example you add an entry 'phi', which should be and IOobject named "psi", and depending on what you want with regard to BCs/ICs you can make it READ_IF_PRESENT, NO_READ, AUTO_WRITE, etc... Don't forget to declare psi in the phaseModel.H file!

After having done this you can add your scalar to one or both of the phases in twoPhaseEulerFoam, look at all the examples in the beginning of the createFields.H file. For example, you can call your psi for phase 1 with the command phase1.psi. Now you should be able to use this psi in your scalar transport equation.

I hope this makes sense to you, if you have any questions let me know.

Kind regards,
Ramon
RjwV is offline   Reply With Quote

Old   February 9, 2015, 11:08
Default
  #3
Member
 
A. Bernath
Join Date: Jun 2011
Location: Karlsruhe, Germany
Posts: 39
Rep Power: 14
derkermit is on a distinguished road
Hi Ramon,
thanks for the help.
Does it make a difference if I declare psi like in my first post compared to the implementation through phaseModel.C/H as you mentioned?
I'm trying that approach at the moment but I'm curious what could be the difference.

Thanks, Alex

EDIT: Actually, the results are the same :-/
derkermit is offline   Reply With Quote

Old   February 9, 2015, 12:49
Default
  #4
New Member
 
Ramon
Join Date: Feb 2014
Location: Eindhoven
Posts: 25
Rep Power: 12
RjwV is on a distinguished road
You could first of all try leaving out any turbulence terms and other complexities, try just a simple convection/diffusion equation, and build it up from there. I know its a little trivial but it that way you can find the part of the equation that is destroying your solver.

K.R.
Ramon
RjwV is offline   Reply With Quote

Old   February 13, 2015, 10:22
Default
  #5
Member
 
Mattia de\' Michieli Vitturi
Join Date: Mar 2009
Posts: 50
Rep Power: 17
demichie is on a distinguished road
Hi,
may I ask why are you using alphaEff coefficient for diffusion?

fvc::interpolate(thermo1.alphaEff(phase1.turbulenc e().mut()))

This coefficient is the thermal diffusion rate (ratio between kinematic visocsity and Prandtl number). Is this what you really want?

Can you also try do add before the solve the following line?

fvOptions.constrain(psiEqn);


Ciao
Mattia
demichie is offline   Reply With Quote

Old   February 16, 2015, 10:51
Default
  #6
Senior Member
 
Gerhard Holzinger
Join Date: Feb 2012
Location: Austria
Posts: 335
Rep Power: 28
GerhardHolzinger will become famous soon enoughGerhardHolzinger will become famous soon enough
Try the PBiCG solver with none as preconditioner.

From the solver output I figure that there are zones in your domain with alpha=0. In these zones the terms of your transport equation are multiplied by zero. Solving the linear equation system involves dividing terms in some cases.

Since the floating point exception happens during the smoothing operation, my guess is, that you run into numerical trouble due to the zero-alpha values in parts of your domain.
GerhardHolzinger is offline   Reply With Quote

Old   November 25, 2015, 13:40
Default
  #7
Member
 
Victor Koppejan
Join Date: May 2015
Posts: 40
Rep Power: 10
vkoppejan is on a distinguished road
I found a much easier option using function objects, that doensn't require you to modify the solver.

please find it in this thread

http://www.cfd-online.com/Forums/ope...tml#post574980
RjwV likes this.
vkoppejan is offline   Reply With Quote

Old   March 7, 2016, 17:38
Default
  #8
New Member
 
Ramon
Join Date: Feb 2014
Location: Eindhoven
Posts: 25
Rep Power: 12
RjwV is on a distinguished road
Great post vkoppejan, this could be very useful for obtaining some quick results. I will definitely add the post to my ever growing list of OpenFOAM tips & tricks.

Keep in mind that this only works for constant isotropic diffusion type problems. So for simple cases this could be the perfect solution, for more advances problems this seems to lack flexibility.

Kind regards,
Ramon
RjwV is offline   Reply With Quote

Old   March 8, 2016, 10:36
Default
  #9
Member
 
Victor Koppejan
Join Date: May 2015
Posts: 40
Rep Power: 10
vkoppejan is on a distinguished road
Hi Ramon,

Your right, it's a very simple solution but the function object can easily adapted to cope with other types of diffusion problems. You just need to make sure the proper files are included and that the makefiles are adapted accordingly.

Glad to have helped you.

Cheers,

Victor
vkoppejan is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Problem diverges: Scalar transport equation, variable: Species transport, material shahrbanoo AVL FIRE 0 July 24, 2014 09:52
turbulent diffusion term in transport equation for additional variables Raijin Thunderkeg CFX 2 May 17, 2014 23:53
Adding temperature equation to twoPhaseEulerFoam nikhil5490 OpenFOAM Programming & Development 0 April 30, 2014 14:28
One transport equation, two user-defined scalar, can it be solved? sharonyue FLUENT 0 April 1, 2014 23:18
Implement species transport equation Tobi OpenFOAM Programming & Development 0 June 2, 2012 14:26


All times are GMT -4. The time now is 06:53.