Nek5000 Case¶
Background¶
This is an example that illustrates a two-dimensional fluidized bed. This example also illustrates linking ppiclF to the incompressible, spectral element, fluid solver Nek5000.
The core particle equations being solved in this case are
where, for each particle we have the position vector \(\mathbf{X}\), the velocity vector \(\mathbf{V}\), the drag force \(\mathbf{F}_{qs}\), the collision force \(\mathbf{F}_{c}\), and the weight force \(\mathbf{F}_{b}\), defined by
and the same weight and collision forces as in the DEM Packing 3D example are used. For the drag force, the Gidaspow drag model has been used, which is
where \(\mathbf{U}\) is the fluid velocity vector interpolated to the particle’s coordinates, \(V_p\) is each particle’s volume, and \(\beta\) is a drag coefficient, computed by
where \(C^*_D\) is another drag coefficient, \(Re_p^* = \phi_f |\mathbf{U}-\mathbf{V}|D_p/ \nu_f\) is the volume weighted Reynolds number, and \(\nu_f = \mu_f/\rho_f\) is the fluid kinematic viscosity. The equation for \(C^*_D\) is
The fluid equations that Nek5000 solves in this example are
where \(\mathbf{u}\) is the fluid velocity, \(\rho_f\) is the fluid density, \(\mathbf{\sigma}_f\) is the Navier-Stokes fluid stress tensor, \(\phi_f\) is the fluid volume fraction, and \(\phi_p\) is the particle volume fraction (\(\phi_f + \phi_p = 1\)), and \(\mathbf{f}_{pf}\) is a particle-fluid coupling force.
The solution to these equations reside on a mesh within Nek5000. Since the particles are solved in the Lagrangian reference frame, their contributions on the mesh must be accounted for. Most notably, the explicit particle contributions from ppiclF to Nek5000 are \(\phi_p\) (and as a result \(\phi_f\)) and \(\mathbf{f}_{pf}\). The method by which these fields are obtained is Projection.
User Interface¶
The H-File for this case (Nek5000 H-File) is given below and corresponds to the equations being solved and the property being stored for each particle. Note that since \(g\) is constant, we do not included in in the list of properties.
#define PPICLF_LRS 4
#define PPICLF_LRP 6
#define PPICLF_LEE 1000
#define PPICLF_LEX 6
#define PPICLF_LEY 6
#define PPICLF_LRP_INT 3
#define PPICLF_LRP_PRO 3
#define PPICLF_JX 1
#define PPICLF_JY 2
#define PPICLF_JVX 3
#define PPICLF_JVY 4
#define PPICLF_R_JRHOP 1
#define PPICLF_R_JDP 2
#define PPICLF_R_JVOLP 3
#define PPICLF_R_JPHIP 4
#define PPICLF_R_JUX 5
#define PPICLF_R_JUY 6
#define PPICLF_P_JPHIP 1
#define PPICLF_P_JFX 2
#define PPICLF_P_JFY 3
The two blocks of lines denote the pre-defined and user-only directives. The pre-defined directives are in the top block and are the number of equations, the number of properties, the sizes of the overlap mesh, the number of interpolated fields, and the number of projected fields. The user-only directives are in the bottom block.
The F-File for this case (Nek5000 F-File) has meaningful information in every routine. The routine ppiclf_user_SetYdot is nearly the same as the DEM Packing 3D example but with an added drag model evaluation that is slightly more complicated than the Stokes 2D example. Also, the ppiclf_user_EvalNearestNeighbor routine is similar to the DEM Packing 3D example. The new addition is the mapping of particle properties to be projected in ppiclf_user_MapProjPart. With some study, it can be found that the three fields being projected in 2D are:
Projected Field (\(a(\mathbf{x})\)) |
Particle Property (\(A^{(i)}\)) |
|---|---|
\(\phi_p(\mathbf{x})\) |
\(V_p/D_p\) |
\(f_{pf,x}(\mathbf{x})\) |
\(-F_{qs,x}/D_p\) |
\(f_{pf,y}(\mathbf{x})\) |
\(-F_{qs,y}/D_p\) |
where \(D_p\) has been used to normalize the values in 2D. Note that the negative signs of the components of \(\mathbf{F}_{qs}\) were added when the forces were stored in the storage array ppiclf_ydotc at the end of the routine ppiclf_user_SetYdot.
The External Interface calls for this example occur within the user initialization Nek5000 routine usrdat2 in the file uniform.usr with the minimum number of initialization and solve subroutines called. In this case:
ppiclf_comm_InitMPI is called to initialize the communication,
ppiclf_comm_InitParticle is called with initial properites and conditions for the particles,
ppiclf_solve_InitGaussianFilter is called to initialize the fitler for projection to the overlap mesh,
ppiclf_comm_InitOverlapMesh is called to initialize the overlap mesh from Nek5000,
ppiclf_solve_InitNeighborBin is called with minimum interaction distance of the largest particle size,
ppiclf_solve_InitWall is called which sets a wall for the particles at the bottom of the domain,
ppiclf_solve_InitPeriodicX is called which sets periodicity in the x dimension along the domain.
Additionally, the solve routines are called every time step in the same file in various Nek5000 user routines. In this example,
ppiclf_solve_InterpFieldUser is called three times to interpolate the fields \(\phi_p\), \(u_x\), and \(u_y\) into the property array,
ppiclf_solve_IntegrateParticle is called to integrate the system at the current time step,
ppiclf_solve_GetProFldIJKEF is called to access the projected fields and use them locally in Nek5000 (force coupling and volume fraction effects).
Also, note that ppiclF has been linked with Nek5000 in the Nek5000 makenek compilation file through the following lines:
SOURCE_ROOT_PPICLF=$HOME/libraries/ppiclF/source
FFLAGS=" -I$SOURCE_ROOT_PPICLF"
USR_LFLAGS+=" -L$SOURCE_ROOT_PPICLF -lppiclF"
Compiling and Running¶
This example can be tested with Nek5000 by issuing the following commands:
cd ~
git clone https://github.com/dpzwick/ppiclF.git # clone ppiclF
git clone https://github.com/Nek5000/Nek5000.git # clone Nek5000
mkdir TestCase # make test directory
cd TestCase
cp -r ../ppiclF/examples/Nek5000/* . # copy example files to test case
cd ../ppiclF # go to ppiclF code
cp ../TestCase/user_routines/* source/ # copy ppiclf_user.f and PPICLF_USER.h to source
make # build ppiclF
cd ../TestCase
./makenek uniform # build Nek5000 and link with ppiclF
echo uniform > SESSION.NAME && echo `pwd`/ >> SESSION.NAME # create Nek5000 necessary file
mpirun -np 4 nek5000 # run case with 4 processors