Section author: Václav Šmilauer <eu@doxos.eu>
Crossanisotropy model with discrete elements¶
Summary¶
The task was to prepare theoretical foundation for modelling rock with crossanisotropic elasticity (schistosity) by using modified DEM formulation.
The following paragraphs present derivation of local stiffnesses of geometrically isotropic linearelastic lattice leading to globally crossanisotropic 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 orientationdependent 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:
pbcpreparesample.py
to prepare monomodal periodic packing with porosity equal to 0.5; alternatively, downloadpbcspheres_N=2556,r=0.03,rRelFuzz=0.txt
which is the result of running this scriptpbcdeformsample.py
to run a single element test on the whole periodic latticepbcdeformsample.table.xls
as input to batch for running element tests on a range of the \(I_R\) parameters by using:woobatch pbcdeformsample.table.xls pbcdeformsample.py
pbcprocesslogs.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 crossanisotropic stiffness computation, woo.utils.muStiffnessScaling
implements computation of the \(\mu\) proportionality coefficient, and woo.dem.WeirdTriaxControl
was used for stresscontrolled element tests.
Acknowledgement
This work was done during postdoc 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¶
Crossanisotropic (also transverseisotropic) materials exhibit rotational symmetry in mechanical response about the crossanisotropy 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
with five independent components.
Lattice¶
The objective of this paper is to model a linearelastic crossanisotropic 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 inbetween 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
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.
The crossanisotropic 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 crossanisotropic nature is introduced by supposing dependency of those moduli on the respective microplane orientation such that symmetries of a crossanisotropic 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 isotropicallyoriented 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 crossanisotropic, 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}'\) 
inplane microplane/lattice moduli 
\(E_*^b\), \({E_*^b}'\) 
outofplane 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
with the transformation Jacobian equal to \(r^2 {\rm s}\theta\):
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 crossanisotropy axis; this lets us write microplane moduli as \(E_N=E_N(\theta)\) and \(E_T=E_T(\theta)\), independent of the azimuth \(\phi\).
where \(\tens{\mathcal{N}}\), \(\tens{\mathcal{T}}\) are projection tensors; in index notation, they read
Stiffness tensor components are subsequently written as
Plugging \(\vec{n}\) in (38) into (40), and reversing integration order, we obtain
As the inner integrals over \(\phi\) do not depend on \(\theta\), they can be evaluated directly, giving
Similarly for other components,
Note that only 5 of the stiffness tensor components are independent, which is consistent with the symmetries of the crossanisotropic 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.
Crossanisotropy¶
\(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 inplane moduli \(E_*^a\) and outofplanemoduli \(E_*^b\), via some weight function \(w(\theta)\):
We require the function \(w\) to be
interpolating between inplane and outofplane 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
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 outofplane 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,
Writing out the independent components separately, we obtain
or in a compact form
with the dependent component
Since the \(\mat{A}\) matrix is invertible, we can determine \(\mathcal{E}\) for given values of \(\mathcal{C}\) as
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 orientationdependent as in (44), with different inplane and outofplane moduli \(\mathcal{E}'\).
Without loss of generality, we consider that the lattice is constructed from spherical packing and that intranodal 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
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)):
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
where \((2r')\) was kept inside the integral, as it is not a constant. Plugging in stiffnesses from (50) yields
where the \(r^2r'\) term appearing inside the integral was written in the front as an average, because it is on average orientationindependent.
Microplanelattice 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\)
We can plug this relation into (39) obtaining
By putting (53) and (55) equal, we solve
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
with \(A=\hat\pi r^2\) and \(l\) being is the crosssection 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
and (47), (49) become respectively
and
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). Crossanisotropy 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:
from prescribed \(\mathcal{C}\), additionally with the dependent component \(C_{44}\) from (48);
from the current lattice stiffnesses using (51);
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.
Evaluating the stiffness tensor¶
Homogeneous periodic lattice is loaded in order to obtain parameters of the crossanisotropic 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 crossanisotropic material. The redundant values can be used to check correctness.
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 crossanisotropic material is a special case:
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
and use measured response \(T_{xx}\), \(\eps_{yy}\), \(\eps_{zz}\) to compute
Shear moduli¶
Shear test is purely deformationcontrolled. The periodic cell is prescribed pure shear \(\eps_{ij}=\eps_{ji}=\hat\eps\) (zero elsewhere). Shear moduli are found from shear stiffness equations
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
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:
while three pure shear simulations gave
Dependent values can be written (we used averages for symmetric couples):
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% 
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 orientationdependent stiffness for geometrically isotropic lattice with global crossanisotropic 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 orientationdependent lattice stiffness was written as an interpolation function between inplane and outofplane moduli, reducing the number of elastic parameters to four (instead of five for a general crossanisotropic 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 bestfit deformation description [LCYC97] to determine global lattice stiffness to obtain lower estimate of the stiffness for crossanisotropic medium, instead of the upper bound obtained with the Voigt hypothesis.
Videos¶
Circular hole in DEMdiscretized medium is subjected to radial hydrostatic pressure applied in the cavity. Contact stiffness are crossanisotropic, 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 crossanisotropy axis not aligned with the cavity axis (it goes downwards under the angle \(\beta\)).
Tip
Report issues or inclarities to github.