4.1. Generalizing to an elasto-plastic framework

To accurately implement a model with the components described – allowing for the possibility of stopped material with an apparent yield stress, dense flow obeying a non-local rheology, and separation/reconsolidation of material – is a very challenging task. We have chosen to address these challenges by representing the model in an elasto-plastic framework, with details provided next. This approach, which we have used previously for a purely local material model (Dunatunga & Kamrin Reference Dunatunga and Kamrin2015), essentially turns our proposed rheology into a complex Maxwell fluid by putting the flow rule in series with an elastic element (a ‘spring’) whose stiffness vanishes when the material density drops below $\rho _c$. As long as the spring is stiff in the dense (positive pressure) regime, its effect on the flow is negligible. However, its presence not only assists in implementing the separation/reconsolidation rule, but allows us to model the onset of a true yield stress without the need for viscous regularization.

4.2. Summary of equations

We use the standard notation for continuum mechanics as defined in Gurtin, Fried & Anand (Reference Gurtin, Fried and Anand2010). The trace of a tensor $\boldsymbol {A}$ is given by $\operatorname {tr} \boldsymbol {A}$, the transpose by $\boldsymbol {A}^{\rm T}$, the inverse by $\boldsymbol {A}^{-1}$ and the deviator by $\boldsymbol {A}_0 = \boldsymbol {A} - \frac {1}{3} (\operatorname {tr} \boldsymbol {A}) \boldsymbol {I}$ in three dimensions. Scalars are represented by lowercase text, points and vectors by lowercase bold text and tensors are represented by uppercase bold text (except the Cauchy stress is $\boldsymbol {\sigma }$). The dot product over vectors is represented with a single centred dot, and simply multiplies the vectors together componentwise and takes the sum, resulting in a scalar (e.g. $\boldsymbol {a} \boldsymbol {\cdot } \boldsymbol {b} = \sum _{i} a_i b_i$). Tensorial contraction is represented with a colon, and multiplies tensors together componentwise and takes the sum, also resulting in a scalar (e.g. $\boldsymbol {A} \colon \boldsymbol {B} = \sum _{i,j} A_{ij} B_{ij}$). The spatial gradient and spatial divergence operators are given by $\operatorname {grad}$ and $\operatorname {div}$ respectively. Overdots indicate a material time rate of change of the underlying variable.

Define the strain-rate tensor $\boldsymbol {D}=(\boldsymbol {\nabla }\boldsymbol {v}+\boldsymbol {\nabla }\boldsymbol {v}^{\rm T})/2$ and the spin (vorticity) tensor $\boldsymbol {W}=(\boldsymbol {\nabla }\boldsymbol {v}-\boldsymbol {\nabla }\boldsymbol {v}^{\rm T})/2$ from the velocity field $\boldsymbol {v}$. From the stress tensor, define the equivalent shear stress $\bar {\tau }=\sqrt {\sum _{i,j}\sigma _{0 ij}\sigma _{0 ij}/2}$ and hydrostatic pressure $p=-(1/3)\sum _i\sigma _{ii}$, which are used to define the stress ratio $\mu =\bar {\tau }/p$. Let the specific body force (here from gravity) be denoted $\boldsymbol {b}$, the plastic part of the strain rate be denoted $\boldsymbol {D}^{p}$ (and elastic part $\boldsymbol {D}^{e}\equiv \boldsymbol {D}-\boldsymbol {D}^{p}$) and the shear and bulk elastic moduli be denoted $G$ and $K$.

With these definitions in hand, let $\varOmega$ be the material domain, let $\varOmega _d$ represent the subdomain where material is dense ($\rho >\rho _c$) and let $\varOmega _s$ represent the subdomain where material is separated ($\rho <\rho _c$). The continuum system we solve is summarized as follows:

For all $\boldsymbol {x}\in \varOmega$

1. (i) Balance of linear momentum, angular momentum and mass

   (4.1a–c)\begin{equation} \operatorname{div} \boldsymbol{\sigma} + \rho \boldsymbol{b} = \rho \dot{\boldsymbol{v}},\quad \boldsymbol{\sigma} = \boldsymbol{\sigma}^{\rm T},\quad \dot{\rho}={-}\rho\operatorname{div} \boldsymbol{v}. \end{equation}

For all $\boldsymbol {x}\in \varOmega _d$

1. (i) Elasticity relation

   (4.2)\begin{equation} \dot{\boldsymbol{\sigma}}-\boldsymbol{W}\boldsymbol{\sigma}+ \boldsymbol{\sigma}\boldsymbol{W} =K \operatorname{tr} (\boldsymbol{D}-\boldsymbol{D}^{p}) \boldsymbol{1}+2G(\boldsymbol{D}-\boldsymbol{D}^{p})_0. \end{equation}

2. (ii) Plastic flow rule

   (4.3)\begin{equation} \dot{\bar{\gamma}}^{p}\equiv \sqrt{\sum_{i,j}2D_{ij}^{p}D_{ij}^{p}}=g\mu,\quad \boldsymbol{D}^{p}=\dot{\bar{\gamma}}^{p}\frac{\boldsymbol{\sigma}_0}{2\tau}. \end{equation}

3. (iii) Non-local fluidity PDE

   (4.4)\begin{equation} t_0\dot{g}=A^{2}d^{2}\nabla^{2}g-\varDelta_\mu \left( \frac{\mu_s-\mu}{\mu_2-\mu} \right)g -b\sqrt{\frac{\rho_s d^{2}}{p}}\mu g^{2}. \end{equation}

For all $\boldsymbol {x}\in \varOmega _s$

1. (i) Open material condition:

   (4.5)\begin{equation} \boldsymbol{\sigma}=\boldsymbol{0}. \end{equation}

The above equations are solved simultaneously with the aforementioned (re)consolidation rule for material passing through $\rho =\rho _c$ and appropriate boundary conditions for $g$, tractions and velocities on walls and free surfaces. Note that, although the mathematical model is fully three-dimensional, our simulations employ 2-D plane-strain assumptions (that is, $D_{zz} = 0$, but $\sigma _{zz} \ne 0$) in our implementation to reduce unnecessary computations.

