User Defined Particle-Wall Interaction via wall film: how many max. fragments?

- October 3, 2024 at 8:51 pm
  
  nilotpalc
  Subscriber
  
  #question Hi everyone, I have a very simple question. I am writing a user defined Fortran subroutine for particle-wall interaction via a wall film, and was referring to the following code documented in the page (8.6.3.1.3. User Defined Particle-Wall Interaction) of CFX Solver Modelling Guide. The question is about a single line in this code: REAL FACT, ANGLE, VEL_PT(3), DIAM_PT, NRATE_PT,
  & BRKUP_IND(3), CHILD_DIAM(3), CHILD_NRATE(3),
  & CHILD_COEF_N(3), CHILD_COEF_P(3), CHILD_MODE(3), DENS_FL
  Is increasing the size of BRKUP_IND array allowed in CFX? I have tried doing BRKUP_IND(10), which means a particle can generate 10 fragments. The simulation ran without any error and a single particle did end up generating 10 fragments. However, the CFX engineers put this line on the page containing this code: "Note that a maximum of 4 child droplets may be generated with the child generation model." What do they mean here? Maximum of 4 fragments by this code or in general? Although I tested with 1 particle, I am concerned it might not work for say a million particles and I might get an "overflow" error. Could anyone of you throw some light on this? I didn't post on CFD Online because some smart people there deviate from answering the question and ask to run a simulation and check. In my case, I already did that, so I think AIS is the right place to ask. Thank you!
  
  #include "cfx5ext.h"
  dllexport(pt_wall_splash)
  SUBROUTINE PT_WALL_SPLASH (NLOC,NRET,NARG,RET,ARG,CRESLT,
  & CZ,DZ,IZ,LZ,RZ)
  CC
  CD User routine: template for particle user routine
  CC
  CC --------------------
  CC Input
  CC --------------------
  CC
  CC NLOC - number of entities
  CC NRET - length of return stack
  CC NARG - length of argument stack
  CC ARG - argument values
  CC
  CC --------------------
  CC Modified
  CC --------------------
  CC
  CC --------------------
  CC Output
  CC --------------------
  CC
  CC RET - return values
  CC
  CC --------------------
  CC Details
  CC --------------------
  CC
  CC======================================================================
  C
  C ------------------------------
  C Preprocessor includes
  C ------------------------------
  C
  #include "cfd_sysdep.h"
  #include "cfd_constants.h"
  C
  C ------------------------------
  C Argument list
  C ------------------------------
  C
  INTEGER NARG, NRET, NLOC
  C
  CHARACTER*(4) CRESLT
  C
  REAL ARG(NLOC,NARG), RET(NLOC,NRET)
  C
  INTEGER IZ(*)
  CHARACTER CZ(*)*(1)
  DOUBLE PRECISION DZ(*)
  LOGICAL LZ(*)
  REAL RZ(*)
  C
  C=======================================================================
  C
  C ---------------------------
  C Executable Statements
  C ---------------------------
  C
  C=======================================================================
  C
  C Return variables:
  C -----------------
  C
  C Child droplet breakup indicator : RET(1,1)
  C Child droplet diameter : RET(1,5)
  C Child droplet number rate : RET(1,9)
  C Child droplet normal coef of restitution : RET(1,13)
  C Child droplet parallel coef of restitution : RET(1,17)
  C Child droplet mode : RET(1,21)
  C
  C Argument variables
  C -------------------
  C
  C Particle impact angle : ARG(1,1)
  C Particle velocity : ARG(1,2)
  C Particle diameter : ARG(1,5)
  C Particle number rate : ARG(1,6)
  C Fluid density : ARG(1,7)
  C
  C We know that NLOC is 1 for the particle user source routines!!!!
  C=======================================================================
  C
  C-----------------------------------------------------------------------
  C Calculate the return variables
  C-----------------------------------------------------------------------
  C
  CALL SPLASH (RET(1,1),RET(1,5),RET(1,9),RET(1,13),RET(1,17),
  & RET(1,21),ARG(1,1),ARG(1,2),ARG(1,5),ARG(1,6),
  & ARG(1,7))
  C
  END
  C
  SUBROUTINE SPLASH (BRKUP_IND,CHILD_DIAM,CHILD_NRATE,
  & CHILD_COEF_N,CHILD_COEF_P,CHILD_MODE,
  & ANGLE,VEL_PT,DIAM_PT,NRATE_PT,DENS_FL)
  C
  C ---------------------------
  C Preprocessor includes
  C ---------------------------
  C
  #include "cfd_sysdep.h"
  #include "cfd_constants.h"
  C
  C---- Do not change these defines, they have to be consistent with what
  C the solver uses
  C
  #define __regular_particle__ 1
  #define __wall_particle__ 20
  #define __pt_uf_get_data__ 1
  #define __pt_uf_save_data__ 2
  #define __pt_uf_iseed__ 1
  #define __pt_uf_pmode__ 2
  C
  C ------------------------------
  C Argument list
  C ------------------------------
  C
  REAL FACT, ANGLE, VEL_PT(3), DIAM_PT, NRATE_PT,
  & BRKUP_IND(3), CHILD_DIAM(3), CHILD_NRATE(3),
  & CHILD_COEF_N(3), CHILD_COEF_P(3), CHILD_MODE(3), DENS_FL
  C
  C ------------------------------
  C Local variables
  C ------------------------------
  C
  REAL SIGMA, VISC_PT, R, SPEED_PT, ANGLE_DEG, ALFA, BETA,
  & GAMMA, VEL_NORM, OH, RE, K, M1M0
  INTEGER ICHILD, ISEED, ACTION
  C
  C ------------------------------
  C Executable statements
  C ------------------------------
  C
  C=======================================================================
  C Prologue
  C=======================================================================
  C
  C---- Initialization
  C
  ICHILD = 0
  C
  C---- Get particle seed info (needed for random number generator)
  C
  ACTION = __pt_uf_save_data__
  CALL USER_PARTICLE_INFO(ACTION,__pt_uf_iseed__,ISEED)
  C
  C---- Some material properties (water)
  C --> Surface tension coefficient
  C --> Particle viscosity
  C
  SIGMA = 0.0201
  VISC_PT = 0.001
  C
  C---- Random number (required for deviation angle)
  C
  CALL GET_RANDOM(R,1,ISEED)
  C
  C=======================================================================
  C Computation Section
  C=======================================================================
  C
  C---- Particle velocity magnitude and impact angle
  C
  SPEED_PT = SQRT(VEL_PT(1)**2 + VEL_PT(2)**2 + VEL_PT(3)**2)
  ANGLE_DEG = ANGLE*180./PI
  C
  C---- Ohnesorge and Reynolds number
  C --> Re, based on wall normal impact velocity!
  C
  VEL_NORM = SIN(ANGLE)*SPEED_PT
  OH = VISC_PT/SQRT(DENS_FL*DIAM_PT*SIGMA)
  RE = DENS_FL*DIAM_PT*VEL_NORM/VISC_PT
  C
  C-----------------------------------------------------------------------
  C Precompute reflection angle
  C-----------------------------------------------------------------------
  C
  C---- Impact angle, ALFA, average reflection angle, BETA. Add random
  C deviation angle, GAMMA, if ALFA < 15 deg.
  C
  ALFA = 90. - ANGLE_DEG
  BETA = 61.88 + 0.326*ALFA
  GAMMA = 17.6 - 0.18*ALFA
  C
  IF (ALFA .LT. 15.) BETA = BETA + (R-0.5)*GAMMA
  C
  C---- Compute K-value
  C --> Clip K to avoid overflow during m1/m0 calculation
  C (221^9.2133 = 3.9773E+21)
  C
  K = MIN(OH*RE**1.25,221.)
  C
  C-----------------------------------------------------------------------
  C K <= 57 -> Particle deposited at wall
  C-----------------------------------------------------------------------
  C
  IF (K .LE. 57) THEN
  C
  ICHILD = ICHILD + 1
  BRKUP_IND(ICHILD) = 1.
  CHILD_DIAM(ICHILD) = DIAM_PT
  CHILD_NRATE(ICHILD) = NRATE_PT
  CHILD_COEF_N(ICHILD) = 0.0
  CHILD_COEF_P(ICHILD) = 0.0
  CHILD_MODE(ICHILD) = __wall_particle__
  C
  C-----------------------------------------------------------------------
  C ... else partially reflected
  C-----------------------------------------------------------------------
  C
  ELSE
  C
  C---- Mass fraction across reflection
  C
  M1M0 = 3.9869E-21*K**9.2133
  C
  C---- Reflected portion of particle
  C
  ICHILD = ICHILD + 1
  BRKUP_IND(ICHILD) = 1.
  CHILD_DIAM(ICHILD) = DIAM_PT*M1M0
  CHILD_NRATE(ICHILD) = NRATE_PT
  CHILD_COEF_N(ICHILD) = COS(BETA*PI/180.)*SPEED_PT
  CHILD_COEF_P(ICHILD) = SIN(BETA*PI/180.)*SPEED_PT
  CHILD_MODE(ICHILD) = __regular_particle__
  C
  C---- Deposited portion of particle
  C
  ICHILD = ICHILD + 1
  BRKUP_IND(ICHILD) = 1.
  CHILD_DIAM(ICHILD) = DIAM_PT*(1. - M1M0)
  CHILD_NRATE(ICHILD) = NRATE_PT
  CHILD_COEF_N(ICHILD) = 0.0
  CHILD_COEF_P(ICHILD) = 0.0
  CHILD_MODE(ICHILD) = __wall_particle__
  C
  ENDIF
  C
  C=======================================================================
  C Epilogue
  C=======================================================================
  C
  C---- Update ISEED in particle database with value computed in
  C this routine
  C
  ACTION = __pt_uf_get_data__
  CALL USER_PARTICLE_INFO(ACTION,__pt_uf_iseed__,ISEED)
  C
  END
