Section author: Václav Šmilauer <eu@doxos.eu>

Cross-anisotropy model with discrete elements

Summary

The task was to prepare theoretical foundation for modelling rock with cross-anisotropic elasticity (schistosity) by using modified DEM formulation.

The following paragraphs present derivation of local stiffnesses of geometrically isotropic linear-elastic lattice leading to globally cross-anisotropic behavior. Lattice elements are assumed to have only normal and shear stiffnesses, which are made dependent on the element orientation. The derivation is based on the microplane theory regarded as continuous limit of a dense lattice, and the Voigt hypothesis of uniform displacement of lattice nodes. Microplane moduli are made orientation-dependent in the same way as for the lattice, via a function which allows us to carry out analytical integration of the stiffness tensor.

Geometrically isotropic lattices can be found within the Discrete Element Method (DEM) with spherical particles, where the mechanical structure is a transient lattice. Due to the local Voigt hypothesis, only packings with relatively high number of contacts per particle will fulfill our assumptions; this can be achieved by introducing distant contacts between particles which are not geometrically touching, a common practice for modeling cohesive materials.

Element tests on a periodic lattice structure were done, and simulation results were compared with stiffness tensor based on which the lattice stiffnesses were assigned.

Behavior of circular cavity under radial pressure was examined using this model; this is shown at the very end.

Note

The results presented here can be reproduced by downloading the following scripts:

• pbc-prepare-sample.py to prepare monomodal periodic packing with porosity equal to 0.5; alternatively, download pbc-spheres_N=2556,r=0.03,rRelFuzz=0.txt which is the result of running this script

• pbc-deform-sample.py to run a single element test on the whole periodic lattice

• pbc-deform-sample.table.xls as input to batch for running element tests on a range of the \(I_R\) parameters by using:

  woo-batch pbc-deform-sample.table.xls pbc-deform-sample.py
  

• pbc-process-logs.py to extract resulting tables (such as those presented in Contact radius I_R=1.8) from log files, once the batch is completed.

woo.dem.Cp2_FrictMat_FrictPhys_CrossAnisotropy implements the cross-anisotropic stiffness computation, woo.utils.muStiffnessScaling implements computation of the \(\mu\) proportionality coefficient, and woo.dem.WeirdTriaxControl was used for stress-controlled element tests.

Acknowledgement

This work was done during post-doc position at Arbeitsbereich für Geotechnik und Tunnelbau at the University of Innsbruck; limited part of this work was presented at the EURO:TUN 2013 conference, but the full version including all derivations never made it to a proper publication.

Introduction

Cross-anisotropic (also transverse-isotropic) materials exhibit rotational symmetry in mechanical response about the cross-anisotropy axis \(\vec{a}\): the material behaves isotropically perpendicular to \(\vec{a}\). Such behavior is often found in layered materials. For the case where \(\vec{a}\) coincides with the global \(z\)-axis, this symmetry is captured by the stiffness tensor

\[\begin{split}\tens{C}= \begin{pmatrix} C_{11} & C_{12} & C_{13} & 0 & 0 & 0 \\ C_{12} & C_{11} & C_{13} & 0 & 0 & 0 \\ C_{13} & C_{13} & C_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{C_{11}-C_{12}}{2} \end{pmatrix}\end{split}\]

with five independent components.

Lattice

The objective of this paper is to model a linear-elastic cross-anisotropic medium using a discrete model. Under discrete models we understand the Discrete Element Method (DEM) with spherical particles, though we are only interested in the static solution in the linear domain. Therefore, an arrangement of DEM particles with contacts in-between them can be seen as a lattice structure, where the contacts (lattice elements) describe interactions with neighboring particles (lattice nodes). We will thus speak of lattice elements, and only turn to the DEM terminology when required.

Constructing the lattice from an arrangement of spherical particles is shown in this figure; the procedure is of geometrical nature, where spheres with positions \(\vec{x}_i\) and radii \(r_i\) are in contact if

(37)\[\lvert\vec{x}_1-\vec{x}_2\rvert≤I_R(r_1+r_2)\]

where \(I_R\) is contact radius [SJS10], dimensionless parameter determining density of the lattice. I_R is equal to 1 if only geometrically adjacent spheres are to be considered in contact, while greater values will allow “contact” between spheres which some distance between them.

Fig. 46 Dependence of lattice density on contact radius \(I_R\): (a) sphere packing for constructing the lattice, (b) lattice with \(I_R=1\), (b) lattice with bigger \(I_R\).

The cross-anisotropic behavior is described by the stiffness tensor. With continuum models, the constitutive law is written in terms of the stiffness tensor itself; however, discrete models are not homogeneous and a distinction between local (microscopic) and global (macroscopic) levels of description is to be made. In the local sense, each lattice element is be characterized by stiffness values. Globally, though, the behavior is given by complicated interplay between numerous lattice elements.

