Is the particle mass in the SQP method calculated based on the boundary sphere diameter of the non-spherical particles or on other parameters?

Hello, everyone!
I found a problem. In the SQP method, the total mass of the non-spherical particles is set to 3kg, and the number of particles displayed is 1856, but in the actual calculation, the number of particles of 3kg is 5161. I want to ask, which particle do we need?
image
image
The particle is 1856, which is calculated based on the boundary sphere diameter of non-spherical particles, but the particle is 5161, which is calculated based on the equivalent diameter of non-spherical particles ? The radius is calculated to be 0.00905 m, but the equivalent diameter of the non-spherical particles is 0.00941 m.
New0.15.zip (41.2 MB)
Thank you!

You are right, the “Estimated # of particles” displayed in the initial conditions does not take non-spherical particles into account. It uses d_p0 for the calculation, which for spheres is the diameter, but for SQP is the bounding diameter.

Your particle is a fairly flat disk, so it takes up a lot less space than the bounding sphere.

So this leads to a low estimate. Note that this number is only “Estimated”, but it should not be so far off! I will fix this for the next MFiX release.

But, I’m not sure where you get the numbers 0.00905 and 0.00941 from.

The bounding diameter is 0.01271 and so the bounding radius is half of that, 0.0063591.

For the equivalent radius, I refer to formula 2.60 in this text, where the volume of a superquadric is computed as:

V = 2a₁a₂a₃ε₁ε₂ B( ε₁/2 + 1, ε₁)B( ε₂/2 , ε₂/2)

Here B is the mathematical beta function, and the parameters are related to our the a,b,c,m,n used in MFiX as follows:

The aᵢ’s are the semiaxes, which we denote a,b,c
ε₁ = 2/n
ε₂ = 2/m

We compute in Python, using the ‘scipy’ library for the beta function (you may not have this installed, but it’s in conda-forge so it’s easy to install: conda -c conda-forge install scipy)

$ python
>>> from scipy.special import beta
>>> m, n = 2, 4
>>> e1, e2 = 2/n, 2/m
>>> a1 = a2 = 0.00635
>>> a3 = 0.00175
# Computing the volume
>>> V = 2*a1*a2*a3*e1*e2*beta( (e1/2) + 1,  e1)*beta(e2/2, e2/2)
>>> V
3.8751305057030724e-07
# Computing the equivalent radius
>>> from math import pi
>>> r3 = 3*V / (4*pi)
>>> r3 ** (1/3)
0.004522715206706299
# Computing the number of SQP particles in 3.0kg, given a density of 1500kg/m³
>>> 3.0 / (V*1500)
5161.116501899943

So the 5161 number is correct, and the incorrect estimated 1856 comes from using the bounding radius to compute the volume, instead of the correct formula as above.

Here’s what we get using the bounding diameter to compute a volume and particle count for 3kg:

>>> d0=0.012718275266250623
>>> r0=d0/2
>>> V0 = pi * r0**3 * (4/3)
>>> V0
1.0771676037098241e-06
>>> 3/(1500*V0)
1856.7212689203525

Again, this number is only correct for the case of spherical particles.

Thanks for the report.

– Charles

PS: rather than installing scipy it’s probably easier to do

>>> from math import gamma
>>> def beta(a,b): return gamma(a)*gamma(b)/gamma(a+b)

Thank you for your reply! @cgw
I get it! In MFIX, the equivalent diameter of the flat disk is 0.00905 m, but in fact it is 0.00941 m, which may be related to the value of m, n. Thank you again for your help.

I calculated an equivalent radius of 0.0045227, multiplying this by 2 yields 0.0090454, I’d like to understand where the 0.00941 number is coming from.

Hello @cgw
The equivalent diameter is calculated according to the volume of the cylinder is equal to the volume of the sphere. The radius of the cylinder is 0.00635 m and the length is 0.0035 m.

Ok, I see. For a cylinder, with r=0.00635, h=0.0035, the volume is 4.43369e-07, and the equivalent radius would be 0.009461.

But your particle is not a cylinder, the sides are not flat, so the volume is a little bit less and we get the equivalent radius 0.009045, which I believe is correct.

Yes, you are right! @cgw
In the SQP method, if the n value is set very large, the calculation will be wrong for a period of time, so the n value can only be reduced. Thank you for your reply!

Sharp corners lead to extremely high values for derivatives which sometimes leads to numerical overflow, loss of accuracy, etc. I’m not sure what can be done about this, other than to avoid creating shapes with sharp corners.