Elastic properties of isotropic discrete systems: Connections between geometric structure and Poisson's ratio

TL;DR

This paper analytically investigates the Poisson's ratio in isotropic discrete systems, establishing theoretical limits based on geometric configurations and comparing them with numerical results, highlighting deviations due to complex geometries.

Contribution

It provides the first analytical derivation of Poisson's ratio limits for isotropic discrete models considering arbitrary geometries.

Findings

Theoretical limits depend on contact normal directions and their alignment with connecting vectors.

Deviations from parallelism reduce the Poisson's ratio limits.

Numerical results show smaller actual restrictions than theoretical predictions.

Abstract

The use of discrete material representation in numerical models is advantageous due to the straightforward way it takes into account material heterogeneity and randomness, and the discrete and orientated nature of cracks. Unfortunately, it also restricts the macroscopic Poisson's ratio and therefore narrows its applicability. The paper studies the Poisson's ratio of a discrete model analytically. It derives theoretical limits for cases where the geometry of the model is completely arbitrary, but isotropic in the statistical sense. It is shown that the widest limits are obtained for models where normal directions of contacts between discrete units are parallel with the vectors connecting these units. Any deviation from parallelism causes the limits to shrink. A comparison of the derived equations to the results of the actual numerical model is presented. It shows relatively large deviations from the theory because the fundamental assumptions behind the theoretical derivations are largely violated in systems with complex geometry. The real shrinking of the Poisson's ratio limit is less severe compared to that which is theoretically derived.

Figures

Equations

$$\displaystyle{\bm{n}}=\begin{cases}\left(\begin{array}[]{cc}\cos\xi&\sin\xi\end{array}\right)&\mathrm{in\ 2D}\\ \left(\begin{array}[]{ccc}\cos\xi\sin\zeta&\sin\xi\sin\zeta&\cos\zeta\end{array}\right)&\mathrm{in\ 3D}\end{cases}$$

