Description: At Time = 0.02
 

Description:

At Time = 0.0278082 in the default damBreak4phaseFine case, we get this message in the stderr output. (Full log file is attached as well: http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif log.tar.bz2 )

#0 Foam::error::printStack(Foam:http://www.cfd-online.com/OpenFOAM_D...part/proud.gifstream&) in "/home/madhavan/OpenFOAM/OpenFOAM-1.4.1/lib/linux64GccDPOpt/libOpenFOAM.so"
#1 Foam::sigFpe::sigFpeHandler(int) in "/home/madhavan/OpenFOAM/OpenFOAM-1.4.1/lib/linux64GccDPOpt/libOpenFOAM.so"
#2 __restore_rt in "/lib64/tls/libc.so.6"
#3 Foam::GAMGSolver::scalingFactor(Foam::Field<double >&, Foam::Field<double> const&, Foam::Field<double> const&, Foam::Field
<double> const&) const in "/home/madhavan/OpenFOAM/OpenFOAM-1.4.1/lib/linux64GccDPOpt/libOpenFOAM.so"
#4 Foam::GAMGSolver::scalingFactor(Foam::Field<double >&, Foam::lduMatrix const&, Foam::Field<double>&, Foam::FieldField<Foam
::Field, double> const&, Foam::UPtrList<foam::lduinterfacefield> const&, Foam::Field<double> const&, unsigned char) con
st in "/home/madhavan/OpenFOAM/OpenFOAM-1.4.1/lib/linux64GccDPOpt/libOpenFOAM.so"
#5 Foam::GAMGSolver::Vcycle(Foam::PtrList<foam::lduma trix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const
&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<foam::field<double> >&, Foam::PtrList<Foam
::Field<double> >&, unsigned char) const in "/home/madhavan/OpenFOAM/OpenFOAM-1.4.1/lib/linux64GccDPOpt/libOpenFOAM.so"
#6 Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/home/madhavan/OpenFOA
M/OpenFOAM-1.4.1/lib/linux64GccDPOpt/libOpenFOAM.so"
#7 Foam::fvMatrix<double>::solve(Foam::Istream&) in "/home/madhavan/OpenFOAM/OpenFOAM-1.4.1/lib/linux64GccDPOpt/libfiniteVol
ume.so"
#8 main in "/home/madhavan/OpenFOAM/OpenFOAM-1.4.1/applications/bin/linux64GccDPOpt/multiphaseInterFoam"
#9 __libc_start_main in "/lib64/tls/libc.so.6"
#10 Foam::regIOobject::readIfModified() in "/home/madhavan/OpenFOAM/OpenFOAM-1.4.1/applications/bin/linux64GccDPOpt/multipha
seInterFoam"

Solver/Application:
multiphaseInterFoam (OpenFOAM Version 1.4.1)

Source file:
Not sure. Check the diff output (posted below) from multiphaseInterFoam source (OF v 1.4) versus that from OF v 1.4.1.

Testcase:
damBreak4phaseFine. Interestingly the damBreak4phase case seems to run OK in OF 1.4.1.

Platform:
So far linux64.

Version:
(OpenFOAM version 1.4.1) OF Version 1.4 seems to be unaffected. The case runs without problems there.

Notes:

Related discussion thread: http://www.cfd-online.com/OpenFOAM_D...ges/1/315.html

For some reason, Nobody reported this bug http://www.cfd-online.com/OpenFOAM_D...lipart/sad.gif

Source code diff output:

[madhavan@tohirdhurd ~]$ diff multiphaseInterFoam.OF1.4/createFields.H multiphaseInterFoam.OF1.4.1/createFields.H

[madhavan@tohirdhurd ~]$ diff multiphaseInterFoam.OF1.4/multiphaseInterFoam.C multiphaseInterFoam.OF1.4.1/multiphaseInterFoam.C

[madhavan@tohirdhurd ~]$ diff multiphaseInterFoam.OF1.4/pEqn.H multiphaseInterFoam.OF1.4.1/pEqn.H

[madhavan@tohirdhurd ~]$ diff multiphaseInterFoam.OF1.4/UEqn.H multiphaseInterFoam.OF1.4.1/UEqn.H

[madhavan@tohirdhurd ~]$ diff multiphaseInterFoam.OF1.4/multiphaseMixture/multiphaseMixture.C multiphaseInterFoam.OF1.4.1/multiphaseMixture/multiphaseMixture.C
344c344
< fvPatchVectorFieldField& nHatb
---
> surfaceVectorField::GeometricBoundaryField& nHatb
347c347
< const fvPatchScalarFieldField& gbf = refPhase_.boundaryField();
---
> const volScalarField::GeometricBoundaryField& gbf = refPhase_.boundaryField();

[madhavan@tohirdhurd ~]$ diff multiphaseInterFoam.OF1.4/multiphaseMixture/multiphaseMixture.H multiphaseInterFoam.OF1.4.1/multiphaseMixture/multiphaseMixture.H
51d50
< #include "fvPatchFieldFields.H"
205c204
< fvPatchVectorFieldField& nHatb
---
> surfaceVectorField::GeometricBoundaryField& nHatb

[madhavan@tohirdhurd ~]$ diff multiphaseInterFoam.OF1.4/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.C multiphaseInterFoam.OF1.4.1/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.C
79c79
< const scalarField& iF
---
> const DimensionedField<scalar,>& iF
90c90
< const scalarField& iF,
---
> const DimensionedField<scalar,>& iF,
102c102
< const scalarField& iF,
---
> const DimensionedField<scalar,>& iF,
116c116
< const scalarField& iF
---
> const DimensionedField<scalar,>& iF
126d125
< // Write

