
[Sponsors] 
February 7, 2015, 09:04 
Adding transport equation to twoPhaseEulerFoam (OF231)

#1 
Member
A. Bernath
Join Date: Jun 2011
Location: Karlsruhe, Germany
Posts: 39
Rep Power: 6 
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; } 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" ); 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.75344e14, 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_64linuxgnu/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_64linuxgnu/libc.so.6" #11 at ??:? Would be very glad to get any help on this! Thanks in advance. Regards, Alex 

February 9, 2015, 10:12 

#2 
New Member
Ramon
Join Date: Feb 2014
Location: Eindhoven
Posts: 24
Rep Power: 3 
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 

February 9, 2015, 11:08 

#3 
Member
A. Bernath
Join Date: Jun 2011
Location: Karlsruhe, Germany
Posts: 39
Rep Power: 6 
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 :/ 

February 9, 2015, 12:49 

#4 
New Member
Ramon
Join Date: Feb 2014
Location: Eindhoven
Posts: 24
Rep Power: 3 
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 

February 13, 2015, 10:22 

#5 
New Member
Mattia de\' Michieli Vitturi
Join Date: Mar 2009
Posts: 23
Rep Power: 8 
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 

February 16, 2015, 10:51 

#6 
Senior Member
Gerhard Holzinger
Join Date: Feb 2012
Location: Austria
Posts: 170
Rep Power: 14 
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 zeroalpha values in parts of your domain. 

Thread Tools  
Display Modes  


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 08:52 
turbulent diffusion term in transport equation for additional variables  Raijin Thunderkeg  CFX  2  May 17, 2014 22:53 
Adding temperature equation to twoPhaseEulerFoam  nikhil5490  OpenFOAM Programming & Development  0  April 30, 2014 13:28 
One transport equation, two userdefined scalar, can it be solved?  sharonyue  FLUENT  0  April 1, 2014 22:18 
Implement species transport equation  Tobi  OpenFOAM Programming & Development  0  June 2, 2012 13:26 