- October 4, 2024 at 11:15 am
  
  Mark Owens
  Ansys Employee
  
  Hi, you might get an overflow or memory error. The reason the help says "Note that a maximum of 4 child droplets may be generated with the child generation model." is becuase PT_WALL_SPLASH chops the single return array into chuncks of 4 as RET(1,1), RET(1,5), etc. You could chop it into larger chunks but it is not clear to me whether sufficient memory will get allocated. I am waiting for the developers to get back to me.
- October 4, 2024 at 1:40 pm
  
  Mark Owens
  Ansys Employee
  
  You can set the maximum number of children created at breakup using the expert parameter:
  
  pt real droplet breakup max children
  
  It doesn't appear in the user interface as an available expert parameter, but you can insert an expert parameter object into the tree in CFX-Pre and with a right mouse click on it select the option to edit in the command editor. Then you can insert
  
  pt real droplet breakup max children = 10
  
  You should then have in PT_WALL_SPLASH
  
  C Return variables:
  C -----------------
  C
  C Child droplet breakup indicator :RET(1,1)
  C Child droplet diameter :RET(1,11)
  C Child droplet number rate :RET(1,21)
  C Child droplet normal coef of restitution :RET(1,31)
  C Child droplet parallel coef of restitution :RET(1,41)
  C Child droplet mode :RET(1,51)
  
  and
  
  CALL SPLASH (RET(1,1),RET(1,11),RET(1,21),RET(1,31),RET(1,41),
  & RET(1,51),ARG(1,1),ARG(1,2),ARG(1,5),ARG(1,6),
  & ARG(1,7))
  
  and in SPLASH
  
  REAL BRKUP_IND(10),CHILD_DIAM(10),CHILD_NRATE(10),CHILD_COEF_N(10),
  & CHILD_COEF_P(10),CHILD_MODE(10)
  REAL ANGLE, VEL_PT(3), DIAM_PT, NRATE_PT, DENS_FL
- October 12, 2024 at 12:03 pm
  
  nilotpalc
  Subscriber
  
  Hi Mark,
  This is very helpful, thanks so much for digging this out!
- January 27, 2025 at 5:27 pm
  
  nilotpalc
  Subscriber
  
  Hi Mark and everyone,
  I have a question related to this code. This code is called for each particle track when it interacts with a wall. Now if I want to write the values of intermediate variables (say the random number R from the line CALL GET_RANDOM(R,1,ISEED))to an external file or just print it on the console everytime the user routine is called, how do I do that? I can write a similar code that generates a random number, say in Python, but am afraid I will end up with a different random number. So, it’d be great if I had access to the random number or any other intermediate value calculated in the particle user routine itself.
  Thank you!