[madhavan@tohirdhurd ~]$ diff multiphaseInterFoam.OF1.4/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.H multiphaseInterFoam.OF1.4.1/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.H
145c145
< const scalarField&
---
> const DimensionedField<scalar,>&
152c152
< const scalarField&,
---
> const DimensionedField<scalar,>&,
162c162
< const scalarField&,
---
> const DimensionedField<scalar,>&,
179c179
< const scalarField&
---
> const DimensionedField<scalar,>&
183c183,186
<virtual> clone(const scalarField& iF) const
---
> virtual tmp<fvpatchscalarfield> clone
> (
> const DimensionedField<scalar,>& iF
> ) const

[madhavan@tohirdhurd ~]$ diff multiphaseInterFoam.OF1.4/multiphaseMixture/phase/phase.C multiphaseInterFoam.OF1.4.1/multiphaseMixture/phase/phase.C
25c25
< \*----------------------------------------------------------------------------*/
---
> \*---------------------------------------------------------------------------*/

[madhavan@tohirdhurd ~]$ diff multiphaseInterFoam.OF1.4/multiphaseMixture/phase/phase.H multiphaseInterFoam.OF1.4.1/multiphaseMixture/phase/phase.H

Addendum: The case of the phas
 

Addendum: The case of the phases needs to be changed in setFieldsDict before running the case. For example, alphaAir needs to be changes to alphaair, alphaMercury to alphamercury and so on. But this isn't the bug. This is just a minor typo correction.

Addendum2: The terminal output
 

Addendum2: The terminal output also stated the following when the solver stopped:

"Floating point exception"

I have found that the fine-mes
 

I have found that the fine-mesh case does not run reliably in either 1.4 or 1.4.1 whereas the coarse mesh case runs fine in both. The problem relates to the phase fraction unboundedness possible with 2nd-order spatial schemes which if larger than the density ratio causes the code to fail. I have solved this problem in the 2-phase VOF codes by developing the MULES solver which guarantees boundedness for well-posed problems but have not yet completed the multi-phase version of this; I hope to have done so for the next release.

In the meantime you can get the fine mesh case to run if you reduce the compressive effect of the scheme used for the "anti-diffusion" term slightly. This can be done by switching from the quartic to the quadratic scheme in OpenFOAM-1.4.1/src/transportModels/interfaceProperties/interfaceCompression/inte rfaceCompression.H
and recompiling the transportModels library. If switching between these schemes is found to be generally useful I will make it run-time selectable in the next release.

Thanks Henry, your suggestion
 

Thanks Henry, your suggestion seems to work. No problems so far. I do have one question however. How do you set the alphas field (i.e. all volume fractions in 1 file) in the 0/ directory?

You don't need to have an alph
 

You don't need to have an alphas file as it is generated automatically by the code and used for post-processing only.

Thanks Henry.
 

Thanks Henry.

Does this problem only relates
 

Does this problem only relates to multiphaseInterFoam solver or does it also affects other interFoam solvers including rasInterFoam?

Krystian

I have solved this problem in
 

I have solved this problem in the 2-phase VOF codes including rasInterFoam by developing the MULES solver which guarantees boundedness for well-posed problems but I have not yet completed the multi-phase version of this; I hope to have done so for the next release.

Henry - do you have a validati
 

Henry - do you have a validation test case that you are using to judge the performance of the multiphaseInterFoam solver, and if so could you please share the reference?

I am guessing that the tutorial test case (an extension of the dam break problem, with three liquids) is an arbitrary problem, chosen simply as a demonstration.

I have not run any rigorous va
 

I have not run any rigorous validation cases as this development was more of a demonstration and was not part of a funded project. Further development is needed, in particular the extension of MULES for multiphase problems and then further validation. Do you have any suggestions for validation cases?

Henry

Understood. I wonder what the
 

Understood. I wonder what the commercial code vendors use to validate their models ... I'll investigate.

I've met the same problem when
 

I've met the same problem when I try to run a custom case set by myself. Unfortunately although I changed the scheme from quartic to quadratic as Henry recommended, the problem still exists.

My case is a 2-D wave tank which has approximately 60000 cells, I just want to employ a viscoid phase in order to absorb waves at the far end of the tank, so I modified the damBreak4phase case, changed blockMeshDict and re-blockMesh and setFields. The program exited in no more than 3 time steps and reported the same error as Srinath met.

by the way, my wave tank based on interFoam works quite well, so is it possible to set up an absorption zone where water has large viscosity in it?

In "multiphaseMixture/multipha
 

In "multiphaseMixture/multiphaseMixture.C" I believe there is an initialization
bug in

Foam::tmp<foam::volscalarfield> Foam::multiphaseMixture::rho()


The line:

tmp<volscalarfield> trho = iter()*iter().rho();

should be replaced by something like

tmp<volscalarfield> trho = 0.0*iter()*iter().rho();

to ensure correct initialization. Currently the density is approximately double
what is intended.

Regards,

RGK.

The plan was to initialize trh
 

The plan was to initialize trho with the density of the first phase and than add the other phases but it looks like I forgot to increment the iterator after the initialization. Try

tmp<volscalarfield> trho = iter()*iter().rho();
++iter;

Henry

or slightly neater tmp<
 

or slightly neater

tmp<volscalarfield> trho = iter()*iter().rho();

for(++iter; iter != phases_.end(); ++iter)

and the same goes for the mu() and muf() functions.

H

Many thanks for the clarificat
 

Many thanks for the clarification and elegant solution.

I only spotted the issue with rho() on account of the extreme
densities I was using, mu() and muf() slipped under the radar.


RGK