PFC3D Fluid-Particle Interaction Formulation

The following table gives symbols and a consistent set of units for typical fluid properties. Any set of units can be used as long as they are consistent with the units used in the PFC mechanical calculation.

Any set of units can be used as long as they are consistent with the units used in the mechanical calculation.

Table 1: System of SI units for fluid flow
Property Unit Symbol
Length m l
Fluid Density kg/m3 \(\rho_f\)
Time s t
Fluid Velocity m/s \(\vec{v}\)
Particle Velocity m/s \(\vec{u}\)
Porosity   \(\epsilon\)
Dynamic Viscosity Pa·s \(\mu\)
Drag Coefficient   \(C_d\)
Reynolds Number   \(Re\)
Fluid Pressure Pa p
Fluid Pressure Gradient Pa/m \(\vec{\nabla}p\)
Fluid Kinematic Pressure m2/s 2 P
Kinematic Viscosity m2/s \(\nu\)

In porous flow there are two definitions of velocity. In the CFD module the fluid velocity \(\vec{v}\) is the macroscopic or Darcy velocity. The macroscopic velocity is the volumetric flow rate per unit cross-sectional area. This is a nonphysical velocity since it assumes that the flow occurs across the entire cross-section, where fluid-flow actually only occurs in the pore space. The interstitial velocity is the actual velocity a parcel of fluid has as it moves through the pore space. To convert from macroscopic velocity to interstitial velocity divide by porosity.

A brief summary of the equations describing the fluid-particle interaction is given here. The equations of motion for the PFC3D particles are given by the standard PFC3D equations described in the “PFC Model Formulation” section, with an additional forcing term to account for interaction with the fluid:

(1)\[\frac{\partial \vec{u}}{\partial t} = \frac{\vec{f}_{mech} + \vec{f}_{fluid}}{m} + \vec{g}\]
(2)\[\frac{\partial \vec{\omega}}{\partial t} = \frac{\vec{M}}{I}\]

where \(\vec{u}\) is the particle velocity, \(m\) is the particle mass, \(\vec{f}_{fluid}\) is the total force applied by the fluid on the particle, \(\vec{f}_{mech}\) is the sum of additional forces (externally applied forces and contact forces) acting on the particle, \(\vec{g}\) is the acceleration due to gravity, \(\vec{\omega}\) is the particle angular velocity, \(I\) is the moment of inertia, and \(M\) is the moment acting on the particle. The force applied by the fluid to the particles (the fluid-particle interaction force), \(\vec{f}_{fluid}\), is made up of two parts: a drag force, and a force due to the fluid pressure gradient.

The drag force the fluid exerts on the particles is defined individually for each particle, based on the conditions in the fluid element that contains the particle. Depending on how the porosity is calculated, a particle may intersect more than one fluid element. In this case, the forces are distributed based on the fractional overlap between the particle and the fluid element. Note that the fluid-particle interaction force is always applied at the particle centroid, and no rotational moment is applied to the particle. For the purpose of calculating a fluid-particle interaction force, PFC clumps are treated as a single sphere with an equivalent radius calculated from the clump volume. The fluid-particle interaction force is applied to the clump centroid and no moment is applied. As with PFC balls, PFC clumps may have fluid-particle interaction forces applied from more than one fluid element. The drag force \(\vec{f}_{drag}\) is defined as

(3)\[\vec{f}_{drag} = \vec{f}_0 \epsilon ^{-\chi}.\]

\(\vec{f}_0\) is the single particle drag force and \(\epsilon\) is the porosity of the fluid element in which the particle resides. The \(\epsilon^{-\chi}\) term is an empirical factor to account for the local porosity. This correction term makes the force applicable to both high- and low-porosity systems and for a large range of Reynolds numbers ([DiFelice1994], [Xu1997]).

The single particle drag force is defined as

(4)\[\vec{f}_0 = \left(\frac{1}{2} C_d \rho_f \pi r^2 |\vec{u}-\vec{v}|(\vec{u}-\vec{v}) \right)\]

where \(C_d\) is a drag coefficient, \(\rho_f\) is the fluid density, \(r\) is the particle radius, \(\vec{v}\) is the fluid velocity, and \(\vec{u}\) is the particle velocity.