Our task is therefore to investigate the relationship between local and global lattice behavior: providing a mathematical description of how the global behavior (expressed as a stiffness tensor) depends on local characteristics (stiffnesses of lattice elements) and how are they playing together (lattice geometry). Subsequently, we will inverse the previous result, i.e. for some desired stiffness tensor, we will determine local characteristics (stiffnesses) leading to the behavior described by that tensor.

To this end, we make use of the following assumptions:

• elements of the lattice have only normal and shear stiffnesses \(k_n\), \(k_t\);

• the lattice is geometrically isotropic: the orientation of the lattice elements is in average uniformly distributed;

• the displacement of lattice nodes does not deviate from the mean displacement; in another words, the lattice is deformed uniformly, as whole. This assumption is commonly referred to as Voigt hypothesis [LCYC97] and is an important restriction for the lattice behavior. Denser lattices (with greater \(I_R\)) fulfill this assumption by themselves better than loose lattices, as they have more contacts restricting deviation of individual nodes from the surrounding deformation. We therefore expect our results to better describe the real behavior in the case of dense lattices.

Microplane and lattice

The transition between discrete (lattice) [KDAddettaLR01] and continuous description will be done via the microplane theory [CJirasekB04]. This theory describes each material point as infinite number of microplanes oriented uniformly in all possible directions at that material point, each of them characterized by volumetric, deviatoric, normal and shear moduli (in our case, we only use the latter two, \(E_N\) and \(E_T\)). The cross-anisotropic nature is introduced by supposing dependency of those moduli on the respective microplane orientation such that symmetries of a cross-anisotropic medium are satisfied. The stiffness tensor is obtained by integration of the moduli over all microplanes.

The lattice structure has only a finite number of nodes and a finite number of isotropically-oriented lattice elements in each node. The stiffness tensor is obtained by summation of stiffnesses over all elements. We can write the lattice stiffnesses \(k_n\), \(k_t\) as functions of some yet unknown moduli \(E_N'\), \(E_T'\) (which depend on the orientation in the same way as \(E_N\), \(E_T\)) and let the lattice density grow without bounds. After limit transition, we obtain the stiffness tensor of the infinitely dense lattice by integration of the unknown moduli over all lattice elements.

By imposing equality of microplane and lattice stiffness tensors, we obtain the values of the unknown lattice moduli \(E_N'\), \(E_T'\) (proportional to the microplane moduli via a yet unknown factor \(\mu\)); using those moduli when constructing the lattice ensures that the stiffness tensor of the discrete lattice will be equal to the stiffness tensor of the microplane model. Consequently, for a given stiffness tensor, we can compute lattice moduli which will lead to the response characterized by that stiffness tensor.

The function expressing how moduli depend on the orientation is the same for the microplane and the lattice model; in our case, we chose a function with only 4 independent parameters. Thus, the resulting stiffness tensor is cross-anisotropic, but has only four linearly independent components.

Tests

The stiffness tensor of a lattice is obtained from element tests (in our case, uniaxial unconfined compression and simple shear tests) on a periodic lattice. It is subsequently compared with the prescribed values, and the dependence between accuracy and lattice density is shown.

Notation

