DEM-Simulation of thin elastic membranes interacting with a granulate

The discrete element method (DEM) solves Newton’s equation for the position $\boldsymbol{r}_{i}$ and the angular orientation $\boldsymbol{\varphi}_{i}$ of each particle, $i$, of mass $m_{i}$ and tensorial moment of inertia, $\hat{J}_{i}$:

Here, $\boldsymbol{F}_{i}^{\text{ext}}$ is an external force, e.g. gravity, and $\boldsymbol{F}_{ij}$ and $\boldsymbol{M}_{ij}$ are the force and torque acting on particle $i$ due to contacts with particles $j$. There are several models for $\boldsymbol{F}_{ij}$ and $\boldsymbol{M}_{ij}$ as functions of the relative position, orientation, velocity and angular velocity of the involved particles, $i$ and $j$, for an extended discussion see, e.g., shaferForceSchemesSimulations1996 ; poschelComputationalGranularDynamics2005 ; kruggel-emdenReviewExtensionNormal2007 ; kruggel-emdenStudyTangentialForce2008 ; matuttisUnderstandingDiscreteElement2014 .

In the current paper, we focus on the description of an ambient membrane and its interaction with the granular particles. In our model, the elastically deformable membrane is modeled by mass-carrying particles that are connected by viscoelastic springs (mass-spring system, MSS). The topology of the membrane is given by a mathematical graph whose vertices and edges are represented by particles and springs, respectively.

The positions of the membrane particles (here called vertex particles), thus, describe the shape of the membrane (Fig. 1).

They are subject of Newton’s equations where the forces acting on the vertex membrane particles originate from three contributions: (a) viscoelastic stretching of the membrane, (b) moments due to bending of the membrane, and (c) interaction of granular particles with the membrane. We shall discuss these contributions in Secs. 2.3-2.5. The total force acting on the vertex particles is the sum of these three contributions.

Given two adjacent vertex particles $i,j$ at positions $\boldsymbol{\rho}_{i}$, $\boldsymbol{\rho}_{j}$ and velocities $\dot{\boldsymbol{\rho}}_{i}$, $\dot{\boldsymbol{\rho}}_{j}$, we define the relative quantities

and the unit vector

The interaction between particles $i$ and $j$ is due to linear elastic spring. Particle $i$ feels the force

which is the force of a damped harmonic oscillator, with equilibrium length $\rho_{ij}^{0}$, effective mass $m_{ij}^{\text{eff}}=m_{i}m_{j}/(m_{i}+m_{j})$, damping coefficient $\gamma$, and spring constant $k$.

To relate the spring constant, $k$, to material characteristics, we notice that each realistic membrane has a finite width, $d$, and the elasticity of the membrane material is characterized by its elastic modulus $E$. For our idealized two-dimensional membrane of vanishing thickness, one obtains the elastic constant kotMassSpringModels2017 ; ostoja-starzewskiLatticeModelsMicromechanics2002 .

To explain the description of membrane flexibility, we consider four adjacent vertex particles at positions $\boldsymbol{\rho}_{1}$, $\boldsymbol{\rho}_{2}$, $\boldsymbol{\rho}_{3}$, $\boldsymbol{\rho}_{4}$ bridsonSimulationClothingFolds2003 , see Fig. 2.

The vertex particles span two triangles with normal vectors

The corresponding angle $\theta_{12}$

characterizes the flection of the triangles, with respect to their common edge $\boldsymbol{\rho}_{43}=\boldsymbol{\rho}_{4}-\boldsymbol{\rho}_{3}$. The restoring torque counteracting the flection can be expressed by elastic and dissipative forces, $\boldsymbol{F}_{i}^{\text{el}}$, acting on the involved vertex particles, $i\in\{1,2,3,4\}$ bridsonSimulationClothingFolds2003 ,

where $k^{\text{el}}$ and $k^{\text{diss}}$ are material parameters.

The directions of the forces are linear combinations of the triangles’ normal vectors given by

Each vertex particle is involved in 6 different pairs of triangles, see Fig. 1. The total force acting on a vertex particle is, thus, the sum of the 6 corresponding forces given by Eqs. (11, 12).

For the description of the interaction between the membrane and the confined granular particles, we assume triangular patches spanned between the time dependent momentary positions of adjacent vertex particles, Fig. 3.

The interaction of the granular particles with the membrane is then described by contacts between the granular particles and the patches. This assures that the patches are always impenetrable disregarding of the sizes of the particles and the deformation of the membrane.

Contacts between a patch and a granular particle are classified as vertex, edge, or face contact:

In case a granular particle contacts two neighboring patches, $A$ and $B$, we chose the adequate contacts according to Tab. 1 in dependence whether these contacts are face, edge, or vertex contacts huNewAlgorithmContact2013

These selection rules do only apply if the patches $A$ and $B$ have a common edge. Otherwise, all contacts are handles regularly.

In case a granular particle contacts three neighboring patches at their common vertex, one of these contacts is selected randomly and the others are disregarded.

Contacts of the membrane with itself may be calculated from contacts between vertex particles and triangular patches.

Once the contact point is defined, we compute the force according to the specified contact law. The relative velocity of the granular particle and the membrane at the contact point which enters the force is interpolated from the velocities of the vertex particles using barycentric weights, as sketched in Fig. 3. Similarly, the obtained force is distributed to the involved vertex particles with barycentric weights $a_{1}$, $a_{2}$, $a_{3}$. The positions and velocities of the vertex particles and, thus, the dynamics of the membrane are obtained by numerical integration in the same way as the granular particles.

