Particle shape change

Hi, how can we change shape of SQP particle by editing particle output file written from SQP simulation after particles are settled?

For example, say, prolate spheroid to spherical particle? (Instead of setting up a whole new different spherical particle simulation).

Best,

If you look at the particle_output.dat file saved at the end of the run you should see something like:

================================================================================================================================================================================================================
|       x       |       y       |       z       |    Phase_ID   |    Diameter   |       u       |       v      |       w       |     SQP a     |     SQP b     |     SQP c     |     SQP m     |     SQP n     |
================================================================================================================================================================================================================

Editing the values in the columns SQP a, b, c, m, n should redefine the particle shape, however it is up to you to make sure there are no overlaps. If you only reduce a,b, and c (the semiaxes) you should be able to get away with it. It’s up to you to write a script to modify the file, I’d use Python.

However, when looking into this I found something odd. There is actually no data in the rightmost column (SQP n) so there may be a bug in the code which writes this file. I’ll follow up on this.

– Charles

Hi Team, greetings.

Did you get a chance to look into this?

Thank you,

@jagan1mohan

Unfortunately, this slipped between the cracks, and I apologize. I just tested with the 24.1.1 release and I see the same problem. Here’s a part of a particle_output.dat I generated.

===========================================================================================================================================================================================================================================================
|       x       |       y       |       z       |    Phase_ID   |    Diameter   |       u       |       v      |       w       |     SQP a     |     SQP b     |     SQP c     |     SQP m     |     SQP n     |     SQP q1    |     SQP q2    |     SQP q3
===========================================================================================================================================================================================================================================================
  0.72614829E-02  0.36370217E-02  0.91296748E-02        1         0.78964441E-02 -0.40615761E-01 -0.15655867E+00 -0.21305198E-01  0.30000000E-02  0.30000000E-02  0.30000000E-02  0.40000000E+01  0.10000000E+01  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.15141674E-01  0.38399431E-02  0.88044984E-02        1         0.78964441E-02 -0.59874150E-01 -0.15150379E+00 -0.35527086E-01  0.30000000E-02  0.30000000E-02  0.30000000E-02  0.40000000E+01  0.10000000E+01  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.24014156E-01  0.39261728E-02  0.86376993E-02        1         0.78964441E-02 -0.38318918E-01 -0.14791545E+00 -0.42361332E-01  0.30000000E-02  0.30000000E-02  0.30000000E-02  0.40000000E+01  0.10000000E+01  0.00000000E+00  0.00000000E+00  0.00000000E+00

The quaternion should have 4 components, not 3, and the number that appears in the SQP n column is actually the first component of the quaternion. So the problem persists.

It’s a bit unfortunate that this didn’t get fixed for the maintenance release we just put out. I’ll work on a fix for this and maybe you can just replace a few files instead of waiting for the next release.

Thanks again for reporting this and sorry for the delay,

– Charles

Two small issues in read_particle_input.f (which also writes particle_output):

diff --git a/model/des/read_particle_input.f b/model/des/read_particle_input.f
index 04e5a22de..3da2c04f5 100644
--- a/model/des/read_particle_input.f
+++ b/model/des/read_particle_input.f
@@ -3366,9 +3366,9 @@ CONTAINS
       character(255):: Phase_ID_setting, Diameter_setting, Density_setting, Velocity_setting
       character(255):: Temperature_setting, Species_setting, User_scalar_setting
       character(255):: SQP_semi_axis_setting, SQP_roundness_setting, SQP_quaternion_setting
-      character(255):: column_header
+      character(300):: column_header
       integer :: header_bar_length
-      character(255):: header_bar
+      character(300):: header_bar
       character(32) :: label
 
       DOUBLE PRECISION, ALLOCATABLE :: ltemp_array(:,:)
@@ -3582,7 +3582,7 @@ CONTAINS
          column_header         = trim(column_header)//"     SQP m     |     SQP n     |"
          SQP_roundness_setting = 'SQP_roundness:   T'
          SQP_roundness_start   = n_vars + 1
-         n_vars                = n_vars + 1
+         n_vars                = n_vars + 2
       ENDIF
 
 ! SQP quaternion option

The header line is not long enough for all the fields, and there’s an off-by-one error (roundness is two variables, m,n)

This will be fixed in the next MFiX release, whether that is 24.1.2 or 24.2. However in the meanwhile you can put this modifed read_particle_input.f in your project directory, build a custom solver, and you should be good to go.

read_particle_input.f (162.3 KB)

– Charles

Getting back to your original question: Once you’ve gotten past the read_particle_input.f issue, you can use something like this script to modify data in the file:

#!/usr/bin/env python
import numpy as np

# Load particle_output.dat
fname = 'particle_output.dat'
lines = open(fname).readlines()
header = lines[:38]
data = np.loadtxt(fname, skiprows=38)

# Define column names for convenience
x,y,z,phase,diameter,u,v,w,a,b,c,m,n,q1,q2,q3,q4 = range(17)

# Modify data here
# Example: modify semi-axis 'a' of particles in phase 1
for row in data:
    if row[phase] == 1:
        row[a] = 0.002

# Write out the modified data
output = 'particle_output.new'
with open(output,'w') as f:
    f.write(''.join(header))
    for row in data:
        for val in row:
            f.write('%16f'%val)
        f.write('\n')

This should work to reduce size of particles (say, prolate to oblate) but if you increase anything the particles will overlap and the simulation will not run.

Also note that an SQP case with a=b=c and m=n=2 (spherical case) will run slower than a straight DEM spherical simulation, so if you are going from prolate to sphere, you might be better off switching to the DEM solver for that.

Hope this helps,

– Charles

1 Like

Hi Team, greetings. Thank you for your replies. At my end, this is work in progress and few important points from last week are:

→ Also note that an SQP case with a=b=c and m=n=2 (spherical case) will run slower than a straight DEM spherical simulation, so if you are going from prolate to sphere, you might be better off switching to the DEM solver for that.
I’m trying to compare non-spherical and spherical systems and want to keep everything same - drag law, particle count, flow conditions. Hence, using same solver in MFiX

→ Problem statement is slightly modified now–original question was changing monodisperse sphere to monodisperse SQP. I have to generate prolates from image data as attached below.


All 3 axes have mean, standard deviation from image post-processing. I’d make an attempt and let you know.

Thank you,
Jagan Mohan

Thanks for the update. Please let us know how things work out for you. Also note that the recently-released glued-sphere particle model (GSP) might be of interest and should run faster than SQP.