\(E_N\), \(E_T\)	microplane normal and tangential moduli
\(E_N'\), \(E_T'\)	lattice normal and tangential moduli
\(k_n\), \(k_t\)	lattice normal and tangential stiffnesses
\(E_*^a\), \({E_*^a}'\)	in-plane microplane/lattice moduli
\(E_*^b\), \({E_*^b}'\)	out-of-plane microplane/lattice moduli
\(\tens{E}\)	stiffness tensor
\(C_{ijkl}\)	stiffness tensor, index notation
\(C_{mn}\)	stiffness matrix (in Voigt notation) element
\(\mathcal{E}\)	\((E_N^a, E_N^b,E_T^a,E_T^b)^T\)
\(\mathcal{E}'\)	\(({E_N^a}', {E_N^b}', {E_T^a}',{E_T^b}')^T\)
\(\mathcal{C}\)	\((C_{11}, C_{33}, C_{13}, C_{12})^T\)
\(\phi\in\langle 0,2\pi)\)	azimuth angle in spherical coordinates
\(\theta\in\langle 0,\pi)\)	inclination angle in spherical coordinates, from the pole
\({\rm s}\theta\), \({\rm s}\phi\)	\(\sin\phi\), \(\sin\theta\)
\({\rm c}\theta\), \({\rm c}\phi\)	\(\cos\phi\), \(\cos\theta\)

Cartesian coordinates of a point \(\theta\), \(\phi\) on sphere with radius \(r\) are given (see the image as

(38)\[\begin{split}\vec{n}=\begin{pmatrix}r {\rm s}\theta{\rm c}\phi \\ r {\rm s}\theta{\rm s}\phi \\ r{\rm c}\theta\end{pmatrix}\end{split}\]

with the transformation Jacobian equal to \(r^2 {\rm s}\theta\):

Fig. 47 Polar coordinate system; \(r\), \(\theta\), \(\phi\) are radius, inclination, azimuth respectively.

Microplane stiffness tensor

We consider the microplane model by [CJirasekB04] where all microplanes only have normal and shear moduli \(E_N\), \(E_T\) volumetric and deviatoric moduli \(E_V\), \(E_D\) are zero. Stiffness tensor is obtained by integration of moduli over all possible orientations of microplanes given by unit vector \(\vec{n}\). For our purposes, we will integrate in spherical coordinates with unit radius over angles \(\theta\), \(\phi\) (polar coordinates), \(z\)-axis being coincident with the cross-anisotropy axis; this lets us write microplane moduli as \(E_N=E_N(\theta)\) and \(E_T=E_T(\theta)\), independent of the azimuth \(\phi\).

(39)\[\begin{split}\tens{\mathcal{C}}&=\frac{3}{2\pi}\int\limits_{\Omega} E_N\tens{\mathcal{N}}\otimes\tens{\mathcal{N}} + E_T \tens{\mathcal{T}}^T\cdot\tens{\mathcal{T}}\d\Omega= \\ &=\frac{3}{2\pi}\int\limits_{\theta=0}^{\pi}\int\limits_{\phi=0}^{2\pi} \left(E_N\tens{\mathcal{N}}\otimes\tens{\mathcal{N}} + E_T \tens{\mathcal{T}}^T\cdot\tens{\mathcal{T}}\right){{\mathrm s}\theta\,}\d\phi\d\theta\end{split}\]

where \(\tens{\mathcal{N}}\), \(\tens{\mathcal{T}}\) are projection tensors; in index notation, they read

\[\begin{split}\tens{\mathcal{N}}\otimes\tens{\mathcal{N}}&=\tens{n}\otimes\tens{n}\otimes\tens{n}\otimes\tens{n}=n_i n_j n_k n_l \\ \tens{\mathcal{T}}&=\tens{n}\cdot\tens{\mathcal{I}}^{\rm sym}-\tens{n}\otimes\tens{n}\otimes\tens{n} \\ \tens{\mathcal{T}}^{T}&=\tens{\mathcal{I}}^{\rm sym}\cdot{n}-\tens{n}\otimes\tens{n}\otimes\tens{n} \\ T_{ijk}&=n_l\frac{1}{2}(\delta_{lj}\delta_{ik}+\delta_{lk}\delta_{ij})-n_i n_j n_k \\ T^T_{ijk}&=n_l\frac{1}{2}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})-n_i n_j n_k \\ \tens{\mathcal{T}}^T\cdot \tens{\mathcal{T}} &= T^T_{ijm} T_{mkl} = \\ &=\frac{n_j n_k}{4} \delta_{il}+\frac{n_j n_l}{4}\delta_{ik}+\frac{n_i n_k}{4}\delta_{jl}+\frac{n_i n_l}{4}\delta_{jk}-n_i n_j n_k n_l\end{split}\]

Stiffness tensor components are subsequently written as

(40)\[\begin{split}C_{ijkl}&=\frac{3}{2\pi}\int\limits_0^{2\pi}\int\limits_0^{\pi}\Big[E_N(\theta)(n_i n_j n_k n_l)+E_T(\theta)\Big(\frac{n_j n_k}{4} \delta_{il}+\frac{n_j n_l}{4}\delta_{ik}+ \\ &\quad+\frac{n_i n_k}{4}\delta_{jl}+\frac{n_i n_l}{4}\delta_{jk}-n_i n_j n_k n_l\Big)\Big]{{\mathrm s}\theta\,}\d\phi\d\theta\end{split}\]

Plugging \(\vec{n}\) in (38) into (40), and reversing integration order, we obtain

\[\begin{split}C_{1111}&=\frac{3}{2\pi}\int\limits_{\theta=0}^{\pi}{\int\limits_{\phi=0}^{2\pi}{\left[E_N {{\mathrm s}^4\theta\,} {{\mathrm c}^4\phi\,} + E_T({{\mathrm s}^2\theta\,}{{\mathrm s}^2\phi\,}-{{\mathrm s}^4\theta\,}{{\mathrm c}^4\phi\,})\right]{{\mathrm s}\theta\,}\,{\rm d}\theta}\,{\rm d}\phi}\,{\rm d}\theta = \\ &=\frac{3}{2\pi}\int\limits_{\theta=0}^{\pi}{E_N{{\mathrm s}^5\theta\,}\int\limits_{0}^{2\pi}{{{\mathrm c}^4\phi\,}\,{\rm d}\phi}+E_T{{\mathrm s}^3\theta\,}\int\limits_{0}^{2\pi}{{{\mathrm c}^2\phi\,}\,{\rm d}\phi}-E_T{{\mathrm s}^5\theta\,}\int\limits_{0}^{2\pi}{{{\mathrm c}^4\phi\,}\,{\rm d}\phi}}\,{\rm d}\theta\end{split}\]