We would like to emphasize that we are primarily interested in solutions where the elastic moduli are high (correspondingly, the elastic time scale is small) and the fluidity time scale $t_0$ is low. The fluidity time scale $t_0$ is a real parameter of the model, but at present we have no experimental data to calibrate it, so we choose a low value for convenience. A low value of $t_0$ implies that the steady-state solution for $g$ is achieved much faster than the geometry of the system changes. This ensures that the system quickly determines if $g = 0$ (the clogged case) is the stable solution, or if a flowing $g \ne 0$ solution prevails. The $t_0 \dot {g}$ term is important, even if $t_0$ is small, as it allows the system to determine the stable solution on its own by advancing forward in time – if we remove the term entirely, giving a steady-state form, we would need to perform a perturbation analysis to determine which solution is the stable one (this is difficult for general geometries).

5.1. Introduction

The material point method (MPM), developed in Sulsky, Chen & Schreyer (Reference Sulsky, Chen and Schreyer1994), sits at the intersection of Lagrangian and Eulerian methods. In MPM, the state of a body is tracked in Lagrangian markers termed material points. These material points transfer their state to an ephemeral background mesh, which is then used to solve the continuum equations of motion. The deformation of this background mesh is used to inform each of the material points of their deformation, which is then used to perform the constitutive update on each material point and move the points, completing the step. At this point, the ephemeral mesh is no longer needed and can be reset. The key is that the state of the body is stored entirely on material points, but the solution for the equations of motion are done on a temporary mesh that resets itself each step. Thus, large deformations do not accumulate in the mesh, only on material points, and the material points maintain a consistent and feasible stress state given their individual deformation histories.

Due to these features, MPM has been used in a variety of granular materials problems. MPM has been used to model individual grains (Bardenhagen, Brackbill & Sulsky Reference Bardenhagen, Brackbill and Sulsky2000) and determine the stresses at contacts along force chains. Others have used continuum descriptions of granular material (e.g. from soil mechanics) to capture phenomena such as landslides (Mast Reference Mast2013; Abe, Soga & Bandara Reference Abe, Soga and Bandara2014; Bandara & Soga Reference Bandara and Soga2015), silo discharge (Wiȩckowski Reference Wiȩckowski1999, Reference Wiȩckowski2001, Reference Wiȩckowski2003; Wiȩckowski & Kowalska-Kubsik Reference Wiȩckowski and Kowalska-Kubsik2011) and column collapse (Mast et al. Reference Mast, Arduino, Mackenzie-Helnwein and Miller2015). Our previous work (Dunatunga & Kamrin Reference Dunatunga and Kamrin2015, Reference Dunatunga and Kamrin2017) used MPM in this manner with the local $\mu (I)$ rheology to make quantitative and qualitative predictions for well-known phenomena such as the steady-state velocity profiles in chute flow, aspect-ratio-dependent run-out scaling during column collapse, the power-law outflow scaling in silo discharge and the impact of an elastic intruder into a granular bulk. In this paper, we extend our previous method rather significantly to encompass non-local effects, to solve NGF problems with MPM. See Appendix A for a comparison of MPM to other methodologies.

While the theoretical underpinnings of MPM involve a similar construction to the finite element method (FEM), when writing a computer implementation it may be simpler to conceptualize the steps as mapping and inverse mapping operators. This is because material points in MPM are not constrained to a single element throughout the simulation, while basic FEM typically keeps the carriers of element state collocated with the Gauss points used during quadrature. For details on how this integration is performed and the variants of MPM which arise from those choices, see Appendix B.

An overview of the broad steps is presented in figure 1. As MPM has many variants, the procedure we used to advance an arbitrary time step $t^{n}$ to $t^{n+1}$ ($\Delta t = t^{n+1} - t^{n}$) is described in words in Appendix C. These steps can be applied with any constitutive model; we describe the specialization to the NGF model in further detail below.

Figure 1.

The procedure to advance an MPM step. First, state information such as mass, volume and momentum are mapped from the material points to the nodes via shape functions on an ephemeral background mesh (step 1). The mesh exists for a single step to solve the equations of motion (step 2), then the results are projected back to the material points again via the shape functions (step 3). Finally, the constitutive update is performed with the new deformation information to update material state variables and the stresses, completing the step. As the background mesh is ephemeral, it can be reset for the next step and cannot accumulate much total deformation.

5.2. Constitutive update for the NGF model

In general, the constitutive update takes deformation-related quantities and computes the stress state and any state variables. For our case, we are given the velocity gradient $\boldsymbol {L}$ over the time step, the density of the material point $\rho ^{n+1}$, the Cauchy stress at the beginning of the time step $\boldsymbol {\sigma }^{n}$ and the granular fluidity at the beginning of the time step $g^{n}$. We wish to calculate the Cauchy stress at the end of the time step $\boldsymbol {\sigma }^{n+1}$ and the granular fluidity at the end of the time step $g^{n+1}$.

As implied by the names, at a very high level, the difference between the non-local model and local models is that spatial gradient information is needed to compute the desired quantities at the end of the step. With the NGF model (and assuming an explicit update), this means we need to do something akin to the following:

1. (i) Consistently project the value of $g$ on each material point onto a mesh.

2. (ii) Use the mesh to obtain the discrete Laplacian of $g$ (call it $\delta ^{2} g$) at the elemental level.

3. (iii) Consistently project $\delta ^{2} g$ back to material points.

4. (iv) Use the value of $\delta ^{2} g$, $g$ and $\boldsymbol {\sigma }$ on the material points to compute the value of $g$ at the end of the step. Note that now $\delta ^{2} g$ can be treated computationally as just another state variable, as it is collocated with the other quantities, although it came from non-local information (based on the state of nearby material points).

5. (v) Use the value of $g$ at the end of the step to determine the equivalent plastic strain rate $\dot {\bar {\gamma }}^{p}$ at the end of the step, which is enough to compute the final stress state.

We need to solidify these steps in an actual numerical implementation, as there are several complications and nuances to implementing the relations and projections. For pedagogical purposes we do not consider substepping yet, although those details will be discussed later. A detailed explanation of each step of the constitutive update we use is provided below.