The drag coefficient is defined as

(5)\[C_d = \left(0.63 + \frac{4.8}{\sqrt{Re_p}}\right)^2,\]

where \(Re_p\) is the particle Reynolds number.

The empirical coefficient \(\chi\) is defined as

(6)\[\chi = 3.7-0.65 \exp \left( -\frac{(1.5-\log_{10}Re_p)^2}{2}\right)\]

The particle Reynolds number is

(7)\[Re_p = \frac{2\rho_f r |\vec{u}-\vec{v}|}{\mu_f}\]

where \(\mu_f\) is the dynamic viscosity of the fluid.

The body force per unit volume exerted on the fluid is

(8)\[\vec{f_b}=\frac{\sum_j \vec{f^j}_{drag}}{V}\]

where \(V\) is the volume of the fluid element and the sum is over the particles which overlap the fluid element.

(9)\[\vec{f}_{fluid} = \vec{f}_{drag} + \frac{4}{3}\pi r^3\left(\nabla p - \rho_f\vec{g} \right)\]

Note that the buoyancy is added explicitly because the fluid pressure field may not contain a hydrostatic component. The addition of the explicit buoyancy term can be controlled with the command cfd buoyancy.

Implementation Details

Porosity Calculation

There are two options available for calculating the porosity in a fluid element: using the location of the particle centroid, or a polyhedron representation of the particles. The command cfd porosity changes this mode.

The centroid method considers a particle to be entirely within a fluid element if its centroid is contained within that element. This method is fast, and it is conservative of particle volume. However, this method can result in jumps in porosity as the particles move, which can reduce the smoothness of the solution. As the particle size decreases relative to the size of the fluid elements, this effect decreases.

The polyhedron method represents a particle as a cube with height, length, and width equal to the diameter of the particle. The volume of intersection of this cube with the fluid element is calculated and adjusted to conserve the particle volume. This method is slower, but it has the advantage that the change in porosity is smooth as a particle moves from one fluid element to another.

A porosity of zero occurring in a fluid element results in a singularity in the governing equations of the fluid. To avoid problems due to the singularity, PFC3D limits the porosity to < 0.005. Note that the full volume of the spherical particle is used in the porosity calculation. No volume adjustment is made for particles that overlap.

The porosity contribution of PFC clumps to fluid elements is handled differently from that of the PFC balls. Each pebble of a clump contributes to the porosity of zero or more fluid elements, and a correction to these contributions is made to conserve total clump volume. The volume of overlap between a clump pebble and an element is determined in the same way as for the PFC balls (using either the centroid or polyhedron method). This volume of overlap is then scaled such that each pebble has an original volume equal to the clump volume divided by the number of pebbles in that clump.

Fluid Mesh and Data Format

The CFD module supports meshes composed of all tetrahedral or all hexahedral elements. Hexahedral element faces must be planar to within a small tolerance. PFC expects the following hexahedral node ordering:

       z
3----------0
|\     ^   |\
| \    |   | \
|  \   |   |  \
|   2------+---1
|   |  +-- |-- | -> x
7---+---\--4   |
 \  |    \  \  |
  \ |     \  \ |
   \|      y  \|
    6----------5

Two methods are available for importing a fluid mesh: the itasca.cfdarray.create_mesh Python function or the cfd read nodes and cfd read elements commands.

The following is an example of using Python to import a fluid mesh

import numpy as np
from itasca import cfdarray as ca

nodes = np.array(((0.0,  1.0,  0.0), (1.0,  1.0,  0.0),
                  (0.0,  1.0,  1.0), (0.0,  0.0,  0.0),
                  (0.0,  0.0,  1.0), (1.0,  0.0,  0.0),
                  (1.0,  1.0,  1.0), (1.0,  0.0,  1.0),
                  (0.0,  1.0,  2.0), (1.0,  1.0,  2.0),
                  (0.0,  0.0,  2.0), (1.0,  0.0,  2.0)))

elements = np.array(((7, 6, 1, 5, 4, 2, 0, 3),
                     (11, 9, 6, 7, 10, 8, 2, 4)), dtype=np.int64)