As the inner integrals over \(\phi\) do not depend on \(\theta\), they can be evaluated directly, giving

(41)\[ C_{1111}=\frac{3}{2}\int\limits_{\theta=0}^{\pi}{\frac{3}{4}E_N{{\mathrm s}^5\theta\,}+ E_T{{\mathrm s}^3\theta\,}-\frac{3}{4}E_T{{\mathrm s}^5\theta\,}}\,{\rm d}\theta.\]

Similarly for other components,

(42)\[\begin{split} C_{2222}&=C_{1111} \\ C_{3333}&=\frac{3}{2}\int\limits_{0}^{\pi}{2E_N{{\mathrm c}^4\theta\,}{{\mathrm s}\theta\,}+2E_T{{\mathrm c}^2\theta\,}{{\mathrm s}\theta\,}-2E_T{{\mathrm c}^4\theta\,}{{\mathrm s}\theta\,}}\,{\rm d}\theta \\ C_{2233}&=\frac{3}{2}\int\limits_{0}^{\pi}{E_N{{\mathrm s}^3\theta\,}{{\mathrm c}^2\theta\,}-E_T{{\mathrm s}^3\theta\,}{{\mathrm c}^2\theta\,}}\,{\rm d}\theta \\ C_{1133}&=C_{2233} \\ C_{1122}&=\frac{3}{2}\int\limits_{0}^{\pi}{\frac{1}{4}E_N{{\mathrm s}^5\theta\,}-\frac{1}{4}E_T{{\mathrm s}^5\theta\,}}\,{\rm d}\theta \\ C_{2323}&=\frac{3}{2}\int\limits_{0}^{\pi}{E_N{{\mathrm s}^3\theta\,}{{\mathrm s}^2\theta\,}+E_T\left(\frac{{{\mathrm c}^2\theta\,}{{\mathrm s}\theta\,}}{4\pi}+\frac{{{\mathrm s}^3\theta\,}}{4}-{{\mathrm s}^3\theta\,}{{\mathrm c}^2\theta\,}\right)} \,{\rm d}\theta \\ C_{1313}&=C_{2323} \\ C_{1212}&=\frac{3}{2}\int\limits_{0}^{\pi}{\frac{1}{4}E_N{{\mathrm s}^5\theta\,}+E_T\left(\frac{{{\mathrm s}^3\theta\,}}{2}-\frac{{{\mathrm s}^5\theta\,}}{4}\right)}\,{\rm d}\theta= \\ &=\frac{1}{2}\left(C_{1111}-C_{1122}\right)\end{split}\]

Note that only 5 of the stiffness tensor components are independent, which is consistent with the symmetries of the cross-anisotropic material.

Isotropy

With \(\theta\)-independent \(E_N\) and \(E_N\), \(C_{ijkl}\) can be evaluated explicitly and assumes the values given in literature [KDAddettaLR01], [SJS10], i.e.

(43)\[\begin{split} C_{1111}=C_{2222}=C_{3333}&=\frac{1}{5}(6E_N+4E_T), \\ C_{2233}=C_{1133}=C_{1122}&=\frac{1}{5}(2E_N-2E_T), \\ C_{2323}=C_{1313}=C_{1212}&=\frac{1}{5}(2E_N+3E_T).\end{split}\]

Cross-anisotropy

\(E_N\) and \(E_T\) are functions of inclination \(\theta\in\langle0,\pi\rangle\). To further evaluate the stiffness tensors, we will assume that they can be written as functions interpolating between in-plane moduli \(E_*^a\) and out-of-plane-moduli \(E_*^b\), via some weight function \(w(\theta)\):

(44)\[ \begin{align}\begin{aligned}E_N(\theta)&=w(\theta)E_N^a+(1-w(\theta))E_N^b\\E_T(\theta)&=w(\theta)E_T^a+(1-w(\theta))E_T^b.\end{aligned}\end{align} \]

We require the function \(w\) to be

• interpolating between in-plane and out-of-plane values, i.e. \(w(0)=0\), \(w(\pi/2)=1\);

• continuous in \(\langle 0, \pi/2\rangle\);

• symmetrical around \(\pi/2\), i.e. \(\forall\theta: w(\theta)=w(\pi/2-\theta)\).