1. (i) For each material point, check the pressure at the beginning of the time step $p^{n} = - \frac {1}{3} \operatorname {tr} \boldsymbol {\sigma }^{n}$ and the density at the end of the time step $\rho ^{n+1}$. If either of these indicate disconnected material ($p = 0$ or $\rho \le \rho _c$), the material point does not contribute to the elemental average of $g$ during this step. Otherwise, the material is dense during the step, and we need to add its contribution to the background grid. As we implement volumetric convected particle domain interpolation (see Appendix B), each material point has both a physical volume and a numerical extent over which it is integrated. If the material point is entirely contained within an element, we simply add the mass-weighted value of $g$ to an elemental $g$ accumulator and separately track the elemental mass. When a material point straddles element boundaries, only the fraction present in a given element contributes to these quantities (the remainder is added to the other elements the material point overlaps). After all material points are accounted for, we divide the mass-weighted $g$ quantities by the total mass in the element to obtain an estimate for $g$ in that element, which is then used for our spatial gradient calculations. This splitting of material point extent was chosen to reduce the cell crossing error, but we found that simply choosing to accumulate $g$ to whichever element contained the material point centre produced nearly identical results. We suspect, however, that convergence of results may be affected, similar to how undeformed generalized interpolation material point (uGIMP) and MPM convergence rates can differ.

2. (ii) Once all dense material points have contributed to the elemental $g$ field (obtained by dividing the accumulated mass-weighted value by the mass of contributing material points within that element), we can consider each element centre as a node and use a finite-difference stencil for the Laplacian operator to calculate

   (5.1)\begin{equation} \delta^{2} g_{i,j} = \frac{1}{(\Delta y)^{2}} (g^{n}_{i+1,j} - 2 g^{n}_{i,j} + g^{n}_{i-1,j}) + \frac{1}{(\Delta x)^{2}} (g^{n}_{i,j+1} - 2 g^{n}_{i,j} + g^{n}_{i,j-1}). \end{equation}
   In this case, we use the comma not to indicate differentiation, but simply as an index based on the logical coordinates of the element ($i$ indicating a row index and $j$ indicating a column index). In this part of the routine, we use the value of $g$ at the beginning of the time step when calculating this quantity, so all terms on the right-hand side are known.

3. (iii) Now iterate over all material points again; compute the elastic trial stress, $\boldsymbol {\sigma }^{tr}$, which is what the updated stress would be if $\boldsymbol {D}^{p}=\boldsymbol {0}$ (i.e. a purely elastic update). While plastic flow can modify the stress deviator from the trial value, in the dense regime the pressure is determined directly from the trial stress since our model assumes no plastic dilation. For any points with $\operatorname {tr} \boldsymbol {\sigma }^{tr}>0$, set the stress at the end of the step to $\boldsymbol {\sigma }^{n+1} = \boldsymbol {0}$. These material points are disconnected and not involved in further calculations within this time step, although they may reconsolidate and participate in later ones. Note that this check is necessary to catch all open material only for numerical reasons, as the pressure and density may be slightly mismatched since they evolve separately. For all other points, set $p^{n+1}=-\frac {1}{3}\operatorname {tr}\boldsymbol {\sigma }^{tr}$.

4. (iv) Recall that the evolution equation for $g$ can be written as

   (5.2)\begin{equation} t_0 \dot{g} = A^{2} d^{2} \nabla^{2} g - \varDelta_\mu \left(\frac{\mu_s - \mu}{\mu_2 - \mu} \right)g - \frac{\varDelta_\mu}{I_0} \sqrt{\frac{\rho_s d^{2}}{p}} \mu g^{2}. \end{equation}
   We may create a split numerical update for this equation into the following two parts via an intermediate variable $g^{*}$, resulting in the following expressions:
   (5.3)$$\begin{gather} \frac{g^{*} - g^{n}}{\Delta t} = \frac{1}{t_0} A^{2} d^{2} \delta^{2} g_{i,j} \end{gather}$$
   (5.4)$$\begin{gather}\frac{g^{n+1} - g^{*}}{\Delta t} = \frac{1}{t_0} \left( - \varDelta_\mu \left( \frac{\mu_s - \mu^{n+1}}{\mu_2 - \mu^{n+1}} \right)g^{n+1} - \frac{\varDelta_\mu}{I_0} \sqrt{\frac{\rho_s d^{2}}{p^{n+1}}} \mu^{n} g^{n} g^{n+1} \right). \end{gather}$$
   Note that the first expression is easily solved for the value of $g^{*}$, as the discrete Laplacian was obtained from $g$ earlier via (5.1).

5. (v) Recalling the plastic flow rule, we can perform an implicit update by setting $g^{n+1} \mu ^{n+1} = (\dot {\bar {\gamma }}^{p})^{n+1}$. The plastic shearing $(\dot {\bar {\gamma }}^{p})^{n+1}$ is used to correct the elastic trial stress, giving the equivalent shear stress at the end of the step via

   (5.5)\begin{equation} \bar{\tau}^{n+1} = \frac{\bar{\tau}^{tr} p^{n+1}}{G g^{n+1} \Delta t + p^{n+1}}. \end{equation}

