View Issue Details

IDProjectCategoryView StatusLast Update
0002589OpenFOAM[All Projects] Bugpublic2017-07-06 14:51
Reporteroertel59Assigned Tohenry 
Status closedResolutionfixed 
PlatformGNU/LinuxOSUbuntuOS Version14.04
Product Versiondev 
Fixed in Versiondev 
Summary0002589: The use of the IshiiZuber dragModel in the bubbleColumn tutorial for reactingMultiphaseEulerFoam results in a sigFpe
DescriptionThe sigFpe occurs within IshiiZuber.C, line 87. Apparently the sqrt() function gets a slightly negative value. This can be mitigated by changing the line from

    volScalarField F((muc/muMix)*sqrt(1 - pair_.dispersed()));
    volScalarField F((muc/muMix)*sqrt(1 - min(pair_.dispersed(), 1.0)));

as done in the attached patch.

The error does not occur for the analogue situation in reactingTwoPhaseEulerFoam.



2017-06-22 18:33


IshiiZuber.C (2,997 bytes)


2017-06-22 18:34


bubbleColumn.tar.gz (3,859 bytes)


2017-06-26 09:18

reporter   ~0008252

The root cause is that pair_.dispersed() > 1, which is unphysical. It might be better to ensure that phase fractions are bounded to 0-1 in the multiPhaseSystem, instead of fixing the symptoms of non-physical values one by one in various submodels.


2017-06-26 09:29

manager   ~0008254

MULES can guarantee boundedness only to within round-off error and the accuracy of the flux fields. To absolutely guarantee boundedness to avoid floating-point exception in sqrt etc. clipping would need to be applied. If this is applied to the phase-fraction field following solution there is the possibility of accumulating a conservation error which would naturally correct itself out if a modest level of unboundedness were permitted. The alternatives are to clip the phase-fraction field everywhere that absolute boundedness is required or store a clipped set of phase-fraction fields and use those where necessary.

I am not particularly happy with any of the posible solutions to this problem.


2017-06-26 09:48

reporter   ~0008255

To me this is nothing else than catching a division by zero by max(value, SMALL) in the denominator and I'd be happy with a "fixing the symptoms" solution.

Following Henry's last solution, what about a getter function for the phase fraction with a switch for the clipping that defaults to none, either supplying a stored clipped field or calculating it when the switch evaluates to true?


2017-06-26 10:39

reporter   ~0008256

Last time this was discussed with reactingTwoPhaseEulerFoam we ended up clipping the phase fractions. Current reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C implementation:

        // Ensure the phase-fractions are bounded


2017-06-26 10:48

manager   ~0008257

Right, but this is not an ideal solution. Have you found any conservation issues following this change?


2017-06-26 11:15

reporter   ~0008259

Not related to clipping, but I do try to set up my simulations in a way that significant clipping doesn't occur.

I believe the previous conclusion was to implement the simple solution (clipping) to see if a better solution is needed.

It might be a good idea to add a warning for the user when clipping occurs and print out integral of the correction so that the user can evaluate its significance. However, this was considered unnecessary in the previous discussion.

(The "max(1.0 - alpha1[celli], 1e-4)" dgdt source limiter can cause significant conservation issues, but that is a different subject.)


2017-06-26 11:56

manager   ~0008261

Clipping to at least round-off error will probably occur every time-step and unboundedness is already reported. It would be possible to provide a clipping error measure but no mechanism to avoid it could be recommended and it would not be clear if removing the clipping would improve conversation or adversely affect stability. It would be possible to implement the clipping on a switch but unless clipping were also added to the models which require absolute boundedness they would not run without it.

For now I will implement clipping into reactingMultiphaseEulerFoam for consistency with reactingTwoPhaseEulerFoam and we can evaluate the consequence and review as required.

> (The "max(1.0 - alpha1[celli], 1e-4)" dgdt source limiter can cause significant conservation issues, but that is a different subject.)

Yes but this relates to a mass rather than a phase-fraction error. It is not clear how this can be avoided other than reducing the "1e-4" which will also affect stability. This could certainly be made an input.


2017-06-26 12:16

reporter   ~0008262

> Yes but this relates to a mass rather than a phase-fraction error. It is not clear how this can be avoided other than reducing the "1e-4" which will also affect stability. This could certainly be made an input.

In CFBs we've successfully used 1e-6 to get rid of the conservation errors, but that cannot be recommended for general use. It would be good to make it a user input.


2017-06-26 16:27

manager   ~0008266

Resolved by commit 5caadae42b73724784934142024612c8f1cdc9e7


2017-07-05 09:32

reporter   ~0008335

Could you include the fix in 4.x as well? Thank you very much.


2017-07-05 09:40

manager   ~0008336

Because this fix may affect conservation and have other side-effects I didn't want to include it in the patch version of the current release, at least not until it has been sufficiently tested. Have you run a range of cases with this change? Have you seen any problems?


2017-07-06 14:26

reporter   ~0008351

No I didn't run a range of cases so far. I understand your concern and will report on this if I experience trouble.

Issue History

Date Modified Username Field Change
2017-06-22 18:33 oertel59 New Issue
2017-06-22 18:33 oertel59 File Added: IshiiZuber.C
2017-06-22 18:33 oertel59 Tag Attached: reactingMultiphaseEulerFoam
2017-06-22 18:34 oertel59 File Added: bubbleColumn.tar.gz
2017-06-26 09:18 Juho Note Added: 0008252
2017-06-26 09:29 henry Note Added: 0008254
2017-06-26 09:48 oertel59 Note Added: 0008255
2017-06-26 10:39 Juho Note Added: 0008256
2017-06-26 10:48 henry Note Added: 0008257
2017-06-26 11:15 Juho Note Added: 0008259
2017-06-26 11:56 henry Note Added: 0008261
2017-06-26 12:16 Juho Note Added: 0008262
2017-06-26 16:27 henry Assigned To => henry
2017-06-26 16:27 henry Status new => resolved
2017-06-26 16:27 henry Resolution open => fixed
2017-06-26 16:27 henry Fixed in Version => dev
2017-06-26 16:27 henry Note Added: 0008266
2017-07-05 09:32 oertel59 Status resolved => feedback
2017-07-05 09:32 oertel59 Resolution fixed => reopened
2017-07-05 09:32 oertel59 Note Added: 0008335
2017-07-05 09:40 henry Note Added: 0008336
2017-07-06 14:26 oertel59 Note Added: 0008351
2017-07-06 14:26 oertel59 Status feedback => assigned
2017-07-06 14:51 henry Status assigned => closed
2017-07-06 14:51 henry Resolution reopened => fixed