Coupling PFC3D with 3rd Party Fluid Flow Software

PFC3D can be coupled with an any fluid flow solver using the same coarse grid coupling method described above. In the two previous examples the fluid-particle coupling problem is solved entirely within PFC3D. It may be necessary for PFC to synchronize with and share data with a separate process possibly running on a different computer.

This section gives a simple example of connecting PFC3D to a third-party flow solver using TCP sockets. In this case, the flow solver will simply return a velocity and pressure of zero for every time step. The file pfc_coupling.py demonstrates the PFC3D side of the coupling and cfd_coupling.py demonstrates the flow solver side of the coupling.

The Python classes itasca.util.p2pLinkServer and itasca.util.p2pLinkClient are used to establish this coupling. These classes are available for use in Python programs outside of the PFC Python environment. This package can be installed into any Python environment with the command pip install itasca (see https://pypi.python.org/pypi/itasca/2016.06.09).

The file pfc_coupling.py is listed below. This file should be run first from within PFC3D.

import itasca as it
from itasca import cfdarray as ca
from itasca.util import p2pLinkServer

import numpy as np

with p2pLinkServer() as cfd_link:
    cfd_link.start()

    nodes = cfd_link.read_data()
    elements = cfd_link.read_data()
    fluid_density = cfd_link.read_data()
    fluid_viscosity = cfd_link.read_data()
    print fluid_density, fluid_viscosity
    nmin, nmax = np.amin(nodes,axis=0), np.amax(nodes,axis=0)
    diag = np.linalg.norm(nmin-nmax)
    dmin, dmax = nmin -0.1*diag, nmax+0.1*diag
    print dmin, dmax

    it.command("""
    new
    cmat default model linear
    domain extent {} {} {} {} {} {}
    """.format(dmin[0], dmax[0],
               dmin[1], dmax[1],
               dmin[2], dmax[2]))
    ca.create_mesh(nodes, elements)
    it.command("""
    configure cfd
    set timestep max 1e-5
    element cfd attribute density {}
    element cfd attribute viscosity {}
    cfd porosity polyhedron
    cfd buoyancy on
    ball create radius 0.005 x 0.5 y 0.5 z 0.5
    ball attribute density 2500
    ball property kn 1e2 ks 1e2 fric 0.25
    set gravity 0 0 -9.81
    def fluid_time
      global fluid_time = mech.age
    end
    ball history id 1 zvelocity id 1
    history add id 2 fish @fluid_time
    plot clear
    plot add hist 1 vs 2
    plot add cfdelement shape arrow colorby vectorattribute "velocity"
    """.format(fluid_density, fluid_viscosity))

    element_volume = ca.volume()
    dt = 0.005

    for i in range(100):
        it.command("solve age {}".format(it.mech_age()+dt))
        print "sending solve time"
        cfd_link.send_data(dt) # solve interval
        cfd_link.send_data(ca.porosity())
        cfd_link.send_data((ca.drag().T/element_volume).T/fluid_density)
        print " cfd solve started"
        ca.set_pressure(cfd_link.read_data())
        ca.set_pressure_gradient(cfd_link.read_data())
        ca.set_velocity(cfd_link.read_data())
        print " cfd solve ended"

    cfd_link.send_data(0.0) # solve interval


    print "ball z velocity", it.ball.find(1).vel_z()

The file cfd_coupling.py is listed below. This files should be run in an instance of Python outside of PFC3D.

from itasca import p2pLinkClient
import numpy as np

rho = 1000.0 # Fluid Density
mu = 1e-3    # Dynamic viscosity

# fluid mesh
nodes = np.array([[ 0.0,  0.0,  0.0], [ 0.5,  0.0,  0.0], [ 1.0,  0.0,  0.0],
                  [ 0.0,  0.5,  0.0], [ 0.5,  0.5,  0.0], [ 1.0,  0.5,  0.0],
                  [ 0.0,  1.0,  0.0], [ 0.5,  1.0,  0.0], [ 1.0,  1.0,  0.0],
                  [ 0.0,  0.0,  0.5], [ 0.5,  0.0,  0.5], [ 1.0,  0.0,  0.5],
                  [ 0.0,  0.5,  0.5], [ 0.5,  0.5,  0.5], [ 1.0,  0.5,  0.5],
                  [ 0.0,  1.0,  0.5], [ 0.5,  1.0,  0.5], [ 1.0,  1.0,  0.5],
                  [ 0.0,  0.0,  1.0], [ 0.5,  0.0,  1.0], [ 1.0,  0.0,  1.0],
                  [ 0.0,  0.5,  1.0], [ 0.5,  0.5,  1.0], [ 1.0,  0.5,  1.0],
                  [ 0.0,  1.0,  1.0], [ 0.5,  1.0,  1.0], [ 1.0,  1.0,  1.0]])

elements = np.array([[ 1, 10, 13,  4,  0,  9, 12,  3],
                     [ 4,  5, 14, 13,  1,  2, 11, 10],
                     [ 4, 13, 16,  7,  3, 12, 15,  6],
                     [13, 16, 17, 14,  4,  7,  8,  5],
                     [10, 19, 22, 13,  9, 18, 21, 12],
                     [13, 14, 23, 22, 10, 11, 20, 19],
                     [13, 22, 25, 16, 12, 21, 24, 15],
                     [22, 25, 26, 23, 13, 16, 17, 14]], dtype=np.int64)

nele = len(elements)

with p2pLinkClient() as pfc_link:  # open connection to PFC
    pfc_link.connect("localhost")
    pfc_link.send_data(nodes)
    pfc_link.send_data(elements)
    pfc_link.send_data(rho)
    pfc_link.send_data(mu)

    while True:
        print "waiting for run time"
        deltat = pfc_link.read_data()
        if deltat == 0.0:
            print "solve finished"
            break
        print "got run time", deltat
        porosity = pfc_link.read_data()
        body_force = pfc_link.read_data()
        print "got runtime and data"

        print "sending data to pfc"
        pressure = np.zeros(nele)
        gradp = np.zeros((nele,3))
        velocity = np.zeros((nele,3))

        pfc_link.send_data(pressure)
        pfc_link.send_data(gradp)
        pfc_link.send_data(velocity)
        print "send finished"

Coupling PFC3D with OpenFOAM®

A community-driven project demonstrates using this methodology to couple PFC3D with OpenFOAM® CFD solvers. See https://github.com/jkfurtney/PFC3D_OpenFOAM for more details.

• This offering is not approved or endorsed by OpenCFD Limited, producer and distributor of the OpenFOAM software and owner of the OpenFOAM® and OpenCFD® trade marks.

• OpenFOAM® is a registered trade mark of OpenCFD Limited, producer and distributor of the OpenFOAM software.

• Documentation for OpenFOAM® can be found here: http://cfd.direct/openfoam/documentation/.

• Note that Itasca cannot provide support for OpenFOAM®.