6. (vi) After obtaining $g^{*}$, note that (5.4) is a semi-implicit update for the value of $g^{n+1}$. Splitting the $g^{2}$ term as shown into a beginning-of-step and end-of-step value allows us to obtain an equation that is quadratic in $g^{n+1}$. After substitution of $\mu ^{n+1}$ obtained from (5.5), the quadratic we obtain is

   (5.6)\begin{equation} \left.\begin{gathered} 0 = \bar{a} (g^{n+1})^{2} + \bar{b} g^{n+1} + \bar{c},\quad \text{where} \\ \bar{a} = (\mu_2 - \mu_s f) \bar{G} \sqrt{p^{n+1}} + \bar{G} \mu_2 f \bar{\xi}, \\ \bar{b} = (S_2 - \bar{\tau}^{\mathrm{tr}} - \bar{G} \mu_2 g^{*} + S_0 f - \bar{\tau}^{\mathrm{tr}} f) \sqrt{p^{n+1}} + (S_2 - \bar{\tau}^{\mathrm{tr}}) f \bar{\xi}, \\ \bar{c} = (\bar{\tau}^{\mathrm{tr}} - S_2) g^{*} \sqrt{p^{n+1}}, \end{gathered}\right\} \end{equation}
   where
   (5.7a–e)\begin{align} \bar{G} = G \Delta t,\quad S_0 = p^{n+1} \mu_s,\quad S_2 = p^{n+1} \mu_2,\quad f = \Delta t \varDelta_\mu,\quad \bar{\xi} = \frac{\dot{\overline{\gamma_n}}^{p}}{I_0} \sqrt{\rho_s d^{2}}. \end{align}

   Picking the correct root of the quadratic involves some care. We have the physical constraints that $\mu \leq \mu ^{tr}$ and $\mu < \mu _2$. The form of the equation implies that there should be only one solution for $g$ which meets both these constraints, so we solve for both using the familiar expression for the roots of a quadratic equation

   (5.8)\begin{equation} g^{n+1}_{{\pm}} = \frac{-\bar{b} \pm \sqrt{\bar{b}^{2} - 4 \bar{a} \bar{c}}}{2 \bar{a}}. \end{equation}
   As in Dunatunga & Kamrin (Reference Dunatunga and Kamrin2015), we may choose to evaluate the roots in mathematically equivalent alternative forms (depending on the signs of $\bar {a}$, $\bar {b}$ and $\bar {c}$) for numerical stability reasons. We then back substitute into (5.5) to obtain $\bar {\tau }^{n+1}$, which in turn is used to compute $\mu ^{n+1}$ and select the correct root.

7. (vii) Once the correct $g^{n+1}$ is selected on the material points, use the corresponding equivalent shear stress at the end of the step to scale the deviator of the trial stress via

   (5.9)\begin{equation} \boldsymbol{\sigma}^{n+1} = \frac{\bar{\tau}^{n+1}}{\bar{\tau}^{tr}}\boldsymbol{\sigma}^{tr}_0 - p^{n+1} \boldsymbol{I}, \end{equation}
   and obtain the stress at the end of the step $\boldsymbol {\sigma }^{n+1}$.

