Random-Graph Models and Characterization of Granular Networks

TL;DR

This paper evaluates various network measures to characterize and distinguish between different spatial random-graph models and empirical granular networks, revealing key differences across dimensions and model plausibility.

Contribution

It identifies network measures capable of differentiating physically plausible models from unphysical ones and highlights dimensional differences in network measure distributions.

Findings

Certain network measures distinguish plausible from unphysical models.

Significant differences in measure distributions between 2D and 3D networks.

Implications for understanding experimental granular networks.

Abstract

Various approaches and measures from network analysis have been applied to granular and particulate networks to gain insights into their structural, transport, failure-propagation and other systems-level properties. In this article, we examine a variety of common network measures and study their ability to characterize various two-dimensional and three-dimensional spatial random-graph models and empirical two-dimensional granular networks. We identify network measures that are able to distinguish between physically plausible and unphysical spatial network models. Our results also suggest that there are significant differences in the distributions of certain network measures in two and three dimensions, hinting at important differences that we also expect to arise in experimental granular networks.

Tables

Table 1. Table ​​ 1: Network diagnostics of the experimental data, the two-dimensional modified-lattice model, and the two-dimensional RGGs. We average the results over 17 realizations for the experimental data and over 1023 realizations for the random-graph models.

	Experimental data		Modified lattice
Diagnostic	Mean	Standard deviation	Mean	Standard deviation
Edge density	0.0038		0.0042
Geodesic node betweenness	0.0163	0.0002	0.0168	0.0001
Geodesic edge betweenness	0.0080	0.0002	0.0083
Clustering coefficient	0.2613	0.0100	0.3019	0.0068
Transitivity	0.2714	0.0070	0.2994	0.0042
Degree assortativity	0.1377	0.0236	0.0648	0.0262
Global efficiency	0.0772	0.0005	0.0812	0.0003
Local efficiency	0.3301	0.0156	0.3886	0.0089
Mean shortest-path length in the LCC	19.1634	0.1069	18.1703	0.0585
Weighted mean shortest-path length	19.1483	0.1046	18.1456	0.0680
Maximized modularity	0.8712	0.0045	0.8487	0.0028
Subgraph centrality	7.6236	0.2958	8.0579	0.0678
Communicability	0.0790	0.0046	0.0937	0.0013

Table 2. Table ​​ 2: Network diagnostics of the 2D proximity-modified RGG and the 2D force-modified RGG. We average the results over 1023 1023 1023 realizations.

	Proximity-modified RGG (mono)		Proximity-modified RGG (bi)
Diagnostic	Mean	Standard deviation	Mean	Standard deviation
Edge density	0.0038		0.0038
Geodesic node betweenness	0.0171	0.0002	0.0166	0.0002
Geodesic edge betweenness	0.0085	0.0001	0.0082	0.0001
Clustering coefficient	0.3540	0.0070	0.3356	0.0066
Transitivity	0.3584	0.0044	0.3389	0.0040
Degree assortativity	0.2799	0.0270	0.2179	0.0268
Global efficiency	0.0734	0.0006	0.0754	0.0006
Local efficiency	0.4663	0.0100	0.4430	0.0099
Mean shortest-path length in the LCC	20.1977	0.1987	19.5844	0.1852
Weighted mean shortest-path length	20.1808	0.2012	19.5716	0.1858
Maximized modularity	0.8754	0.0031	0.8670	0.0032
Subgraph centrality	8.6839	0.2133	8.5141	0.2257
Communicability	0.0874	0.0036	0.0879	0.0040

Equations

$$\rho =\frac{|E|}{|V|(|V|-1)},$$

$$c_B(v)=\sum _{⟨s,t⟩\in V}\frac{{\sigma }_{st}(v)}{{\sigma }_{st}}\underset{\text{normalization factor}}{\underbrace{\frac{2}{(|V|-1)(|V|-2)}}}.$$

$$C_B(G)=\frac{1}{|V|}\sum _{v\in V}{c}_B(v)$$

$$c_B(e)=\sum _{⟨s,t⟩\in V}\frac{{\sigma }_{st}(e)}{{\sigma }_{st}}\underset{\text{normalization factor}}{\underbrace{\frac{2}{|V|(|V|-1)}}}.$$

$$C_B(G)=\frac{1}{|E|}\sum _{e\in E}{c}_B(e).$$

$$c(v)=\frac{2T(G;v)}{\text{deg}(v)(\text{deg}(v)-1)}.$$

$$c(G)=\frac{1}{|V|}\sum _{v\in V}c(v).$$

$$\text{Transitivity}=\frac{T(G)}{\text{Number of possible triangles in }G}.$$

$$E(u,v)=\frac{1}{d_{\text{min}}(u,v)}.$$

$$E_{\text{glob}}(G)=\frac{1}{|\{⟨u,v⟩\in V\}|}\sum _{⟨u,v⟩\in V}E(u,v).$$

$$E_{\text{loc}}(v)=\frac{1}{|⟨u,w⟩\in \Gamma (v)|}\sum _{⟨u,w⟩\in \Gamma (v)}E(u,w),$$

$$E_{\text{loc}}(G)=\frac{1}{|V|}\sum _{v\in V}{E}_{\text{loc}}(v),$$

