View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0001576OpenFOAMBugpublic2015-03-19 13:532015-03-19 21:46
ReportercraigWhite Assigned Tohenry  
PriorityhighSeveritymajorReproducibilityalways
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSUbuntuOS Version12.04
Summary0001576: Bug in DSMC LarsenBorgnakkeVariableHardSphere class
DescriptionIn the source code for the LarsenBorgnakkeVariableHardSphere class, there is a small bug that makes the results incorrect for more than 2 degrees of rotational freedom.

  230 if (iDofP > 0)
  231 {
  232 if (inverseCollisionNumber > rndGen.scalar01())
  233 {
  234 availableEnergy += preCollisionEiP;
  235
  236 scalar ChiA = 0.5*iDofP;
  237
  238 EiP = energyRatio(ChiA, ChiB)*availableEnergy;
  239
  240 availableEnergy -= EiP;
  241 }
  242 }
  243
  244 if (iDofQ > 0)
  245 {
  246 if (inverseCollisionNumber > rndGen.scalar01())
  247 {
  248 availableEnergy += preCollisionEiQ;
  249
  250 // Change to general LB ratio calculation
  251 scalar energyRatio = 1.0 - pow(rndGen.scalar01(),(1.0/ChiB));
  252
  253 EiQ = energyRatio*availableEnergy;
  254
  255 availableEnergy -= EiQ;
  256 }
  257 }

The loop for particle P (lines 230-242) is ok. But you can see it's different for particle Q, and in fact a comment has been left in to change it to the general ratio calculation, but it hasn't been performed. This means it treats particles as if they always have 2 degrees of freedom even if they have more.

I think it just needs to be changed to:

  230 if (iDofP > 0)
  231 {
  232 if (inverseCollisionNumber > rndGen.scalar01())
  233 {
  234 availableEnergy += preCollisionEiP;
  235
  236 scalar ChiA = 0.5*iDofP;
  237
  238 EiP = energyRatio(ChiA, ChiB)*availableEnergy;
  239
  240 availableEnergy -= EiP;
  241 }
  242 }
  243
  244 if (iDofQ > 0)
  245 {
  246 if (inverseCollisionNumber > rndGen.scalar01())
  247 {
  248 availableEnergy += preCollisionEiQ;
  249
  250 scalar ChiA = 0.5*iDofQ;
  251
  252 EiQ = energyRatio(ChiA, ChiB)*availableEnergy;
  253
  254 availableEnergy -= EiQ;
  255 }
  256 }
Additional InformationYou could even add a statement so that if iDofP or iDofQ is 2, then it will use the:

scalar energyRatio = 1.0 - pow(rndGen.scalar01(),(1.0/ChiB));

