CMT-nek Case

Background

This is an example that illustrates a two-dimensional multiphase shock tube. This example also illustrates linking ppiclF to the compressible Discontinuous Galerkin solver CMT-nek.

The core particle equations being solved in this case are

dXdt=V,MpdVdt=Fqs+Fpg+Fam+Fc,CpMpdTdt=Qqs

where, for each particle we have the position vector X, the velocity vector V, the particle temperature T, the particle mass Mp, the particle specific heat Cp, the viscous drag force Fqs, the pressure gradient force Fpg, the added-mass force Fam, the collision force Fc, and the quasi-steady heat transfer Qqs. The position and velocity vectors are given by

X=[XY],V=[VxVy].

The collision force in the same as in the DEM Packing 3D example. The drag force is complicated so it is not detailed here. It is given by the model found in Parmar et al. (2010).

The pressure gradient force is given by

Fpg=Vpp

where Vp is the particle volume and p is the fluid pressure at each particle’s position. The added-mass force is given by

Fam=VpCM(pd(ρfV)dt)

where CM is an added mass coefficient given by the Mach number corrections of Parmar et al. (2008) and the volume fraction corrections of Zuber (1964). For simplicity, we assume that the density of the fluid at each particle’s location ρf is relatively constant in time so that it can be pulled out of the time derivative.

The quasi-steady heat transfer is given by

Qqs=2πκDp(TfT)Φqs

where κ=Cp,f/Pr, Cp,f is the specific heat at constant pressure of the fluid at the particle’s location, Pr is the Prandtl number of the fluid at the particle’s position, Tf is the fluid temperature at the particle’s position, and Φqs is a correction factor for the heat transfer that is a function of the particle Reynolds number and the Prandtl number.

With some manipulation, the final form of the particle equations which are being solved are

dXdt=V,(CMVpρf+Mp)dVdt=FqsVp(1+CM)p+Fc,CpMpdTdt=Qqs

As a compressible fluid solver, CMT-nek without particles solves the Euler equations given by

ρft+(ρfu)=0,(ρfu)t+(ρfuu+pI)=0,(ρfE)t+(ρfuE+up)=0,

where ρf is the fluid density, u is the fluid velocity, E is the total energy of the fluid, and I is the identity tensor.

Similarly, the governing multiphase equations are given by Ling et al. (2016) and are

(ϕfρf)t+(ρfϕfu)=0,ϕfρf(ut+uu)+p=fpf,ϕfρf(Et+uE)+(ϕfup+ϕpvp)=gpf+epf.

Here, ϕf is the fluid volume fraction, fpf is the total hydrodynamic force from the particles on the fluid, gpf is the energy removed from the fluid by the work done by that force, and epf is heat transfered directly from the particles to the fluid. Following a similar process as Ling et al. (2016), these equations can be rearranged to yield

ρft+(ρfu)=Aϕf,(ρfu)t+(ρfuu+pI)=uAϕf+fpfϕf,(ρfE)t+(ρfuE+up)=EAϕfpϕfuϕfpϕf(ϕpv)+gpfϕf+epfϕf,

where the asterisk on the terms fpf and gpf indicate that they pressure gradient force is no longer explicitly included in the hydrodynamic coupling as it is accounted for implicitly. Additionally, we have

A=ρf(ϕft+uϕf).

Note in the present application, the term ϕf/t isn’t included, but may be incorporated with minimal effort. In the above form, it is clear that the original CMT-nek solver can be used with minimal modification when including multiphase coupling from ppiclF, since every particle contribution can be included in the formulation through volumetric source terms only.

User Interface

The H-File for this case (CMT-nek H-File) is given below and corresponds to the equations being solved and the property being stored for each particle.

#define PPICLF_LPART 1500
#define PPICLF_LRS 5
#define PPICLF_LRP 12
#define PPICLF_LEE 2700
#define PPICLF_LEX 5
#define PPICLF_LEY 5
#define PPICLF_LRP_INT 8
#define PPICLF_LRP_PRO 6

#define PPICLF_JX 1
#define PPICLF_JY 2
#define PPICLF_JVX 3
#define PPICLF_JVY 4
#define PPICLF_JT 5