Elliptic interpolation

One candidate for \(w\), used in the rest of this paper, is the function

(45)\[ \begin{align}\begin{aligned}w(\theta)&={{\mathrm s}^2\theta\,},\\1-w(\theta)&={{\mathrm c}^2\theta\,}.\end{aligned}\end{align} \]

Fig. 48 Elliptic interpolation of stiffnesses.

The advantage is that the interpolation function has simple elliptical shape and that we are able to evaluate stiffness tensor components in a closed form. The disadvantage is that \(w\) contains no additional parameter, therefore only 4 stiffness tensor components will be independent; however, such an approximation can be in many cases accurate enough; we will choose the out-of-plane shear modulus \(C_{44}\) to be the dependent component, for it is difficult to measure experimentally [Bli11], but any other component could be selected instead.

Plugging (45) into (44), and then into (41) − (42), we obtain all stiffness tensor components as function of the lattice moduli,

(46)\[\begin{split}\begin{pmatrix}C_{1111} \\ C_{3333} \\ C_{2233} \\ C_{1122} \\ C_{2323} \end{pmatrix} = \begin{pmatrix}C_{11}\\C_{33}\\C_{13}\\C_{12}\\C_{44}\end{pmatrix}= \frac{1}{35} \begin{pmatrix} 36 & 6 & 20 & 8 \\ % E_1111 12 & 30 & 16 & 12 \\ % E_3333 8 & 6 & -8 & -6 \\ % E_2233 12 & 2 & -12 & -2 \\ % E_1122 8 & 6 & 13 & 8 \\ % E_2323 \end{pmatrix} \begin{pmatrix}E_N^a \\ E_N^b \\ E_T^a \\ E_T^b \end{pmatrix}.\end{split}\]

Writing out the independent components separately, we obtain

\[\begin{split}\begin{pmatrix}C_{11}\\C_{33}\\C_{13}\\C_{12}\end{pmatrix}&= \frac{1}{35} \begin{pmatrix} 36 & 6 & 20 & 8 \\ % E_1111 12 & 30 & 16 & 12 \\ % E_3333 8 & 6 & -8 & -6 \\ % E_2233 12 & 2 & -12 & -2 \\ % E_1122 \end{pmatrix} \begin{pmatrix}E_N^a \\ E_N^b \\ E_T^a \\ E_T^b \end{pmatrix}\end{split}\]

or in a compact form

(47)\[\mathcal{C}=\mat{A}\mathcal{E}\]

with the dependent component

(48)\[C_{44}=\frac{1}{35}\begin{pmatrix}8& 6& 13& 8\end{pmatrix}\mathcal{E}=\begin{pmatrix}\frac{1}{4}& \frac{1}{4}& -\frac{1}{2}& 0\end{pmatrix} \mathcal{C}.\]

Since the \(\mat{A}\) matrix is invertible, we can determine \(\mathcal{E}\) for given values of \(\mathcal{C}\) as

(49)\[\mathcal{E}=\mat{A}^{-1}\mathcal{C}.\]

Lattice stiffness tensor

We start with a lattice composed of elements characterized by normal and shear stiffnesses \(k_n\), \(k_t\). Since we want to compare lattice with the microplane model, we will write lattice element stiffnesses as functions of some lattice moduli \(E_N'\), \(E_T'\); we will further suppose that \(E_N'\), \(E_T'\) are orientation-dependent as in (44), with different in-plane and out-of-plane moduli \(\mathcal{E}'\).

Without loss of generality, we consider that the lattice is constructed from spherical packing and that intra-nodal stiffnesses are given as \(k_n=E_N'\pi r^2 /(2r')\), where \(\pi r^2\) is a fictitious contact area divided by the contact length \(2r'\); depending on the algorithm, the actual stiffness can be different up to multiplicative constant, so we will write

\begin{align*} k_n &=E_N'\frac{\hat\pi r^2}{2r'}, & k_t&=E_T'\frac{\hat\pi r^2}{2r'}. \end{align*}

In the case we mentioned, \(\hat\pi=\pi\) but we assume no particular value in the following.

Stiffness tensor of a lattice is, assuming Voigt hypothesis of uniform node displacements, written as [KDAddettaLR01] (eq. (35)):

(51)\[\tens{\mathcal{C}}=\frac{1}{V}\sum|l|^2\left[k_n \tens{\mathcal{N}}\otimes\tens{\mathcal{N}} +k_t\tens{\mathcal{T}}^T\cdot\tens{\mathcal{T}}\right],\]