relation, which is correct for 2 degrees of freedom (see Bird's 1994 book, equation 5.45 and surrounding text). Otherwise it should go to the general relation as particle P currently does, but Q does not.
TagsNo tags attached.

Activities

henry

2015-03-19 14:01

manager   ~0004150

I agree this is odd. I am happy to implement either of the correction you propose. Is there any particular advantage in using the specifically 2D model for 2D or is the general model a good choice for anyD?

craigWhite

2015-03-19 14:13

reporter   ~0004151

For diatomic molecules (which have 2 degrees of rotational freedom), it might mean less floating point operations, but I suppose it does mean performing more IF operations, so it perhaps doesn't really matter too much. The general model should work fine for any number of rotational degrees of freedom, i.e. diatomics, triatomics, polyatomics. I have a version of the DSMC code that I have heavily modified, but this is an example of what I have in mine and I can confirm that it works:

    if (rotationalDofP > VSMALL)
    {
        if (inverseRotationalCollisionNumber > rndGen.scalar01())
        {
            scalar EcP = translationalEnergy + preCollisionERotP;
            
            scalar energyRatio = 0.0;
            
            if(rotationalDofP == 2.0)
            {
                energyRatio = 1.0 - pow(rndGen.scalar01(),(1.0/ChiB));
            }
            else
            {
                scalar ChiA = 0.5*rotationalDofP;
                
                energyRatio = cloud_.energyRatio(ChiA, ChiB);
            }

            ERotP = energyRatio*EcP;
        
            translationalEnergy = EcP - ERotP;
        }
    }

I've known about this for a while (I found it when benchmarking my code for polyatomic species), should have reported it earlier, sorry. It was preventing the rotational mode from reaching equilibrium.

henry

2015-03-19 14:24

manager   ~0004152

Why is rotationalDofP a floating-point number? In you modelling do you consider degrees of freedom a continuous function? If so the test

if(rotationalDofP == 2.0)

is a bit risky as rotationalDofP may never be exactly 2.0.

craigWhite

2015-03-19 14:31

reporter   ~0004153

It's a constant property of the gas species and is defined in the dsmcProperties dictionary before the simulation begins. A diatomic molecule always has exactly 2 degrees of rotational freedom, it can only spin around its bond in two planes.

craigWhite

2015-03-19 14:34

reporter   ~0004154

scalar iDofP = cloud.constProps(typeIdP).internalDegreesOfFreedom();

henry

2015-03-19 14:40

manager   ~0004155

I note that currently in OpenFOAM internalDegreesOfFreedom_ is also a scalar, is this necessary? Is there a need to be able to set internalDegreesOfFreedom_ to something other than 0, 1, 2 or 3? In order to implement a reliable test for
internalDegreesOfFreedom_ == 2 it would be MUCH better if internalDegreesOfFreedom_ were an integer.

craigWhite

2015-03-19 14:52

reporter   ~0004156

My understanding is that no matter their size, a linear molecule will always have 2 rotational degrees of freedom, a non-linear molecule will have 3 rotational degrees of freedom and an atom obviously has zero. So, really, this value should only ever be 0, 2, or 3 and so it could easily be an integer rather than a scalar value.

Here's a chemistry textbook explaining it (page 234): https://books.google.co.uk/books?id=zzxLTIljQB4C&pg=PA234&lpg=PA234&dq=5+rotational+degrees+of+freedom&source=bl&ots=wbaPHKmZjM&sig=NX-dR2asRc_dVG3kM_KSGPH6D4A&hl=en&sa=X&ei=5OEKVZGyLcTi7AbttIDAAw&ved=0CD0Q6AEwBw#v=onepage&q=5%20rotational%20degrees%20of%20freedom&f=false

henry

2015-03-19 14:57

manager   ~0004158

This is also my understanding, I just wanted to check with you that there isn't some special DSMC modeling in which the degrees of freedom is non-integer.

I will change the code so that internalDegreesOfFreedom_ is an integer and add the 2D condition.

henry

2015-03-19 21:46

manager   ~0004162

Resolved by commit ec68d3384e695a4f3ad55a6fa0367ce3f341d45c in OpenFOAM-dev

If you find any issues with the fix to the LarsenBorgnakkeVariableHardSphere or the change of internalDegreesOfFreedom to an integer type please reopen this report.

Issue History

Date Modified Username Field Change
2015-03-19 13:53 craigWhite New Issue
2015-03-19 14:01 henry Note Added: 0004150
2015-03-19 14:13 craigWhite Note Added: 0004151
2015-03-19 14:24 henry Note Added: 0004152
2015-03-19 14:31 craigWhite Note Added: 0004153
2015-03-19 14:34 craigWhite Note Added: 0004154
2015-03-19 14:40 henry Note Added: 0004155
2015-03-19 14:52 craigWhite Note Added: 0004156
2015-03-19 14:57 henry Note Added: 0004158
2015-03-19 21:46 henry Note Added: 0004162
2015-03-19 21:46 henry Status new => resolved
2015-03-19 21:46 henry Resolution open => fixed
2015-03-19 21:46 henry Assigned To => henry
2015-03-24 00:17 liuhuafei Issue cloned: 0001599