phorgue / porousMultiphaseFoam

A porous multiphase toolbox for OpenFOAM
Other
124 stars 64 forks source link

Bug in `CoatsNo.H` in anisotropic version #7

Closed FoamScience closed 5 years ago

FoamScience commented 5 years ago

I believe we need to divide on mesh.V() only once here:

https://github.com/phorgue/porousMultiphaseFoam/blob/bd371c936d05cffd640845e3329b0c5639a753cc/solvers/anisoImpesFoam/CoatsNo.H#L18-L25

This should be done in this order (With line 18 postponed after capillarity is taken into account)

// - capillarity part of CFL 
if(activateCapillarity) 
{ 
     CFLCoats += (runTime.deltaT()/eps)*2*mag(pcModel->dpcdS())*fvc::surfaceSum(Kf*mesh.magSf()/mag(mesh.delta()))*(kra*krb/(mub*kra+mua*krb)); 
}
CFLCoats.ref() /= mesh.V(); 

Let me if I'm missing anything...

Also, from what I understand, the lines https://github.com/phorgue/porousMultiphaseFoam/blob/bd371c936d05cffd640845e3329b0c5639a753cc/solvers/anisoImpesFoam/CoatsNo.H#L28-L40 are meant to find maximal CFL value in all tensor directions. Yeah, this can be retrieved with

// I'm assuming CFLCoats is a tensorField, didn't test with it being volTensorField but I imagine .ref() would then work :)
// gMax returns a tensor of component-wise maximum values for CFLCoats in all cells
// cmptMax then returns maximal algebraic value (not in magnitude) in the tensor
// which should represent the global max in all tensor directions.
CFLUse = cmptMax(gMax(CFLCoats));
phorgue commented 5 years ago

Thanks for finding and correct the bug.