where \(|l|=(2r')\) is the length of a lattice element. After limit transition, for some representative volume \(V\) and large number of lattice elements \(N\), the tensor is written as an integral over all orientations

(52)\[\tens{\mathcal{C}}=\frac{N}{4V\pi}\int_\Omega (2r')^2 \left( k_n \tens{\mathcal{N}}\otimes\tens{\mathcal{N}} +k_t\tens{\mathcal{T}}^T\cdot\tens{\mathcal{T}}\right) \d\Omega,\]

where \((2r')\) was kept inside the integral, as it is not a constant. Plugging in stiffnesses from (50) yields

(53)\[\tens{\mathcal{C}}=\frac{N}{2V\pi}\hat\pi \overline{r^2r'}\int_{\Omega}E_N'\tens{\mathcal{N}}\otimes\tens{\mathcal{N}}+E_T'\tens{\mathcal{T}}^T\cdot\tens{\mathcal{T}} \d\Omega.\]

where the \(r^2r'\) term appearing inside the integral was written in the front as an average, because it is on average orientation-independent.

Microplane-lattice moduli proportion

Microplane moduli \(E_N\), \(E_T\) and lattice moduli \(E_N'\), \(E_T'\) have the same dimensions and both depend on the orientation via (44); they are therefore related via an unknown dimensionless parameter \(\mu\)

(54)\[ \begin{align}\begin{aligned}E_N&=\mu E_N',\\E_T&=\mu E_T'.\end{aligned}\end{align} \]

We can plug this relation into (39) obtaining

(55)\[\tens{\mathcal{C}}=\frac{3}{2\pi}\mu\int_{\Omega} E_N'\tens{\mathcal{N}}\otimes\tens{\mathcal{N}} + E_T' \tens{\mathcal{T}}^T\cdot\tens{\mathcal{T}}\d\Omega.\]

By putting (53) and (55) equal, we solve

(56)\[\mu=\frac{E_N}{E_N'}=\frac{E_T}{E_T'}=\frac{1}{3}\frac{N\overline{r'r^2}\hat\pi}{V}\]

If we consider the special case of \(r\) and \(r'\) being constant, we can omit the average and write \(\mu\) in the form better revealing its geometrical meaning as

\[\mu=\frac{1}{6}\frac{N (\hat\pi r^2)(2r')}{V}=\frac{1}{6}\frac{N A l}{V}\]

with \(A=\hat\pi r^2\) and \(l\) being is the cross-section area and length of one lattice element; the fraction is therefore proportion between the stiff volume of lattice elements \(NAl\) to the overall volume \(V\).

Stiffness tensor from lattice moduli

Since stiffness tensor components \(C_{ij}\) are linear with respect to microplane moduli \(E_*\), we can write \(C_{ij}(E_*)=C_{ij}(\mu E_*')=\mu C_{ij}(E_*')\). Thus, for instance, (43) becomes

\[C_{1111}=\frac{1}{5}(6E_N+4E_T)=\mu\frac{1}{5}(6E_N'+4E_T')\]

and (47), (49) become respectively

(57)\[ \mathcal{C}=\mu\mat{A}\mathcal{E}'\]

and

(58)\[ \mathcal{E}'=\frac{1}{\mu}\mat{A}^{-1}\mathcal{C}.\]

Element test

Random dense (porosity equal to 0.5) periodic packing of spheres with equal radius is considered. The lattice structure is created by finding contacts between particles with varying contact radius \(I_R\) according to (37). Local lattice moduli are computed using (44) with the elliptic interpolation (45). Cross-anisotropy axis is aligned with the global \(z\)-axis. Given some input values of the \(\mathcal{C}\) vector (i.e. \(C_{11}\), \(C_{33}\), \(C_{13}\), \(C_{12}\)), we use the current packing geometry to determine \(\mu\) via (56) and obtain \(\mathcal{E}'\) via (58).

The goal is to compare stiffness tensor from 3 sources:

1. from prescribed \(\mathcal{C}\), additionally with the dependent component \(C_{44}\) from (48);

2. from the current lattice stiffnesses using (51);

3. from simulated lattice response during element tests (described below), with stiffnesses (60), (61), (62).

The test will be run for a range of contact radius \(I_R\); higher values of \(I_R\) lead to stabilization of the packing, better satisfying the Voigt hypothesis of uniform deformation, as mentioned in the introduction. It is therefore expected that lower values of \(I_R\) will overestimate the microplane/lattice stiffnesses with respect to the simulated response and that this difference will be smaller with growing \(I_R\). The minimum value of \(I_R\) is 1.2, since smaller values yielded unstable lattices. A rendering of the lattice for is shown in fig. Periodic lattice.

Fig. 49 Periodic lattice of the element test for \(I_R=1.4\), shown in the \(yz\) plane. The right picture varies color depending on stiffness, showing the preference for more compliant behavior (dark) along the \(z\) axis.

Evaluating the stiffness tensor

Homogeneous periodic lattice is loaded in order to obtain parameters of the cross-anisotropic material. For each axis, unconfined uniaxial compression is performed to obtain the normal modulus and Poisson’s ratios, and shear test is performed to obtain the shear modulus (<element test figure). In total, six tests (for each value of \(I_R\)) are run.

Within the 12 values (3 normal moduli, 3 shear moduli, 3 Poisson’s ratios, each measured twice) obtained, there are 5 symmetric couples (\(E_{xy}=E_{yx}\), \(G_{yz}=G_{zx}\), \(\nu_{xy}=\nu_{yx}\), \(\nu_{yz}=\nu_{xz}\), \(\nu_{zx}=\nu_{zy}\)) and 2 dependent values (e.g. \(\nu_{xz}=\nu_{zx}E_x/E_z\) and \(G_{xy}=E_x/2(1+\nu_{xy})\)). This results in five independent material constants of the cross-anisotropic material. The redundant values can be used to check correctness.

Fig. 50 Configuration of (a) unconfined uniaxial compression test in the \(y\)-direction and (b) shear test in the \(xy\) plane.

Normal moduli and Poisson’s ratios

To find normal moduli from simulations, we make use of the normal compliance relationship for an orthotropic material, of which cross-anisotropic material is a special case:

(59)\[\begin{split}\begin{pmatrix}\eps_{xx} \\ \eps_{yy} \\ \eps_{zz}\end{pmatrix}=\begin{pmatrix}\frac{1}{E_x} & -\frac{\nu_{yx}}{E_y} & -\frac{\nu_{zx}}{E_z} \\ -\frac{\nu_{xy}}{E_x} & \frac{1}{E_y} & -\frac{\nu_{zy}}{E_z} \\ -\frac{\nu_{xz}}{E_x} & -\frac{\nu_{yz}}{E_y} & \frac{1}{E_z} \end{pmatrix}\begin{pmatrix}T_{xx} \\ T_{yy} \\ T_{zz} \end{pmatrix}\end{split}\]

and prescribe extension \(\hat{\eps}\) (the hat signifies that the value is prescribed) in one direction with zero lateral stresses and evaluate the response; since stress cannot be prescribed directly in DEM, we adjust lateral strains until stress is smaller than some relative tolerance \(t_r\) times axial stress. Because the behavior is linear elastic, the way lateral strains are adjusted has no influence on the result. For instance, for \(x\) extension we prescribe

\[\begin{split}\eps_{xx}&=\hat{\eps}_{xx}=\hat\eps, \\ T_{yy}=T_{zz}&=0,\end{split}\]

and use measured response \(T_{xx}\), \(\eps_{yy}\), \(\eps_{zz}\) to compute

(60)\[\begin{split}E_x&=\frac{T_{xx}}{\hat{\eps}_{xx}} \\ \nu_{xy}&=-\eps_{yy}\frac{E_x}{T_{xx}}=-\frac{\eps_{yy}}{\hat{\eps}_{xx}} \\ \nu_{xz}&=-\frac{\eps_{zz}}{\hat{\eps}_{xx}}.\end{split}\]

Shear moduli

Shear test is purely deformation-controlled. The periodic cell is prescribed pure shear \(\eps_{ij}=\eps_{ji}=\hat\eps\) (zero elsewhere). Shear moduli are found from shear stiffness equations

(61)\[G_{ij}=\frac{T_{ij}}{2\hat\eps_{ij}}=\frac{T_{ij}}{\gamma_{ij}}.\]

---

Components of the stiffness tensor are found by inversion of orthotropic compliance matrix (59) using symmetries \(E_y=E_x\) and \(\nu_{xz}=\nu_{yz}\) and abbreviating \(e=E_x/E_z\), \(m=1-\nu_{xy}-2e\nu_{xz}^2\) as

(62)\[\begin{split}\begin{pmatrix} C_{11} \\ C_{12} \\ C_{13} \\ C_{33} \\ C_{44}\end{pmatrix}&=\begin{pmatrix}E_x\frac{1-e\nu_{xz}^2}{(1+\nu_{xy})m} \\ E_x\frac{\nu_{xy}+e\nu_{xz}^2}{(1+\nu_{xy})m} \\ E_x\frac{\nu_{xy}}{m} \\ E_z\frac{1-\nu_{xy}}{m} \\ G_{yz} \end{pmatrix}.\end{split}\]

Numerical results

Contact radius \(I_R=1.8\)

For the particular value of \(I_R=1.8\), the input values are as follows:

number of spheres	2556
sphere radius	0.03 m
interaction radius \(I_r\)	1.8
\(\hat\pi\) parameter	\(\pi\)
\(C_{11}\)	130 MPa
\(C_{33}\)	55 MPa
\(C_{13}\)	28 MPa
\(C_{12}\)	40 MPa
rel. stress tolerance \(t_r\)	\(10^{-3}\)
prescribed axial strain \(\hat\eps\)	-1%

These are the derived values, computed from geometry of the actual lattice:

average number of contacts	22
lattice density scaling \(\mu\) (56)	2.02
\({E_N^a}'\)	59.7 MPa
\({E_N^b}'\)	5.94 MPa
\({E_T^a}'\)	2.84 MPa
\({E_T^b}'\)	0.989 MPa

Simulation of normal loading yielded the following results; fives symmetric couples are marked with symbols (\(\dagger\), \(\ddagger\), \(*\), \(\circ\), \(+\)), each line corresponds to one axis of compression:

\begin{align*} E_{x}&=103.8\,{\rm MPa}^{\circ}, & \nu_{xy}&=0.2244^{\ddagger}, & \nu_{xz}&=0.3962^{*},\\ E_{y}&=103.1\,{\rm MPa}^{\circ}, & \nu_{yz}&=0.3916^{*}, & \nu_{yx}&=0.223^{\ddagger},\\ E_{z}&=43.93\,{\rm MPa}^{}, & \nu_{zx}&=0.1669^{\dagger}, & \nu_{zy}&=0.1664^{\dagger}, \end{align*}

while three pure shear simulations gave

\begin{align*} G_{yz}&=30.67\,{\rm MPa}^{+} , & G_{zx}&=30.88\,{\rm MPa}^{+} , & G_{xy}&=42.91\,{\rm MPa}^{} . \end{align*}

Dependent values can be written (we used averages for symmetric couples):

\[\begin{split}\nu_{xz}=\nu_{yz}=\frac{E_x}{E_z}\nu_{zx}&=0.392, \\ G_{xy}=\frac{E_x}{2(1+\nu_{xy})}&=42.3\,{\rm MPa}.\end{split}\]

We see that both symmetric couples and dependent values are consistent, with some error due to imperfect isotropicity of the lattice geometry.

These are resulting stiffness tensor components (in MPa) from different sources; bold values were prescribed:

	\(C_{11}\)	\(C_{33}\)	\(C_{13}\)	\(C_{12}\)	\(C_{44}\)
microplane (57)	130	55	28	40	32.2
lattice (51)	130	55.1	28.1	40	32.1
error	0.2%	0.2%	0.3%	0.09%	-0.4%
simulations (60), (61), (62)	122	52.8	26.7	37.9	30.8
error	-5%	-3%	-4%	-5%	-4%

All element tests

The procedure described above was repeated for a range of contact radii \(I_R\in\{1.2, \dots, 2.8\}\); values lower than \(1.2\) (for spheres of equal radii) lead to unstable packing, whereas high values bring no additional information.

This figure shows how the actual sample response (only \(C_{11}\) is shown, figures for other values are similar) approaches the theoretical value for higher \(I_R\) due to better stabilization of the packing, leading to more uniform deformation (the Voigt hypothesis).

Conclusions

An analytical derivation of orientation-dependent stiffness for geometrically isotropic lattice with global cross-anisotropic behavior was presented. The derivation was done withing the microplane framework, related to lattice behavior under the assumption of the Voigt hypothesis (uniform displacements). The orientation-dependent lattice stiffness was written as an interpolation function between in-plane and out-of-plane moduli, reducing the number of elastic parameters to four (instead of five for a general cross-anisotropic material). The \(\mu\) coefficient relates lattice local moduli to global moduli using geometrical configuration of the lattice. Dense geometrically isotropic lattices can be obtained from the Discrete Element Method (DEM) for larger values of contact radius \(I_R\); this procedure was used to perform element tests on periodic lattice structures. It was shown that larger values of \(I_R\) give better agreement between analytical and simulated stiffness tensor, due to stabilization leading to more uniform deformation.

Future work could use the best-fit deformation description [LCYC97] to determine global lattice stiffness to obtain lower estimate of the stiffness for cross-anisotropic medium, instead of the upper bound obtained with the Voigt hypothesis.

Videos

Circular hole in DEM-discretized medium is subjected to radial hydrostatic pressure applied in the cavity. Contact stiffness are cross-anisotropic, with the strike angle \(\alpha=270^{\circ}\) and the dip angle \(\beta=25^{\circ}\). This video shows behavior of the medium when the pressure is applied.

It can be seen that:

• the sound wave propagates faster in the sense of higher stiffness (horizontally)

• the deformation is bigger in the sense of lower stiffness (vertically)

• radial loading leads to axial deformation due to cross-anisotropy axis not aligned with the cavity axis (it goes downwards under the angle \(\beta\)).

Tip

Report issues or inclarities to github.