View Issue Details
|ID||Project||Category||View Status||Date Submitted||Last Update|
|0002589||OpenFOAM||[All Projects] Bug||public||2017-06-22 18:33||2017-07-06 14:51|
|Fixed in Version||dev|
|Summary||0002589: The use of the IshiiZuber dragModel in the bubbleColumn tutorial for reactingMultiphaseEulerFoam results in a sigFpe|
|Description||The 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.
IshiiZuber.C (2,997 bytes)
bubbleColumn.tar.gz (3,859 bytes)
||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.|
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.
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?
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
||Right, but this is not an ideal solution. Have you found any conservation issues following this change?|
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.)
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.
> 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.
||Resolved by commit 5caadae42b73724784934142024612c8f1cdc9e7|
||Could you include the fix in 4.x as well? Thank you very much.|
||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?|
||No I didn't run a range of cases so far. I understand your concern and will report on this if I experience trouble.|
|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|