$${\bar{d}}_{\text{min}}=\frac{1}{|\{⟨u,v⟩\in C_{\text{max}}\}|}\sum _{⟨u,v⟩\in C_{\text{max}}}{d}_{\text{min}}(u,v).$$

$${\bar{d}}_{\text{min}}^w=\frac{1}{|V|}\sum _{G_C\in G}{\bar{d}}_{\text{min}}(G_C)|G_C|,$$

$$Q=\frac{1}{2|E|}\sum _{⟨u,v⟩\in G}\left[A_{uv}-\frac{\text{deg}(u)\text{deg}(v)}{2|E|}\right]\delta (g_u,g_v),$$

$$S(v)=\sum _{k=0}^{\infty }\frac{{\left(A^k\right)}_{vv}}{k!}$$

$$S(G)=\frac{1}{|V|}\sum _{v\in V}S(v).$$

$$C_{\mathrm{com}}(u,v)=\sum _{k=0}^{\infty }\frac{(A^k{)}_{uv}}{k!},$$

$$C_{\mathrm{com}}(G)=\frac{1}{|\{⟨u,v⟩\in V\}|}\sum _{⟨u,v⟩\in V}{C}_{\mathrm{com}}(u,v).$$

$$P(p_k)=\{\begin{array}{ll}0,& \text{for}{p}_k<2\alpha r_k,\\ 1,& \text{for}{p}_k\ge 2\alpha r_k,\end{array}$$

$${\mathbf{f}}_k=$$

$$+{\left[r_k-x_k\right]}_{+}^{\beta }\hat{\mathbf{x}}-{\left[r_k+x_k-L_x\right]}_{+}^{\beta }\hat{\mathbf{x}}+{\left[r_k-y_k\right]}_{+}^{\beta }\hat{\mathbf{y}}-{\left[r_k+y_k-L_y\right]}_{+}^{\beta }\hat{\mathbf{y}},$$

$$[x{]}_{+}=\{\begin{array}{ll}x,& \text{for }x>0,\\ 0,& \text{for }x\le 0.\end{array}$$

$$V=$$

$$+{\left[r_k-x_k\right]}_{+}^{\beta +1}-{\left[r_k+x_k-L_x\right]}_{+}^{\beta +1}+{\left[r_k-y_k\right]}_{+}^{\beta +1}-{\left[r_k+y_k-L_y\right]}_{+}^{\beta +1}.$$

$${\mathbf{f}}_k=$$

$$+{\left[r_k-x_k\right]}_{+}^{\beta }\hat{\mathbf{x}}-{\left[r_k+x_k-L_x\right]}_{+}^{\beta }\hat{\mathbf{x}}+{\left[r_k-y_k\right]}_{+}^{\beta }\hat{\mathbf{y}}-{\left[r_k+y_k-L_y\right]}_{+}^{\beta }\hat{\mathbf{y}}$$

$$+{\left[r_k-z_k\right]}_{+}^{\beta }\hat{\mathbf{z}}-{\left[r_k+z_k-L_z\right]}_{+}^{\beta }\hat{\mathbf{z}},$$

$$V=$$

$$+{\left[r_k-x_k\right]}_{+}^{\beta +1}-{\left[r_k+x_k-L_x\right]}_{+}^{\beta +1}+{\left[r_k-y_k\right]}_{+}^{\beta +1}-{\left[r_k+y_k-L_y\right]}_{+}^{\beta +1}$$

$$+{\left[r_k-z_k\right]}_{+}^{\beta +1}-{\left[r_k+z_k-L_z\right]}_{+}^{\beta +1}.$$

$$V_{\text{mono}}^{\ast }$$

$$=\sum _{k}6\times {\left(\frac{r_k}{200}\right)}^{\beta +1}=1122\times 6\times {\left(\frac{4.5}{200}\right)}^{1.5}$$

$$\approx 22.7$$

Full text

I Introduction

Various tools from network analysis have yielded insights into transport, failure mechanisms and other system-level properties of granular networks, which are typically constructed by interpreting the particles of the underlying packings as nodes and interpreting their contacts as edges smart2007effects ; smart2008granular ; berthier2018forecasting ; papadopoulos2018network . Network analysis may help lead to better understanding and the ability to control fracture processes, design new materials, and assess structural degradation in various engineering applications herrmann2014statistical ; anderson2017fracture ; heisser2018controlling . Researchers have employed continuum anderson2017fracture and particle-level herrmann2014statistical ; cates1999jamming descriptions to study granular materials jaeger1996 ; andreotti2013 , but neither framework is designed to study intermediate-scale organization, which is important for understanding and characterizing granular packings bassett2015extraction ; papadopoulos2018network .

In studies of granular and particulate networks, it is not clear which network measures are most suitable for distinguishing between physical granular networks and randomly embedded and overlapping particle systems. To take advantage of methods of network analysis, it is necessary to identify appropriate network measures to best highlight the distinctions between different structures barthelemy2018complex . Building on ideas from Ref. porterstudentthesis , we compute a variety of common network diagnostics and apply them to several spatial network models and empirical granular networks in two and three dimensions. We identify measures that are able to distinguish between the network characteristics of physical and unphysical spatial networks that we construct from random particle systems. We test the convergence properties of our results by considering large ensembles of granular networks. Our results are useful both for the study of granular networks and because they help provide a roadmap for the development of random-graph models for studying spatially-embedded networks (such as transportation, communication and vascular systems). As has also been demonstrated in applications like sensor systems coon2012impact ; dettmann2016 , incorporating a small amount of physics into models such as random geometric graphs can provide important insights into the structure of such systems barthelemy2018complex .

Our paper proceeds as follows. In Sec. II, we describe the network measures that we use in our investigation of spatial random graphs and granular networks. In Sec. III, we present several models of both unphysical and physical granular networks. The network structures range from random geometric graphs (RGGs) to physically plausible contact networks. To ensure comparability of our results with properties of real granular networks, we also apply our network measures to empirical granular network structures, and we generate our simulated networks according to empirically observed edge densities. We thereby examine if certain network measures are able to distinguish between clearly unphysical (i.e., overlapping) granular packings and ones that are potentially physically plausible. In Sec. IV, we compare the distributions of the network diagnostics on the various networks. We also examine the convergence properties of the considered network diagnostics by computing them for ensembles of large networks. In Sec. V, we conclude our study and discuss our results.

II Network diagnostics

The structure of a network depends fundamentally on the processes and interactions that define it newman2018networks . Many diagnostics have been developed to study and characterize networks, and we seek to use some of them to study granular networks and examine connections between network and particulate systems.

Consider a graph (i.e., network) $G(V,E)$, where $V$ is a set of nodes and $E$ is a set of edges. For each node $v\in V$, the degree $\text{deg}(v)$ is the number of edges that are attached to $v$. For each node pair $\langle s,t\rangle\in V$, let $d_{\text{min}}(s,t)$ denote a shortest path between $s$ and $t$. We use $\sigma_{st}$ to denote the number of shortest paths in the graph from $s$ to $t$. The number of shortest paths from $s$ to $t$ that traverse $v\in V$ is $\sigma_{st}(v)$, and $\sigma_{st}(e)$ denotes the number of shortest paths from $s$ to $t$ that traverse an edge $e\in E$. We let $T(G)$ denote the number of triangles (i.e., a set of three nodes with an edge between each pair of nodes) present in the graph. A triangle is a set of three nodes that are all pairwise connected by an edge. We denote the number of triangles that include a node $v\in V$ by $T(G;v)$. We count each triangle three times when computing the number of triangles in a graph $G$. All of the networks that we consider are unweighted and undirected, and we discuss several ways for characterizing them in the following paragraphs. To compute the values of these diagnostics, we use the Python library NetworkX. See Ref. networkx_doc for a detailed documentation.

Edge density.

The edge density (or sometimes just “density”) of an undirected network is newman2018networks

[TABLE]

where $|E|$ is the numbers of edges and $|V|$ is the number of nodes. We use edge density as a reference measure for the construction of our networks. To compare the results of different network models, we first make sure that their edge densities are the same. From this density and the number of nodes, we derive the number of edges, which we use in our generation of random graphs.

Geodesic node betweenness.

The geodesic node betweenness centrality $c_{B}$ for a node $v$ is the sum (over all node pairs $\langle s,t\rangle$ in a network) of the number of shortest paths from $s$ to $t$ that traverse $v$ divided by the total number of shortest paths from $s$ to $t$ newman2018networks . Including a normalization factor, we write geodesic node betweenness centrality as

[TABLE]

Note that $\sigma_{st}=1$ if $s=t$ and that $\sigma_{st}(v)=0$ if $v\in\{s,t\}$. We define a network’s geodesic node betweenness

[TABLE]

as the mean of geodesic node betweenness centrality over all nodes in a network.

A large geodesic betweenness centrality $c_{B}(v)$ suggests that many shortest paths traverse a particular node. This may in turn suggest that such a node is relevant for transportation processes or failure-propagation mechanisms in granular networks smart2007effects ; smart2008granular ; berthier2018forecasting ; papadopoulos2018network .

Geodesic edge betweenness.

We can also calculate a geodesic betweenness centrality for edges girvan2002community :

[TABLE]

Note that $\sigma_{st}(e)=0$ if $s=t$, and we again take $\sigma_{st}=1$ when $s=t$. Based on Eq. (4), the geodesic edge betweenness of a network is

[TABLE]

Clustering coefficient.

The local clustering coefficient $c(v)$ of a node $v$ with degree $\text{deg}(v)\geq 2$ is the number of triangles that include that node divided by all possible triangles that could include $v$ (i.e., the connected triples) newman2018networks . That is,

[TABLE]

If the degree of a node is [math] or $1$, the corresponding clustering coefficient of that node is [math]. The mean local clustering coefficient of a network is

[TABLE]

Transitivity.

Transitivity is the number of triangles in a network divided by the number of possible triangles newman2018networks :

[TABLE]

Based on the definition of clustering coefficient and transitivity, we see that the two measures are closely related.

Degree assortativity.

Degree assortativity is a measure to quantify the correlation between the degree of a certain node and the degrees of its neighbors newman2018networks . Larger values of degree assortativity entail greater similarity in the degree of neighboring nodes. We compute degree assortativity by calculating a Pearson correlation coefficient, whose values lie in the interval $[-1,1]$ Newman2003mixing . A value of $1$ indicates a network with perfect degree assortativity, a value of [math] corresponds to a nonassortative network, and a value of $-1$ indicates perfect degree disassortativity.

Global efficiency.

The efficiency $E(u,v)$ of a node pair $\langle u,v\rangle$ is the inverse of the shortest-path length between these two nodes latora2001efficient . That is,

[TABLE]

Global efficiency is the mean of the efficiencies over all pairs of nodes in a network:

[TABLE]

Local efficiency.

The local efficiency latora2001efficient of a node $v$ is

[TABLE]

where $\Gamma(v)$ denotes the neighborhood of $v$. The local efficiency of a network is

[TABLE]

the mean of the local efficiencies over all nodes of the network.

Mean shortest-path length in the largest connected component (LCC).

The mean shortest-path length over all node pairs $\langle u,v\rangle$ of the largest connected component (LCC) $C_{\text{max}}$ of a network is

[TABLE]

Weighted mean shortest-path length.

The weighted mean shortest-path length $\bar{d}_{\text{min}}^{w}$ is the weighted arithmetic mean of the mean shortest-path lengths $\bar{d}_{\text{min}}(G_{C})$ that we compute for connected components $G_{C}$ of a graph. We use the number of nodes in each component as a weight, so

[TABLE]

where $|G_{C}|$ denotes the number of nodes in the considered connected component.

Maximized modularity.

Subgraphs of a network $G$ are often called “communities” in the context of dense sets of nodes that are connected sparsely to other dense sets of nodes porter2009 ; fortunato2016 . Given a set of communities of a network $G$, the modularity $Q\in[-1,1]$ measures the density of edges (adjusted for edge weights) inside communities compared with that between communities Blondel2008communities . For an unweighted graph, the modularity $Q$ of a partition of a graph $G$ is

[TABLE]

where $A$ is the adjacency matrix of the network, $\delta$ is the Kronecker delta, and $g_{u}$ denotes the community that includes node $u$.

We maximize modularity by applying the locally greedy Louvain method for community detection Blondel2008communities . Maximized modularity measures how well one can partition a network into disjoint communities porter2009 .

Subgraph centrality.

For a node $v$ and an adjacency matrix $A$ associated with a graph $G$, the subgraph centrality EstradaRodriguez

[TABLE]

is the sum of weighted closed walks (of all lengths) that start and end at node $v$. Instead of weighting each term with $1/k!$, one can also make other choices for the weights. The $(u,v)$th element of $A^{k}$ have the following interpretation: each element $(u,v)$ represents the number of paths of length $k$ from node $u$ to node $v$. The mean subgraph centrality of a network $G$ is

[TABLE]

Communicability.

The communicability estrada2008communicability

[TABLE]

between nodes $u$ and $v$ measures the number of closed walks between nodes $u$ and $v$.

The communicability of a network $G$ is the mean of the communicabilities between all pairs of nodes in the network:

[TABLE]

III Random-graph models of granular networks

As baselines to compare to empirical granular networks, we consider various models of spatial networks in two and three dimensions. We start with a standard RGG, which is an unphysical model, and then we examine some more realistic models (such as force-modified RGGs). After discussing these models in two dimensions, we also examine some models in three dimensions. The reason for studying different models, ranging from unphysical ones to more physical ones, is that we aim to identify diagnostics that are able to capture properties of realistic granular networks, which one constructs from non-overlapping particles that interact with each other porterstudentthesis . The model networks also allow us to analyse the convergence properties of network measures, because we can generate large numbers of different network configurations and examine how they behave when they have progressively more nodes. We also apply our network measures to two-dimensional empirical granular networks and use the empirically determined edge densities as input for our granular network models in two dimensions. Given a specific edge density, we perform all of the computations for the same configuration as in the experiments (see Sec. III.5). Specifically, the empirical granular systems are confined by a two-dimensional box with dimensions $290~{}\text{mm}\times 380~{}\text{mm}$ and have a mean of $1122$ particles. To study models of three-dimensional materials, we examine model networks with $1122$ particles in a box with dimensions $102~{}\text{mm}\times 102~{}\text{mm}\times 102~{}\text{mm}$.

III.1 Random geometric graphs

As a first, overly simplistic model of a granular system, we consider a monodisperse random geometric graph (RGG). In this model, we place $n$ particles (represented by nodes) of the same radius $r$ uniformly at random in the box. The nodes are adjacent if the distance between them is smaller than $2R$, where we choose the parameter $R$ based on the edge density that matches an empirical granular network that we will study later (see Sec. III.5). Because of the possibility of multiple overlapping particles, this monodisperse RGG is not a good model of a physical granular system. We also consider a bidisperse RGG. In this model, the particles can have two different radii, $r_{1}$ and $r_{2}$; nodes $i$ and $j$ are adjacent to each other if their distance is smaller than $R\left(r_{i}+r_{j}\right)$. For each particle, we generate a uniformly distributed random number $\epsilon\sim\mathcal{U}(0,1)$. If $\epsilon<0.5$, we set the particle radius to $r_{1}$; otherwise, we set it to $r_{2}$. We use the same procedure for all bidisperse models in this article. We want 50%, on average, of particles to have radius $r_{1}$ to facilitate our comparison with the employed experimental granular network data (see Sec. III.5). As in the monodisperse case, we place particles uniformly at random. We again need to choose $R$ to match a desired edge density. We show examples of both a monodisperse and a bidisperse RGG in Fig. 1.

The parameters in our computations are $n=1122$ nodes, a radius of $r=4.5$, and a distance parameter of $R\approx 5.69$ for the monodisperse RGG and $n=1122$, $r_{1}=4.5$, $r_{2}=5.5$, and $R\approx 1.13$ for the bidisperse RGG.

III.2 Proximity-modified RGGs

To compensate for some of the overlaps between particles in an RGG, we now define a proximity measure $p$ that describes the distance between each point $(x,y)$ and its nearest particle(s). Initially, the system is empty. We then place each particle as follows. For each integer $k\in\{1,\dots,n\}$, we place particle $k$, which has a radius of $r_{k}$, in the box. To do this, we choose a position $(x_{k},y_{k})$ in the box uniformly at random, and we place the particle with probability

[TABLE]

where $p_{k}$ is the proximity measure for the position $(x_{k},y_{k})$ and $\alpha\in[0,1]$ is a parameter that enforces a maximum particle density in a given region. We use this procedure for both monodisperse ($r_{k}=r$) and bidisperse ($r_{k}\in\{r_{1},r_{2}\}$) particle configurations. We generate the edges in the same way as in the bidisperse RGG model.

In Fig. 2, we show examples of monodisperse and bidisperse proximity-modified RGGs. Based on the depicted configurations, we see that the proximity-modified RGGs have fewer overlaps than the standard RGGs. However, overlaps are still present, so the model does not describe a physical granular network.

For our computations, we use the parameters $n=1122$, $r=4.5$, $\alpha=0.75$, and $R\approx 1.39$ for the monodisperse case and $n=1122$, $r_{1}=4.5$, $r_{2}=5.5$, $\alpha=0.75$, and $R\approx 1.24$ for the bidisperse case. In the bidisperse case, we generate a uniformly distributed random number $\epsilon\sim\mathcal{U}(0,1)$ for each particle. If $\epsilon<0.5$, we set the particle radius to $r_{1}$; otherwise, we set it to $r_{2}$.

III.3 Force-modified RGGs

To obtain a network whose characteristics better resemble those of a physical granular packing, we now apply a different modification to an RGG. Initially, we start from an RGG configuration. Subsequently, for each particle $k$, we suppose that particles interact via a Hertzian-like contact owens2011sound and compute the force $\mathbf{f}_{k}$ that acts on particle $k$ from overlaps with other particles. After computing $\mathbf{f}_{k}$ for all $k$, we update the corresponding particle locations $\mathbf{x}_{k}$ according to $\mathbf{x}_{k}\rightarrow\mathbf{x}_{k}+\mathbf{d}_{k}$, with $\mathbf{d}_{k}=\varepsilon\mathbf{f}_{k}$. We repeat this process until the particles no longer overlap (see Appendix A). One can use the parameter $\varepsilon$ to rescale the displacements $\mathbf{d}_{k}$, which may be too large and potentially lead to new overlaps between particles if we do not reduce them. For a box with dimensions $L_{x}\times L_{y}$ and according to Hertzian contact theory, the force on particle $k$ is owens2011sound ; chong2017nonlinear

[TABLE]

where the exponent $\beta$ (which is $\beta=3/2$ in classical Hertzian theory, but was estimated to be $\beta=5/4$ for the 2D granular material in Ref. owens2011sound ) describes the interaction strength between particles. The vectors $\hat{\mathbf{x}}$ and $\hat{\mathbf{y}}$ are the unit vectors in the $x$ and $y$ directions, respectively. Forces occur only between overlapping particles or between a particle and the box boundary if both are in contact. Therefore, the bracket $[x]_{+}$ in Eq. (21) is

[TABLE]

The potential energy that corresponds to Eq. (21) is

[TABLE]

The first term of Eqs. (21) and (23) accounts for overlaps between particles, and the other terms take into account the left, right, bottom and top boundaries of the box. We generate the edges in this model as we did for bidisperse RGGs. We show examples of monodisperse and bidisperse force-modified RGGs in Fig. 3. In this figure, it appears (based on visual inspection) that particles do not overlap anymore. At least visually, they more closely resemble real granular networks than is the case for the standard RGGs and the proximity-modified RGGs.

To attempt to ensure that our final configurations are free of overlaps, we monitor their energy. If the energy is 0, there are no overlaps, so we stop updating locations when we cannot distinguish the energy from 0. However, in some cases, the energy does not approach 0, because some particles are locally jammed behringer2018physics . In such cases, we accept a configuration if its energy is smaller than the energy of a system in which every particle $k$ has an overlap distance of no more than $r_{k}/200$ with all of its neighbors. In Fig. 4, we show the monitored energy for two cases of the monodisperse force-modified RGG. A situation in which many particles overlap with their neighbours is unlikely to occur. Instead, there are often some particles that overlap with one or two neighbours with a large fraction of their areas. To avoid large overlaps, we also check for every configuration that each particle $k$ does not have an overlap distance with a neighbour that is larger than $r_{k}/15$. See Appendix A for details.

The parameters in our computations are $n=1122$, $r=4.5$, $\beta=0.5$, $L_{x}=290$, $L_{y}=380$, and $R\approx 1.2900$ for the monodisperse case. For the bidisperse case, they are $n=1122$, $r_{1}=4.5$, $r_{2}=5.5$, $\beta=0.5$, $L_{x}=290$, $L_{y}=380$, and $R\approx 1.0850$. To generate these configurations, we perform $100$ iterations for the monodisperse case and $400$ iterations for the bidisperse one. As in the RGG model, 50% of the particles, on average, in the bidisperse case have a radius of $r_{1}$; and the others have a radius of $r_{2}$.

III.4 Modified lattice

We also consider a modified lattice. In this model, we modify a hexagonal lattice of particles of the same radii by removing a fraction $\mu$ of edges uniformly at random to match a certain edge density. We show an example of a modified lattice in Fig. 5. In our computations, we use a total of 1020 particles with $r=4.5$ and $\mu\approx 0.26$.

III.5 Experimental data

To compare our model networks with real granular networks, we extract granular networks from experimental data bassett2012influence . Henceforth, when we write “experimental data”, we mean data from a bidisperse granular system with radii $r_{1}=4.5~{}\text{mm}$ and $r_{2}=5.5~{}\text{mm}$ in a two-dimensional box with side lengths $L_{1}=290~{}\text{mm}$ and $L_{2}=380~{}\text{mm}$. The employed experimental data set consists of $17$ different configurations of granular networks. We compute the diagnostics of Sec. II for these $17$ cases, and we present the corresponding diagnostic values in Appendix B. On average, the number of particles in the experimental configurations is $1122$, and $49.2$% of the particles have a radius of $r_{1}$. Therefore, we use $1122$ particles in our model networks to facilitate comparisons between our models and the experimental data. For the bidisperse network models, we generate a uniformly distributed random number $\epsilon\sim\mathcal{U}(0,1)$. If $\epsilon<0.5$, we set the particle radius to $r_{1}$; otherwise, we set it to $r_{2}$. We consider bidisperse packings, as we wish to avoid the crystallization of monodisperse packings. We show an example of the experimental setup in Fig. 6. To obtain such a configuration, Ref. bassett2012influence placed about $1000$ particles (which were cut from Vishay PSM-4 photoelastic material) into a container with an open top, such that the particles were confined only by gravity.

III.6 Models in 3D

As with the 2D examples, we calculate the network diagnostics of Sec. II for 3D granular networks. In this section, we briefly discuss extensions of the previously-discussed models to 3D. The standard RGGs do not require any special modification for the generation of the configurations, although the underlying box is now in 3D. We show examples of 3D monodisperse and bidisperse RGGs in the top panels of Fig. 7. As in the 2D RGGs, we observe that the particles overlap, so the model is not physical. Our parameters for these models are $n=1122$, $r=4.5$, and $R\approx 6.0750$ for the monodisperse case and $n=1122$, $r_{1}=4.5$, $r_{2}=5.5$, and $R\approx 1.1950$ for the bidisperse case. As in the 2D models, 50% of the particles in our 3D models have radius $r_{1}$ on average and the others have radius $r_{2}$.

It is also straightforward to extend the proximity-modified RGGs to 3D. We show examples of 3D monodisperse and bidisperse proximity-modified RGGs in the bottom panels of Fig. 7. We see that the particles overlap less than for the unmodified 3D RGGs, but they still overlap. To obtain these configurations, we use the parameters $n=1122$, $r=4.5$, $\alpha=0.75$, and $R\approx 1.4200$ for the monodisperse case and $n=1122$, $r_{1}=4.5$, $r_{2}=5.5$, $\alpha=0.75$, and $R\approx 1.4385$ for the bidisperse case.

For the force-modified RGG, we have to adapt Eqs. (21) and (23) to 3D. Equation (21) becomes

[TABLE]

where $\hat{\mathbf{z}}$ denotes the unit vector in the $z$ direction. Equation (23) becomes

[TABLE]

We show examples of 3D force-modified RGGs in the top panels of Fig. 8. As in the 2D case, we conclude based on visual inspection that they resemble real granular networks more closely than the other model. The parameters in our computations are $n=1122$, $r=4.5$, $\beta=0.5$, and $R\approx 1.4080$ for the monodisperse case and $n=1122$, $r_{1}=4.5$, $r_{2}=5.5$, $\beta=0.5$, and $R\approx 1.2600$ for the bidisperse case. In both examples, we set $L_{x}=L_{y}=L_{z}=102$. To generate these configurations, we perform $100$ iterations for the monodisperse case and $400$ iterations for the bidisperse one. We describe the energy thresholds for the 3D force-modified RGGs in Appendix A.

We extend the modified lattice to 3D by considering spheres in a hexagonal close packing (HCP). An HCP consists of a series of two alternating layers of spherical particles, where each layer has a hexagonal arrangement. We create the edges in the same way as in the 2D modified lattice. In the bottom panel of Fig. 8, we show an example of a 3D modified lattice that consists of $1100$ particles. The particles have a radius of $r=4.5$, and we remove a fraction $\mu\approx 0.4950$ of the edges uniformly at random from the network.

IV Results

IV.1 Diagnostics for 2D granular networks

We first focus on the results of calculating the network diagnostics for the 2D models from Sec. III. We describe the convergence characteristics of all diagnostics and then discuss the distributions of several of them.

IV.1.1 Convergence characteristics

To monitor the convergence of the diagnostics, we generate $10$ ensembles with different numbers of granular-network configurations. Specifically, ensemble $m$ has $2^{m}$ networks (where $m\in\{0,\dots,9\}$), so there are $1023$ networks in total. For each $m$, we compute the mean of the diagnostic values and normalize each value by the mean that we calculate for $m=9$. Therefore, each of the mean diagnostic values approaches $1$ for $m=9$ (see Fig. 9). This procedure allows us to identify a reasonable number of networks such that the deviations of the mean diagnostic values are sufficiently small in comparison to the values in the largest ensemble. We thereby study the convergence properties of the mean diagnostic values with respect to the values in the largest ensemble. Based on the data in Fig. 9, we see that the diagnostics behave differently in the different network models.

In the unmodified RGG models, some of the normalized mean diagnostics deviate more from $1$ than in the other models. We observe this behavior for geodesic node betweenness, geodesic edge betweenness, subgraph centrality, and communicability. Because we place the particles in the RGGs uniformly at random, different configurations can differ substantially from each other. For example, in our computations, we observe that some particle configurations have local accumulations of particles, and such outliers strongly influence the mean over ensembles with a small number of networks. Such outliers exert less influence for ensembles with many networks, so they have a smaller effect on those ensembles. For the aforementioned diagnostics, such particle accumulations in some samples lead to the clearly visible zigzag behavior in Fig. 9. In the other network models, this issue is less prominent, because local particle accumulations are suppressed by the proximity and force modifications. The larger deviations of the mean diagnostics in the unmodified RGG models is also reflected by their distributions (see Sec. IV.1.3). For the proximity-modified RGGs and the force-modified RGGs, the mean diagnostic values approach the values that we observe in the largest ensemble faster than in the unmodified RGG models. For the modified-lattice model, we observe that almost all mean diagnostic curves are constant; the only exception is degree assortativity. In this model, the only process that has any randomness is deleting edges from the initially-complete hexagonal structure (in 2D) or HCP packing (in 3D).

IV.1.2 Numerical values

We summarize the diagnostic values for the models and experimental granular networks in Appendix B. For each model, our results are means over all $1023$ configurations that we generated to study the convergence characteristics of the diagnostics in Sec. IV.1.1. For the experimental data, we take means over the $17$ network realizations.

Based on the numbers in Tables 1 and 2 in Appendix B, we note several differences between our models. From the results, it is clear that the standard RGGs are the most unphysical of these model; we already made this point when visualizing an associated particle system. The proximity-modified RGGs have diagnostic values that resemble the experimental values more closely than the standard RGGs, and the force-modified RGGs generate the most realistic configurations. The diagnostic values for the modified-lattice model agree partly with those in the experimental data, but they are not better than the values for the force-modified RGGs, and the values of degree assortativity in the modified-lattice model differ significantly from those in the experimental configurations. The modified lattices are also visually rather different from the experimental packings, as noted in Ref. porterstudentthesis . A closer look at the diagnostics of the bidisperse force-modified RGG reveals that it is the model that best matches the experimental data. This is also apparent from the distributions in Sec. IV.1.3 and Appendix C.

IV.1.3 Distributions of network characteristics

We now discuss the distributions of the network characteristics in ensembles of networks for the 2D models and the experimental data. We show the distributions of the clustering coefficient, transitivity, mean shortest-path length in the LCC, subgraph centrality, geodesic edge betweenness, and communicability in Fig. 10. We show the distributions of the other measures in Appendix C. In each plot, the black curve is the distribution of the experimental data.

The top panels of Fig. 10 suggest that the mean local clustering coefficient and transitivity are useful measures for distinguishing between the unphysical and physical models. Specifically, the distributions of the standard RGGs do not overlap with the distributions from the other models, so we can clearly distinguish between the RGG distributions and those of the other models. A closer look at the distributions of mean local clustering coefficient and transitivity also reveals that the bidisperse force-modified RGG distributions and the empirical granular networks have characteristics that are distinct from those of the other models. Based on these results, it is possible to distinguish these two network types from the proximity-modified RGGs and the monodisperse force-modified RGG. The modified-lattice distribution overlaps slightly with the one from the experimental data, although visual inspection of the networks indicates that these two types of networks are rather distinct from each other.

The distributions of the mean shortest-path length in the LCC in the experimental granular networks are similar to those that we obtain for the bidisperse force-modified RGG. By contrast, the distributions of the mean local clustering coefficient and transitivity in the experimental networks differ from those in all of the models. This difference is extremely large for the unphysical models, but there is some overlap between experimental and model distributions for other models (e.g., more physical ones). Therefore, computing these quantities may not yield a clear distinction between experimental and model distributions when sampling from data.

In the remaining three panels of Fig. 10 we show the distributions of subgraph centrality, geodesic edge betweenness and communicability. In the corresponding insets, we show the distributions for the standard monodisperse and bidisperse RGG models, and we conclude that they differ qualitatively from the distributions of the corresponding diagnostics in the other models. We can connect this observation to the convergence characteristics of subgraph centrality, geodesic edge betweenness and communica- bility. To study these convergence characteristics, we compare the mean values of all diagnostics for different numbers of network realizations to the corresponding mean values in the largest ensemble (see Sec. IV.1.1). For subgraph centrality, geodesic edge betweenness and communicability in the unmodi- fied RGGs, we observe larger deviations of the mean diagnostic values in the smaller ensembles from the mean values in the largest ensemble than is the case for the other diagnostics (see Fig. 9). This observation is also reflected by the qualitatively different distributions of those diagnostics (see Fig. 10).

IV.2 Diagnostics for 3D granular networks

We now briefly discuss the behavior of the diagnostic distributions for the 3D network models. We again focus on mean local clustering coefficient, transitivity, mean shortest-path length in the LCC, subgraph centrality, geodesic edge betweenness and communicability. We show the corresponding distributions in Fig. 11 and the remaining ones in Appendix C. For each model, we base the distributions on $256$ network realizations.

In the top panels of Fig. 11, we show the distributions of mean local clustering coefficient and transitivity. As in our 2D models, these two diagnostics distinguish successfully between physical and unphysical models of granular networks. The only major difference between our results in 2D and 3D are the distributions (which we show in the insets in some panels) in the modified-lattice model. In 3D, the mean local clustering coefficient and transitivity values for this model are to the left of the bidisperse force-modified RGG model.

For mean shortest-path length in the LCC, we observe that, in contrast to our observations in 2D, the distribution of the bidisperse proximity-modified RGG is to the left of the bidisperse force-modified RGG. Moreover, the distributions of the monodisperse and bidisperse force-modified RGGs and the monodisperse proximity-modified RGG overlap more in 3D than in 2D. We also observe these differences between 2D and 3D in the distributions of geodesic edge betweenness. The distributions of subgraph centrality and communicability in 3D are qualitatively more similar to the ones in 2D than is the case for the other network diagnostics.

V Conclusions and discussion

We described and computed several common network diagnostics for a variety of 2D and 3D models of granular networks, ranging from unphysical ones with overlapping particles to physical ones that are free (or at least mostly free) of overlaps. We studied the convergence properties of the diagnostics to identify a reasonable number of network realizations for the computation of the diagnostic distributions

We examined the ability of the various network diagnostics to distinguish between physical and unphysical models of granular networks in 2D and 3D. Our results suggest that mean local clustering coefficient and transitivity, which are related measures of triadic closure, are appropriate diagnostics to distinguish between the physical and unphysical models. The power of these two measures lies in the fact that their distributions for the different models are readily distinguishable from each other. In some cases, the diagnostics distributions of the 2D and 3D networks appear to satisfy markedly different statistical distributions, a phenomenon that we also expect to observe in comparisons of 2D and 3D granular networks from experiments.

Our results advance spatial random-graph models of granular networks and illustrate that simple modifications of random geometric graphs that incorporate a minimal amount of physics are useful models for gaining insights into granular and particulate networks. In future studies, it will be useful to investigate the effects of shape variations, such as by examining the effects of different rectangular box geometries estrada2015random and irregular boundaries dettmann2015more on network diagnostics in spatial random-graph models.

Acknowledgements.

We thank Karen E. Daniels for providing the experimental granular network data bassett2012influence . We also thank David Kammer for the use of his cluster computing infrastructure and Kornel Kovacs for the cluster maintenance.

Appendix A Force-modified configurations

To obtain the potential-energy thresholds $V^{\ast}$ (below which we no longer update particle positions) in Sec. III.3, we consider the first term of Eq. 23 and set the distance between two neighboring particles, $k$ and $l$, to be $r_{k}/200$ (which depends only on one of the particles). In a monodisperse packing, we use $r_{k}=r$, whereas $r_{k}\in{r_{1},r_{2}}$ in bidisperse packings. For simplicity, we consider the case in which every particle has $6$ neighbors. In a monodisperse packing, we use $1022$ particles of radius $r=4.5~{}\text{mm}$ to determine a value of the potential-energy threshold. In a bidisperse packing, we consider a system with $561$ particles of radius $r_{1}=4.5~{}\text{mm}$ and $561$ with radius $r_{2}=5.5~{}\text{mm}$.

Supposing that the radii are unitless, the potential-energy thresholds are

[TABLE]

for the monodisperse packing and

[TABLE]

for the bidisperse packing. In both cases, we set $\beta=1/2$, which is the value that we used for our computations in Sec. III. These two values of the potential-energy threshold provide an estimate to assess the quality of network configurations that we generated using the force-modified RGG model. For 3D models, we take the number of neighbors of a particle to be $12$. Therefore, the thresholds in 3D are $V^{\ast}_{\text{mono}}\approx 45.4$ and $V^{\ast}_{\text{bi}}\approx 53.4$.

We also examine the maximum overlap of our configurations. Based on visual inspection of Fig. 12, we conclude that an overlap of $r_{k}/200$ is sufficiently small. However, overlap distances of $r_{k}/20$ and $r_{k}/15$ are visible.

Appendix B Diagnostic values

We present the numerical results for the diagnostics of the various models in Tables 1 and 2.

Appendix C Additional diagnostic distributions

In Figs. 13 and 14, we show the distributions of the diagnostics that we did not show in Sec. IV.