在粒子束的颗粒云上定义用户编码的体积力

在这个简单的例子中，用户程序比场函数更长且更难，但在更复杂的模拟中，用户程序可能更易于编写和调试。

在这两种方法中：

• 粒子束的颗粒云上的体积力与从原点到粒子束位置的距离的立方成反比。
• 为消除奇点，最小径向距离将剪切为 1 cm。
• 颗粒以 -10 m/s 的速度从顶面喷射 10 ms。
• 在喷射后的时间步，体积力开始起作用，将颗粒吸引到域中心的原点。
• 250 ms 过后，粒子束将聚集在原点（实际上并非如此）。

场函数方法和用户程序方法将产生相同的结果。

对于场函数和用户程序，喷射质量流率在喷射定义中设为场函数，用于在 10 ms 后关闭喷射（即，使用场函数：($Time < 0.01) ? 0.001 : 0.0）。

场函数方法

定义了两个场函数来表示体积力：

• rCubed，定义如下：

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

• UserParticleBodyForce，定义如下：

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

用户程序方法

定义相同体积力的用户程序为：

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

该用户程序在粒子束的颗粒云上指定矢量力（在结果矢量中返回）。 粒子束形心的尺寸为 3*CoordReal，将作为矢量场函数传递到用户程序。 注册此函数和参数列表的 uflib.f 函数为：

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

在这种情况下，例程中需要“ParcelCentroid(粒子束形心)”和时间场函数。