View Issue Details
|ID||Project||Category||View Status||Date Submitted||Last Update|
|0002544||OpenFOAM||[All Projects] Bug||public||2017-05-08 11:56||2017-05-09 10:06|
|Fixed in Version|
|Summary||0002544: reactingTwoPhaseEulerFoam/RAS/bubbleColumn tutorial crashes immediately if the Burns turbulent dispersion is used|
|Description||Copying the reactingTwoPhaseEulerFoam/RAS/bubbleColumn tutorial and just adding the Burns turbulent dispersion force will result in a crash at the second timestep. This happens independently of the choice of dragModels and the blending setup for the turbulent dispersion.|
The problem appears to be the phase inversion. A large phaseFraction gradient appears at the interface. The value calculated in
//- Turbulent diffusivity
// multiplying the gradient of the phase-fraction
virtual tmp<volScalarField> D() const;
is very high in those regions. The case chrashes with a sigFpe in the dragModel choosen for (water in air).
The case be run for 10s without turbulent dispersion, giving a properly initialized field. If the turbulent dispersion is switched on starting from 10s, it does not chrash.
The proposed fix is to pull out the dispersed phase fraction from the expression within the brackets back to the original formulation of Burns et al. (Equation 27) as attached (see Burns.C)
Writing out D in reactingTwoPhaseEulerFoam/pU/pEqn (line 39) using the original and the new implementation gives identical values (see bubbleColumn.origBurns and bubbleColumn.newBurns).
With the "new" implementation however, the case runs fine starting from 0s.
bubbleColumn_Burns.tar.gz (1,529,615 bytes)
After reading my own report again, I realized that I should have been more clear about what I changed by adding the relevant code directly in the report description. The function returning the turbulent diffusivity within the Burns model was simply rewritten as follows
+ 1.0/max(pair_.continuous(), residualAlpha_)
*(1.0 + pair_.dispersed()/max(pair_.continuous(), residualAlpha_));
which is mathematically equivalent, but more stable at startup for cases with phase inversion.
Thank you for the report. I can confirm that this change makes the case in question stable from time zero. It is also mathematically equivalent, so I have no problem with it going in.
I find it a bit bizarre that it works, though. I distinctly remember writing some models in the (1 + alpha_d/alpha_c) form rather than the alpha_d*(1/alpha_d + 1/alpha_c) form in order to prevent an additional division and hopefully improve stability. It reality, the latter form is more stable. I find this counter-intuitive. If anyone can come up with an explanation as to why this works I would be interested to hear it.
||The modified version is zero when pair_dispersed() is zero. The original is 1 when pair_dispersed() is zero.|
||Resolved by commit 80c101b14f71d0ddd140fe7f9832d4b27524df1e|
||Could you include the fix in 4.x as well?|
||@Juho So it is. The modification is preventing a finite force being applied to a (near) zero field, then. Thanks.|
||@oertel59 Yes, this is OK for 4.x. Commit cf5b0901eadc288336ae7e7e66ae7d3e69c8c4ce.|
|2017-05-08 11:56||oertel59||New Issue|
|2017-05-08 11:56||oertel59||File Added: bubbleColumn_Burns.tar.gz|
|2017-05-08 11:56||oertel59||Tag Attached: reactingtwoPhaseEulerFoam|
|2017-05-09 08:25||oertel59||Note Added: 0008125|
|2017-05-09 08:52||will||Note Added: 0008126|
|2017-05-09 09:24||Juho||Note Added: 0008128|
|2017-05-09 09:25||will||Assigned To||=> will|
|2017-05-09 09:25||will||Status||new => resolved|
|2017-05-09 09:25||will||Resolution||open => fixed|
|2017-05-09 09:25||will||Note Added: 0008129|
|2017-05-09 09:31||oertel59||Status||resolved => feedback|
|2017-05-09 09:31||oertel59||Resolution||fixed => reopened|
|2017-05-09 09:31||oertel59||Note Added: 0008130|
|2017-05-09 09:39||will||Note Added: 0008131|
|2017-05-09 10:06||will||Status||feedback => resolved|
|2017-05-09 10:06||will||Resolution||reopened => fixed|
|2017-05-09 10:06||will||Note Added: 0008134|