Defining a User-Coded Body Force on a Parcel Cloud

In this simple case, the user code is longer and more difficult than the field functions, but in more elaborate situations, user code can be easier to write and debug.

In both methods:

• The body force on the parcel cloud is inversely proportional to the cube of the distance from the origin to the parcel position.
• The minimum radial distance is clipped to 1 cm to eliminate the singularity.
• Particles are injected for 10 ms from the top surface with a velocity of -10 m/s.
• On the time-step after injection, the body force begins to act, attracting the particles to the origin at the center of the domain.
• After 250 ms, the parcels are (unphysically) gathered at the origin.

Both the field function method and the user code method produce the same results.

For both field functions and user code, the injection mass flow rate is set to a field function that turns the injection off after 10 ms (that is, using a field function of: ($Time < 0.01) ? 0.001 : 0.0) in the injection definition.

Field Function Method

Two field functions are defined to describe the body force:

• rCubed, defined as:

  max(0.01,pow($$ParcelCentroid.mag(),3))

• UserParticleBodyForce, defined as:

  ($Time > 0.01) ? -9.81e5*($$ParcelCentroid)/$rCubed : [0,0,0]

User Code Method

The user code defining the same body force is:

subroutine parcelBodyForce(result,size,centroid,Id,Time)
 c parcel body force example
   use StarRealMod
   integer, intent(in) :: size
   real(StarReal), intent(out), dimension(3,size) :: result
   real(CoordReal), intent(in), dimension(3,*) :: centroid
   real(StarInt), intent(in), dimension(*) :: Id
   real(CoordReal), intent(in) :: Time
 c body force is proportional to 1/r^2 after 0.01 seconds
   real(StarReal), parameter:: aScale = -9.81e5
   real(StarReal), parameter:: center(3) = (/0.0, 0.0, 0.0/)
   real(StarInt), parameter :: NULL_INDEX = -1
   real(StarReal) radius(3), rmag
   integer i
 c set the magnitude and direction for valid parcels
   do i = 1, size
     if ((Id(i) .gt. NULL_INDEX) .and. (Time > 0.01)) then
       radius(:) = centroid(:,i)-center(:)
       rmag = max(0.01,sqrt(sum(radius**2))**3)
       result(:,i) = aScale*(radius)/rmag
     endif
   end do
   return
   end

The user code specifies a vector force on the parcel cloud, returned in the result vector. The parcel centroid is passed into the user code as a vector field function with a size of 3*CoordReal. The uflib.f function that registers this function and argument list is:

subroutine uflib()
    use StarRealMod
 c Register user functions
    external parcelBodyForce
    call uffunc(parcelBodyForce,
    &    "ParcelProfile",
    &    "Parcel Body Force")
    call ufarg(parcelBodyForce,
    &    "Parcel",
    &    "$$ParcelCentroid",
    &    3*CoordRealSize)
    call ufarg(parcelBodyForce,
    &    "Parcel",
    &    "ParcelId",
    &    StarInt)
    call ufarg(parcelBodyForce,
    &    "Parcel",
    &    "$Time",
    &    CoordRealSize)
    return
    end

In this case, both the ParcelCentroid and Time field functions are required in the routine.