$$f_{\xi }(\xi )=\{\begin{array}{ll}\frac{1}{2\pi }& \mathrm{for}\xi \in [0,2\pi ]\\ 0& \mathrm{otherwise}\end{array}{f}_{\zeta }(\zeta )=\{\begin{array}{ll}\frac{\sin \zeta }{2}& \mathrm{for}\zeta \in [0,\pi ]\\ 0& \mathrm{otherwise}\end{array}$$

$$f_{\theta }(\theta )$$

$$\mathrm{in}2D:f_{\chi }(\chi )$$

$$\mathbf{t}=\mathbf{\rho }\cdot \mathbf{n}$$

$$\displaystyle{\bm{\rho}}(\chi)=\left[\begin{array}[]{cc}\cos\chi&-\sin\chi\\ \sin\chi&\cos\chi\end{array}\right]$$

$$\mathbf{n}$$

$$\mathbf{\rho }(\xi ,\zeta ,\chi ,\theta )={\mathbf{\rho }}_z(\xi )\cdot {\mathbf{\rho }}_y(\zeta )\cdot {\mathbf{\rho }}_z(\theta )\cdot {\mathbf{\rho }}_y(\chi )\cdot {\mathbf{\rho }}_y^T(\zeta )\cdot {\mathbf{\rho }}_z^T(\xi )$$

$$\cos \chi =\mathbf{n}\cdot \mathbf{t}=\mathbf{n}\cdot \mathbf{\rho }\cdot \mathbf{n}=\mathbf{\rho }:\left(\mathbf{n}\otimes \mathbf{n}\right)=\mathbf{\rho }:\mathbf{\nu }$$

$$V=\sum _e{V}_e=\sum _e\cos {\chi }_e\frac{A_e{l}_e}{N_{\dim }}=\sum _e{\mathbf{\rho }}_e:{\mathbf{\nu }}_e\frac{A_e{l}_e}{N_{\dim }}$$

$$\phi$$

$$\mathbf{\Delta }$$

$$e_N$$

$$s_N$$

$$\delta W$$

$${\mathbf{A}}_{\dots i\dots j\dots }^{T_{ij}}={\mathbf{A}}_{\dots j\dots i\dots }$$

$${\mathcal{I}}^{\mathrm{vol}}$$

$$\mathbf{T}$$

$$e_N$$

$$\delta W$$

$$=Al{E}_0\left(\left[\left(\mathbf{\nu }\cdot {\mathbf{\rho }}^T\right):\mathbf{\epsilon }\right]\left[\left(\mathbf{\nu }\cdot {\mathbf{\rho }}^T\right):\delta \mathbf{\epsilon }\right]+\alpha \left[\left(\mathbf{T}\cdot {\mathbf{\rho }}^T\right):\mathbf{\epsilon }\right]\cdot \left[\left(\mathbf{T}\cdot {\mathbf{\rho }}^T\right):\delta \mathbf{\epsilon }\right]\right)$$

$$=Al{E}_0\left[\mathbf{\epsilon }:{\left(\mathbf{\rho }\cdot \mathbf{\nu }\otimes \mathbf{\nu }\cdot {\mathbf{\rho }}^T\right)}^{T_{12}}:\delta \mathbf{\epsilon }+\alpha \mathbf{\epsilon }:\left(\mathbf{\rho }\cdot {\mathbf{T}}^{T_{13}}\cdot \mathbf{T}\cdot {\mathbf{\rho }}^T\right):\delta \mathbf{\epsilon }\right]$$

$$=Al{E}_0\mathbf{\epsilon }:\left({\left(\mathbf{\rho }\cdot \mathbf{\nu }\otimes \mathbf{\nu }\cdot {\mathbf{\rho }}^T\right)}^{T_{12}}+\alpha \mathbf{\rho }\cdot {\mathbf{T}}^{T_{13}}\cdot \mathbf{T}\cdot {\mathbf{\rho }}^T\right):\delta \mathbf{\epsilon }$$

$$=Al{E}_0\mathbf{\epsilon }:\left(\mathcal{N}+\alpha \mathcal{T}\right):\delta \mathbf{\epsilon }$$

$$\mathcal{N}$$

$$\delta W^{\mathrm{dis}}=\sum _e\delta W_e=\sum _e{A}_e{l}_e{E}_0\mathbf{\epsilon }:\left({\mathcal{N}}_e+\alpha {\mathcal{T}}_e\right):\delta \mathbf{\epsilon }$$

$$\delta W^{\mathrm{con}}=V\mathbf{\sigma }:\delta \mathbf{\epsilon }=V\mathbf{\epsilon }:\mathcal{D}:\delta \mathbf{\epsilon }$$

$$\delta W^{\mathrm{dis}}=\delta W^{\mathrm{con}}$$

$$\mathcal{D}={⟨\frac{1}{V}\sum _e{A}_e{l}_e{E}_0\left({\mathcal{N}}_e+\alpha {\mathcal{T}}_e\right)⟩}^{\mathrm{SYM}}$$

$$\mathrm{major}\mathrm{symmetry}:{\mathcal{D}}_{ijkl}$$

$${⟨\bullet ⟩}^{\mathrm{SYM}}=\frac{\bullet +{\bullet }^{T_{34}}}{2}$$

$$\mathcal{D}=\frac{E_0}{V}{⟨\mathrm{E}\left[\mathcal{N}\right]+\alpha \mathrm{E}\left[\mathcal{T}\right]⟩}^{\mathrm{SYM}}\sum _e{A}_e{l}_e$$

$$\mathrm{E}\left[\bullet (\mathbf{x})\right]=\underset{-\infty }{\overset{\infty }{\int }}\cdots \underset{-\infty }{\overset{\infty }{\int }}\bullet (\mathbf{x})f_{\mathbf{X}}(\mathbf{x})\text{mbox}d\mathbf{x}$$

Full text

††journal: International Journal of Solids and Structures

1 Introduction

Discrete modeling is a well established technique in mechanics. It allows the explanation or prediction of the complex behavior of heterogeneous, cohesive or granular materials. The main advantages it offers are the straightforward representation of material random heterogeneity, the simple formulation of constitutive equations in vectorial form and also the direct consideration of discrete and oriented cracks. In contrast, the elastic behavior of these models still poses open challenges. Besides the minor issue of the inevitable boundary layer with different elastic properties [11], the most serious problem lies in the inability of such models to exhibit Poisson’s ratios greater than 1/3 for plane stress simplification and 1/4 for plane strain and three dimensional models ([5, 19, 11] or see Eq. 53).

Recently, four remedies providing the full range of Poisson’s ratio in discrete systems were presented. The first one [2, 3] introduces artificial auxiliary stresses within an iterative loop to achieve an elastically homogeneous system with an arbitrary Poisson’s ratio. The other three methods are similar as all of them are based on estimation of tensorial stresses or strains. The tensorial quantity is always computed nonlocally, in some neighborhood of the contact or body. These methods take into account the lateral stresses (confinement effect). The second remedy [9] proposes constitutive model as a function of the volumetric and deviatoric strain split. The third remedy [8] adds into the standard vectorial constitutive model terms accounting for the lateral stress. The fourth remedy [24] modifies distance between particles by integration tensorial strain over the body and evaluating body deformation. The stress oscillations caused by the heterogeneity of the material are unfortunately full or partially smeared out. Therefore, such models do not seem to be convenient for studying elastic behavior of highly heterogeneous structures at the mesoscale.

This paper is motivated by the author’s long belief that the Poisson’s ratio of discrete systems can be increased by changing the model geometry. Most of the papers published in this field use contact faces between discrete model units perpendicular to contact vectors [14, 6]. Examples of models with skewed normals are mostly from the field of granular materials when non-spherical particle shapes are used [25, 13, 17] but also static homogeneous models can be found [26]. The assumption of perpendicularity is abandoned here allowing model of completely arbitrary geometry to be constructed. Poisson’s ratio is then analyzed using strong assumptions about rotations and translations (Voigt’s hypothesis) in the model adopted according to [18]. It is proven here that abandoning the perpendicularity leads only to shrinking of the interval of achievable Poisson’s ratios.

The studied discrete system fills space continuously (without gaps or overlapping) with rigid bodies that possess translational ($\bm{u}$) and rotational ($\bm{\varphi}$) degrees of freedom. It is assumed that the system is isotropic – arbitrary rotation of the domain does not change its geometrical properties in the statistical sense. The rigid bodies interact via contacts defined at their boundaries. The normal and tangential displacement discontinuities $\bm{\Delta}$ at these boundaries are dictated by rigid body kinematics and give rise to normal and tangential forces linearly dependent on the corresponding component of $\bm{\Delta}$. A critical parameter governing the macroscopic Poisson’s ratio is the ratio between tangential and normal contact stiffness, hereinafter denoted by $\alpha$. The parameter $\alpha$ must be non-negative, otherwise the system would exhibit negative stiffness and instability. One can find several examples of these models in literature [23, 12, 16, 1].

The analytical derivation utilizes the equivalence of virtual work arising in the discrete system and continuum when they are subjected to equal straining. The Boltzmann continuum is used, and therefore the stress tensor must be symmetric (Boltzmann axiom). However, the discrete system yields a non-symmetric stress tensor, as it is a discrete instance of polar (Cosserat) continua instead [22]. The virtual work equivalence is therefore accomplished with the help of the symmetrization of the tensor of elastic constants from the discrete model.

2 Normal and contact vector, volume

The domain is divided into rigid bodies, each of which has degrees of freedom associated with the translations and rotations of its governing node, $\bm{x}_{a}$. The contact between two nodes $\bm{x}_{a}$ and $\bm{x}_{b}$ is provided by a mechanical element with contact area $A$, length $l=\|\bm{x}_{b}-\bm{x}_{a}\|$, unit normal vector ${\bm{n}}$ and contact vector ${\bm{t}}=(\bm{x}_{b}-\bm{x}_{a})/l$. The situation is depicted in Fig. 1a in two dimensions.

We assume that the system has no directional bias, and that therefore all normal directions share the same probability of occurrence. The vector $\bm{n}$ is here defined in the Cartesian coordinate system by two angles, $\xi$ and $\zeta$

[TABLE]

In two dimensions (2D), $\xi$ represents the angle between the $x$ axis and the normal, and uniformly covers the solid angle. In three dimensions (3D), $\zeta$ is the angle between the $z$ axis and the normal, and $\xi$ is the rotation of ${\bm{n}}$ around the $z$ axis - see Fig. 1b. For 3D isotropic systems, $\xi$ must also be uniform over the interval from 0 to $2\pi$ and $\zeta$ has the following probability density function

[TABLE]

The second fundamental vector governing the behavior of the contact is the contact vector ${\bm{t}}$. It is defined relative to the normal vector ${\bm{n}}$ by angles $\chi$ and $\theta$ - see Fig. 1b. The requirement of isotropicity dictates that (i) in 2D $\chi$ must have probability density function symmetric around zero and (ii) in 3D $\theta$ must be uniformly distributed over the 0–2$\pi$ interval.

[TABLE]

The probability distribution $f_{\chi}$ can be arbitrary (but symmetric around zero in 2D). For the sake of simplicity, it will now be assumed that the maximum angle between ${\bm{n}}$ and ${\bm{t}}$ is $\gamma\in\left[0,\,\pi\right]$ and that all directions within this range are equally probable.

[TABLE]

This assumption will later be removed in Sec. 6.

Let us define a second order tensor (rotation matrix) that provides the following relation between ${\bm{n}}$ and ${\bm{t}}$

[TABLE]

In two dimensions, this tensor is the standard rotation matrix

[TABLE]

A more complex situation is in 3D. One can imagine the construction of ${\bm{n}}$ by taking the vector \left(\begin{array}[]{ccc}0&0&1\end{array}\right), rotating it along the $y$ axis by angle $\zeta$ and then along the $z$ axis by angle $\xi$ (Fig. 1b). In the same way, the construction of ${\bm{t}}$ is achieved via four successive rotations along axes $y$, $z$, $y$ and $z$ by angles $\chi$, $\theta$, $\zeta$ and $\xi$, respectively.

[TABLE]

The rotation matrix from Eq. (5) is therefore

[TABLE]

The cosine of angle $\chi$ between ${\bm{n}}$ and ${\bm{t}}$ can be calculated using Eq. (5)

[TABLE]

where the second order tensor ${\bm{\nu}}$ is defined according to Kuhl et al. [18] as ${\bm{\nu}}={\bm{n}}\otimes{\bm{n}}$. Based on the assumption of no gaps or overlapping between the rigid bodies of the model, the volume of the domain is a summation over the volumes of individual mechanical elements

[TABLE]

where the number of dimensions ${N_{\mathrm{dim}}}$ is 2 or 3. Note that the volume of an individual element is negative if $|\chi|>\pi/2$.

3 Equivalence of virtual work

The fundamental assumption about system degrees of freedom is taken from [18]. It is known as Voigt’s hypothesis and widely used as a homogenization technique (e.g. [21]). It is assumed that when a discrete system is subjected to constant strain $\bm{\varepsilon}$, all the rotations are zero and differences in translations are dictated by differences in position

[TABLE]

The displacement jump at the contact between cells $a$ and $b$ is, based on the previous assumption, given by the rigid body kinematics of bodies without rotation

[TABLE]

where $l$ and $\bm{t}$ are respectively the length and contact vector belonging to the element connecting bodies $a$ and $b$. The normal and shear strain directly follow

[TABLE]

The stresses read

[TABLE]

where $E_{0}$ is the normal stiffness coefficient and $\alpha$ is the tangential/normal stiffness ratio. Both of the material parameters are considered constant throughout the whole domain.

The virtual work done by a single element is obtained via integration of the product of the stress and the displacement jump (both constant) over the contact face

[TABLE]

The integration is simple because the assumed zero rotations imply constant stresses and the displacement jump over the whole contact face. The total virtual work in the discrete system is the summation of these individual contributions.

To simplify the notation, we introduce the transpose operation $T_{ij}$ on an arbitrary tensor $\bm{A}$ of sufficient order by swapping indices $i$ and $j$.

[TABLE]

Let us now define two additional tensors according to Kuhl et al. [18]: the fourth order tensor ${\boldsymbol{\mathscr{I}}^{\mathrm{vol}}}$ and the third order tensor ${\bm{T}}$.

[TABLE]

where $\bm{1}$ is the identity matrix of size $N_{\mathrm{dim}}$. Note that ${\bm{T}}$ is different from the definition in [18, 11] because the symmetry implied by equality ${\bm{t}}={\bm{n}}$ is no longer present. The transposition $T_{13}$ means that dimensions 1 and 3 are swapped.

With the previously defined tensors ${\bm{\nu}}$ and ${\bm{T}}$, Eq. (17) can be rewritten as

[TABLE]

using the transposition $T$ of the second order tensor swapping its two dimensions $T=T_{12}$. The virtual work of a single element (Eq. (19)) can be rewritten as well.

[TABLE]

where

[TABLE]

The total virtual work of the discrete assembly is

[TABLE]

The virtual work of an equally strained elastic isotropic Boltzmann continuum occupying the same volume $V$ is

[TABLE]

with the constitutive equation $\bm{\sigma}={\bm{\mathcal{D}}}:\bm{\varepsilon}$ where ${\bm{\mathcal{D}}}$ is fourth order tensor of elastic constants.

The equivalence of the discrete and continuous system implies the equality of virtual works (Hill-Mandel condition [15])

[TABLE]

Substituting Eqs. (26) and (27) into Eq. (28), the following expression for the tensor of elastic constants is derived

[TABLE]

The symmetrization is needed because the tensors ${\boldsymbol{\mathscr{N}}}$ and ${\boldsymbol{\mathscr{T}}}$ do not possess the symmetries required for the Boltzmann continuum, which are major symmetry (derived from the equivalence of mixed derivatives of elastic potential) and minor symmetry (derived from symmetry of stress and strain tensors $\sigma_{ij}=\sigma_{ji}$, $\varepsilon_{ij}=\varepsilon_{ji}$).

[TABLE]

Because of the non-symmetric stress tensor in the discrete system, the minor symmetry is violated. The symmetric part can be easily obtained using transposition $T_{34}$.

[TABLE]

Thanks to the assumed statistical independence between the normal and contact vector and the elemental area and length, the summation in Eq. (29) can be broken into the following expression

[TABLE]

where $\mathrm{E}\left[\bullet(\bm{x})\right]$ is the mean value of function $\bullet$, which is dependent on vector $\bm{x}$ with the distribution function $f_{\bm{X}}(\bm{x})$

[TABLE]

Substituting $V$ from Eq. (14) and utilizing the statistical independence again, one obtains

[TABLE]

4 Calculation of expectations

4.1 Two dimensional case

In two dimensions, the rotation matrix depends only on angle $\chi$, while the normal ${\bm{n}}$ depends only on angle $\xi$. The calculation of the mean value can be separated into two steps. Let us first calculate all the quantities dependent solely on ${\bm{n}}$.

[TABLE]

where the fourth order tensor ${\boldsymbol{\mathscr{I}}}=\mathscr{I}_{ijkl}=(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})/2$ delta is employed; $\delta_{ij}\equiv\bm{1}$ is the the Kronecker delta.

In the second step, theses quantities are used in the calculation of the mean values of terms involving both ${\bm{\rho}}$ and ${\bm{n}}$.

[TABLE]

Only the symmetric parts of these expectations are needed. The following identities, which are valid for both the 2D and the 3D model, are derived from Eq. (31)

[TABLE]

and used to obtain the final symmetric expectations of tensors.

[TABLE]

4.2 Three dimensional case

In three dimensions, integration is substantially more complex. It is performed over four independent angles and cannot be separated since rotation matrix ${\bm{\rho}}$ depends on all four angles. Calculation by hand is extremely tedious; it was performed by computer instead. The following three integrations were delivered with the help of the Python library for symbolic mathematics, SymPy [20].

[TABLE]

Using identity (41), the symmetric part yields

[TABLE]

5 Relation between the elastic parameters of discrete system and continuum

The mechanical behavior of a linearly elastic isotropic solid is determined by two constants. Here we choose the elastic modulus ($E$) and Poisson’s ratio ($\nu$). The tensor of the elastic constants is expressed using these variables

[TABLE]

Equation (34), along with symmetrized expectations (38), (42) and (43) in two dimensions as well as (44), (47) and (48) in three dimensions, provides

[TABLE]

Equality between meso and macroscopic elastic tensors (50) and (49) requires equality between the respective scalar multipliers of tensors ${\boldsymbol{\mathscr{I}}^{\mathrm{vol}}}$ and ${\boldsymbol{\mathscr{I}}}$. One can modify these algebraic equations into relations between macroscopic parameters $E$ and $\nu$ and mesoscopic parameters $E_{0}$, $\alpha$ and $\gamma$.

[TABLE]

These equations are plotted in Fig. 2 for the range $\alpha\in[0,\,3]$.

Decreasing $\gamma$ towards zero must yield relations for a discrete system with ${\bm{n}}={\bm{t}}$.

[TABLE]

Indeed, the calculation of limits provides correct expressions (derived for example in [5, 19, 11] under the assumption of the perpendicularity of the contact vector and contact face). They are also identical to those from microplane theory [7].

What are the maximum and minimum Poisson’s ratios that can be achieved? One can differentiate the expression with respect to $\gamma$ and search for a stationary point (leaving out the degenerative case $\alpha=1$). In 2D, such an analysis reveals a local extreme at points $\gamma=0$ and $\gamma\approx 2.24670$ (solution of $2\gamma=\tan 2\gamma$). In 3D, the stationary points are $\gamma=0$, $\pi$ and $\approx 2.09440$ (exactly $\arccos(-0.5)$). Plotting the Poisson’s ratio with respect to the limit angle $\gamma$ (Fig. 3) shows that the maximum range of $\nu$ is obtained for $\gamma=0$, i.e. when the contact vector equals the normal vector. This is the classic solution stated in Eq. (53). Increasing $\gamma$ up to $\pi/2$ causes the interval of achievable Poisson’s ratios to shrink to zero. Then, the interval opens again with opposite signs; its width maximizes at $\gamma=2.24670$ (2D) or $\gamma=2.09440$ (3D). The interval of possible values of $\nu$ for these $\gamma$ reads: $\left[-0.122,\,0.098\right]$ for 2D plane stress, $\left[-0.139,\,0.089\right]$ for 2D plane strain and $\left[-0.091,\,0.034\right]$ for 3D, respectively. These maximum and minimum value of Poisson’s ratio occur always for $\alpha=0$ or $\alpha\rightarrow\infty$.

We have proven that under assumption (4), one cannot increase the Poisson’s ratio limits beyond that which is provided by the model with ${\bm{n}}={\bm{t}}$ in Eq. (53).

6 Arbitrary distribution $f_{\chi}(\chi)$

This section proves that the same conclusion unfortunately holds for the arbitrary distribution of angle $\chi$. The only restriction on the probability density of $\chi$ applied here comes from the isotropy requirement in 2D which demands that $f_{\chi}(\chi)$ be symmetric around zero. In 3D, $f_{\chi}(\chi)$ can be arbitrary.

Let us denote the following auxiliary integrals

[TABLE]

$\Omega_{\chi}$ is the domain of the probability distribution function $f_{\chi}(\chi)$ which is the interval $[-\pi,\,\pi]$ in 2D and $[0,\,\pi]$ in 3D. With the help of the Python symbolic mathematics library SymPy [20], the expectations needed in Eq. (34) were computed. Using identity $\int_{\Omega_{\chi}}f_{\chi}(\chi)\mbox{\,d}{\chi}=1$ resulting from the definition of probability density; and in 2D $\int_{\Omega_{\chi}}\sin(2\chi)f_{\chi}(\chi)\mbox{\,d}{\chi}=0$ derived from the symmetry of $f_{\chi}$, the expectations are

[TABLE]

Substituting these expressions into Eq. (34), requiring equality with Eq. (49) and solving for unknown $E$ and $\nu$, provides

[TABLE]

The Poisson’s ratio predicted by Eq. (57) is plotted in Fig. 4.

Let us focus again on the theoretical limits of Poisson’s ratio provided by Eq. (57). The meaningful values of $\alpha$ range from 0 to $\infty$. The integral $I_{2}$ is limited to an interval from -1 to 1 because these are the maximum and minimum values of the continuous function $\cos(2\chi)$, which in the integral is “weighted” by an arbitrary non-negative function with the unit integral over its domain. The only way $I_{2}=1$ can be obtained is when $f_{\chi}(\chi)$ is a discrete distribution with zero probability everywhere except points $\chi=\pi$, 0 or $-\pi$ (the last of which is possible only in 2D). All of these angles imply that ${\bm{t}}$ is parallel to ${\bm{n}}$. The second limit case of $I_{2}=-1$ is only possible when $f_{\chi}(\chi)$ is zero everywhere except points $\chi=\pm\pi/2$, i.e. when ${\bm{t}}$ is perpendicular to ${\bm{n}}$.

Differentiating Eq. (57) with respect to $I_{2}$ reveals that, except in the degenerative case of $\alpha=1$, there is no stationary point within the investigated domain for both 2D and 3D models. The extreme values therefore lie on the boundaries. One can see in Figure 4 that the maximum value of $\nu$ is reached by minimizing $\alpha$ and maximizing $I_{2}$ or by maximizing $\alpha$ and minimizing $I_{2}$; the opposite is true for the opposite goal of reaching the minimum Poisson’s ratio.

In 2D, the same Poisson’s ratio limits are obtained for both $I_{2}=1$ or $I_{2}=-1$, and these limits are the same as those of the model where ${\bm{n}}={\bm{t}}$. Indeed, the equality ${\bm{n}}={\bm{t}}$ implies $I_{2}=1$, and substituting that into Eq. (57) yields Eq. (53). In 3D, the case $I_{2}=1$ obtained for parallel ${\bm{n}}$ and ${\bm{t}}$ again yields Eq. (53) and the widest range of Poisson’s ratio. The case with perpendicular ${\bm{n}}$ and ${\bm{t}}$ (when $I_{2}=-1$), leads in 3D to a narrower interval as well as any other distribution of $\chi$.

We have proven that an isotropic discrete structure with an arbitrary relation between ${\bm{n}}$ and ${\bm{t}}$ cannot increase the Poisson’s ratio limits beyond the limits obtained by the model with ${\bm{n}}={\bm{t}}$ in Eq. (53).

7 Macroscopic elastic properties of actual discrete systems

Let us now observe the behavior of actual discrete systems and compare it to our analytical formulas. We divide a domain into (generally non-convex) polygons or polyhedrons; this tessellation defines the shapes of rigid bodies. Contact forces ${\bm{f}}=A(s_{N}{\bm{n}}+\bm{s}_{T})$ arise at the centroids of the edges (or faces) of polygons (or polyhedrons, respectively). We search for displacements and rotations fulfilling linear and angular momentum balance equations. These equations, for a single discrete body without external load, read

[TABLE]

where $e$ runs over all contacts with body neighbors, $\bm{r}$ is a position vector of the contact force ${\bm{f}}$ and ${\bm{\mathcal{E}}}$ is the Levi-Civita permutation symbol. In 3D, the contraction of ${\bm{\mathcal{E}}}$ from both sides gives the vector product $\bm{b}\cdot{\bm{\mathcal{E}}}\cdot\bm{a}=\bm{a}\times\bm{b}$.

We limit this section to 2D models, as it is expected that 3D models would yield similar results. Four different 2D tessellation types are considered, namely the Voronoi, randomized Voronoi, random and centered random tessellation.

7.1 Voronoi tessellation

The first tessellation type, referred to as Voronoi hereinafter, has parallel normal and contact vectors. It can be obtained via Voronoi tessellation, which is widely used in discrete modeling, or via Power tessellation, which is capable of taking into account the size of the inclusions (mineral aggregates in concrete) associated with rigid bodies [10].

The Voronoi model is created here by placing points randomly into a domain in sequence and accepting only those with a minimum distance from previously placed points that is greater than the length parameter $l_{\min}$. The sequential placement process is terminated after no point is accepted for a sufficiently large number of trials. The random points serve as nuclei for clipped Voronoi tessellation and their translations and rotations constitute model degrees of freedom.

The macroscopic elastic behavior of the Voronoi model is described by Eq. (53). Numerical verification is performed in [10] for both the 2D and 3D models, revealing increasing deviation from Eq. (53) with the increasing distance of parameter $\alpha$ from 1. The deviation is caused by the violation of assumption (15), which is exactly met only for $\alpha=1$. Indeed, assuming $\alpha=1$, applying Eqs. (16) and (17) that were derived from assumption (15), and using ${\bm{t}}={\bm{n}}$, which is valid for Voronoi tessellation, the equilibrium equations become

[TABLE]

The sum over the contacts is transformed into integration over the enclosed surface $\Gamma$ of the rigid body. The first integral is the zero vector, while the second integral is the identity matrix multiplied by the rigid body volume. This is derived via component-wise integration with the help of the divergence theorem and unit standard Cartesian basis vector $\bm{j}$.

[TABLE]

Since the first integration is the zero vector, the right-hand side of Eq. (61) is zero and Eq. (59) is exactly satisfied. Substituting the result of the second integration, the right-hand side of Eq. (64) becomes $VE_{0}{\bm{\mathcal{E}}}:\bm{\varepsilon}$, which is always zero thanks to the symmetry of the strain tensor and the antisymmetry of the Levi-Civita symbol, and Eq. (60) is exactly satisfied as well.

We have shown that for Voronoi tessellation with $\alpha=1$, assumption (15) holds. However, for different $\alpha$ or non-parallel ${\bm{n}}$ and ${\bm{t}}$, assumption (15) is incorrect. The actual system is more compliant and has a higher Poisson’s ratio than that predicted by assumption (15) due to the removed constraint on rotations and translations.

7.2 Randomized Voronoi tessellation

Randomized Voronoi tessellation is generated by modifying the Voronoi model. We simply move each Voronoi vertex in a random direction by a random distance. The normal vectors, $\bm{n}$, are randomly rotated by such random movements, while the contact vectors, $\bm{t}$, remain intact. The random distance is generated from the uniform distribution in the interval $(0,\,k)$, where $k$ is half of the minimum distance to the closest vertex. This upper limit on $k$ is introduced to prevent overlapping of the rigid bodies. The resulting tessellation has generally non-convex body shapes and non-parallel $\bm{n}$ and $\bm{t}$.

By generating 50 discrete structures of size $150l_{\min}\times 150l_{\min}$, the distribution function of angle $\chi$ is obtained. Numerical integration of Eq. (54) then provides constants $I_{1}\approx 0.97724$ and $I_{2}\approx 0.91372$. An example of randomized Voronoi tessellation as well as the probability distribution function of $\chi$ is shown in Figs. 5 and 6.

7.3 Random tessellation

The randomization of the Voronoi structure defined in section 7.2 is still rather limited and strongly resembles the original Voronoi model. When searching for a tessellation with more deviation from the parallelism of ${\bm{t}}$ and ${\bm{n}}$, the following process was found to be quite effective. The random model is based on two sets of nodes: basic nodes, which are generated by the same sequential process as described above with the minimum distance $l_{\min}$, and auxiliary nodes (vertices), which are created independently in the same way with the minimum distance $l_{\min}/2$. Delaunay triangulation is performed on the vertices, and triangles containing a basic node are directly assigned to it. Then, an iterative loop is performed over all unassigned triangles. Whenever a neighboring triangle already belonging to a basic node is found, the unassigned triangle is assigned to the same basic node. The rigid body associated with a given basic node is then a union of all the triangles assigned to that node. The obtained shapes are highly non-convex and typically several contact faces are created between two neighboring bodies.

Again, the distribution function of angle $\chi$ is estimated from 50 discrete structures of the size $150l_{\min}\times 150l_{\min}$. This is shown along with an example of random tessellation in Figs. 5 and 6. Numerical integration of Eq. (54) provides $I_{1}\approx 0.73516$ and $I_{2}\approx 0.28884$.

7.4 Centered random tessellation

The last tessellation type has bodies which are identical to those in the random model. The difference is that for centered random tessellation, a final step is performed in which the governing nodes bearing the rigid body degrees of freedom are moved into the centroids of the generated bodies. The situation is depicted in Fig. 5, where the hollow circles (initial tessellation nuclei) are replaced by solid circles (centroids). Though the shape of the bodies is unchanged, the geometrical characteristics of the tessellation changes because the vectors ${\bm{t}}$ and lengths $l$ are updated. The statistical evaluation of angle $\chi$ on 50 discrete structures with a size of $150l_{\min}\times 150l_{\min}$ provides the probability distribution function shown in Fig. 6 that, when numerically integrated in Eq. (54), provides $I_{1}\approx 0.78830$ and $I_{2}\approx 0.38688$.

One would assume that since the rigid body shapes are equal for the random and centered random model, these models would behave equally. We will now prove the opposite by contradiction.

Let us assume two structures with equal shapes and connectivity of the rigid bodies but two different sets of governing nodes (with coordinates $\bm{x}_{i}$ and $\hat{\bm{x}}_{i}$, respectively) which bear their degrees of freedom. The normals ${\bm{n}}$ and contact areas $A$ remain the same because these are dictated solely by rigid body shape. However, the contact lengths and contact vectors are different, $l\neq\hat{l}$ and ${\bm{t}}\neq\hat{{\bm{t}}}$. The starting point of the proof by contradiction is that under equal load the movements of the bodies are equal as well, and therefore the displacement jumps at the boundaries are identical, $\bm{\Delta}=\hat{\bm{\Delta}}$. According to the first parts of Eqs. (17) combined with Eq. (18), the stresses differ by factor $l/\hat{l}$

[TABLE]

The same factor also holds for contact forces as these are only stresses multiplied by areas: $\hat{{\bm{f}}}=l/\hat{l}{\bm{f}}$

The first set of forces (${\bm{f}}$) satisfies equilibrium equations (59) and (60) because we assume it originates from the true solution with the first set of nuclei. We will show now that the second set ($\hat{{\bm{f}}}$) violates equilibrium equations because the factor $l/\hat{l}$ differs for each element

[TABLE]

The assumption of the equality of rigid body movements for different governing nodes leads to a contradiction. Therefore, the movements are generally different and the macroscopic elastic properties of the random and centered random model differ as well.

7.5 Comparison of analytical formulas to actual characteristics

A two dimensional model with a size of $150l_{\min}\times 150l_{\min}$ was generated for all four tessellation types. The structure was loaded by the prescribed translations and rotations along the whole boundary, $\Gamma$: $u^{\Gamma}_{1}=px_{1}$, $u^{\Gamma}_{2}=qx_{2}$ and $\varphi^{\Gamma}=0$. The resulting macroscopic strain components are $\varepsilon_{11}=p$, $\varepsilon_{22}=q$ and $\varepsilon_{12}=\varepsilon_{21}=0$, respectively. Alternatively, the strain tensor can be obtained via the linear regression of location-translation dependence, see [10]. The stress tensor was estimated from all inner contacts with a distance from the boundary greater than $3l_{\min}$ [4]

[TABLE]

where $\bm{c}$ is the centroid of the contact face and the symmetrization of the second order tensor reads $\langle\bullet\rangle^{\mathrm{SYM}}=(\bullet+\bullet^{T})/2$.

Since the diagonal strain term $\varepsilon_{12}$ is zero, the macroscopic elastic properties are easily evaluated

[TABLE]

The model was assembled and solved for each tessellation type and various $\alpha$ ratios. The macroscopic characteristics obtained by Eqs. (69) and (70) are plotted in Fig. 7 along with the derived analytical predictions (57) and (58). There is reasonable correspondence for Voronoi and random Voronoi tessellation when $\alpha>0.3$. In all the other cases the predictions severely departed from numerical solution due to the unfulfilled fundamental assumption (15). However, the elastic tensor evaluated using Eq. (29) on the generated structure correspond in all the cases with one from Eq. (34) when analytical expectations (56) and (56) are used.

Though the correspondence is rather weak, there is a clear trend seen in both the analytical and the numerical results. The Poisson’s ratio shrinks into a narrower interval as the integral $I_{2}$ grows due to the loss of parallelism between ${\bm{n}}$ and ${\bm{t}}$; simultaneously, the elastic modulus increases. We consider this result to be a verification of our conclusion that Poisson’s ratio limits are maximized for Voronoi (or Power) tessellation. Any departure from perpendicularity between ${\bm{t}}$ and the face plane causes the narrowing of these limits.

8 Conclusions

Analytical estimations are derived for the macroscopic elastic behavior of (i) discrete isotropic assemblies with (ii) normal and shear force linearly dependent on normal and shear displacement discontinuity and (iii) a continuously filled domain. The derivation takes advantage of a strong assumption about the rotations (assumed zero everywhere) and translations (assumed uniformly distributed over domain) of discrete bodies. Comparison with numerical results reveals that the assumption is only reasonable for a limited range of model types, though the overall trend emerging from the derived equations is confirmed.

It is proven that the widest limits of Poisson’s ratio are obtained when contact and normal vectors are parallel, such as when Voronoi or Power tessellation is used to generate rigid body shapes.

• 2.

Any other model geometry (where normal and contact vectors are not strictly parallel) shrinks the interval of achievable Poisson’s ratio. It is not possible to extend the limits via geometrical manipulations. These limits are $\nu\in\left[-1,\,\nicefrac{{1}}{{3}}\right]$ for 2D plane stress, $\nu\in\left(-\infty,\,\nicefrac{{1}}{{4}}\right]$ for 2D plane strain and $\nu\in\left[-1,\,\nicefrac{{1}}{{4}}\right]$ for 3D models, respectively.

• 3.

The mechanical parameters $E_{0}$ and $\alpha$ are considered constant throughout the whole volume. Also, no overlapping or gaps between rigid bodies are allowed. The abandonment of these conditions seems to be one of the possible ways to proceed further in searching for other methods of extending the Poisson’s ratio interval of discrete models.

• 4.

It is shown that the position of the governing node affects model results, even though the actual rigid bodies are identical. It is questionable whether the body centroids, Voronoi sites (generators of the Voronoi diagram) or other points should be used as the governing nodes. It would be advantageous to find a model modification that would remove this dependence.

The discrete models are often used for analyzing inelastic phenomena (mostly fracture). The connection between the angle formed by the normal and contact vectors and the local and global inelastic behavior of the model is an interesting open topic deserving further investigation.

Acknowledgement

Financial support provided by the Czech Science Foundation under project No. GA19-12197S is gratefully acknowledged.