Discrete element simulations of micro-penetration in cohesive granular materials

(1)

Master thesis

Discrete element simulations of

micro-penetration in cohesive granular

materials

Author:

Mirte van der Eyden

1,2

Student number:

10358307

Institute:

1

_{WSL Institute for Snow and}

Avalanche Research SLF

2

_{University of Amsterdam}

Daily supervisor:

Dr. Henning L¨

owe

1

Supervisor:

Dr. Edan Lerner

2

Second corrector:

Dr. Corentin Coulais

2

August 19, 2017

(2)

Abstract

The SnowMicroPenetrometer (SMP) is an instrument that was developed to characterize snow stratigraphy for various applications. It penetrates the snow and measures the pene-tration force signal with high spatial resolution. While some structural parameters can be estimated from the measured force signal within a simple, stochastic penetration model, their relation to the true microstructural and micromechanical properties of snow is yet poorly understood. To improve the understanding of the penetration process we have analyzed Dis-crete Element Simulations of a small cone penetrating into cohesive sphere assemblies. We have used initial sphere configurations covering a variety of packing fractions and (cohesive) coordination numbers and conducted simulations for different (cohesive) bond strengths and cone sizes. Our main results show that the mean penetration force is primarily controlled by the bond failure rate, with a super-linear dependence indicating a feedback of failed material. A secondary influence stems from particle jamming in front of the cone. These influences are incorporated in a proposal towards a new model that describes the penetration force in terms of simulation parameters and failure processes around the cone.