The above selection rule above in combination with the barycentric partition of the force leads to smooth and physically plausible forces acting on the vertex particles. huNewAlgorithmContact2013

The triaxial test is commonly used to investigate the mechanical properties of a granular sample. To that aim, the sample is placed between two parallel platens and wrapped by a cylindrical membrane, see Fig. 4.

A confining stress $\sigma_{\text{c}}$ in radial direction is applied through the membrane. By controlling the initial distance $h=h_{0}$ between the platens, an initial stress $\sigma_{\text{p}}=\sigma_{\text{c}}$ is applied in axial direction. Then, the platens are displaced at constant relative velocity $\boldsymbol{v_{p}}$ to apply a strain $\epsilon=\log(h_{0}/h)$. We record the corresponding deviatoric stress $\sigma_{\text{d}}=\sigma_{\text{p}}-\sigma_{\text{c}}$.

In the simulation, we represent the platens by rigid walls and the membrane by an MSS. We place small particles at random positions inside the membrane such that they do not contact one another. To generate the initial conditions, we apply the Lubachevsky–Stillinger algorithm lubachevskyGeometricPropertiesRandom1990 . We then gradually apply an initial stress, $\sigma_{p}$, to the platens and a confining stress, $\sigma_{\text{c}}$, to the membrane. To this end, each triangular patch of area $A^{w}_{i}$ and normal unit vector $\hat{\boldsymbol{n}}^{\,w}_{i}$ is loaded with a force

where $\hat{\boldsymbol{n}}^{\,w}_{i}$ is defined such that force acts from outside the membrane to the granulate located inside. After this initialization, we displace the platens at relative velocity $\boldsymbol{v}_{\text{p}}$ and record the deviatoric stress, $\sigma_{\text{d}}$. Figure 4 shows the initial and final states of such a simulation.

To demonstrate the performance of our model, we perform simulations for two different cases:

1. 1.

   We describe the membrane by a MSS where the confining stress of the membrane is provided by contacts between vertex particles and granular particles. This approach was used before, e.g., in quDiscreteElementModelling2019 . We term this setup $\text{MSS}_{\text{particle}}$.

2. 2.

   We describe the membrane as described in Sec. 2. Here the confining stress of the membrane is provided by contacts between the granular particles and the triangular patches. We term this setup $\text{MSS}_{\text{patch}}$.

For $\text{MSS}_{\text{particle}}$, we place the vertex particles of the membrane spaced by $s=1.33\,\text{mm}$, such that the radius of a vertex particle is about $1/3$ the radius of a granular particle. This value is a trade-off between keeping the computational cost low and having a smooth membrane that prevents penetration of granular particles debonoDiscreteElementModelling2012 ; quDiscreteElementModelling2019 . For $\text{MSS}_{\text{patch}}$, we use $s=4.0\,\text{mm}$ because the membrane has a smooth surface and is impenetrable by design.

We perform the triaxial test at velocity $v_{\text{p}}=0.05\,\text{m}/\text{s}$ and confining pressure $\sigma_{\text{c}}=100\,\text{kPa}$. The material parameters are given in Tab. 2. Furthermore, we choose the parameters $\gamma\approx 0.15$, $k^{\text{el}}=10^{-3}\,\text{N/m}$ and $k^{\text{diss}}=10^{-4}\,\text{Ns/m}$. Using the initial areas of the triangles, we use Eqs. (11-16) to compute the confining forces.

We consider two different systems: case 1 – a cylinder of initial height $100\,\text{mm}$ and radius $25\,\text{mm}$, and case 2 - a cylinder of initial height $140\,\text{mm}$ and radius $35\,\text{mm}$. The numbers of particles used for the membrane and the granulate are given in Tab. 3.

Figure 5

shows the deviatoric stress, $\sigma_{\text{d}}$, as a function of the axial strain, $\epsilon$. The different membrane representations, $\text{MSS}_{\text{particle}}$ and $\text{MSS}_{\text{patch}}$, do not lead to significant differences of the stress-strain behavior.

Table 4 compares the computer time used to simulate a real time of $20\,\text{ms}$.

The reduced number of particles in $\text{MSS}_{\text{patch}}$ accelerates the simulations by about the factor 5, compared to simulations using $\text{MSS}_{\text{particle}}$.

A correct representation of frictional forces at contacts with the membrane is important for many applications. For instance, one of the mechanisms allowing a granular gripper to grasp an object relies on frictional forces brownUniversalRoboticGripper2010 . By construction, an object gliding on a membrane modeled by particles cannot feel a constant (or, at least, smnooth) force.

We demonstrate the smoothness of a membrane modeled by the here described model, and the resulting frictional forces by means of a simple sliding test: In $\text{MSS}_{\text{particle}}$, a membrane is modeled by vertex particles with spacing $1.33\,\text{mm}$. In $\text{MSS}_{\text{patch}}$, the membrane is modeled by patches of side length $1.33\,\text{mm}$.

For the test, we place a spherical particle of radius $2\,\text{mm}$ on a membrane, the free motion of this particle is restricted to the vertical coordinate, perpendicular to the membrane. Its horizontal motion at constant velocity, $0.05\,\text{m/s}$, is enforced externally. Figure 6 shows the particle’s vertical position and the frictional force in sliding direction.

For $\text{MSS}_{\text{particle}}$ we see oscillations in both horizontal force and vertical position. For $\text{MSS}_{\text{patch}}$, we instead observe constant vertical position and constant friction force, as expected when a particle slides on an even plane.