View Issue Details

IDProjectCategoryView StatusLast Update
0002544OpenFOAM[All Projects] Bugpublic2017-05-09 10:06
Reporteroertel59Assigned Towill 
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSUbuntuOS Version14.04
Product Versiondev 
Fixed in Version 
Summary0002544: reactingTwoPhaseEulerFoam/RAS/bubbleColumn tutorial crashes immediately if the Burns turbulent dispersion is used
DescriptionCopying 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.
TagsreactingtwoPhaseEulerFoam

Activities

oertel59

2017-05-08 11:56

reporter  

bubbleColumn_Burns.tar.gz (1,529,615 bytes)

oertel59

2017-05-09 08:25

reporter   ~0008125

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

Proposal:

    return
        0.75
       *Ctd_
       *drag.CdRe()
       *pair_.continuous().nu()
       *pair_.continuous().turbulence().nut()
       /(
            sigma_
           *sqr(pair_.dispersed().d())
        )
       *pair_.continuous().rho()
       *pair_.dispersed()
       *(
            1.0/max(pair_.dispersed(), residualAlpha_)
          + 1.0/max(pair_.continuous(), residualAlpha_)
        );

Original version:

    return
        0.75
       *Ctd_
       *drag.CdRe()
       *pair_.continuous().nu()
       *pair_.continuous().turbulence().nut()
       /(
            sigma_
           *sqr(pair_.dispersed().d())
        )
       *pair_.continuous().rho()
       *(1.0 + pair_.dispersed()/max(pair_.continuous(), residualAlpha_));

which is mathematically equivalent, but more stable at startup for cases with phase inversion.

will

2017-05-09 08:52

manager   ~0008126

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.

Juho

2017-05-09 09:24

reporter   ~0008128

The modified version is zero when pair_dispersed() is zero. The original is 1 when pair_dispersed() is zero.

will

2017-05-09 09:25

manager   ~0008129

Resolved by commit 80c101b14f71d0ddd140fe7f9832d4b27524df1e

oertel59

2017-05-09 09:31

reporter   ~0008130

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

will

2017-05-09 09:39

manager   ~0008131

@Juho So it is. The modification is preventing a finite force being applied to a (near) zero field, then. Thanks.

will

2017-05-09 10:06

manager   ~0008134

@oertel59 Yes, this is OK for 4.x. Commit cf5b0901eadc288336ae7e7e66ae7d3e69c8c4ce.

Issue History

Date Modified Username Field Change
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