The above explanation does not include substepping, which is used to avoid unnecessarily small time step sizes due to small values of $t_0$ compared with the elastic time scales. To implement substepping, we simply use a smaller time step for the constitutive update alone and repeat the integration until the end of the overall MPM step. The velocity gradient $\boldsymbol {L}$ is held constant throughout the step, but the stress $\boldsymbol {\sigma }$ does update (and as the constitutive update's time step $\Delta t_{{substep}}$ is smaller, the trial stress differs from the case without substepping). The computational savings compared with having an overall smaller $\Delta t$ occur mainly because nodal calculations are elided and because the mapping operations between material points and mesh nodes are reduced to only the initial and final substeps, whereas new mapping matrices must be constructed and used for every full MPM step.

7.1. Silo geometry and numerical parameters

Our implementation is in 2-D plane strain, corresponding to a thick silo in the out-of-plane direction. When choosing the geometry of the mesh, note that the silo needs to have enough elements across the orifice to initiate and resolve the flow (in a flowing case), but to reduce computation time we would like to keep this number manageable. We chose 8 units across the orifice, with each unit being a 0.0175 by 0.0175 m square. Although the background mesh is ephemeral in MPM, since we solve the equations of motion over this mesh it does have an influence on the space of possible solutions on the material points. As with the inclined chutes, each square unit contained 16 material points (four per spatial dimension). In order to reasonably match the geometry of an arch, each unit consists of two triangular elements with the diagonal going from bottom left to top right for the left half of the silo and vice versa for the right half (see figure 3b). We show later that results are insensitive to the exact number of elements as long as there are enough to resolve flow features around the orifice. A full convergence study is not attempted in the present work, as the computational time and space required quickly becomes prohibitive due to the explicit parts of the scheme and the number of auxiliary variables respectively.

Figure 3.

The geometry of the silo is shown in (a). We kept the height $H$, width $L$ and orifice width $W$ consistent throughout the simulations and changed the material parameters, mainly grain size $d$ and non-local amplitude $A$. The background mesh used by the MPM scheme for coarse cases is shown in (b) to aid in understanding the apparent shape of the arch in static cases.

Material parameters are presented in table 1, which are largely consistent with values commonly used for a 3-D packing of glass beads. Of note, however, is the large value of $\mu _2$, explained as follows. At a jammed arch in a plane-strain silo, the stress state must be compressive along the arch, but vanish normal to it because it is a free surface. Thus, the stress state at the arch must attain $\mu = 1$; this fact is true regardless of the constitutive model, and we have observed this in our own results. Note that this does not apply to flowing silos. In the flowing case, a local $\mu (I)$ model would give that if any shear is present on the boundary between dense and disconnected states then $\mu = \mu _2$, and our simulation outputs do indicate $\mu \approx \mu _2$, although we are unsure if a non-local model must also attain exactly $\mu _2$ or simply a value greater than one. To properly capture the clogged stress state, it follows that $\mu _2$, the feasible upper limit for $\mu$ in models like ours that extend the $\mu (I)$ rheology, cannot be less than $1$. The constraint is slightly weaker for fully 3-D silos with a hole orifice, where $\mu _2$ must exceed $\sqrt {3}/2\approx 0.866$ to admit a jammed arch state (and an overall clogged silo). Both these values are higher than, for example, $\mu _2=0.6453$ reported in Jop et al. (Reference Jop, Forterre and Pouliquen2005), but recall from the section on inclined chutes that there is ambiguity on the value of the upper limit of $\mu$ and little experimental data in that regime. We observed in our simulations that above $\mu _2\sim 10$, the solution becomes insensitive to the value of $\mu _2$ so we simply set $\mu _2$ to an asymptotically high limit, giving effectively a linear $\mu (I)$ fit in the local limit.

Table 1.

Common material parameters used for the silo simulations. Our simulations vary $A$, hence it is not present in this table, but instead will be specified on a case-by-case basis.

Each simulation is run with the same orifice size $W$, given by 0.014 m (measured from fixed node to fixed node), but the grain diameter $d$ is varied via the material properties to set the ratio $W / d$. The total width of the silo $L$ is fixed at 0.595 m (34 units) and, importantly, the same background mesh and material point density is used for all simulations so as to control the possible influence of numerical length scales. The initial fill height $H$ is also fixed at 0.595 m. Initially, the material point configuration is assigned a lithostatic pressure; the density of each material point is set to be commiserate with this pressure by computing $\rho = \rho _c K / (K - p)$, where $K$ is the bulk modulus and $p$ is the pressure of that material point. The density is initially set by changing the mass of the material point – all material points have the same initial volume in these simulations. However, as in Dunatunga & Kamrin (Reference Dunatunga and Kamrin2015), the true volume of material points evolves over time (although the numerical extent they are integrated over does not, see Appendix B).

Similar to the inclined chutes, we set the bottom of the silo to have Dirichlet $g = 0$ and $\boldsymbol {u} = \boldsymbol {0}$ (no-slip) boundary conditions – note these are only applied along the solid portion of the boundary (i.e. they do not apply directly over the orifice, which instead allows the granular material to act as a free surface with $\partial _n g = 0$ and a traction-free condition). The sidewalls allow friction-free sliding (perfect slip), but do not allow material points to go through the boundary; that is $v_n = 0$ while the transverse direction is unhindered, and $\partial _n g = 0$.

No motion is imparted to the material points as an initial condition. Instead, the retaining wall at the orifice is removed upon the first step of the simulation and flow is allowed to develop ‘naturally’. Recall from the section on inclined chutes that we initialize $g$ by running the local model for a short period of time. As the local model will always produce a flow field regardless of the ratio of orifice to grain size, this procedure will always result in a non-zero $g$ field. After initializing with the local model, the numerical procedure detailed above for the non-local model is used to determine if the solution decays towards a stable $g = 0$ state corresponding to a clogged silo or evolves to obtain some other steady flowing state.

We found the results largely insensitive to the specifics of how the seeding procedure was performed; changing from a local model $\mu (I)$ to a simple rate-independent Drucker–Prager material only slightly affected transient states after switching to the non-local procedure, and decreasing or increasing the local model run time similarly had almost no discernible effect on the final state. There were two notable exceptions to this. First, if the material points deviate too much from their initial positions, finding a non-flowing solution became difficult, at least at the spatial and temporal sampling parameters we had used. This is likely because a strong arch of material points needs to develop simultaneously without any holes, and as points move further from their initially ideal placement, numerical noise in MPM makes this unlikely. Refinement would likely mitigate this issue, as the crossing noise can be reduced and the additional material points are more likely to sample a good configuration for numerical purposes, but we did not explore this further. Secondly, if the initialization process were done too quickly, the flow field (and therefore $g$) may not have enough time to establish through the system and ends up too small in magnitude to seed the simulation in light of numerical damping.

Outside of initialization, elastic waves through the system may serve to jostle a clogged silo away from that state. These elastic waves may cause fluctuations in $\mu$, which can interfere with the evolution of $g$. While some of these effects are physical, real systems also include elastic-wave damping effects not included in our model; moreover, we are most interested in the elastically rigid limit where elastic waves would not occur anyway (the motion would simply be unsupported by the system). With a high (but finite) wave speed, the time scale for inertial variation is far larger than that of elastic waves, so even a small amount of dissipation should cause elastic waves to dissipate quickly relative to the amount of time it takes for macroscopic granular flow features to appear or change. Adding bulk viscosity could be one approach to adding a controlled dissipative term, but after observing preliminary results with this enabled, we felt it interfered too much with the solution. Instead, we implemented an elastic-wave absorbing boundary by zeroing out the speed of any material point moving upward that is located at $y$-position 0.49 m (28 units) or higher. Upward motion is entirely caused by elastic waves in this problem since gravity points down. This absorbing condition is located sufficiently far from the orifice that we believe it should not affect any true dynamics near the opening when the system clogs or in flowing cases while serving to remove unwanted elastic waves. Results were largely insensitive to the position or strength of the penalty (e.g. moving the dividing line down by 8 cm did not change the results of our analysis, nor did implementing the penalty as a halving of upward velocity instead of zeroing).

7.2. Silo clogging phase diagram

Even if it may be obvious to the eye, determining a precise criterion for when a silo configuration is clogged in a numerical simulation is a non-trivial exercise. It can be comparable to studying avalanche effects and the presence of open/falling material can trick kinetic-energy-based criteria for solidification (Aranson & Tsimring Reference Aranson and Tsimring2001; Martin et al. Reference Martin, Dubois, Monerie and Radjai2009). Instead, we use a heuristic to mark states as ‘likely static’, ‘likely flowing’ and ‘unsure’. If the average velocity of the dense material is small and at most exhibits decaying oscillatory waves, we suspect that the material is in a static configuration. Additionally, if the mass flow rate out of the orifice corresponds to only a few grains per second, we can be more sure that the configuration is static. Conversely, if the average velocity points downward, is large, and does not change size much, we suspect that the material is flowing; we can be even more sure if a large number of grains are leaving per unit time. To make this more precise, in our simulations we check if the most recent zero crossing of the vertical component of velocity happened within the last ten per cent of the simulation; if so, we give a point in favour of a static configuration. If the last zero crossing happened earlier than half-way through the simulation, we give a point to the flowing state; in between, we are unsure. We then check the number of material points leaving between each saved frame of the simulation; if this is less than or equal to one on average in the last tenth of the simulation (with our parameters, this is a mass flow rate of 28.7 g s$^{-1}$ per unit width), that is a point in favour of a static configuration. If this is greater than ten per frame on average, that is a point in favour of a flowing simulation; in between, we are unsure. Both tests must agree (i.e. two points are required for us to mark the simulation as static or flowing) – any other state is marked as unsure. Visual inspection of the last few frames of a selection of simulations was also done for confirmation.

Simulations took approximately 2 hours on an Intel Silver 4112, with faster-flowing simulations finishing earlier as they had fewer material points over time. In all silo simulations, we used a substepping time of $\Delta t_{{substep}} = \Delta t / 150$, which we estimate from smaller cases reduced the wall time by a factor of two. Using the criteria outlined above, we plotted these data on a lattice with simulations having varying $A$ and $W/d$ ratios (as $W$ is fixed, we simply changed the material grain size $d$ to achieve the desired ratios). The value of $A$ is either $0.1$, $0.5$, $1.0$, $1.5$, $2.0$, $2.5$ or $3.0$ while $W/d$ must be either $1$, $2$, $3$, $4$ or $5$. The results are shown in figure 4(a). Note that we see a separation into three regions (flowing, static and unsure) which can be well described simply by two lines of different slope; the ratio of $W/d$ to $A$ seems sufficient to predict if a simulation will flow or reach a static state.

Figure 4.

(a) Silo phase diagram varying non-local amplitude $A$ and the relative opening size $W/d$. A blue circle means the simulation flows, a grey square means the simulation attains a static configuration and an orange triangle indicates that we were unsure how to classify the simulation – see the main text for details on the classification criteria. Note that the geometric properties of the simulation were held constant and that only the indicated material properties of the NGF model (that is, the grain size $d$ and $A$) were varied here. All other material properties were held constant for these tests. Note that the diagram can be partitioned into the three regions with only two lines which go through the origin. (b) Still frames of the arch formation process for the $A = 1.0$, $W/d = 1.0$ case, which is a static simulation in our classification.

This is not too surprising if we consider the form of the NGF PDE and the set-up we used. Note that the term $d$ almost always appears with $A$ – it only shows up alone in the $g^{2}$ term for describing the local rheology, which we believe only has a minor effect on the stability of the solution. Since $W$ is fixed in these simulations, lines of constant $Ad$ should therefore exhibit very nearly the same stability behaviour, as the only difference would appear in the $g^{2}$ term. In our graph, starting at the origin, we can fan out to form these lines of constant $Ad$, and we do see a physically reasonable boundary between flowing and non-flowing cases. Keep in mind that the ‘unsure’ cases and region correspond to error bars in our classification – we believe the ‘true’ boundary is indeed nearly a straight line within this uncertainty cone relating $Ad$ to $W$. As a benchmark against experimental data, note that at $A = 0.5$ (an approximate value for 3-D beads) the flowing/static boundary occurs somewhere between $W/d = 1$ and $W/d = 1.5$. This corresponds well with experimentally observed values for the critical $W/d$. In the flat-bottomed quasi-2-D silo, which we are attempting to model, $W/d\approx 1.0$ was obtained for beads in Choi et al. (Reference Choi, Kudrolli and Bazant2005). This is similar to the range of 1.0–1.16 that was obtained in Mankoc et al. (Reference Mankoc, Janda, Arevalo, Pastor, Zuriguel, Garcimartín and Maza2007) for beads, albeit for a circular orifice. As before, we note that since we have a continuum solver, we achieve the static configuration when the ensemble average behaviour would be to form a static arch. In a single experiment, it may be possible to observe orifice openings of larger size that nevertheless still can find a static configuration, but these become rare as the ratio $W/d$ increases. Figure 4(b) shows an example of an arch forming in one of our continuum simulations (where $\Delta t_f = 1/240$ s) when the silo opening is sub-critical, with material points below the arch falling out leaving the rest of the silo material intact as the fluidity decays to zero. The numerical procedure models the outline of the material domain from edges of the background mesh geometry, which leads to a ‘polygonal’ arch shape; as the mesh refines more, a rounder arch shape emerges (see § 7.3 and Supplemental Movie available at https://doi.org/10.1017/jfm.2022.241). Next, in figure 5(a), we observe the pressure distributions in the silos for three configurations $W/d = 1$, $2$ and $5$ with a common $A = 1$. The fields are plotted near the end of the simulation and use a spatial and temporal average over frames from $350 \Delta t_f$ to $360 \Delta t_f$. The spatial smoothing was done as in Andersen & Andersen (Reference Andersen and Andersen2009) to project point-wise pressure data onto the background mesh – note that this is a purely visual post-processing technique and does not affect the fidelity of the simulation itself. In the clogged configurations ($W/d = 1$ and $2$), we see that the pressure immediately above the arch is low, but does not go to zero. The simulation which is closer to the static/flowing boundary ($W/d=2$) experiences a lower pressure in the region immediately above the arch, although some distance away from the opening the pressure looks essentially the same as the other static case. In contrast, the flowing simulation ($W/d = 5$) contains many materials points in the region above the orifice where material points are in free fall, and the pressure indeed goes to zero when approaching this region since the material points are stress free. The peak pressure is also smaller, although the fill height has already decreased substantially by the time this snapshot was taken due to the flow. Figure 5(b) uses the same data, but shows the pressure along the vertical centreline of the silos, confirming our observations from the patch plots.

Figure 5.

(a) Visualization of pressure in silos with $A = 1$ and $W/d = 1, 2$ and $5$ (left, centre and right, respectively). These correspond to data points in the $A=1$ column of figure 4(a). The pressure shown here is a temporal average projected onto the mesh over several neighbouring frames from the view of grid elements (from $350 \Delta t_f$ to $360 \Delta t_f$ inclusive). White corresponds to a completely empty element. (b) Plots of the pressure along the centreline for silos with $A=1.0$ and $W/d$ ratios of 1.0, 2.0 and 5.0 (pink, green and blue, respectively). While the pressure decreases drastically directly above the orifice in the clogged configurations, it does not reach zero. In contrast, the flowing case has much lower pressure and is actually zero in the region above the opening.

Now we will turn our attention more towards the flowing simulations. In figure 6 we compare two silo flows which have the same material parameters but different values of $A$. Snapshots are taken at slightly different times to account for the difference in outflow rates and were chosen to keep a similar amount of mass remaining in the silo; the snapshot for the $A = 0.1$ case is taken at $70 \Delta t_f$ and for the $A = 1.0$ case is taken at $86 \Delta t_f$, where again $\Delta t_f = 1/240$ s. In the top half of figure 6, we see the velocity fields plotted over material points and velocity contours. The spacing of the contour lines indicates that shear zones in the $A = 0.1$ case (left) are thinner than in the $A = 1.0$ case (right). This is confirmed by plotting the vertical velocity across horizontal slices at heights $1.5W$, $2.0W$ and $2.5W$ in the lower half of the figure – the flow is more plug like in the low $A$ case, and is noticeably similar to a local $\mu (I)$ model solution, which would predict plug flow down the centre with a thin shear band before dropping to zero vertical velocity close to the sidewalls where $\mu < \mu _s$ (Kamrin Reference Kamrin2010; Staron et al. Reference Staron, Lagrée and Popinet2012; Staron, Lagrée & Popinet Reference Staron, Lagrée and Popinet2014). In contrast, the higher $A$ value results in a spreading of the velocity field as expected due to the non-local effect, which will imply a length scale for shear bands. As both cases do actually have non-zero $A$, we do observe a small amount of flow even where $\mu < \mu _s$ as expected, and less flow under $\mu _s$ is observed for the smaller $A$. The decay length scale depends on the size of $A$, and we find that larger $A$ diffuses out the velocity more readily, also as expected.

Figure 6.

Two silos with the same parameters and initial conditions except for the value of $A$. Both simulations have nearly the same amount of mass remaining in the silo when imaged ($70 \Delta t_f$ left, $86 \Delta t_f$ right). The velocity cuts off more sharply in the low $A$-value case as expected. Iso-speed lines (drawn at 0.05, 0.10, 0.15, 0.20, 0.25, 0.30 and 0.35 m s$^{-1}$) are drawn over a scatterplot of the speed confirm that the velocity is diffusing more when $A=1.0$ compared with when $A=0.1$. Vertical velocity profiles are drawn at the heights indicated as horizontal lines in the first row of images located at heights $1.5W, 2.0W$ and $2.5W$ above the orifice. The grain size is given by $d = 0.028$ m, and the orifice size is $W = 0.14$ m resulting in $W/d = 5$.

We also look at the form of the vertical velocity profiles obtained. In order to compare with existing data, which uses glass beads primarily, we choose $A = 0.5$ and set $W/d = 5$. We waited for steady-state flow, and then looked at vertical velocity slices at the same three heights as before ($1.5W$, $2.0W$ and $2.5W$); results are plotted in figure 7. Many have noticed in experiments that the downward velocity profile, in a region not too far from the opening, appears to spread diffusively but with height playing the role of ‘time’ — that is, the downward velocity appears Gaussian with $v_y=C \exp {{-x^{2}}/{4By}}$ where $B$ controls the width of the Gaussian and $y$ is the height above the orifice (Tüzün & Nedderman Reference Tüzün and Nedderman1979; Medina et al. Reference Medina, Cordova, Luna and Trevino1998; Samadani, Pradhan & Kudrolli Reference Samadani, Pradhan and Kudrolli1999; Choi et al. Reference Choi, Kudrolli and Bazant2005; Bazant Reference Bazant2006; Rycroft et al. Reference Rycroft, Bazant, Grest and Landry2006; Zuriguel et al. Reference Zuriguel, Maza, Janda, Hidalgo and Garcimartín2019). The apparently diffusive character of the flow field was central to some of first silo flow models (Mullins Reference Mullins1972; Nedderman & Tüzün Reference Nedderman and Tüzün1979). Interestingly, data from our simulations also match this observation, see figure 7, with two caveats. First, we have an additive constant at each height due to our walls having no friction (and thus allowing downward motion at the walls); in experiments, wall friction ensures velocity decays to zero at the sidewalls. Secondly, the width parameter that we obtain, $B = 0.6d$, is a smaller than the typical values observed in the literature ($1d$–$5d$ in Kamrin & Bazant Reference Kamrin and Bazant2007). We believe this is due to the relatively large orifice size we used compared with the geometry of our system (unlike the experiments which used a much smaller opening relative to the silo width) as well as the aforementioned lack of sidewall friction. Recent studies (Pongó et al. Reference Pongó, Stiga, Török, Lévay, Szabó, Stannarius, Hidalgo and Börzsönyi2021) have shown the flow rates are also a function of particle stiffness. Possibly related effects such as layering near boundaries are not presently incorporated in the continuum model, so exploration of boundary conditions and edge effects are natural areas for future work. Finally, we plot the $g$ field in these same simulations, shown in figure 8 (note that these cases correspond to the top left of figure 4(a) – all of them flow, even $A = 1.0$). The plots show averages over both elements and temporal frames ($280 \Delta t_f$ to $300 \Delta t_f$) to understand what the numerical scheme regularly observes when computing the spatial gradients. Note that the lower values of $A$ show less spreading of the $g$ field and a slightly larger flow rate compared with higher values of $A$ (indicated by the height of the free surface) as expected. Considering $g$ as a surrogate for equivalent plastic shear strain rate, we can also observe the distinctive ‘arms’ of high shearing emanating from the orifice. The region just above the orifice where grains are in free fall is also clearly visible in these plots, indicated by a dark brown colour ($g$ is set to zero in the disconnected state, but those material points do not contribute to the solution of the NGF equation).

Figure 7.

Cross-sections of vertical velocity taken at horizontal slices at height $1.5W$, $2.0W$ and $2.5W$ (pink, green and blue colours, respectively) in a silo with the same parameters as figure 6 except $A = 0.5$ (chosen to match experimental material). Simulation data are plotted with dots, while a diffusing Gaussian fit (with offsets) is plotted as the lines. The value of $B$, which controls the Gaussian's width at each height, is given by $0.6d$ ($d$ grain diameter).

Figure 8.

Plots of $g$ (logarithmic scale) of a silo with $W/d = 5$ and the same material parameters as previously used except in the value of $A$, going from $0.1$ to $0.5$ to $1.0$ from left to right. Our averaging procedure assigns the value of $g = 0$ to open material; the region immediately above the orifice appears to have low $g$ value, although it may be more correct to say there is no value of $g$ at those points, as they are stress free in the disconnected state in our model. All plots are taken over an average of frames from $280 \Delta t_f$ to $300 \Delta t_f$ and use the same mesh-wise projection technique as was used to display the pressure in figure 5. Note that the lower values of $A$ result in less spreading of the $g$ field, as anticipated.

7.3. Convergence test and caveats

Convergence testing is difficult due to the amount of wall time each simulation takes and the total amount of data generated. Performing a (linearized) von Neumann stability analysis on the scheme used for evolving the NGF PDE shows that the stable $\Delta t$ is proportional to $\Delta x^{2}$ due to the Laplacian term (Dunatunga Reference Dunatunga2017). While the substepping is designed to mitigate the impact of this harsh Courant–Friedrichs–Lewy (CFL) condition, running a refined simulation still takes substantially more CPU and wall time than a coarse one.

To at least show qualitative agreement, we performed another set of simulations at twice the spatial resolution of the previous cases (and correspondingly a time step of one quarter the original size). As with the other simulations, the mesh consists of triangles stacked to create squares, and which diagonal used depends on which half of the domain the element is in. Numerical parameters are summarized in table 2. We choose $A = 1$ and other material parameters are given as before in table 1. In summary, we consider a total of four simulations in this section; we look at all combinations of two different refinement levels and two different grain diameters. We wish to observe agreement in the silo geometry as well as fields such as the velocity and pressure.

Table 2.

Numerical parameters used for the refined and coarse simulation cases used in the convergence study. The only differences are that the element size is halved and the refined time step size is correspondingly cut down to a quarter of the coarse value. The material properties are given in table 1.

Our coarse simulations use the same data as before, with a grid size of $\Delta x = 0.0175\,\textrm {m}$, a time step size of $\Delta t = 1\times 10^{-4}$ s and take around 2 hours to complete on an Intel Silver 4112. In contrast, the refined simulations use a grid size of $\Delta x = 0.00875$ m, a time step size of $\Delta t = 2.5\times 10^{-5}$ s and take approximately 30 hours to complete on the same machine; both refined simulations were run at the same time, with the only difference between them being the grain diameter $d = 0.014$ and $d = 0.0028$ m (resulting in $W/d = 1$ and 5). See Supplemental Movie for side-by-side video of these simulations. The flowing simulation (with a smaller grain diameter) finishes slightly earlier, as removed material points do not require much computational work to handle. As with the coarse silos, we use substepping with $\Delta t_{{substep}} = \Delta t / 150$.

We compared the velocity profile in this refined case with the coarse cases run previously with the results displayed in figure 9. Both the static configuration at $W/d = 1$ (snapshots taken at $350 \Delta t_f$, where $\Delta t_f = 1/240$ s as before) and the flowing configuration at $W/d = 5$ (snapshots taken at $105 \Delta t_f$) show good agreement when refined. The increased resolution allows the arch to form with a slightly more natural circular geometry in the static case. In the flowing case, the refined simulation has a slightly higher velocity throughout the bulk, however, this is expected and is also due to geometric refinement; as we measure the orifice size between the fixed points, decreasing the element size allows a larger fraction of material to leave without sensing the orifice edge. This is a minor effect, but does serve to increase the flow rate a little as observed. We expect subsequent refinements to show increasingly smaller differences in velocity fields since the affected elements will get smaller.

Figure 9.

Speed contours of coarse simulations (a,b) and refined simulations (c,d). The $W/d = 1$ cases are taken at time $350 \Delta t_f$ and the $W/d = 5$ cases at $105 \Delta t_f$ (recall $\Delta t_f = 1/240\,\textrm {s}$). The refinement allows a more accurate representation of the orifice size across all simulations and the arch formed in the $W/d = 1$ case. Velocities are plotted as given on each material point.

Similarly, we compared the pressure throughout both sets of silos in figure 10. As before, good agreement is shown between the refined and coarse cases. The static case shows a few differences due to the more circular geometry of the arch, however, the key features of small-but-non-zero pressure above the arch and a largely lithostatic pressure distribution away from the orifice remain. In the flowing case, the refined simulation has a noticeably smoother pressure field and resolves the free surface with more fidelity as expected, but otherwise qualitatively looks similar to the coarse case. We can clearly see the free-fall region above the orifice in both cases, where the pressure of the material points goes to zero as they enter the disconnected state.

Figure 10.

Pressure distributions of coarse simulations (a,b) and refined simulations (c,d). All plots are taken over a three-frame average centred at $350 \Delta t_f$. As with the previous pressure plots, the technique presented in Andersen & Andersen (Reference Andersen and Andersen2009) for visualization is used.

To understand the shape of the arch, we turn our attention to the $g$ field, as $g$ controls how much plastic flow occurs in the dense phase. Since entire elements contribute to the $g$ field and are not considered at a sub-grid scale, this in turn restricts the shapes that the boundary can take; the effect is most apparent in simulations such as a clogged silo since the arch formed must conform to the background mesh geometry. Sufficient mismatch between the true solution and the one forced by the grid geometry might result in flow where none was expected or vice versa. We only observed this in boundary cases (e.g. in the silo, ones in which we were unsure were flowing or static), but the underlying mechanism may persist in pathological grid topologies. We also note that computing the solution is noticeably more expensive in time than a local version of the model – the extra operations performed do contribute to this, but the increase in number of synchronization points also should not be underestimated. Naively coding the algorithm presented can easily result in a method which is tens to hundreds of times slower than the local-only $\mu (I)$ model; while none of our implementations have been heavily optimized, many of the simulations in this work took approximately two to five times longer than our implementation in Dunatunga & Kamrin (Reference Dunatunga and Kamrin2015) of a local-only model implementing the $\mu (I)$ rheology from Jop et al. (Reference Jop, Forterre and Pouliquen2006) (caching the results where possible during substepping results in significant time savings). It may be possible to mitigate much of the computational expense by switching to an implicit solver, but this must be done carefully due to the nature of the material tangents. We would like to explore these solvers in future work as they would improve the applicability of the method to engineering problems.