Images front page: SLF (http://www.slf.ch/ueber/organisation/schnee_permafrost/schneephysik/ SnowMicroPen)

(3)

Acknowledgements

I would first like to thank my daily supervisor Henning L¨owe of the WSL Institute for Snow and Avalanche Research SLF. Discussions with him always gave me new insights and motivation. He made a lot of time for me and was very concerned with the project. While writing this final thesis he helped extensively and his feedback will be the basis for all the

work I will write in the future.

I would also like to thank Johan Gaume, of the EPFL, École polytechnique fédérale de Lausanne, for doing a big part of the simulations for this project and helping me to use the

simulation software myself. He was always available to think along about the results and discuss the possible interpretations. In addition, I want to thank him for organising the SLAM3 _{(SLab Avalanche Multiscale Mechanical Modelling) workshop during my stay at} the SLF institute. This was a great oppurtunity to get a broad view of the research in this

field and to present our work.

Finally, I would like to thank Quirine Krol for getting me in contact with the SLF institute and with Henning. It was a lot of fun to share the office (and ski) with you.

(4)

1 Introduction

For various practical reasons such as snow stability and avalanche formation [30], but also more fundamental questions about snow metamorphism, people are interested in snow stratigraphy and microstructure. One way of measuring the stratigraphy and microstructure is by digging a snow pit and identify the layers. Experienced snow-researchers can classify the different crystal types acoording to an international system for snow classification [6]. They classify and measure the size of the snow grains and measure the density at centimeters resolution in the snowpack. With this information, a rough picture of the stratigraphy can be obtained. This gives a lot of information about the stability of the snowpack. However, this method is very time-consuming and subjective [23]. To overcome these difficulties and make objective, fast measurements in the field, various field methods and instruments were developed in the last two decades (cf. references in [24]), among which the SnowMicroPen [29] is one of them.

1.1 SnowMicroPen

The SMP (Fig. 1) is an instrument (cf. Fig. 1, left) that penetrates the snow with a constant velocity and measures the penetration resistance. It consists of a motor that drives a long rod with a conical tip of 0.5 cm diameter into the snow and a sensor measures the force on the cone every 4 µm. A typical example where the measured data is aligned with the snow microstructure is shown in Fig. 1, as well as a micro-CT image of the conical SMP tip penetrating snow. The advantage of the SMP is that measurements are quick and have high spatial resolution. The instrument is therefore used in various studies such as the recen extensive time series of the snow stratigraphy for a complete season by Calonne et al. [5] or the snow hardness measurements [32]. The force signal gives information about the hardness of the snow. The hardness differs in various layers in the snowpack, so the SMP is helpful to distinguish these layers. The difficulty however, lies in the interpretation of the force signal and to find the relation to key snow microstructural parameters such as snow density and grain size[15, 14].

1.2 Signal Interpretation

There are various models to extract micromechanical properties of the snowlayers from the force signal [15, 22, 28]. The Poisson shot noise model [21] combines previous approaches and describes the force signal as a superposition of linear-elastic, brittle breaking of randomly distributed, uncorrelated elements. Due to these assumptions, this theory allows certain structural parameters of snow to be calculated analytically. Based on this model, a statistical model to extract snow density and grain size from the SMP measurement data is developed by Proksch et al [24]. This model matches structural parameters computed from the micro-CT images of snow with the SMP force signal of the same samples. A regression model can

(7)

Figure 1: (left) SnowMicroPen (SMP). (middle) Example of SMP data next to the snow microstruc-ture. Layers of different density can be observed in both the penetration force signal and the snow microstructure from CT. (right) CT image of the conical SMP tip penetrating snow. The SMP measures the force on the red tip during penetration.

be developed to extract the density and grain size from the penetration force data. The analysis is based on the shot noise model [21] that has some limitations. A first limitation of this theory is the lack of spatial correlation between the elements. Snow is an anisotropic material [20] because of its layering and this is not implemented in the theory. Furthermore, snow grains are expected to compactify under the cone [17]. Therefore the material that is measured is not the same as the initial snow. Bonds between the snow grains are expected to break during penetration, but the force is, in addition, influenced by the friction with the remaining snow grains in the cone region. The shot noise model does not include any interaction bewtween elements. Finally, it is unclear how the abstract ’elements’ can be related to actual snow microstructural entities.

As a remedy, by viewing snow as a foam of ice, as proposed by Kirchner in [16], we can identify the abstract elements with beams of a cellular material. Existing theories on the penetration of cellular foams [9] can be used to describe the penetration force in terms of the material properties of ice. In [14], the penetration of a foam is combined with the idea of local compaction of broken elements under the cone. This theory is however not applicable to snow yet because it uses many unknown geometrical parameters. All theories mentioned above do not give a satisfying understanding of the relation between the SMP force signal and the snow microstructure yet. It is therefore appealing to use discrete element methods to investigate further details of the penetration process.

1.3 Discrete Element Modelling

Discrete Element Modelling (DEM) is an important tool to compute properties of cohesive granular materials [25]. There are various attempts to apply DEM to snow microstructure. The difficulty lies in the representation of the complex microstructure of snow in terms of discrete (clumps of) elements. The first attempt to model deformation mechanisms in snow

(8)

using DEM was by Johnson et al. [13] who employs structures made of generic hard-core objects. More recently Hagenmuller et al. [11] suggested a faithful deterministic reconstruc-tion of snow made of grains. In the present study a different route is pursued and the initial states for penetration are created with Baxter’s method of sticky hard spheres (SHS) [3] because this method allows for a stochastic reconstruction of snow [19]. The inital states are the same as the ones used in [18] for uniaxal compression tests. The mechanical response from these tests can be related in a definite way to the system parameters which serves as a reference for more complicated cone penetration simulations. The initial states cover a variety of packing fractions and contact densities. The spheres are equipped with a cohesive contact law and penetrated with a cone while the penetration force is measured to generate a signal similar to the SMP force signal. DEM enable us to investigate the relation between statistical properties of the force on the cone and the initial parameters, the breaking of bonds and the formation of new contacts between particles in the system.

1.4 Goal and outline

The goal of this study is to make a step forward in relating the SMP force signal and mi-crostructure parameters by analyzing DEM simulations of penetration in cohesive granular materials. Starting from the ideas of the shot noise model, we investigate the relation be-tween the penetration force and the bond failure rate in the system. We will do sensitivity studies on various simulation parameters to gain information about their influence on the force. By varying the size of the penetration cone, we specifically investigate indentation size effects, as observed when penetrating foams [9]. As a next step we look for signatures of local compaction in the region under the cone. The goal is to combine the results from these tests and suggest key ingredients for a model that describes the statistical properties of the force signal for all simulations.

The thesis is organized as follows. In Sec. 2 the current theoretical models of snow micro-penetration, namely the Poisson shot noise model [21], the statistical model [24] and a micromechanical theory [14], will be explained. The latter theory is based on the penetration of cellular materials [9], of which a summary is given as well. Sec. 3 will summarize the DEM method and the sticky hard spheres and descibes the performed simulations and data analysis. Sec. 4 shows the results from the simulations, that gave the ideas for a starting point of a new theoretical model. The description of this model can be found in Sec. 5. Finally we discuss the results, limitations and ideas for future research in the discussion, Sec. 6, and conclude in Sec. 7.

(9)

2 Theoretical background

2.1 Poisson shot noise model

The interpretation of SMP data has been a challenge over the last two decades. Based on previous work [22, 28], the current theoretical model of deriving the micromechanical parameters of snow from the SMP data is the Poisson shot noise model [21]. The goal of the theory is to obtain structural parameters of snow directly from the statistical properties of the SMP force signal. We will summarize this model by stating the assumptions and the most important results.

2.1.1 Model description

The main underlying assumption of the Poisson shot noise model is that the force signal is a superposition of spatially uncorrelated, elastic-brittle failures of elements (Fig. 2). There are three parameters that can be derived from the signal. ` is the average distance between two elements that can be interpreted as the typical size of one element, δ is the deflection before the element fails and f the rupture force.

A penetrometer with a cone of projected area A that moves over a distance H would pene-trate a volume V = AH (Fig. 3). This volume consists of elements of size ` so the average number of elements in this volume is ¯N = AH/`3_{. The second assumption of the theory is} that the elements in this volume are distributed without spatial correlations. This random distribution is obtained by drawing the positions from a three dimensional, homogeneous Poisson point prosses (Fig. 3. The number density or intensity λ3d = ¯N /V characterizes the point process uniquely.

Figure 2: Parameters in the elastic brittle force behaviour: ` is the structural element size, δ is the deflection before rupture and f the rupture force. Image taken from [21].

Figure 3: The cone with projected area A penetrates a volume with a elements that are distributed randomly by draw-ing from a Poisson point process. Image taken from [21].

(10)

The third assumption is that each element contributes to the response force just by its deflection in the z-direction. The force characteristic of an element follows an ideal elastic-brittle behavior, that can be described by:

F (z) =        f δz, 0 < z < δ 0, else. (1)

The total force signal is a superposition of these responses and therefore only depends on the z-component of the positions of the elements in V . Therefore we can reduce the problem to a one-dimensional Poisson point process with intensity λ = ¯N /H = λ3dA. Since the structural element size ` is defined by ¯N `3 _{= AH, the expression of the one-dimensional} intensity is given by

λ = A

`3. (2)

The total force is then expressed as a superposition of N single element contributions F (z) that behave like Eq. (1):

FT(z) = N X n=1

F (z − zn), (3)

where the positions zn are drawn from a 1-d point process with intensity λ. A stochastic process of form (3) is referred to as shot noise process [21]. Due to the absence of spatial correlations many properties of (3) can be computed analytically.

2.1.2 Force statistics

The basic model parameters `, δ and f can be computed analytically by evaluating the cumulants κn of the force signal. We assume that the rupture force f is the same for all elements, with value f0. The first and second cumulant are the mean SMP force κ1= FT(z) and the variance κ2= (FT(z) − κ1)2. When the rupture force is fluctuating, the n−th order cumulants κn can be expressed in terms of the moments of the rupture force fn. Since we work with a uniform rupture force, fn _{= ¯}_f

0 n and κn = λδ ¯f0 n n + 1. (4)

From this we see that the cumulants all depend on the product λδ so we will have to look at higher order correlation functions to be able to compute them seperately. The force covariance is defined by:

(11)

This can be computed exactly to get the correlation function of the shot noise model: C(r) = λδf 2 3 1 − 3 2 |r| δ +1 2 |r| δ 3! θ(δ − |r|), (6)

where θ(r) is the Heaviside step function, defined by θ(r) = 1 for r > 0 and θ(r) = 0 otherwise.

The structural parameters can now be estimated from equations (4) and (6). First, δ can be inferred from the slope of the correlation function in the origin:

δ = −3C(0)

2C0₍₀₎. (7)

We can solve the first two equations in (4) by inserting hfn_{i = f}n

0, yielding: λδ =4κ 2 1 3κ2 (8) f0= 3κ2 2κ1 . (9)

Combining this with Eq. (7) allows to eliminate δ.

2.1.3 Practical aspects of parameter estimation

The shot noise analysis is based on a spatially uniform intensity of the point process λ and on a uniform rupture strength f0. Since snow is spatially inhomogenious due to layering, this analysis has to be applied using a moving window on the force data. For every window the cumulants and correlation function can be computed from the force signal in the window to get the distribution of the structural parameters. For a full practical recipe and application to a simulated point process and snow, see [21].

2.1.4 Limitations

The main interesting parameters of a snowpack are the density and grain size in different layers. It is clear that this model is not suited to obtain these parameters, since they do not explicitly enter the model. The model is rather based on elements as abstract entities which do not have an immediate interpretation in terms of actual microstructural properties. An improved model should explicitly include density or packing fraction φ. For low density, failed grains may move immediately out of the zone of influence of the cone. For higher density however, the failed particles create a compaction zone around the cone which contributes to the penetration force [17]. To account for this ’jamming effect’, the model must explicitly contain φ. The grain size must enter the model too, because the size of snow grains ranges from micrometers to millimeters. Indentation size effects, due to the variation in the ratio of cone size and grain size, dicsussed in detail in Sec. 2.3, may have influence on the penetration force.

(12)

2.2 Statistical model for density and grain size

Despite the mentioned limitations, the shot noise model still provides means to infer some structural parameters from the force signal. Proksch et all. [24] developed a statistical method that is based on the shot noise theory to extract density, correlation length and specific surface area (SSA) from the SMP force signal . This method is used extensively in practice, because it allows to gain the main information about the snow without digging extensive snow pits.

This method is based on the comparison between the micro-CT images and SMP force signals of the same snow samples. The micro-CT gives very detailed information of the 3-d snow microstructure, as shown in Fig. 4 for different kinds of seasonal snow. From this image the density ρ , SSA and correlation length `ex can be computed and related to the SMP force signal on a statistical basis.

Using regression, an estimate for the different parameters can be found in terms of the average force ¯F and the structural element size L, which is derived from the signal using the shot noise model.

ρ(smp)= a1+ a2ln( ¯F ) + a3ln( ¯F )L + a4L (10) The pcoefficients ai were computed by least square methods. Similar models can be written down for the SSA and correlation length. An example of a density profile computed from SMP measurements for both Arctic and Alpine snow is shown in Fig. 5. Here they are compared to density computed from the micro-CT image made in the SLF cold lab. The SMP density shows reasonable agreement, but is not able to predict the detailed peaks measured with micro-CT. Although this method is based on the shot noise model that has some clear limitations, it is performing reasonably well, indicating some significance of the element properties implemented in the shot noise picture. To improve it, however, a more physical description of the penetration process is needed.

2.3 Snow as a foam of ice

Kirchner proposed to describe snow as a cellular material, a foam of ice [16]. Since the density of ice (ρ = 917kg/m3) and its strength and elastic properties are better constrained than those of snow, this could be a way to gain information about the same mechanical properties for snow, that has densities between 50 kg/m3and 800 kg/m3. Gibson and Ashby (1988) have extensively described the mechanical and structural properties of cellular solids in their textbook [9]. Relevant for the present study are the indentation tests they performed on brittle foams, described in chapter 5.3(f) of [9]. In this section we will describe the results of these tests, summarize the theoretical description and compare it to the shot noise model.

(13)

Figure 4: (top) SMP force signals for typical snow types, compared to the micro-CT 3-d reconstruction of the snow microstucture (bottom). The corresponding element size ` (in the figure denoted as L), average penetration force F , density ρ, specific surface area SSA and correlation length lex are given. (a) shows low-density new snow, (b) shows a dense refrozen structure with a higher mean penetration force and more ”noise’, (c) shows medium density depth hoar with an SMP signal that shows high peaks when touching a big grain and drops to zero in the pore space between the crystals. Images taken from [24].

(14)

Figure 5: The density profile for (left) Antarctic snow and (right) Alpine snow. Both snow blocks are measured in the micro-CT and with the SMP. The results of the density computed with (10) is shown together with the density measured from the micro-CT image. In the left figure, two SMP measurements (blue and red), with mean value (green) are compared to the CT density (black). In the right figure, the density measured with a 100 cm3 _{density cutter (red) is compared to the one measured with micro-CT} (black), SMP (green). Images taken from [24]

(15)

(a)

(b)

(c)

Figure 6: (a) Results for the dependence of indentation stress σion indenter area A for a brittle cellular foam made of zirconia [9]. The penetration size effect is observed, σi is proportional to A1/2 (Data from [1], figure taken from [9]). (b) An indenter in contact with the cell edges of a cellular material. The material is idealized as consisting of rods of length ` and thickness t. Figure taken from [9]. (c) Schematic representation of a beam bending over a distance δ under an applied force f0. The beam has length ` and thickness t.

2.3.1 Size effects

Penetration/indentation tests on brittle foams provide a force signal that looks similar ap-parently to the SMP force signal. Tests of indentation with varying indenter sizes showed that the indentation stress dependeds on the area A of the indenter as shown in Fig. 6. Such a behavior is called a size effect and it can be seen that the indentation stress σ := Fm/A is inversely proportional to A1/2. Fm is defined as the mean of the highest 10 maxima in the force signal. Indentation size effects have also been investigated for snow in [12].

2.3.2 Theoretical description

The observed size effect can be understood as follows. An open-cell foam can be idealized as a network of rods with length ` and thickness t (Fig. 6). The rods are made of a linear-elastic material with Young’s modulus Es which fails when the tensile stress locally exceeds the strength σfs. In the process of indentation, rods deform and break. Broken fragments are assumed to drop and disappear from the contact area. Note that this assumption is of great interest for our application to snow penetration. If the fragments would not clear out of the

(16)

way, the structure and density of the foam would locally change during indentation and this will result in changes in the force signal.

When we consider the behavior of a single rod being in contact with the indentor, it will first elastically deflect over a certain distance δ and then break. Using beam theory (cf. Fig. 6.c) the rupture force and deflection can be expressed in terms of the geometrical and material properties of the the rod [9]:

f = σfst3/6` (11)

δ = σfs`2/Est

In a brittle material the ratio σfs/Es is small, resulting in small ratio values for the δ/`. As a consequence the probability of a cell edge being in contact with the indentor is:

P = δ

`. (12)

and the number of cells within A is given by

N = A

`2. (13)

Then, over the total area A the probability of a given number of contacts is given by a binomial distribution with the mean number of active contacts N P and standard deviation N P (1 − P )1/2. If P is small, the standard deviation can be approximated by (N P )1/2_. An estimation of this maximum number of active contacts is

Nmax= (mean) + m(standard deviation) = N P + m(N P )1/2, (14) where m is of the order of 3. The maximum force is then the maximum number of active contacts multiplied by the force at each contact, given by Eq. (12).

Fm= f Nmax= f

N P + 3√N P (15)

This leads to an indentation pressure of:

σi= Fm A ∼ σfs t ` 3δ ` " 1 + ml 3 Aδ 1/2 # . (16)

(17)

2.3.3 Comparison with shot noise model

This model of foam penetration can be compared to the poisson shot noise model that also consists of the elastic-brittle breaking of elements. When we consider Eq. (16) we can write this in terms of the parameters of shot noise by identifying δ → δ, ` → ` and noting that f = σf st3/6` → f0. Using equations (2), (12) and λ = λ3DA the expression for indentation pressure yields σi∼ f0δλ3d " 1 + m 1 Aδλ3d 1/2 # . (17)

In the shot noise model, we can use Eq. (4) to find an expression for the mean force and the variance: κ1= ¯F = 1 2λf0δ (18) κ2= (standard deviation)2= 1 3λf 2 0δ. (19)

Applying the same assumption that the maximum pressure is given by (mean)+m(standard deviation) we conclude that the maximum indentation pressure in the shot noise model can estimated

by: σi= κ1+ m √ κ2= 1 2f0δλ3d " 1 + 2m 1 Aδλ3d 1/2 # . (20)

This is, up to some factors, the equivalent of the expression for foams (20) which shows the conceptual similarity of the shot noise model to foam indentation modeling.

2.4 Other theoretical approaches

In 2003 Johnson developed a statistical, micromechanical theory of cone penetration in gran-ular materials [14]. The aim was to write an accurate physically based theory of penetration that allows for direct calculation of the material properties. Since then there were mostly statistical attempts to compare measurements to laboratory data, like the Proksch method described in Sec. 2.2. The work from [14] contains interesting aspects relevant for the present thesis which are worth to recap. The theory has not been applied yet to obtain parameters of interest for snow, but some ideas can be used towards the development of a new SMP model. In the following we will briefly describe the theory and point out ideas that are of interest for the present study.

2.4.1 Model description

The model is again based on the elastic-brittle failure of rods in a geometry similar to the one shown in Fig. 6. In this case however the indentor is conical, like the SMP. The most

(18)

Figure 7: Schematic picture of the geometry of the penetration area for the model by Johnson [14]. The cone penetrates a cellular material and forms a compaction zone of broken rods that pile up around the cone. Together with the cone they form the Penetrometer effective surface (PES). Figure taken from [14].

important difference to Sec. 2.3 is that a ’clearing effect’ as in the theory of Ashby and Gibson, where broken elements are not accounted for anymore, is dropped. In this theory, fragments from failed elements form a compaction zone around the cone as observed in [17]. Fragments pile up in this zone until they lock up due to a certain maximal density. The cone, together with this compacted area form a Penetrometer effective surface (PES). Only elements in contact with the PES can fail.

This model involves many geometrical factors, as shown in Fig. 7, but the theory starts from the same ideas as the foam penetration theory. The penetration force at any instant is given by the number of active contacts multiplied by the average resisting force of a single element (Eq. (14)). However, the number of active contacts is influenced by the local density in the compaction zone. The compaction is quantified by means of a volumetric strain that describes the current density ρ of the material compared to its inital density ρ0 . It is given by β = 1 −ρ0 ρ . (21)

When a critical density ρcr is reached, elements lock up and the compaction coefficient is:

(19)

βcr= 1 − ρ0 ρcr . (22)

The expression for the number of active contacts (14) is again used, but the probability of an element being in contact that was P = δ/l is modified to account for the . Before any initial failure of elements near the cone, P = P1, where P1 is the maximum probability of contact for elements in the inital material. After the first rupture of an element, fragments are compacting and the probability of contact is modified by the increasing density:

P = P1

β(1 − βcr) βcr(1 − β)

. (23)

This improved probability of contact P can now be used to calculate the maximum penetration force, which is given by eq (15). P is modified according to Eq. (23). The number of active contacts N that is given by Eq. (13) is modified too because the PES grows due to compaction and therefore N is inversely proportional to β. Johnson implements more improvements relative to (16) by using the geometry of a cone, but we will not describe this though. It is difficult to directly test this theory using our DEM simulations, but it gives an idea of how local compaction can be implemented in the theory.

(20)

3 Methods

3.1 Discrete Element Methods

The Discrete Element Method (DEM) is a tool to compute the motion of a large number of particles by solving Newton’s equations of motion [25]. The fundamental idea is that the simulated material consists of particles, for which the positions and velocities are computed for every timestep. Particles interact via pair-interactions which are defined in a so-called contact model. Depending on the aim of the study, different contact models can be used. In our simulations the main focus lies on the impact of initial microstructures. To this end, initial states from [18] are used that are generated with the method of sticky hard spheres. Subsequently, we use these initial particle positions with an cohesive contact law and simulate the penetration of a cone into the assembly. Both steps of this method will be described in this section.

3.2 Initial states - Sticky hard spheres

The initial sphere configurations we used are exactly the same as those used in [18] for uniax-ial compression. For the present study these packings are used to do cone penetration tests. The initial positions are created with the sticky hard sphere (SHS) model as an auxiliary tool, which is explained in the following. Baxter’s model of SHS was used extensively to study the complex dynamical behavior of colloidal suspensions with short-range interactions [3]. The model consists of pairwise interactions with a hard-core repulsion and a short range attractive well. The strength of this adhesion is controlled by a stickiness parameter τ . The interaction between spheres with diameter d can be descriped by the pair potential:

V (r) kBT =          ∞, r < d ln12τ δ/(d + δ), d < r < d + δ 0, r > d + δ. (24)

SHS are a way to create sphere configurations with identical packing fraction φ but consid-erably different microstructures. For the present study most important is that the average coordination number zc increases with decreasing τ . Thermodynamic and configurational properties of the system can be computed using the Percus-Yevick (PY) approximation [2]. This approximation allows to compute the correlation function from the Ornstein-Zernike equation that can be solved analytically for hard shperes. In the PY approximation, zc is dependent on φ and τ according to

zc(φ, τ ) = 2φλ(φ, τ ), (25)

(21)

Figure 8: Schematic of the phase diagram of the SHS in the φ, τ plane. The results of the PY approximation for the derivation of coordination number zc is included. For all the chosen simulation points shown we have 3 realizations. Image taken from [18].

φ 12λ 2₋_{τ +} φ 1 − φ λ + 1 + φ/2 (1 − φ)2 = 0. (26)

In [18] this approximation is tested for the used SHS smaples and confirmed to be relevant. Figure (8) shows a contour plot of the average coardination number for varying φ and τ . The black dots in the figure indicate the simulation points, that will be used in our study. For each point 3 realizations of SHS configurations were created by a Monte-Carlo scheme [18]. Afterwards all spheres in contact are equipped with a new pair interaction, the cohesive contact model.

3.3 Contact model

For the DEM simulations we use the software PFC3D v5 by Itasca. We used the PFC parallel bond model as a cohesive contact law that acts in parallel to the classical linear contact model [8]. This is schematically shown in Fig. 9.

For the unbonded part, the normal component of the force is given by the sum of a linear elastic term and a viscous term, called a spring-dashpot model. It can be visualised as a spring with a damper that is compressed by the amount of overlap of the particles. The shear force is also linear elastic with a Coulombian friction term, implying energy dissipation due to friction between the particles. The parameters are the contact elastic modulus Eu and Poisson’s ratio νu of the particles and the restitution coefficient eu, that sets the viscosity

(22)

and the friction coefficient µ. The value of the Euwas chosen in such a way that the normal interpenetration of the particles is small. Hence, we work in the quasi-rigid grain limit.

The bonds in the simulation can be considered as cylindrical beams between spheres that are in contact. The elastic material properties of the beam are the elastic modulus Eband Poisson’s ratio νb. The beams can bend and rotate within certain stress threshold values set by the shear and tensile strength σsand σt. The beam has radius rband area A = πrb2. The moment of inertia of the bond cross-section is I = πr4

b/4 and the polar moment of inertia J = πr4

b/2. Using this, the maximum shear and tensile stresses σs,max and σt,max can be computed using beam theory.

σs,max= − Nb A + ||Mb,1|| rb I (27) σt,max= ||Sb|| A + ||Mb,2|| rb J , (28)

where Nband Sbare the normal and shear forces in the bond, Mb,1 is the bending moment and Mb,2 the twisting moment. The bond fails if either σtor σsexceeds its maximal value. This is also illlustrated in Fig. 9. In our model, there are no new bonds formed during the simulation. Every new contact that forms will interact only according to the classical unbonded part of the linear contact model, described above.

The parameters that are used in the model and their values are given in table 1.

3.3.1 Results for uniaxial compression

For later purposes we briefly state the main results for uniaxial compression from [18]. The elastic modulus E and compressive strength σc of the cohesive sphere assembly both grow with coordination number zc. The dependence of E on zc and the initial contact density νc,0= φ0zc,0 is shown in Fig. 10.a. By fitting the data, the regression

Efit(νc) Eb

= aν_cµ, a = 1.92 · 10−3, µ = 4.91 (29)

was obtained, where Eb is the bond elastic modulus. The compressive strength of the material is shown in 10.b. Fitting this data the obtained regression is

σc,fit(νc) σt

= bν_cα, b = 4.1 · 10−3, α = 3.04, (30)

with σt the bond tensile strength.

Some sensitivity studies on the bond radius rb and the bond tensile strength σfs were conducted. The results are shown in Fig. 11. The macroscopic compressive strenght in-creases almost linearly with σfs and is proportional to rb3. These results will be relevant for comparison with our own simulation results.

(23)

Figure 9: Representation of the PFC parallel bond model used in the simulations. (a) The bonded part (black) and the unbonded part (grey) are represented with their different components. (b) The bond normal force Nb as a function of normal interpenetration/overlap of the particles δnscaled by the bond radius rb. (c) Bond shear force ||Sb|| as a function of tangential interpenetration δs scaled by rb. (d) Bond bending moment ||Mb,1|| as a function of bending rotation θ1 scaled by rb. (e) Bond twisting moment ||Mb,2|| as a funtion of twist rotation θ2scaled by rb. Bonds break when the shown threshholds are exceeded. Image taken from [18].

unbonded part (linear model)

E

u

(MPa)

ν

u

µ

u

e

u

10

0.3

0.2

0.1 bonded part

E

b

(MPa)

ν

b

σ

s

(MPa)

σ

t

(MPa)

r

b

10

0.3

1

1 0.5 r

Table 1: Parameters for the parallel bond model used in the simulations. Eu: contact elastic modulus, νu: contact Poisson’s ratio, µu contact friction, eu: contact normal restitution coefficient, Eb: bond elastic modulus, νb: bond Poisson’s ratio, σs: bond shear strength, σt: bond tensile strength, rb: bond radius, r: particle radius.

(24)

Figure 10: a) Elastic modulus E vs in initial contact density νc = φ0zc,0. The smaller inset graph shows E vs initial coordination number zc,0. b) Compressive strength σc vs νc. The smaller inset graph shows σc vs zc,0. Figure taken from [18].

Figure 11: Sensitivity studies from [18] show the dependence of the macroscopic compressive strength σc on (a) the bond tensile strength σt which has a linear dependence and (b) the bond relative radius rb/rp which has a cubic dependence. Figure taken from [18].

(25)

3.4 Cone penetration simulations

In the present study we simulate the penetration of a cone into a cohesive granular assembly of particles in a cubic box with sizes of length L0 (Fig. 12). There is no gravity in the simulations. The initial particle positions are taken from the SHS states for different values of volume fraction φ and stickiness τ . We simulate 2048 particles, which will be shown to be sufficient to avoid finite system size effects (Sec. 4.4). The radius of the particles is determined by the number of particles N , the size of the box and the volume fraction φ via:

rp = 3 4 φL3₀ N π 1/3 . (31)

The penetrator consists of a cylinder with a conical tip of 60◦_{, which is the same cone} geometry as the SMP.

Figure 12: Image of the DEM simulation of the penetration of a sphere packing in a cubic box.

Figure 13: Schematic of the penetration path of the cone penetrating the sphere system of size L0. The first and last part equal to the cone length of the data is omitted in the analysis because of boundary effects. The analysis of the data is thus limited to the length L.

Simulation parameters

We have initial packings with φ and τ for the values indicated in Fig. 8. For every combina-tion, we have three different realizations. For every realizacombina-tion, simulations are conducted with three different cone sizes γ defined by

γ := cone diameter sphere radius.

(26)

Parameter Symbol Values

packing fraction φ 0.20 − 0.25 − 0.30 − 0.35

stickiness τ 0.07 − 0.09 − 0.11 − 0.13 − 0.17 − 0.21 − 0.35

ratio cone diameter to sphere radius γ 2 − 4 − 6

Elastic modulus E (N/m2₎ Reference: 10 7

Scaled: 109

Friction co¨efficient µ 0.0 − 4.0

System size L0(m)

Reference: 1.0 Scaled: 0.01

Bond radius rb0.25 - 0.50 - 0.75 - 1.00 rp

Bond tensile strenght σb 105- 5 · 105- 106- 5 · 106- 107

Table 2: System parameters used in the simulations.

All parameters that are varied during simulations and their values are given in Table 2.

Datasets

Our dataset consists of:

189 simulations with elastic modulus 107_{comprising three realizations for all the points} in Fig. 8 with three cone sizes. We refer to it as the ”reference” set. Note that these simulations are done with exactly the same parameters as in [18].

A subset of the reference simulations, consisting of 99 simulations, is repeated with changes in the parameters. To make the properties of the system more similar to snow, E = 109_{and the size of the box is reduced to L}

0= 0.01. We refer to this set of simulations as the ”scaled” set.

Individual (φ,τ) pairs from the scaled simulations are used to test the sensitivity of the force signal on the system size N , the bond radius rb, the friction coefficient µ and the bond strength σb.

3.5 Data analysis

3.5.1 Global quantities

In each simulation, the data for the first and last part of the indentation is omitted, as shown in Fig. 13, to exclude boundary effects. We define L as the total length of the used data. Therefore the average value of some quantity ψ(z) along the penetration path is defined as

(27)

¯ ψ = 1 L Z L 0 dzψ(z). (32)

For all simulations we focused on the following quantities: f(z): the instantaneous position dependent force on the cone,

nc,u(z): the instantaneous number of unbonded contacts in the entire system, ∆nc,b(z): the number of broken bonded contacts in the entire system.

Although the simulations give access to any quantity of interest, these global contact measures provide a convenient starting point. The focus lies on relating the averages of these quantities to each other and to the initial simulation parameters.

3.5.2 Local quantities

As a next step, we also adress the problem more locally in the vicinity of the cone. For a few examples we evaluate the number of particles in contact with the cone ncone_p,c and the local packing fraction around the cone. To measure the local packing fraction in the region under the cone, we defined a measurement region in PFC. The measurement region is a sphere centered around the tip of the cone, with diameter equal to the cone diameter, shown in Fig. 14. Inside the measurement region, the porosity can be measured. The porosity p(z) is defined as the part of space that is not occupied by spheres. From this, the local packing fraction in the sphere can be computed via φ(z) = 1 − p(z).

There are three different types of densities evaluated from these measurements (Fig. 15). In the first measurement the sphere follows the same trajectory as the cone would, but without actually penetrating the sample. This gives the initial density profile of the packing φinit(z) as a function of the penetration distance (Fig. 15.a). Note that the variations in this profile are microscopic density fluctuations at a scale R in an otherwise homogenous sample with packing fraction φ.

During the second measurement the packing is penetrated with the cone while the mea-surement region moves along with the tip of the cone. This gives the local packing fraction in the sphere during penetration (Figure15.b). However, the porosity and therefore also the packing fraction that is derived in the sphere contains the volume of the cone and part of the penetration rod. This leads to an underestimation of the actual packing fraction in the region under the cone. Since this value is therefore not a physical quantity, we call it the auxiliary packing fraction φaux(z). As a third packing fraction we estimated φloc(z), which is the corrected packing fraction in the measurement sphere (Figure15.c). By comparing φauxand φloc an estimate of compaction can be given.

We perform the tests for three realizations with φ = 0.30, τ = 0.13 and τ = 0.21, γ = 6, both for the scaled value of E = 109 _{and E = 10}8_.

(28)

Figure 14: Snapshot of a PFC measurement with measurement region indicated by the large red sphere. The small spheres are colored by velocity here.

Figure 15: Schematic representation of (a) φinit(z), the initial packing fraction in the measurement sphere (b) φaux(z) the packing fraction measured in the total sphere and (c) φloc(z), the corrected packing fraction in the sphere without the volume of the cone.

(29)

4 Simulation results

In this section the results of our simulations are presented. For an overview, a number of examples of the penetration data is shown and discussed first. Then we present the results of our analysis of averaged quantities, namely size effects, the number of broken bonds, some sensitivity studies and the local density under the cone. The results are compared to predictions and ideas of various existing models and it is pointed out where improvements of the models should be made. In the next section, ideas for a new model will be presented.

4.1 Data examples

Fig. 16 shows examples of the simulation data for a selection of simulations. All figures show the position-dependent quantities force f , number of unbonded contacts nc,uand the number of broken contacts ∆nc,b(z) along the cone trajectory (cf. Sec. 3) for different parameter sets. Fig. 16.a shows data for φ = 0.25, τ = 0.15, γ = 2, in the reference set, Fig. 16.b for φ = 0.30, τ = 0.21, γ = 4, in the reference set, Fig. 16.c for φ = 0.30, τ = 0.21, γ = 2, in the scaled set and Fig. 16.d for φ = 0.30, τ = 0.09, γ = 4, in the scaled set. Some interesting features can already be observed from the raw data.

Qualitative observations

The characteristics of the force signal change considerably for different values of E. The reference dataset (Fig. 16.a,b) shows generally a few large force peaks with some apparent correlations with the number of unbonded contacts (blue) while the breaking of bonds (green) is a gradual process throughout the penetration. Fig. 16.a shows that there is an increase in the number of broken bonds when the force drops, but this is not observed for all simulations as can be seen in Fig. 16.b.

The scaled dataset (Fig. 16.c,d) shows different features. Note that the force in the scaled data is an order 104 _{lower because all lengthscales are reduced by a factor 10}2_{. A} higher elastic modulus means that the bonds are stiffer and therefore do not bend as much prior to failure. This gives rise to the noisy appearance of the force data. The number of broken bonds is higher in this case, and peaks in the force signal correspond to ’avalanches’ of breaking events. The unbonded contacts seem to be less relevant.

4.2 Size effects

The first feature we investigated is the size effect coming from the penetration of the system with three different cone-diameters (cf. Sec. 3). For foams and in the shot noise model, the dependence of the indentation pressure σ on the area A of the indentor is given by Eq. (20), while for A → ∞ the strength of the system in uniaxial compression (3.3.1) should be recovered. For uniaxial, the compression strength is only a function of the contact density νc, as described by the equation (30).

(30)

(a) (b)

(c) (d)

Figure 16: Example of data for the simulation with (a) φ = 0.30, τ = 0.09, γ = 4, in the reference set, (b) φ = 0.30, τ = 0.21, γ = 4, in the reference set, (c) φ = 0.30, τ = 0.21, γ = 2, in the scaled set, (d)φ = 0.30, τ = 0.09, γ = 4, in the scaled set.

(31)

Figure 17: Maximum indentation pressure σmax= Fmax/A vs contact density for the reference dataset (E = 107). The different colours correspond to the relative indentor area γ, that is 2 (red), 4 (green) and 6 (blue). Different markers are the different packing fractions φ (circle: 0.20, cross: 0.25, star: 0.30, square: 0.35). The macroscopic strength for the same samples (eq (30)) is given by the black line.

Figures 17 and 18 show the dependence of the indentation pressure on νc for the three different indentor sizes γ, for both values of E. Fig. 17 shows the data from the reference data set which is compared to the limit A → ∞ from the uniaxial compression, for which Eq. (30) is plotted in there. Fig. 18 shows the results for the scaled simulations. In both cases, a size effect due to different values of A is observed. Penetration with a smaller size leads to a higher pressure. Three values of γ are however not sufficient to verify the A−1/2 _{dependence, but in the limit of large A the data is consistent with the macroscopic} compression strength.

All figures in the present section will contain datapoints that correspond to certain values of φ, τ and γ. Unless specified explicitly, we follow a consistent scheme: A datapoint is always the mean over the three realizaztions that share the same parameters. Different colors will represent the different values of γ: 2 (red), 4 (green) and 6 (blue). Different markers will represent the different values of φ: 0.20 (circle), 0.25 (cross), 0.30 (star), 0.35 (square). Since the simulations for the scaled dataset were performed in a smaller box of size L0 = 0.01, the force will always be normalized by a squared lengthscale in the system. Thereby the comparison in magnitude of force can be made between simulations of the reference and the scaled dataset.

(32)

Figure 18: Maximum indentation pressure σmax = Fmax/A vs contact density for the scaled dataset (E = 109). The different colors correspond to the relative indentor area γ, that is 2 (red), 4 (green) and 6 (blue). Different markers are the different packing fractions φ.

4.3 Bond breaking

The various models described in Sec. 2 all express the average penetration force in terms of a number density of elements, or in terms of the number of elements that fail. Therefore it is natural to investigate the dependence of ¯f on the rate of bond breaking during penetration. In the shot noise model, ¯f depends linearly on the intensity of the Poisson point process λ (Eq. (19)). This is equivalent to saying that ¯f depends linearly on the ’rate’ at which bonds break per unit length during penetration since we can identify ∆nc,b/L with λ.

To this end we show a scatter plot between ¯f and ∆nc,b/L in Fig. 19.a, which shows that for all simulations in the scaled and reference set, the average force can still be explained primarily by ∆nc,b/L. However, the log-log plot reveals that the dependence is rather superlinear. In order to understand the origin of the value of ¯f and the dependence on ∆nc,b(L)/L, Eq. (19) is rewritten in terms of the simulation parameters using the beam formulas (12). The beams are the bonds between two connecting particles, so t is identified with 2rb and ` = 2rp. Then, ¯f can be expresesed as

¯ f = 4σ 2 brprb2 3E ∆nc,b L . (33)

In this equation, all parameters are known and we test its validity by plotting ¯f against the right hand side of the equation in Fig. 19.b. The 1:1 line is added as a guideline. It can be seen that for the reference data set, the order of magnitude of the prefactor is estimated

(33)

well with the ingredients of the shot noise model. The scaled simulations are however off by a factor of order 10. In addition, the nonlinear dependence on the bond failure rate is obviously not captured by the model. We will later come back to this point and suggest extensions to the model to introduce a nonlinearity.

As an emperical attempt, a nonlinear dependence on ∆nc,b/L is tested as an empirical extension of Eq. (33). When the power of ∆nc,b/L is changed from 1 to 4/3, almost all simulations collapse on the 1:1 line (Fig. 19.c). Another test is shown in Fig. 19.d wich reveals that the dependence of ¯f on the average number of broken bonds during penetration ∆¯nc,b is rather linear. Of course this property is not intensive because it will grow when the penetration distance grows.

4.4 Sensitivity tests

In order to investigate the dependence of the averaged penetration force on different simu-lation parameters, we conducted various sensitivity tests on system size, friction coefficient, bond radius and bond tensile strength. The dependence on bond radius and bond tensile strength can be compared to the corresponding results for macroscopic compression [18]. All sensitivity studies are done for the three realizations the parameters φ = 0.30 and τ = 0.13. For the study on friction coefficient and bond tensile strength only penetration with a cone of size γ = 4 is performed.

System size

As indicated by Fig. 14, the cone has a certain region of influence around it where the particles move and force chains develop. If the simulation volume is too small, boundary effects would influence the force signal. We varied the system size by varying the number of particles N from 254 to 8096. Technically, the system size remains the same but the size of the particles is adjusted to correspond to a packing fraction φ = 0.30. The mean pressure and pressure-variance for all values of N are shown in Fig. 20. We see that neither shows a trend. This indicates that simulating systems with 2048 particles is reasonable. The zone of influence around the cone is smaller than the size of the whole system.

Friction coefficient

In the system there is friction between the spheres itself, as well as between the spheres and the cone (Sec. 3.3). When this friction is increased, we expect a higher penetration force. The sensitivity analysis for the friction coefficient µ is shown in Fig. 21. For low values of µ (<1) the average force shows a linear dependence. All simulations are done for φ = 0.30, τ = 0.13 and γ = 4.

(34)

(a) (b)

(c) (d)

Figure 19: Both for the reference and scaled data: (a) average indentation force ¯f , relative to the system size, vs the average rate of bond breaking ∆nc,b/L over the total penetration. A linear line is plotted for comparison. (b) ¯f relative to bond radius vs the right hand side of Eq. (33). The 1:1 line shows the exact Eq. (33). (c) Test where the dependence on ∆nc,b(L)/L is made nonlinear. (d) ¯f versus the average number of broken bonds during penetration ∆¯nc,b.

(35)

Figure 20: Average penetration pressure and pressure variance for systems with a increasing number of particles. The simulations were performed for φ = 0.30 and τ = 0.13.

(a) (b)

Figure 21: (a) Average force ¯f and (b) maximum penetration pressure σmax vs friction coefficient µ. The simulations are done for φ = 0.30, τ = 0.13 and γ = 4.

(36)

(a) (b)

Figure 22: (a) Average force ¯f and (b) maximum penetration pressure σmax vs bond radius rb. The simulations are done for φ = 0.30, τ = 0.13 and γ = 4. The datapoints are the average over 3 realizations.

Bond radius

From the shot noise model as well as the foam theory, we have seen that the indentation pressure and the average force are proportional to f0δ. The dependence of f0δ on the bond radius is revealed by Eq. (33)

We checked the dependence of ¯f on the bond radius in the simulations by varying the bond radius but having the particle radius fixed. The results are shown in Fig. 22 and are consistent with the quadratic relation of ¯f on the bond relative radius rb/rpas predicted by eq (33). We also compare this to the sensitivity studies in [18], of which the influence on the macroscopic compression strength is shown in Fig. 11, predicting a cubic denpendence. We compare the maximum penetration pressure shown in Fig. 22 to these macroscopic results and see that they do not show the same trend. σmax has the same dependence on rb as ¯f and we observe a quadratic dependence again. The range and number of varying bond radii is however too small for a robust conclusion on an exponent.

Bond tensile strength

We varied the bond tensile strength σb of the simulations over two orders of magnitude. From Eq. (33) we expect that ¯f will show a quadratic dependence on σb. The results are shown in Fig. 23, where a quadratic and linear line are plotted for comparison. We see that ¯f seems to be linear in σb and not quadratic. As above, we compare our result for σmax in Fig. 23.b to the result for the macroscopic compression strength (Fig. 11.b). The dependence on σb is not linear like for macroscopic compression.

(37)

(a) (b)

Figure 23: (a) Average force ¯f and (b) maximum penetration pressure σmax vs bond tensile strength σb. The simulations are done for φ = 0.30, τ = 0.13 and γ = 4, for only one realization. A linear and quadratic line are plotted for comparison.

4.5 Lenthscales beam theory

It is appealing to also apply the shot noise model to the simulations of SHS. By computing the force correlation functions of the simulation data, δ and f0 are obtained using the shot noise method (cf Sec. 2.1). Inverting Eq. (12), the two lengthscales corresponding to the length ` and width t of a beam can be expressed as:

` = 6f0δ 3_E3 s σ4 b 1/5 , (34) t = δEs(6f0) 2 σ3 b 1/5 . (35)

Computing the values of these lengths for our simulation data helps to assess the interpreta-tion of the abstract ”elements” in the shot noise model. Therefore the values of t and ` are compared to the bond radius rb and the particle radius rp. Identifying the bonds between particles with the rods that break, the relation between t and rb as well as between ` and rp is shown in Fig. 26. For the reference dataset Fig. 24.a shows the computed values of t/rb for all simulations plotted versus νc. The computed values of t are in the order of five to ten times larger than rb. Similarly, l/rp is plotted versus νc as shown in Fig. 24.b. The computed length ` is in the order of ten to fourty times rp. The intensity λ of the shot noise process gives us a lengthscale `, the structural element size (2). No correlation is found between the ` from shot noise theory and the estimate from Eq. (34).

Figures 24.a and 24.b suggest a strong correlation between t and ` which is shown in Fig. 24.c.

(38)

(a) (b)

(c)

Figure 24: Computed lengthscales t (a) and ` (b) relative to respectively the bond radius rb and the particle radius rp, plotted versus the contact density νtextrmc. (c) t/rb vs `/rp shows the correlation between the two lengthscales

(39)

4.6 Local quantities

This section contains preliminary results how to proceed with the analysis in future work. In the previous section, the focus is on ’global’ quantities of the system and relations between averaged quantities or initial parameters and features of the force were established. To further assess deviations of the data (e.g. from prediction (33)) it is necessary to adress the problem more locally in the vicinity of the cone. Therefore we have started to analyze local density fluctuations around the cone to assess aspects like those discussed in Sec. 2.4.

Packing fraction under the cone

Three types of packing fractions are measured in the region under the cone (Fig. 14): φinit(z), φaux(z) and φloc(z) ( cf. Sec. 3.5.2). The tests are performed for three realizations with φ = 0.30, γ = 6, two values of τ : τ = 0.13 and τ = 0.21, both for the scaled value of E = 109_{and E = 10}8_{. We expect the packing fraction under the cone to change due to the} interaction between the cone and the spheres. The local packing fractions are compared to the smoothened force signal fsmooth(z), that is a running average over a windowsize of 5000 timesteps. The results for some tests are shown in Fig. 25. We observe fluctuations in φinit that apparently correlate to fsmooth(z). The local packing fraction φloc, which is the true density in the space around the cone, is consistently higher than φinit. φaux is consistently very similar to φinit. This means that essentially all initial particles that would have been in the volume occupied by the cone, are moved to the space underneath, such that the total amount of spheres in the measurement region is not altered by the penetration.

(40)

(a)

(b)

(c)

(d)

Figure 25: Results of local measurements for simulations with τ = 0.13 in the (a) scaled dataset where E = 109 _{and (b) where E = 10}8 _{and with τ = 0.21 in the (c) scaled dataset where E = 10}9 _{and (d)} where E = 108. The force signal (light-red) and smoothened force signal (red) are compared to the three different kinds of local packing fraction explained in Fig. 14: φinit(z) (black), φaux(z) (blue) and φloc(z) (green) (cf. Sec. 3.5.2).

(41)

5 Towards a new model

The goal is to propose a starting point for a new model that is based on the Poisson shot noise model but includes new relevant features of the penetration process. An expression for the mean force ¯f is written in terms of simulation parameters and local contacts around the cone.

5.1 Theory

The mean force is mainly determined by the rate at which bonds break during penetration (Fig. 19). Furthermore, there is a significant correlation between the number of unbonded contacts at any moment and the force signal (cf. Fig. 16). This motivates us to describe the force in terms of two populations of elements which are in contact with the cone. ncone_b is the number of bonded elements in contact with the cone and ncone

u is the number of unbonded elements in contact with the cone. In terms of these populations, we propose the following expression for the instantaneous total force on the cone:

f (z) = ncone_b (z)f0+nconeb (z) + n cone

u (z)µvsmp. (36)

The first term comes from the deflection of the bonds of all the bonded elements in contact with the cone at position z. Each of them contribute with an average force f0. This term is exactly in accordance with the Poisson shot noise model, if equipped with the right distribution of the positions of the bonded elements. The second part of Eq. (36) describes a frictional force coming from all elements in contact with the cone. During penetration, the force on the cone is now determined by the total number of particles in contact with the cone. This is an implementation of the result that the fluctuating density sets the global penetration force, shown in Sec. 4.6.

The primary goal of this new approach is to relate the penetration force to all the parameters that can be controlled in the simulations. Therefore the expression (36) is averaged to find the average force in terms of the mean number of bonded elements ncone

b and unbonded elements ncone

u in contact with the cone: ¯

f = ncone

b f0+nconeb + nconeu µvsmp. (37) We assume that the variation in the force is due to the spatial arrangement of the particles in the sample. To make the connection to the shot noise model, we assume this arrangement is a Poisson point process with intensity λ = Aλ3d. The force-displacement characteristic of each bond breaking is assumed to be an elastic-brittle responce (Eq. (1)).

Every particle that comes in contact with the cone follows a sequence of states from a bonded to an unbonded to an inactive particle. This means that ncone

u (z) consists of particles of which the bonds are broken, but remain in contact with the cone over a certain penetration distance. Therefore, we can express the rate of change of unbonded contacts as:

(42)

d dzn cone u (z) = n0fail(z) − 1 ξn cone u (38)

The first term n0_fail is the gain term from the number of bonded elements that failed in [z, z + δ]. The second term is the loss which is a certain fraction of the unbonded contacts under the cone. Unbonded contacts stay in contact with the cone over a length ξ. After having moved along with the cone for a verticle distance ξ, the unbonded element will get out of contact with the cone and becomes inactive.

In order to obtain a stationary force signal as observed, the number of unbonded elements under the cone must be stationary as well. This means that dncone

u (z)/dz = 0 in Eq. (38), and thus

ncone

u = ξn0fail. (39)

Similarly, the bonded particles stay in contact with the cone for a distance δ before they break. Therefore we can express the average number of bonded contacts under the cone in terms of the number of the average number of them that fail in [z, z + δ]:

ncone

b = n0failδ. (40)

Assuming that the bonds that break in the system are only the bonds that fail under the cone, the main failure rate n0

fail can be related to the simulation quantity:

n0_fail =∆nc,b

L (41)

where L is the total length of the simulation box and ∆nc,b the total number of broken bonds when the cone has arrived at the bottom of the box. Combining these expressions, the average force yields:

¯ f = ∆nc,b L f0δ 1 + µvsmp f0 + ξµvsmp f0δ . (42)

Beam theory

To relate (42) to the simulations we have to make further assumptions in order to express the force explicitly in terms of bond parameters. As used before, the bonds between particles are beams that bend and break. For a beam with bond-radius rb and length 2rp (the distance between two particles in contact) the situation is schematically shown in Fig. 26. In this specific case, Eq. (12) can be written as:

δ ∼ 2σbr 4 p Er4 b (43) f0 ∼ 2σbr3b 3rp . (44)

(43)

Figure 26: Schematic beam bending

Substituting the latter relations in Eq. (42), all parameters are known from the simula-tions and we could compare the model to existing data. The aim would be that all different simulations where we vary E, µ, rb an σb are explained by this expression.

5.1.1 Jamming effects

The length scale ξ that enters the model mimics already some compaction in front of the cone. Elements that break stay in contact over a distance ξ before they clear out and have no further influence. It is expected that this length ξ is not constant but depends on the number of particles that is already in contact with the cone. The escape from the cone will become increasingly difficult when particles pile when the total number of elements approaches a maximum number nmax. To mimick this jamming effect we assume that the length scale ξ depends on the total number of elements under the cone according to

ξ = ξ0 1 −n cone b + nconeu nmax −1 . (45)

In case that the particles pile up there will be an increasing force due to the friction term in Eq. (36). This setting changes Eq. (38) according to:

d dzn cone u (z) = n0fail(z) − 1 ξ0 ncone_u 1 −n cone b + nconeu nmax . (46)

As before the number of unbonded elements under the cone must be stationary. Aver-aging Eq. (46) and ignoring correlations such that we can replace averages of products with products of averages, we obtain

0 = n0 fail− 1 ξ0 ncone u 1 − n cone b + nconeu nmax . (47)

(44)

Figure 27: ¯f computed via Eq. (49) is plotted as a function of ∆nc,b/L. The relation is super-linear. ncone u = nmax− nconeb 2  1 ± s 1 − 4ξ0 n0 failnmax (nmax− nconeb )2  . (48)

The + solution is asymptotically unstable if (46) is interpreted as a Verhulst-type of growth model. Therefore we proceed with the − solution. Substituting equations (41) and (47), the average force can finally be expressed as

¯ f =1 2 ∆nc,b L δf0+ µvsmp   ∆nc,b L δ + nmax− ∆nc,b L δ 2  1 − v u u t1 − 4ξ₀ ∆nc,b L nmax (nmax− ∆nc,b L δ)2    . (49) Due to the jamming effect, in this model ¯f has a super-linear dependence on ∆nc,b/L. This is shown for a chosen set of parameters in Fig. 27. This demonstrates how observed features (cf. Fig. 19) can be approximately accounted for in phenomenological model ap-proached like Eq. (37). Further analysis is however required to test the validity of the idea from Eq. (45).

(45)

6 Discussion

Main results

Various existing models that describe penetration of snow [21] or cellular materials in general [9, 14] focus on the mean or maximum penetration force and describe it in terms of elements that fail. In the case of the shot noise model [21] the elements are abstract entities and spatially uncorrelated, while the models for cellular materials use the picture of a foam comprising a network of rods [9, 14]. In this case beam theory (Eq. (12)) is used to express the elastic-brittle force to break the rods in terms of material properties. In our simulations, the elements that fail are the bonds between the particles (cf. Sec. 3.3). In view of the existing models as a starting point, it is of key importance to establish a relation between the average force and the rate at which bonds break in the system. It is generally not obvious that the number of bonds that fail in the entire system is well suited to describe the penetration force with, because there is dissipation of energy via particle interactions. In addition, existing models explicitly take into account only elements in contact with the cone [21, 9] or in touch with the penetrometer effective surface [14].

For the present set of DEM simulations (Sec. 3.4) it has been shown that the average force ¯

f can still be described mainly by the average rate ∆nc,b/L at which bonds break in the entire system (Sec. 4.3, Fig.19.a). The relation between both quantities is however non-linear. The first step to discuss the mean force is to try to explain it with the shot noise model. Eq. (33) expresses the mean force ¯f from the shot noise model in terms of simulation variables and data. This equation already gives the right order of magnitude for ¯f for the reference dataset, as shown in Fig. 19.b. The agreement is particularly good for low density (red circles) where the model assumptions are best met. The forces in the scaled dataset are off by approximately a factor of 10. To understand this difference in coefficient we must look beyond the shot noise model. The dependence of the average force on different parameters in Eq. (33) is further tested in some sensitivity studies (Sec. 4.4). It is confirmed that ¯f has a quadratic dependence on the bond radius rb, as predicted by (33), but the dependence on σb disagrees with the model. Eventually the various disagreements constitute the basis for new model ideas.

A neat empirical benchmark for a new model is obtained by imposing a nonlinear relation of ¯f on ∆nc,b/L in Eq. (33): Instead of the linear dependence, an empirical test was done by relating ¯f to ∆nc,b/L4/3. The results are shown in Fig. 19.c. Amost all simulations collapse on the 1:1 line. This implies that the dependence of ¯f on ∆nc,b/L in a new model should be super-linear. A possible explaination for this super-linearity is a post-failure influence of broken bonds/material on the penetration force. Every bond that breaks does not only contribute to the force by its rupture, but adds an additional non linear contribution due to friction. When a bond breaks, the particles are expected to remain in the region of the cone and give an extra contribution to the force. This makes ¯f super-linearly depending on ∆nc,b/L. To investigate this effect further and describe the physical origin of this

Discrete element simulations of micro-penetration in cohesive granular materials

Master thesis