#define PPICLF_R_JRHOP 1
#define PPICLF_R_JRHOF 2
#define PPICLF_R_JDP 3
#define PPICLF_R_JVOLP 4
#define PPICLF_R_JPHIP 5
#define PPICLF_R_JUX 6
#define PPICLF_R_JUY 7
#define PPICLF_R_JDPDX 8
#define PPICLF_R_JDPDY 9
#define PPICLF_R_JCS 10
#define PPICLF_R_JT 11
#define PPICLF_R_JE 12

#define PPICLF_P_JPHIP 1
#define PPICLF_P_JFX 2
#define PPICLF_P_JFY 3
#define PPICLF_P_JPHIPU 4
#define PPICLF_P_JPHIPV 5
#define PPICLF_P_JE 6

The first block of lines denote the pre-defined directives. These directives are the maximum number of particles per processor, 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 remaining blocks of lines indicate the index names of the particle variables being solved for, the property names for each particle, and the projected field names.

The F-File for this case (CMT-nek F-File) has meaningful information in every routine. For this case, the routines are similar to those defined in the Nek5000 Case. The routine ppiclf_user_SetYDot is more complicated due to the more involved force and heat transfer models described previously. Additionally, the routine ppiclf_user_MapProjPart includes six fields being projected, which are given in the table below:

Projection mapping in ppiclf_user_MapProjPart.

Projected Field (a(x))

Particle Property (A(i))

ϕp

Vp/Dp

fpf,x

(Fqs,x+Fam,x)/Dp

fpf,y

(Fqs,y+Fam,y)/Dp

epf+gpf

(Qqs+FqsV+FamU)/Dp

ϕpvx

VpVx/Dp

ϕpvy

VpVy/Dp

As mentioned previously, the projected forces are the added mass and the quasi-steady forces. Since in added mass force isn’t directly computed in the final rearranged form of the particle equations, we save dV/dt before computing the new ˙Y so that the entire added mass force may be computed. Additionally, the work done by the hydrodynamic forces and the heat transfer are stored in the same projected field. Similar to the Nek5000 Case, the projected fields are normalized by Dp in 2D. In 3D, this factor would not be included.

The External Interface calls for this example occur within the user initialization Nek5000 routine usrdat2 in the file mstube.usr with the minimum number of initialization and solve subroutines called. The majority of the solve routines are found in the routine cntchk, which is called every RK stage. The called ppiclF routines are nearly identical to the Nek5000 Case, but with different initial conditions for the particle placement. Additionally, the extra volumetric source terms are computed in cmtchk and added as forcing in the routine userf.

Similar to Nek5000 Case, 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 CMT-nek by issuing the following commands:

cd ~
git clone https://github.com/dpzwick/ppiclF.git           # clone ppiclF
git clone -b jason https://github.com/dpzwick/Nek5000.git # clone CMT-nek
mkdir TestCase                                            # make test directory
cd TestCase
cp -r ../ppiclF/examples/CMT-nek/* .                      # 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 mstube                                          # build CMT-nek and link with ppiclF
echo mstube > SESSION.NAME && echo `pwd`/ >> SESSION.NAME # create CMT-nek necessary file
mpirun -np 4 nek5000                                      # run case with 4 processors

Simulation Output

Previous work has been done with legacy codes and experiments for this setup. In this problem, a Mach 1.66 shock impacts a particle curtain of 2 mm at 20% volume fraction. Previous experiments and simulations (not ppiclF) can be found in Ling et al. (2016). The upstream and downstream fronts of the particle curtain from ppiclF/CMT-nek can be compared to the legacy code and experiments. For ppiclF, the upstream and downstream fronts are taken by simply computing the maximum and minimum particle positions. The figure below shows a good comparison.

../_images/ling_compare.png

Comparison of ppiclF/CMT-nek to previous experiments and a legacy code. ppiclF/CMT-nek is shown in the magenta curve. The background image is reproduced from Figure 14 from Ling, Y., et al. “Interaction of a planar shock wave with a dense particle curtain: Modeling and experiments.” Physics of Fluids 24.11 (2012): 113301, with the permission of AIP Publishing.