ca.create_mesh(nodes, elements)

PFC can read fluid meshes from text files. The mesh information is divided into two files: a node file and an element file.

Example Node File

        12        12
         1  0.00000E+00  1.00000E+00  0.00000E+00
         2  1.00000E+00  1.00000E+00  0.00000E+00
         3  0.00000E+00  1.00000E+00  1.00000E+00
         4  0.00000E+00  0.00000E+00  0.00000E+00
         5  0.00000E+00  0.00000E+00  1.00000E+00
         6  1.00000E+00  0.00000E+00  0.00000E+00
         7  1.00000E+00  1.00000E+00  1.00000E+00
         8  1.00000E+00  0.00000E+00  1.00000E+00
         9  0.00000E+00  1.00000E+00  2.00000E+00
        10  1.00000E+00  1.00000E+00  2.00000E+00
        11  0.00000E+00  0.00000E+00  2.00000E+00
        12  1.00000E+00  0.00000E+00  2.00000E+00

The two integers on the first line are the number of nodes. The subsequent lines have an integer node index (starting with 1) followed by the \(x\)-, \(y\)-, and \(z\)-components of the position of the node.

Example Element File

         2  2
         8
         1  8  7  2  6  5  3  1  4  1.004E-3  1.428E-4  9.98230E+02  1.18E-4 
    0
         8
         2  12  10  7  8  11  9  3  5  1.004E-3 1.428E-4  9.98230E+02  1.18E-4
    0

The two integers on the first line are the number of elements. The elements are then specified with three lines each.

The first line of each element description is a single integer giving the element type code (8 for hexahedral meshes; 6 for tetrahedral meshes). The second line starts with an integer element index (starting with 1) followed by 8 (or 4 for tetrahedral meshes) integer values that correspond to node indices. The final four floating point values on this line are required but are not used by PFC3D. The third line contains a single integer. This line is required but is not used by PFC3D.

Reading Fluid Flow Data into PFC

The CFD module associates fluid properties with each fluid element. There are several ways to set these properties. Fish or Python scripting can be used to set these properties programmatically. The command element cfd attribute can set these properties. Finally, the commands:

can be used to set these properties from a text file. The format for these files is given by the following rules,

  • The first line should be the number of elements.
  • There should be one subsequent line per element.
  • These lines should have a single number for scalar quantities (pressure, viscosity, density) or three numbers of vector quantities (velocity or pressure gradient).
  • If more than one number occurs on a line, the numbers should be separated by a space.

Practical Considerations

This section outlines practical limits on the input parameters. The fluid grid should be sufficiently fine to resolve flow structures. The following inequality should be met:

(10)\[\frac{d_c}{\Delta x_{cfd}} > 5\]

where \(d_c\) is the minimum width of the flow domain and \(\Delta x_{cfd}\) is the fluid element length.

This coupling methodology is designed to describe the average coupling forces occurring within one fluid element. Flow around the particles is not explicitly represented, since the local porosity is assumed to be evenly distributed within one element. For good results, several PFC3D particles should fit inside one CFD element:

(11)\[\frac{\Delta x_{cfd}}{2r} > 3\]

Where the exact particle micro-properties are not critical to the problem, there is some scope to reduce the particle stiffness to increase the allowable PFC3D timestep. However, particles colliding at maximum velocity with the PFC3D walls should not overlap more than half the radius of the particle. The following inequality describes this limit:

(12)\[\frac{2}{3} \pi r^3 \rho_p \left|\vec{u}_{max}\right|^2 < \frac{k_n r^2}{8}\]

where \(\rho_p\) is the particle density and \(k_n\) is the normal contact stiffness of the particle.

Timing

The coupling interval should be small enough to resolve the desired coupling behavior. Coupling information should be exchanged several times as a particle moves across a fluid element. This condition is satisfied when the following inequality is met:

(13)\[\frac{\Delta x_{cfd}}{|\vec{u}|t_c} > 3\]

where \(t_c\) is the coupling interval.

In practice, the PFC3D timestep is often smaller than the fluid timestep. Therefore, several cycles of PFC3D are often needed to meet one fluid step.

Meshing

