CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   CFX (http://www.cfd-online.com/Forums/cfx/)
-   -   Check particle impaction with User Fortran (http://www.cfd-online.com/Forums/cfx/94546-check-particle-impaction-user-fortran.html)

Julian K. November 18, 2011 13:01

Check particle impaction with User Fortran
 
Hi!

I am simulating the impaction of particles on a moving circular cylinder using Lagrangian particle tracking. I found out, that CFX does not consider the particle's radius while checking if a particle has collided with the cylinder's surface. Thus, the particle is captured, when its centre coordiantes of the coincide with the cylinder surface.
However, my case requires, that the particle is captured when the particle surface and the cylinder surface touch. Hence, I need to implement this simple algorithm:

Code:

if (x <= r_part + r_cyl)
then impaction
else no impaction

Where x is the distance between the centre coordiantes of the particle and the cylinder, r_part is the particle's radius and r_cyl is the radius of the cylinder.

I assume, that I need to implement the code via a User FORTRAN. Since I have absolutely no experience with this, I would like you to ask for some help. I found a file called 'pt_termination.F' (CFX\Samples\UserFortran) with the following code:

Code:

IF (DIST.GT.0.75) THEN
  FATE = DEAD
ENDIF

I wonder if this is the right point to start for me?! Do you have any other suggestions for me?

Thanks!

Julian K. November 23, 2011 19:25

In the file <CFXROOT>/etc/VARIABLES I found the variable 'ptpos' with this description:
Code:

  VARIABLE: ptpos
    Option = Definition
    MMS Name  = CRD
    Long Name = Particle Position
    Quantity = Length
    Tensor Type = VECTOR
    Status = M
    User Level = 2
    Output to Postprocessor = No
    Output to Jobfile = No
    General Availability = PTM_UR
    Physical Availability = PARTICLE
    Component Short Names = \
      ptpos x, \
      ptpos y, \
      ptpos z
    Component Long Names = \
      Particle Position X, \
      Particle Position Y, \
      Particle Position Z
    Old Component Long Names = \
      Particle Position x, \
      Particle Position y, \
      Particle Position z
    Component MMS Names = \
      CRD-1, \
      CRD-2, \
      CRD-3
    Variable Scope = PARTICLE
  END

I guess, that I can use this variable to calculate x, the distance between the cylinder's and particle's centre coordiantes?!

Now, what's the name of the variable which gives me the centre coordinates of the cylinder? Btw, the cylinder is a moving rigid body, thus it's position is not constant.

Julian K. January 12, 2012 10:27

By now, I found a way to determine if a particle touches the surface of an oscillating cylinder.

I use this User Fortran Code and put it into the file 'pt_termination.F':
Code:

#include "cfx5ext.h"
dllexport(pt_termination)
      SUBROUTINE PT_TERMINATION (NLOC,NRET,NARG,RET,ARG,CRESLT,
    &                            CZ,DZ,IZ,LZ,RZ)
CC
CD User routine: template for particle user termination 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    Particle fate                                        : RET(1,1)
C
C    Argument variables 
C    -------------------
C
C    Particle mean diameter                              : ARG(1,1)
C    Particle position                                    : ARG(1,2:4)
C    Cylinder x-position                  : ARG(1,5)
C    Cylinder y-position                  : ARG(1,6)
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 USER_FATE (RET(1,1),ARG(1,1),ARG(1,2),ARG(1,5),ARG(1,6))
C
      END
 
      SUBROUTINE USER_FATE (FATE,DIA,PPOS,CPOSX,CPOSY)
C
C ---------------------------
C    Preprocessor includes
C ---------------------------
C
#include "cfd_sysdep.h"
#include "cfd_constants.h"
C
C ------------------------------
C        Argument list
C ------------------------------
C
      REAL              FATE, DIA, PPOS(3), CPOSX, CPOSY
      REAL              CYLRAD, PRAD
      REAL              CONTACT
C
C ------------------------------
C        Local variables
C ------------------------------
C
      REAL              DEAD, ALIVE
C ------------------------------
C    Executable statements
C ------------------------------
C
      DEAD  = 0.
      ALIVE = 1. 
C
      CYLRAD = 0.005 
C
C---- Particle starts as alive
C
      FATE = ALIVE
C
C---- Calc particle radius
      PRAD = ((CPOSX-PPOS(1))**2 + (CPOSY-PPOS(2))**2)**0.5
      CONTACT = CYLRAD - (PRAD - DIA/2.0)
C
C---- Check if user defined criterion was exceeded, stop particle
C    --> Set particle mode to __dead__
C
      IF (CONTACT.GT.0.0) THEN
        FATE = DEAD
      ENDIF
C
      END

In CFX I create a 'Particle User Routine', which uses the output of the F-file compilation. In case of an oscillating cylinder, you need to define the cylinder coordinates as two 'Additional Variables', called 'Xcyl' and 'Ycyl' (Varible Type = Unspecified, Units = [m], Tensor Type = Scalar). Now, these variables can be selected from 'Default Domain > Fluid Models > Additional Variable Models > Additional Variable'. You define both variables with an 'Algebraic Equation' with is for 'Xcyl'
Code:

rbstate(Position X)@Rigid Body Cylinder
and 'Ycyl'
Code:

rbstate(Position Y)@Rigid Body Cylinder
Now comes the crucial part: The in- and output arguments of the User Fortran need to be defined with a 'Particle User Termination Criteria', however, this can not be selected in the 'Particle Control' Panel of the CFX GUI. The 'Particle User Termination Criteria' needs to be defined using two seperate CCL-files, 'pt_term.ccl' and 'pt_extra.ccl'. There is the contents:
pt_term.ccl: (make sure the 'Library Path' is set corretly)
Code:

LIBRARY:
  USER ROUTINE DEFINITIONS:
    USER ROUTINE: termination
      Calling Name = pt_termination
      Library Name = pt_termination
      Library Path = /work/simulations/test/part_capt/oscillation_capture
      Option = Particle User Routine
    END
  END
END
FLOW: 
  SOLVER CONTROL: 
    Turbulence Numerics = First Order
    ADVECTION SCHEME: 
      Option = High Resolution
    END
    CONVERGENCE CRITERIA: 
      Conservation Target = 0.01
      Residual Target = 1.E-4
      Residual Type = MAX
    END
    DYNAMIC MODEL CONTROL: 
      Global Dynamic Model Control = On
    END
    PARTICLE CONTROL: 
      PARTICLE INTEGRATION: 
        Option = Forward Euler
      END
      PARTICLE TERMINATION CONTROL:
        Maximum Tracking Time = 10 [s]
        Maximum Tracking Distance = 10 [m]
        Maximum Number of Integration Steps = 100000
        FLUID: Particle
          PARTICLE USER TERMINATION CRITERIA:
            Particle User Routine = termination
            Argument Variables List = \
              Particle.Mean Particle Diameter,\
              Particle.Particle Position,\
          Xcyl,\
          Ycyl
            Particle User Fate Return Variables List = \
              Particle Fate
          END
        END
      END
    END
  END
END

pt_extra.ccl:
Code:

RULES:
  SINGLETON: PARTICLE TERMINATION CONTROL
    Description = Sets details of the particle tracking integration \
                  algorithm.
    Solver Name = PT_TERM_CTRL
    Optional Parameter List = \
      Maximum Tracking Time, \
      Maximum Tracking Distance, \
      Maximum Number of Integration Steps, \
      Minimum Diameter, \
      Minimum Total Mass
    Editable Option = True
    Optional Child List = \
      MASS FRACTION LIMITS, \
    FLUID
  END
 
 OBJECT: FLUID
  Description = This object encloses fluid-specific models and data.
  #
  # ...Context....  ...Allowed Contents.........
  # DOMAIN          FLUID MODELS, SOLVER CONTROL
  # SUBDOMAIN        SOURCES, BULK SOURCE DISTRIBUTION
  # BOUNDARY        BOUNDARY CONDITIONS
  # INITIALISATION  INITIAL CONDITIONS
  #
  Solver Name = FL
  Object Name Restriction = CEL
  Optional Child List =  \
    ADDITIONAL VARIABLE, \
    COMBUSTION MODEL, \
    COMPONENT, \
    DILATATIONAL STRESS, \
    ELECTROMAGNETIC MODEL, \
    FLUID BUOYANCY MODEL, \
    HEAT TRANSFER MODEL, \
    KINETIC THEORY MODEL, \
    SOLID PRESSURE MODEL, \
    SOLID BULK VISCOSITY, \
    SOLID SHEAR VISCOSITY, \
    THERMAL RADIATION MODEL, \
    TURBULENCE MODEL, \
    TURBULENT WALL FUNCTIONS, \
    WALL SLIP MODEL, \
    EROSION MODEL, \
    NUCLEATION MODEL, \
    INITIAL CONDITIONS, \
    BOUNDARY CONDITIONS, \
    SOURCES, \
    BULK SOURCE DISTRIBUTION, \
    SOLVER CONTROL, \
    INJECTION CONDITIONS, \
    PARTICLE ABSORPTION, \
    PARTICLE ROUGH WALL MODEL, \
    PARTICLE USER TERMINATION CRITERIA
 END
 
 SINGLETON: PARTICLE USER TERMINATION CRITERIA
    Description = Controls user defined particle fate
    Solver Name = USR_TERM_CTRL
#  Context Rule = Option
#  Allowed Option List = \
#    User Defined
#  Essential Parameter List = \
#    Option
#  CONTEXT: User Defined
      Essential Parameter List = \
        Particle User Routine, \
        Argument Variables List, \
        Particle User Fate Return Variables List
#  END
  END
 
  PARAMETER: Particle User Fate Return Variables List
    Parameter Type = String List
    Solver Name = CVAR_RET
    Allowed String List = \
        Particle Fate
  END
END

These need to be added to the CFX solver call. From the command line, you add:
Code:

-ccl pt_term.ccl -ccl pt_extra.ccl
to your regular function call.

Now, after each time step. The solver writes the number of particles, which touched the surface of the cylinder in the 'Particle Diagnostics'.

Julian K. January 12, 2012 10:46

Now, I know, that a particle touch the cylinder. However, I also would like to know on which side, lee- or windward, the particle was captured. I came up with two ideas:
  1. open an external file from within the User Fortran call and write leeward, or windward in it
  2. return the particle and cylinder coordinates in the 'Particle User Fate Return Variable List'
to 1.: I already tried this, using a very simple test code, which I added to the subroutine call of the User Fortran:
Code:

C---- Open and close file again to see if it works
    n = 0.005
      OPEN(UNIT=1, FILE='OUTPUT', STATUS='NEW')
    WRITE (1,'(F4.3)') N
      CLOSE(1)

The result is that the file 'OUTPUT' is being created an the value n=0.005 is written into it. However, eventually, an error occurs and the solver exits with return code 10.

Does anyone know, how why this error occurs?

to 2.:
I haven't tried this, yet. But I suppose it won't work. I'll post my resutls, once I've tried it.


All times are GMT -4. The time now is 11:16.