If the particles are interacting with the fluid boundaries, then PFC3D walls corresponding to the fluid boundaries need to be created within PFC3D. If the geometry is simple, PFC walls can be used to produce the domain boundaries. If the boundary geometry is complex, it is possible to mesh the surface of the fluid domain with triangles and import them into PFC3D as an STL file. A large number of wall segments can result in slow execution in PFC3D. Often, a combination of the two techniques discussed above can be used.

Coupling

When the CFD module is active the fluid-particle interaction forces (\(\vec{f}^{j}_{fluid}\)) are applied to the PFC particles during the PFC cycle sequence. For more information on the PFC cycle sequence see Cycling. The fluid-particle interaction force, \(\vec{f}^{j}_{fluid}\), and the porosity in each fluid element, \(\epsilon^i\), are recalculated at given intervals during PFC3D cycling according to the cfd interval command. The superscript i refers to the fluid element and the superscript j refers to the particle. The command cfd update can be used to force an update these quantities. Use the command cfd coupling to enable or disable the coupling terms.

Two-way hydro-mechanical coupling is realized as a series of data exchanges between a fluid-flow solver and PFC3D. The porosity, \(\epsilon^i\), in each fluid element is determined by PFC3D. The body force per unit volume, \(\vec{f}^{i}_{b}\) , in each fluid element is determined by PFC3D as

(14)\[\vec{f}^{i}_{b} = \frac{\sum_j \vec{f}^{i}_{drag}}{V_i}\]

where the sum is over all the particles in a given fluid element, and \(V_i\) is the volume of that element. The fluid velocity, \(\vec{v}^i\) , fluid pressure, \(p^i\), fluid pressure gradient, \(\vec{\nabla} p^i\) , fluid density, \(\rho^i_f\) , and fluid dynamic viscosity, \(\mu^i\), in each element are determined by the fluid-flow software.

Exchanging information with the fluid-flow solver is not done automatically and must be done explicitly in the PFC data file. The synchronization and data exchanges can be accomplished via TCP socket communication with FISH or Python scripting.

Calculation Cycle

Solving a fluid-particle interaction problem with PFC requires the following steps:

  1. Initialize PFC: create particles and walls.

  2. Initialize fluid flow solver.

  3. Flow solver sends fluid mesh into PFC.

  4. PFC initializes the CFD module and calculates porosity and drag force.

  5. PFC sends porosity and drag force to the flow solver.

  6. Flow solver updates and sends \(\vec{v}^i\), \(p^i\), \(\vec{\nabla} p^i\), and optionally fluid density, \(\rho^i_f\) , and fluid dynamic viscosity, \(\mu^i\) to PFC.

  7. PFC cycles forward a time increment equal to the coupling interval \(t_c\).

    7.1 PFC Equations of Motion

    7.2 If the porosity and fluid-particle interaction force are being updated (when cycle % interval == 0) the following three steps are performed, otherwise jump to step 8.

    7.2.1 PFC calls functions registered at the CFD_BEFORE_UPDATE call point.

    7.2.2 PFC recalculates calculates the porosity field, \(\epsilon^i\), the fluid-particle interaction force, \(\vec{f}^i_{fluid}\) and the fluid drag force, \(\vec{f}^i_b\).

    7.2.3 PFC calls functions registered at the CFD_AFTER_UPDATE call point. This call point is executed after the fluid particle interaction force has been calculated but before this force has been added to the ball or clump unbalanced force. This allows the fluid particle interaction force to be changed via FISH or Python scripting.

    7.3 PFC Force-Displacement Law

    7.4 If the PFC mechanical time is less than \(t_c\) return to step 7.1.

  8. Return to step 5.

Step 7 in this list is done automatically by PFC during cycling. The other steps need to be explicitly done in the PFC datafile.

References

[DiFelice1994]Di Felice, R. D. “The voidage function for fluid-particle interaction systems,” Int. J. Multiphase Flow, Vol. 20, No. 1, 153-159 (1994).
[Xu1997]Xu, B. H., and A. B. Yu. “Numerical simulation of the gas-solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics,” Chemical Engineering Science, Vol. 52, No. 16, 2785-2809 (1997).