Setting up virgin stress conditions in discrete element models.

In the present work, a methodology for setting up virgin stress conditions in discrete element models is proposed. The developed algorithm is applicable to discrete or coupled discrete/continuum modeling of underground excavation employing the discrete element method (DEM). Since the DEM works with contact forces rather than stresses there is a need for the conversion of pre-excavation stresses to contact forces for the DEM model. Different possibilities of setting up virgin stress conditions in the DEM model are reviewed and critically assessed. Finally, a new method to obtain a discrete element model with contact forces equivalent to given macroscopic virgin stresses is proposed. The test examples presented show that good results may be obtained regardless of the shape of the DEM domain.

Introduction

In underground construction, often discontinuous behavior occurs around the excavation (i.e. an excavation disturbed/damaged zone (EDZ) is formed). The prediction of the damage induced by underground excavation is important for the designer. Numerical modeling of the EDZ requires a suitable numerical method [1]. Discontinuous behavior of rocks cannot be adequately modeled by continuum models, i.e. the finite element method (FEM), the boundary element method (BEM) or the finite difference method (FDM). However, the discrete element method (DEM) can be used since it models effectively discontinuous behavior. In one approach, the rock mass is represented by rigid circular or spherical particles that are in contact with each other. As has been shown by [2] the method can be used to model a jointed rock mass. However, if the DEM would be used to model a whole tunnel excavation process, too many particles would be required, leading to unrealistic computational cost. Moreover, the assumption of continuous behavior is appropriate for modeling the far field, pointing to a combination of DEM with other methods as the best approach. In the attempt to model realistic problems a combination of the DEM with the FDM [3], the FEM [4], [5] and the BEM [6] has been suggested.

The problem in using the DEM combined with continuum models is the conversion of the initial stress field, that exists in the ground, to forces between the DEM particles. While the macroscopic stresses can be expressed quite easily in terms of inter-particle forces [7], [8], [9], [10], the inverse problem, the one we treat in this work, is more difficult due to its non-uniqueness [7], [8].

In previous works dealing with the discrete modeling of underground structures different approaches have been used to consider virgin stresses. The stressed configuration in the discrete element model can be obtained by applying gravitational load to all the rock mass surrounding an underground structure [11]. In that work, the excavation was characterized with a relatively shallow overburden, which was included in the geometrical model. In some discrete models, in situ stresses can be introduced as traction boundary conditions [12], [13]. In [12], [13] two dimensional discretely jointed models defined on rectangular domains have been loaded on three sides assuming tractions equal to prescribed virgin stresses.

In the mentioned publications, the discrete models employ polygonal elements. In such models the traction can be applied directly on the element side. Application of forces equivalent to traction boundary conditions is not straightforward when cylindrical or spherical discrete element models are used in an irregular computation domain. Special procedures like that developed in [14] are required.

In [2] a method to transfer stresses to a spherical particle based discrete model by means of kinematic conditions corresponding to the strain field from the continuous solution has been presented. Application of the velocities has been suggested either to the boundary or to the whole subdomain, discretized with discrete elements. However, a rigorous examination of the accuracy and efficiency of the proposed methods is missing in [2]. It should be noted that virgin stresses are a consequence of the geological history of the ground and especially near subduction zones may include tectonic stresses. Virgin stresses should be in a condition of steady state (i.e. on or inside the yield surface of the material). In continuum analyses usually these stresses are assigned directly to Gauss points (FEM) or internal stress points (BEM) without applying a “consolidation” load case as it is used in the DEM.

In the present work, a method for transfering a constant virgin stress field to arbitrary shaped DEM domains is presented. This is of particular interest in cases where DEM is coupled with other methods and the particles have to be embedded in FEM or BEM subdomains. After the motivation for this work, which is presented in the next section, a brief summary of the particle method is provided, highlighting the differences to continuum methods. Next the suggested method of stress installation is introduced and verified with test examples, that follow.

The proposed algorithm has been implemented for 2D model and verified in plane strain examples. Analysis of complex engineering problems in geomechanics requires usually a full three dimensional modeling. Nevertheless two-dimensional models are suitable for some cases and are still popular in engineering practice [15]. Two-dimensional simplification enables faster and easier calculations. This is especially important in the discrete element method in which three-dimensional analysis are time consuming. The proposed approach, however, is general and can be easily extended in future to three-dimensional formulation.

Motivation

The motivation of the present work stems from geotechnical engineering and in particular from tunnel excavation. When a tunnel is excavated (Fig. 1a), non-linear behavior may occur in small portions of the total domain. Although the DEM is capable of treating such behavior effectively, using it for the whole domain would not be efficient. We propose that the DEM will be used only in parts of the total domain, where non-linear behavior is expected to appear (Fig. 1b). When such an area has been identified, a DEM subdomain will be introduced and coupled to the remaining domain, modeled by the FEM or BEM (Fig. 1c). However, before the installation of the DEM subdomain, the virgin stresses have to be transfered. A description of the proposed coupling strategy algorithm is discussed elsewhere [16].

Fig. 1

The problem domain: (a) before excavation, (b) after excavation with non-linear zones, (c) with the DEM mesh installed.

Discrete element method formulation

Basic assumptions

Within the discrete element method (DEM) framework, it is assumed that a material can be represented by an assembly of rigid particles interacting with one another. The shape of the particles can be arbitrary. In this work however, cylindrical elements are employed. A discrete element formulation using spherical or cylindrical particles was first proposed by Cundall and Strack [17], [18], [19]. A similar formulation has been developed and implemented by the first author in the discrete and finite element code Dempack [20], [5], [21]. Simulation results presented in this work have been obtained using Dempack.

Equations of motion

The translational and rotational motion of discrete elements (particles) is described by means of the Newton–Euler equations of rigid body dynamics. For the ith element we havewhere u is the element centroid displacement in a fixed (inertial) coordinate frame X, – the angular velocity, m – the element mass, J – the moment of inertia, F – the resultant force, and T – the resultant moment about the central axes. Vectors F and T are respectively composed of the forces and moments due to the external load ( and ), due to the contact interaction with adjacent particles and due to the external damping :where is the number of elements in contact to the ith discrete element, and is the vector connecting the center of mass of the ith element with the contact point with the jth element (Fig. 2). In the present work, only the force-type contact interaction will be considered, resulting in zero values for the interaction moments of Eq. (4).

Fig. 2

Definition of inter-particle interaction.

The damping terms and used in Eqs. (3), (4) are of non-viscous type:where αt and αr, are the respective damping factors for translational and rotational motion.

Constitutive contact model

The overall behavior of the system is determined by the contact laws assumed for the particle interaction. The contact law can be seen as the material model on the microscopic level. Contact between particles can be considered using theoretical models developed in contact mechanics [22]. Numerical studies [23] show that interparticle contact can also be represented properly using a simple linear spring models, which are used in the present work.

The formulation of the contact model employs the decomposition of the contact force between two elements1 Fc into normal and tangential components, Fn and Ft, respectively:where n is the unit normal vector at the contact point (Fig. 2).

Modeling of rock or other cohesive materials requires contact models with cohesion allowing tensile interaction force between particle [2], [24], [25]. In the present formulation, rock materials are modeled using the elastic perfectly brittle model of contact interaction, in which initial bonding between neighboring particles is assumed. The bonds are established between particles i and j satisfying the condition:where r and r are the particle radii, g0 is the overlap or gap at the contact point, and is a length parameter, being a tolerance in contact verification. The cohesive bonds can be broken under excessive load which allows us to simulate initiation and propagation of material fracture. After decohesion, the contact is treated assuming the standard contact model with Coulomb friction. In the present work, however, we will work in the elastic range, only, in which no decohesion occurs.

The normal and tangential particle interactions are modeled (Fig. 3) by linear springs connected in parallel with dashpots providing an additional mechanism to dissipate contact oscillations. Thus, the normal and tangential contact forces are decomposed into the elastic, Fne and Fte, and damping parts, Fnd and Ftd, respectively:The elastic normal force is calculated from the linear constitutive relationship:where kn is the interface stiffness in the normal direction, and un is the change of the distance between the particles evaluated with respect to the distance when the cohesive bond has been establishedConsidering the distance g0 in the formula for un allows us to eliminate initial stresses in the specimen due to overlaps or gaps between bonded particles produced by a generation procedure. The tangential elastic force is given by the linear equationkt – interface stiffness in the tangential direction, ut – relative displacement at the contact point in tangential direction. The relative tangential displacement ut must be evaluated incrementally, cf. [26], [23]:where is the vector of the relative tangential displacement from the previous time step rotated to the present contact plane and Δu is the incremental relative tangential displacementwith v being the relative tangential velocity at the contact point determined aswhere vrc is the relative velocity at the contact pointand vrn is the normal velocity at the contact pointThe contact damping forces in the normal and tangential directions are given byrespectively, where cn and ct are the damping coefficients. The damping coefficients, cn and ct, can be taken as appropriate fractions, ξn and ξt, of the critical damping in the normal and tangential direction, and , respectively:For the system of two rigid bodies with masses m and m, connected with a spring of stiffness k, the critical damping ccr is given by, cf. [27]:By taking k = kn or kt in Eq. (23) the critical damping and is obtained.

Fig. 3

Rheological scheme of the bonded particle interaction model.

Micro–macro relationships in the discrete element model

Discrete element modeling has all the features of material modeling at a lower scale. Specific macroscopic material properties can be obtained assuming adequate contact laws for the interaction between discrete elements (Fig. 4). Depending on the size of discrete elements we can have models at the mezoscopic or microscopic scale. Here, however, we shall not define a specific size of the elements, referring to the models as microscopic models.

Fig. 4

Transition from the microscopic constitutive model to the macroscopic one.

Micro- and macroscopic parameters and variables

The discrete medium at the microscale is characterized by a set of microscopic parameters and microscopic state variables, defined at characteristic points of microstructure, i.e. particles, particle centers and contact points between elements. Having assumed the elastic-brittle model of interaction described above and limiting the analysis to the elastic range we have the following set of microscopic constitutive parameters: kn, kt, αt, αr, ξn, and ξt, and the following microscopic state variables: displacements un and ut, velocities vrn and vt, and forces Fn, and Ft.

In the macroscale we deal with equivalent continuous medium characterized by fields of macroscopic variables. In the description of the macroscopic state we will introduce the macroscopic constitutive parameters: E – Young’s modulus and ν – Poisson’s ratio, and macroscopic state variables: – strain tensor and – stress tensor.

Data transfer between scales

Transfer of the data between scales is an important issue in multiscale or lower scale modeling. It allows us to establish the relationship between the microscopic model parameters and macroscopic variables. The transfer can be applied to different types of parameters. In this work, special attention is paid to the transfer involving macroscopic stresses and microscopic forces.

Transition between scales can be performed two ways, from the lower scale to the upper one (this is called upscaling) or from the upper scale to the lower one (downscaling). Theoretical bases of upscaling have been developed within the theory of homogenization which was also applied to granular materials modeled with the discrete element method [7], [8]. Some authors use the term continualization to refer to homogenization applied to upscaling of discrete models reserving the term homogenization to heterogenous continuous models [9]. Downscaling, the procedure to determine parameters and variables at lower scale using higher-scale information is more difficult, first of all because this is usually non-unique and requires some additional assumptions.

Setting up initial stress conditions (given in macroscopic stress measures) in the discrete element model is a downscale transfer, however, verification of the stress “installation” procedure will require an inverse transfer, calculation of macroscopic stresses corresponding to forces in the microscale. Therefore in this work both downscaling and upscaling of stresses will be explained. Downscaling procedure will require setting of appropriate constitutive parameters of the microscopic model. We will need a knowledge of the relationship between the microscopic and macroscopic constitutive relationship.

Transition from the microscopic to macroscopic model

Effective macroscopic variables and properties in micromechanical models can be determined by various analytical and numerical homogenization and averaging methods [28], [29], [30], [31], [32], [33]. In this work averaging methods based on the concept of the representative volume element (RVE) will be used [34], [35].

The idea of this approach is shown schematically in Fig. 5. We assume that we can define a continuum body occupying the volume Ω, equivalent to the discrete medium considered. The body occupying the volume Ω can be regarded continuous with respect to a certain physical quantity Q, if this variable can be defined at every point x ∈ Ω. For this purpose, around each point x ∈ Ω we define a representative element having volume V, over which the quantity Q will be averaged according to the following formula:The average 〈Q〉, obtained in this way, will be assigned to the considered point x. Performing the averaging procedure for all points in the domain Ω, we transform the microscopic discrete description into a macroscopic continuum description.

Fig. 5

Transition from the microscopic to macroscopic description by averaging on the representative volume element.

The size of the representative volume element should be much smaller than macroscopic dimensions of the considered body Ω (d ≪ L). The RVE size, however, should be sufficiently large in order to eliminate fluctuations typical for the lower scale (Fig. 6). On the other hand the RVE should not be too large, if we want to treat the values of averaged quantities as local ones (if the field of averaged quantity is non-homogeneous).

Fig. 6

Determination of the size of the representative volume element.

The RVE shape can be arbitrary – Eq. (24) is valid for any shape. In mathematical homogenization methods requiring a solution of the boundary problem, a square (in 2D) or cubic (in 3D) RVE is usually assumed because it is easy to define appropriate boundary conditions for these shapes. In our case circular representative volume elements centered at the points x ∈ Ω is used.

Macroscopic stress tensor

We assume that the discrete elements interact among themselves with contact forces determined according to the model described earlier, being the material model at the microscopic scale. The macroscopic stress tensor will be determined using the concept of two-level averaging procedure presented in [36], [10]. In the first stage, the stress tensor for a single discrete element (particle) is calculated as follows [10]:where V is the element volume, n – number of elements being in contact with the pth element, s – vector, connecting the element center with the contact point, Fc – contact force, and the symbol ⊗ denotes the outer (tensor) product. The contact force Fc and the vector s are shown in Fig. 2.

The formula (25) is derived assuming that the particle is in equilibrium and all the forces acting on the particle are taken into account. In case of a constrained particle, except for contact forces we have also reaction forces (Fig. 7). Including reaction forces into the formula (25) is not straightforward since the reaction forces are not applied to the particle surface. Assuming, that reactions are opposed to the other forces acting on the particle we can take into account the reactions in stress calculations in the following way:where (σ) is a component of the stress tensor calculated according to Eq. (25), r is the particle radius, R is the particle reaction, is a corrected micro-stress at the particle on the boundary. Formula (26) means that by including the reaction forces we increase the absolute value of the stress components calculated according to Eq. (25) keeping its sign. Numerical results presented later confirm that the formula (26) improves macroscopic stress calculation. Analogically we can include a reaction moment M in calculation of shear stress components:

Fig. 7

Forces used in stress calculation for a single particle.

After calculation of stresses for single elements (with possible correction for constrained particles), we will perform averaging over representative volume elements, defined for any point x ∈ Ω. The average stress in the representative volume element can be calculated as:where is the intersection of the pth particle and the RVE, . For practical reasons, averaging will be done for particle centroids, only, and the continuous stress field can be constructed by interpolation. Special treatment should be employed for the points on the boundary or near the boundary. In such cases the averaging will be performed considering the intersection of the RVE with the volume Ω occupied by the body.

The expression (28) for the mean stress of a representative volume can be written in an alternative equivalent form, cf. [36]:in which summation is over all N contacts in the representative volume element and L is the so called branch vector connecting the centroids of two particles, for two particles i and j defined as followswhere and are the position vectors of the two particle centroids.

Dimensionless relationships between micro- and macroscopic constitutive properties

Macroscopic properties for the assembly of discrete elements in the linear elastic range can be estimated theoreatically in terms of parameters defining the microscopic models [37], [10]. The approach, used in this work, and proposed in [38], employs dimensionless relationships obtained on the basis of the Buckingham π theorem [39] and numerical simulation of laboratory tests.

Here, we will search functions defining the macroscopic elastic parameters, Young’s modulus E Poisson’s ratio ν, in terms of the introduced earlier microscopic parameters. The dimensionless functional relationships for Young’s modulus E and Poisson’s ratio ν have been assumed in the following form, cf. [38]:The elastic constants are assumed to be functions of the contact normal and tangential stiffness kn and kt, and a certain function ψ characterizing particle packing, dependent on porosity, particle size distribution, density of contact connections (coordination number), etc.

The specific form of the relationships (31), (32) have been obtained by performing numerical simulations of the uniaxial compression test. The results of simulation are presented in Fig. 8. The failure of the specimen is shown with the distribution of macroscopic stresses (averaged over RVE) in the direction of loading are presented in Fig. 8a and the stress–strain curve is plotted in Fig. 8b. The stress values in the plots in Fig. 8b have been evaluated in two ways, firstly, by taking average contact pressure between the specimen and plates, and secondly, by taking the mean value of averaged stresses. Both methods have given similar results which is shown in Fig. 8b. The results of simulation allow us to determine Young’s modulus E and Poisson’s ratio ν for given microscopic parameters. Performing a series of analyses for different values of the ratio kt/kn the form of the functions (31), (32) are determined. The curves obtained are plotted in Fig. 9. Detailed results of calibration procedure are given in [40]. These curves will be used to estimate microscopic parameters yielding required macroscopic elastic properties in the numerical examples presented in this work.

Fig. 8

Results of simulation of the unconfined compression test: (a) specimen after failure with distribution of stress along the loading direction, (b) stress–strain curve.

Fig. 9

Dimensionless relationships between the microscopic parameters and macroscopic elastic constants: (a) relationship for Young’s modulus, (b) relationship for Poisson’s ratio.

Given the Poisson’s ratio ν, we determine first the ratio kt/kn from the curve in Fig. 9b. In order to have a unique result we assume that kt/kn ∈ (0, 1). For the established value of the ratio kt/kn we determine E/kn from the plot in Fig. 9a and calculate kn using the given value of the Young’s modulus E.

The dimensionless curves have been obtained neglecting the function ψ characterizing the particle assembly, therefore microscopic parameters determined from these curves can be treated as approximate ones requiring a certain correction. This correction will be made employing a tuning procedure based on the sensitivity calculation which is presented later. We assume that the characteristics (i.e. particle size distribution, porosity and density of contacts) of particle assemblies used in numerical examples are not very much different from those used in determination of the dimensionless relationships. Therefore the correction should be relatively small.

The virgin stress installation is thought to be an initial stage followed by a subsequent analysis of stress redistribution and possible failure induced by underground excavation. The elastic properties have a big influence on rock deformation due to excavation therefore micromechanical properties used in virgin stress installation should be carefully determined using realistic values of the Poisson’s ratio and Young’s modulus. Keeping the elastic micromechanical parameters determined for the initial stage, it will be also necessary to estimate other model parameters, namely the bond strength and Coulomb friction coefficient, which influence failure and postcritical behavior of the rock. These parameters should be calibrated taking into account relevant macroscopic properties such compressive and tensile strength and the failure envelope, cf. [2].

Downscaling of the stress tensor

We assume we know the macroscopic stresses in the continuum equivalent to a discrete element assembly. The problem to be solved consists in determination of the corresponding discrete element force state. The solution of this problem in general is non-unique. Some additional assumptions are necessary. Different assumptions lead to different solutions. In [36] a hypothesis has been put forward that the distribution of forces on inter-particle contacts is influenced by the packing structure of particles. The hypothesis has been expressed by the formula for the contact force at the ith contact pointwhere n is the unit normal vector at the ith contact point and A is a a certain tensor related to packing structure. Assuming further that the macroscopic stress is consistent with the average stress given by Eq. (29) the tensor A is determined asContact forces evaluated from formulae (33), (34) do not satisfy the equilibrium conditions for each particle, the equilibrium condition is fulfilled in an average sense only. An alternative to Eq. (33) form of the localization formula has been proposed in [8] for the uniform stress field at the micro-scale. Both methodologies mentioned here do not consider rotational equilibrium of particles, they give solutions for non-rotating particles.

As we can see the downscaling procedure based on homogenization theories have their drawbacks. In this work, we will introduce the downscaling procedure employing constitutive relationships at micro-scale, allowing us to obtain the contact forces satisfying translational and rotational equilibrium of particles. Since the stress transfer will be performed indirectly, via the displacement field. The whole downscaling procedure adapted for the virgin stress installation will be presented in the next section.

Virgin stress installation method

Let us consider subdomain A of a continuum body Ω, A ⊂ Ω. In the present work, we will limit ourselves to a two-dimensional problem, Ω ⊂ E2, where E2 is a two-dimensional Euclidean space. The subdomain A is enclosed by the boundary S (Fig. 10). We assume a constant virgin stress field in the subdomain A under plane strain conditions. The subdomain A will be next discretized with spherical discrete elements interacting with one another by contact forces producing averaged stress field equivalent to the continuous stress field . The problem to be solved requires determination of a particle configuration and appropriate microscopic parameters.

Fig. 10

Idea of virgin stress installation.

The proposed solution method employs a virtual stress-free configuration A′. The stress-free configuration is obtained by applying the displacement field −u to the current configurationwith u being the displacement field producing the strain field corresponding to the initial stress field . The displacement field u is defined in the current configuration which is taken as the reference configuration for strain and stress definition.

In order to obtain the displacement field u(x) the strain field (x) must be determined first. We use an inverse macroscopic constitutive relationship, which in the matrix notation can be written as followswhere D is the elastic compliance matrix. Assuming the isotropic elasticity and plane strain conditions the relationship (36) can be written explicitly in the following form:where ν and E are the Poisson’s ratio and the Young’s modulus, respectively, The displacement field is obtained by integration of the strain tensor whose components assuming small strains are given by:To be integrable, the strains must satisfy the compatibility conditions [41]. For the constant strain field, the compatibility conditions are satisfied trivially. The displacement field, however, can be determined up to rigid modes, rotation and translation. Assuming zero rotation and fixing zero displacements at a certain point x0 , integration of expressions (38), (39), (40) with respect to x1 and x2 gives the following displacements:which can be written in a compact form aswhere [] is a strain matrixUsing Eq. (43) into Eq. (35) we can transform the current configuration A to the reference stress-free configuration A′. Zero displacement can be set at an arbitrary point, since any rigid body translation added to a displacement field does not affect the strain field.

The particle assembly is generated in the stress-free configuration A′, the microscopic parameters are estimated and the particles are pushed to the original configuration A. Considering initial separation or penetration between particles according to Eq. (12) allows us to obtain zero stresses in the discrete model after generation.

Transformation of the discrete element assembly to the original (current) configuration can be performed using the relationship obtained by inverting transformation (35) after substitution Eq. (43)where I is the unit matrix. Kinematic conditions defined by transformation (45) are prescribed either to all the particles or to the particles on the boundary, only. In both cases the kinematic conditions are applied gradually assuming a certain constant velocity at respective time interval. When the displacement/velocity conditions are applied to the boundary the particles inside accommodate their positions in accordance with kinematically controlled loading. When the kinematic conditions are applied to the particles in the whole domain, the internal particles are also released and allowed to accommodate their positions after translation to the current configuration. In both approaches, the rotations of all the particles are free. The two methods will be employed and compared in the numerical examples.

After kinematic loading and stress relaxation, when quasi-static equilibrium is obtained, stress averaging is performed and mean values of averaged stresses are calculated and compared with the target values. If the stresses are determined with insufficient accuracy the microscopic parameters are modified and kinematic loading is applied again until we obtain the solution with sufficient accuracy. In the fine tuning procedure we will use sensitivity analysis which indicates necessary change of the model parameters. This will be explained in detail in the first numerical example.

The proposed algorithm of virgin stress installation can be summarized as follows:

Given: E, ν, , S

Calculation of the strain field corresponding to the initial stress field according to Eq. (37).

Discretization of the boundary S with N points x, i = 1, …, N, x ∈ S.

Set the fixed point x0.

Move the points discretizing the boundary S according to Eq. (35) applying the displacements calculated from Eq. (43). This will produce the discretized boundary S′ enclosing the subdomain A′ in the stress free configuration.

Discretize the subdomain A′ with discrete elements.

Estimate the microscopic parameters. First, from the curve in Fig. 9b determine the ratio kt/kn corresponding to the given value of the Poisson’s ratio ν (we assume kt/kn ∈ (0, 1)). For the established value of the ratio kt/kn determine E/kn from the plot in Fig. 9a and calculate kn using the given value of the Young’s modulus E. Assume appropriate damping parameters.

Solve the DEM problem with applied kinematic loading according to (A) or (B):

Calculate new positions x of all the particle centers in the original configuration according to Eq. (45):Solve the DEM problem under kinematic loading with prescribed velocity , with free rotation of all the particles, t being a certain time interval, i = 1, …, n, n is the number of all the particles.

Calculate new positions x of all the external particles (on the boundary) according to Eq. (45):Solve the DEM problem under kinematic loading with prescribed velocity to the particles on the boundary with free rotation of all the particles, i = 1, …, n, n is the number of particles that reside on the boundary.

Stress relaxation with fixed positions of the points on the boundary only until quasi-static equilibrium is achieved.

Calculate averaged macroscopic stresses according to Eqs. (25), (28). Calculate the mean values of averaged stresses.

Verify the mean values of the averaged stresses comparing them to the target stresses. If the accuracy is sufficient then end, otherwise modify kt and kn and repeat steps 6 to 9. Estimation of new constitutive parameters will be explained in the first numerical example.

Numerical examples

In order to demonstrate the proposed methodology, two numerical examples are solved. The first one deals with a square sample that is brought to a compressed state. The second example consists of an arbitrary-shaped subdomain of rock mass under virgin stress.

Square sample

A 2D square sample with dimensions L × L, L = 5 m, shown in Fig. 11a is considered. For this example, the center of the square is taken as the fixed point x. The macroscopic material properties of the sample are chosen to be E = 16.5 GPa and ν = 0.19, for the Young’s modulus and Poisson’s ratio respectively. The sample is considered to be in an initial stress state with , and , having thus a stress ratio .

Fig. 11

Square sample subjected to initial stress: (a) geometry definition, (b) discretization of the uncopmressed domain with particles.

Discretization of the stress-free geometry with particles was carried out using a high-density particle packing algorithm presented in [42] and implemented in GiD graphical preprocessor [43].

The algorithm is based on the solution of an optimization problem, i.e. to minimize the global distance between particles, that have been generated randomly on a background finite element mesh. The algorithm delivers particle assemblies characterized by low porosity. Specifically for 2D problems, porosity less than 10% can be obtained.

The discretization of the stress-free configuration for the square sample is shown in Fig. 11b. The particle mesh has 4111 discs with the mean radius 〈r〉 = 40.6 mm. The distribution of the particle sizes is shown in Fig. 12a and the overall statistics of the particle assembly are listed in Table 1. Apart from the mean radius, Table 1 contains also the following data, that describe the particle assemblies used in the two numerical examples: standard deviation of particle radii, coefficient of variation (CV) defined as the ratio of the standard deviation and the mean value, minimum radius, maximum radius, porosity and coordination number. The porosity of the particle assembly discretizing the square sample is 0.089, and coordination number 5.91 assuming the contact tolerance defined in Eq. (8) to be 0.2rmin, where rmin is the minimum radius in the assembly. Polar orientation of contact directions is presented in Fig. 12b. It can be seen that particle contact directions is quite isotropic although there are 3 preferential directions corresponding to some regular structures in the particle configuration.

Fig. 12

Characterization of the particle assembly discretizing the square domain (a) distribution of the particle radii, (b) polar distribution of the contact directions.

Table 1

Statistics summary of the particle assemblies.

Description	Square sample	Arbitrary shaped sample
Mean radius (mm)	40.60	5.10
Standard deviation (mm)	10.33	1.26
Coefficient of variation – CV	0.254	0.247
Minimum radius (mm)	12.95	0.77
Maximum radius (mm)	87.86	10.77
Porosity	0.089	0.090
Coordination number	5.91	5.88

The microscopic material properties were evaluated from the curves in Fig. 9. Given the Poisson’s ratio ν = 0.19, the ratio kt/kn = 0.2 is taken from Fig. 9b. Then, from the plot in Fig. 9a, the value of E/kn = 1.03 corresponding to kt/kn = 0.2 is determined. Thus, the values kn = 16.019 GPa and kt = 3.2038 GPa for the normal and tangential stiffness, respectively, were calculated.

Having determined the model parameters, the DEM problem is solved in order to obtain the compressed state of the sample. The two strategies presented in Section 5 are employed. In the first strategy, denoted as (A), the kinematic conditions are applied to all the particles. Subsequently, internal particles are released, thereby allowing stress relaxation leading to a quasi-static equilibrium. In the second strategy, denoted as (B), the kinematic conditions are applied to the boundary particles. In both cases rotational degrees of freedom are not constrained. The kinematic loading is followed by an interval of stress relaxation. In both cases adequate damping parameters should be used.

Microscopic stresses after applying displacement to all the particles calculated for each particle according to formula (25) are shown in Fig. 13. Heterogeneity of the stresses can be clearly seen. It can also be seen that the particles on the boundary have significantly different stress values from those of the internal particles. Employing the correction given by Eq. (26) we obtain microscopic stresses shown in Fig. 14. We can see clearly that the values of stresses σ11 (Fig. 14a) are improved at vertical sides and the values of stresses σ22 at horizontal sides (Fig. 14b). Distribution of microscopic stresses after stress relaxation are shown in Fig. 15. The results shown in Fig. 13, Fig. 14, Fig. 15 have been obtained using the global damping factors αt = αr = 0.4 and contact damping factors ξt = ξn = 0.8. The effect of damping has been studied performing simulation with different values of global and local damping. Fig. 16 shows evolution of the micro-stress at the particle No. 1829 (its location is shown in Fig. 15b) during stress relaxation. It can be seen that after 9 ms a steady solution is obtained and similar stress is obtained for different values of damping. Higher values of damping result in faster energy dissipation but the final solution does not depend on the value of damping. Further simulations will be carried out assuming global damping factors αt = αr = 0.4 and contact damping factors ξt = ξn = 0.8.

Fig. 13

Micro-stress distribution obtained by applying the analytical displacements in the whole domain (strategy A) before stress relaxation, kn = 16.019 GPa, kt/kn = 0.2, αt = αr = 0.4, ξn = ξt = 0.8: (a) (Pa), (b) (Pa) (stresses are calculated without correction considering reactions).

Fig. 14

Micro-stress distribution obtained by applying the analytical displacements in the whole domain (strategy A) before stress relaxation, kn = 16.019 GPa, kt/kn = 0.2, αt = αr = 0.4, ξn = ξt = 0.8: (a) (Pa), (b) (Pa) (stresses are calculated with correction considering reactions).

Fig. 15

Micro-stress distribution obtained by applying the analytical displacements in the whole domain (strategy A) and subsequent stress relaxation, kn = 16.019 GPa, kt/kn = 0.2, αt = αr = 0.4, ξn = ξt = 0.8: (a) (Pa), (b) (Pa), marked location of the node 1829 and detail A (stresses are calculated with correction considering reactions).

Fig. 16

Effect of damping on stress relaxation – change of microstress at node 1829, (a) effect of contact damping at constant global damping αt = αr = 0.4, (b) effect of global damping at constant contact damping ξn = ξt = 0.8.

Having calculated microscopic stresses averaging can be carried out assuming an appropriate RVE size. The RVE size should eliminate local stress fluctuations due to microscopic stress heterogeneities as it is explained in Fig. 6. In order to answer the question on an optimum REV size we will study average stress values evaluated for different RVE sizes at a group of neighboring particles in a zone characterized with strong heterogeneity of microscopic stresses. Location of this zone (detail A) is shown in Fig. 15b and the zoom of it is given in Fig. 17a. Variation of the stresses at the particles in this zone evaluated for different RVE sizes is shown in Fig. 17b. It can be seen that for rRVE > 4 rmax, where rmax is the maximum radius in the particle assembly, we obtain practically equal average stresses at the particles investigated. This indicates that rRVE = 4rmax can be accepted as appropriate allowing to eliminate local fluctuations. Distribution of averaged stresses for rRVE = 4 rmax after stress relaxation is shown in Fig. 18. Although local fluctuations are eliminated we can see certain inhomogeneities in stress distributions. We must remember, however, that heterogeneity is an inherent feature of the discrete element model.

Fig. 17

Effect of RVE size on averaged stresses, (a) detail of the investigated area – microstresses for damping parameters kn = 16.019 GPa, kt/kn = 0.2, αt = αr = 0.4, ξn = ξt = 0.8 after stress relaxation, (b) averaged stresses at selected neighboring nodes for different RVE sizes.

Fig. 18

Distribution of averaged stresses for rRVE = 4rmax obtained by applying the analytical displacements in the whole domain (strategy A) and subsequent stress relaxation, kn = 16.019 GPa, kt/kn = 0.2, αt = αr = 0.4, ξn = ξt = 0.8: (a) (Pa), (b) (Pa).

Distribution of microscopic and average stresses is presented in the form of histograms in Fig. 19. It can be clearly seen that averaging leads to more homogenous distributions than that of microscopic stresses. The averaged stress distributions are shifted with respect to the microscopic stresses which can be explained applying Eq. (28) to homogenous microscopic stress distribution with equal for all the particles in RVE. Then we haveInequality (46) indicates that absolute values of averaged stresses must be smaller than absolute values of microscopic stresses.

Fig. 19

Stress histograms – distribution of micro- and averaged stresses obtained by applying the analytical displacements in the whole domain (strategy A) and subsequent stress relaxation, kn = 16.019 GPa, kt/kn = 0.2, αt = αr = 0.4, ξn = ξt = 0.8: (a) σ11, (b) σ22.

Statistics summary of the stress distributions for the square specimen is given in Table 2, in which stress distributions are characterized with the mean value, standard deviation, coefficient of variation, minimum and maximum value. The data confirm the observations from the histograms, in particular the coefficient of variation (CV), which is a useful measure of heterogeneity, is smaller for distributions of averaged stresses in comparison to the CV calculated for microscopic stresses.

Table 2

Statistics summary of the stress distributions for the square sample obtained by applying the analytical displacements in the whole domain (stress values are given in MPa).

Description	Initial microstresses, kn = 16.019 GPa, kt/kn = 0.2		Final microstresses, kn = 16.019 GPa, kt/kn = 0.2		Averaged final stresses, kn = 16.019 GPa, kt/kn = 0.2, rRVE = 2.5rmax		Averaged final stresses, kn = 16.019 GPa, kt/kn = 0.2, rRVE = 4rmax		Averaged final stresses, kn = 16.019 GPa, kt/kn = 0.18, rRVE = 4rmax		Averaged final stresses, kn = 12.55 GPa, kt/kn = 0.24, rRVE = 4rmax
	σ11	σ22	σ11	σ22	σ11	σ22	σ11	σ22	σ11	σ22	σ11	σ22
Mean	−2.06	−8.59	−1.97	−8.31	−1.75	−7.37	−1.79	−7.57	−1.846	−7.516	−1.33	−6.00
Std. deviation	0.44	0.63	0.40	0.48	0.08	0.10	0.071	0.070	0.073	0.069	0.053	0.056
CV	0.21	0.07	0.20	0.06	0.048	0.013	0.04	0.009	0.04	0.009	0.04	0.009
Minimum	−7.22	−14.97	−6.45	−11.67	−1.98	−7.94	−1.99	−7.90	−2.04	−7.84	−1.48	−6.26
Maximum	0.16	−5.54	0.19	−4.84	−0.51	−2.38	−1.49	−7.27	−1.53	−7.23	−1.10	−5.76
σ11/σ22	0.2398		0.2367		0.2374		0.2369		0.2456		0.2217

The ratio of mean values of stress components σ11/σ22 is calculated for each distribution and given in Table 2. It can be noticed that averaging and RVE size do not affect the stress ratio much.

The mean values of the averaged stresses given in Table 2 for the assumed microscopic parameters, are different from the target values, the stress ratio of the mean values does not agree with the target value, either.

We will perform fine tuning of the model parameters aiming to obtain the initial stress state as close as possible to the target values. In the calibration procedure we will try to get a satisfactory stress ratio first. We assume that the stress ratio depends on the ratio of the tangential and normal contact stiffness kt/kn. In the search of the correct value of kt/kn we will use the sensitivity analysis. We will change kt/kn to 0.18 in our model keeping kn without change and solve the DEM problem. The mean values of the averaged stresses for kt/kn = 0.18 are given in Table 2. We can calculate the stress ratio 0.2456 and its changecorresponding to the change of the ratio Δ(kt/kn) = −0.02. Mean values of averaged stresses are denoted here by 〈·〉m. The sensitivity of the stress ratio with respect to the contact stiffness ratio is the followingFor the required change of the stress ratiowe estimate the change of the stiffness ratiogiving kt/kn = 0.24, which is used in the next simulation, which yields the following resultsThe new stress ratiois very close to the target value. Now, we try to correct the level of stresses. Keeping the stiffness ratio kt/kn = 0.24, we scale the stiffness parameters kt and kn multiplying them by the ratio of the target and present value of stressesFinally, we have the following parameters kn = 12.55 GPa and kt = 3.01 GPa. Stress distributions for these parameters are presented in Fig. 20. Main parameters of the averaged stress distribution are summarized in 2. We can see that obtained mean values of the averaged stresses are in perfect agreement with the target stresses. This proves effectiveness of the method of parameter estimation and good performance of the proposed stress installation procedure.

Fig. 20

Distribution of averaged stresses for rRVE = 4rmax obtained by applying the analytical displacements in the whole domain (strategy A) and subsequent stress relaxation, kn = 12.55 GPa, kt/kn = 0.24, αt = αr = 0.4, ξn = ξt = 0.8: (a) (Pa), (b) (Pa).

The alternative procedure of stress installation based on the application of kinematic conditions on the boundary only will be investigated for two sets of parameters, the initial parameters: kn = 16.019 GPa and kt/kn = 0.2 and the final ones kn = 12.55 GPa and kt/kn = 0.24.

The kinematic conditions have been applied assuming a constant velocity. Translational degrees of internal particles and all the rotational degrees are free so during the kinematically controlled loading particle can accommodate their positions. After the final configuration is obtained relaxation of the particles in the domain is continued with the boundary particles fixed at the final positions. Microscopic stresses obtained using the initial contact parameters kn = 16.019 GPa and kt/kn = 0.2are shown in Fig. 21. Comparison with Fig. 15 shows that the final results obtained applying kinematic conditions on the boundary are very similar to the results obtained by applying kinematic conditions to all the particles. The results shown in Fig. 21 have been obtained using the global damping factors αt = αr = 0.4, contact damping factors ξt = ξn = 0.8, prescribing the kinematic loading during the interval of 1 ms (the maximum loading velocity on the boundary was 0.825 m/s) and allowing the stress relaxation during 9 ms. The influence of damping on the solution is shown in Fig. 22 by curves of microstress evolution at node 1829 for different damping values. It can be seen that in all the cases a quasi-static equilibrium has been achieved and the final stress is the same for different damping. Effect of the loading velocity has also been investigated. Fig. 23 shows the microstress evolution at node 1829 for different loading velocity. It can be seen that the final stresses at equilibrium are the same for different loading velocity. The lower velocity, the smaller are stress oscillations. Oscillations are practically not visible for the lowest velocity. In this case the loading is introduced at quasi-static conditions and no additional stress relaxation after loading is necessary.

Fig. 21

Micro-stress distribution obtained by applying the analytical displacements on the boundary (strategy B) and subsequent stress relaxation, kn = 16.019 GPa, kt/kn = 0.2, αt = αr = 0.4, ξn = ξt = 0.8: (a) (Pa), (b) (Pa).

Fig. 22

Effect of damping on microstress evolution at node 1829 – strategy B, kn = 16.019 GPa, kt/kn = 0.2: (a) effect of contact damping at constant global damping αt = αr = 0.4, (b) effect of global damping at constant contact damping ξn = ξt = 0.8.

Fig. 23

Effect of loading velocity on microstress evolution at node 1829 – strategy B (kn = 16.019 GPa, kt/kn = 0.2, αt = αr = 0.4, ξn = ξt = 0.8).

Taking the RVE size rRVE = 4rmax average stresses have been calculated and shown in Fig. 24. Stress distribution in Fig. 24 is practically identical with that in Fig. 18. Similar equivalence can be seen between Fig. 25, Fig. 20 showing the results for the set of final parameters kn = 12.55 GPa and kt/kn = 0.24. Equivalence of the results obtained with the two compared methods is confirmed by the similar values of mean stresses. Statistics summary of the stress distributions obtained applying kinematic conditions at the boundary is given in Table 3. The respective values are identical or very close to the respective results given in Table 2.

Fig. 24

Distribution of averaged stresses for rRVE = 4rmax obtained by applying the analytical displacements on the boundary (v = 0.00825 m/s) and subsequent stress relaxation (strategy B), kn = 16.019 GPa, kt/kn = 0.2, αt = αr = 0.4, ξn = ξt = 0.8: (a) (Pa), (b) (Pa).

Fig. 25

Distribution of averaged stresses for rRVE = 4rmax obtained by applying the analytical displacements on the boundary (v = 0.00825 m/s) and subsequent stress relaxation (strategy B), kn = 12.55 GPa, kt/kn = 0.24, αt = αr = 0.4, ξn = ξt = 0.8: (a) (Pa), (b) (Pa).

Table 3

Statistics summary of the stress distributions for the square sample obtained by applying the analytical displacements on the boundary – strategy B (stress values are given in MPa).

Description	Final microstresses, kn = 16.019 GPa, kt/kn = 0.2		Averaged final stresses, kn = 16.019 GPa, kt/kn = 0.2, rRVE = 4rmax		Averaged final stresses, kn = 12.55 GPa, kt/kn = 0.2.4, rRVE = 4rmax
	σ11	σ22	σ11	σ22	σ11	σ22
Mean	−1.97	−8.34	−1.80	−7.60	−1.34	−6.02
Std. deviation	0.40	0.49	0.072	0.070	0.053	0.056
CV	0.20	0.06	0.040	0.009	0.04	0.009
Minimum	−6.48	−11.72	−2.00	−7.93	−1.49	−6.29
Maximum	0.20	−4.85	−1.49	−7.30	−1.11	−5.78
σ11/σ22	0.237		0.237		0.2220

Similarities between the results obtained with the two methods of application of kinematic conditions are also seen in Fig. 26, Fig. 27, Fig. 28 showing the histograms of microscopic and averaged stress distributions. Fig. 28 demonstrates that the macroscopic stresses obtained by our procedure are very close to the desired values of initial stresses.

Fig. 26

Stress histograms – comparison of micro-stress distributions obtained by applying the analytical displacements in the whole domain (strategy A) and on the boundary (strategy B) after stress relaxation, kn = 16.019 GPa, kt/kn = 0.2, αt = αr = 0.4, ξn = ξt = 0.8: (a) , (b) .

Fig. 27

Stress histograms – comparison of distributions of averaged stresses obtained by applying the analytical displacements in the whole domain (strategy A) and on the boundary (strategy B) after stress relaxation, kn = 16.019 GPa, kt/kn = 0.2, αt = αr = 0.4, ξn = ξt = 0.8: (a) , (b) .

Fig. 28

Stress histograms – comparison of distributions of averaged stresses obtained by applying the analytical displacements in the whole domain (strategy A) and on the boundary (strategy B) after stress relaxation, kn = 12.55 GPa, kt/kn = 0.24, αt = αr = 0.4, ξn = ξt = 0.8: (a) , (b) .

Arbitrary-shaped domain

The second example is defined on an arbitrary-shaped domain (Fig. 29a). Such shapes occur as yielded zones during tunnel excavation. The width W, height H and inner circle radius R of the domain are chosen to be W = 0.95 m, H = 1.7 m and R = 1 m, respectively. The material properties, the Young’s modulus and Poisson’s ratio, are the same as in the previous example, i.e. E = 16.5 GPa and ν = 0.19. The initial stress ratio was assumed to be with initial stress components, and and −2.75 MPa, respectively.

Fig. 29

Arbitrary-shaped domain: (a) definition of main dimensions, (b) discretization of the domain with particles, (c) compressed and stress-free uncompressed domains (deformation is magnified 500 times).

The domain boundary is described by a closed polyline, i.e. a sequence of points connected with each other by straight lines. Subtracting the displacements given by the analytical expression from the points on the boundary of the domain, the shape of the uncompressed domain can be obtained. The zero point P0 used in the displacement analytical expression lies this time outside the domain, and is located at the center of the virtual circle whose part is the arc being the interior boundary of the domain (Fig. 29a). The uncompressed shape meshed with particles is shown in Fig. 29b. The compressed and uncompressed domains are compared in Fig. 29c.

The particle assembly discretizing the domain has 8513 particles with an average radius of 〈r〉 = 0.51 × 10−2 m. Particle size distribution is illustrated with a histogram in Fig. 30a. Statistics summary of the particle distribution is given in Table 1. The particle assembly is characterized with dense packing, porosity is 0.090.

Fig. 30

Characterization of the particle assembly discretizing the arbitrary-shaped domain (a) distribution of the particle radii, (b) polar distribution of the contact directions.

Similar values of the coefficient of variation for both meshes given in Table 1 indicate that the particle assemblies have similar degree of heterogeneity. Having assumed the contact tolerance , the coordination number 5.88 is obtained. Uniform polar orientation of contact directions presented in Fig. 30b shows isotropy of the particle assembly.

Since the coefficient of variation and coordination number for the particle mesh is similar as in the first example we can expect that the dimensionless parameters used in the previous example can be adopted in this example as well. Thus, the same microscopic parameters kn = 12.55 GPa and kt/kn = 0.24 (the final set of microscopic parameters) are taken. Similar damping defined by the global damping factors αt = αr = 0.4 and contact damping factors ξt = ξn = 0.8 has been assumed.

The problem has been analyzed employing the two strategies investigated – in the first one the kinematic conditions have been applied in the whole discretized domain and in the second one on the boundary only. Spatial distribution of microscopic stresses after applying displacement to all the particles is shown in Fig. 31 and after subsequent stress relaxation is shown in Fig. 32. Distributions of corresponding averaged stresses are presented in Fig. 33, Fig. 34. It can be seen that the average stresses distributions are more homogenous than those of microscopic stresses. Fig. 35 shows distribution of average stresses obtained by applying the displacement analytical expression as a boundary condition. Comparison of Fig. 34, Fig. 35 confirms equivalence of the two methods of application of kinematical loading. This equivalence is also demonstrated by the histograms in Fig. 36. It can be also easily seen in Fig. 36 that the peaks in the histograms are very close to the required stress values. Quantitative agreement between the numerical results obtained using the two methods and reaching the target stress values in numerical solution are proved by the numerical results reported in Table 4. For the target stresses and , the respective mean values of averaged stresses and are −1.74 MPa and −2.74 MPa both for kinematic conditions applied in the whole domain and on the boundary. This confirms correct results of our analysis and good performance of the algorithm proposed for setting-up initial stress conditions in discrete element models.

Fig. 31

Distribution of microscopic stresses obtained by applying the displacement condition in the whole domain (strategy A) before stress relaxation, kn = 12.55 GPa, kt/kn = 0.24, αt = αr = 0.4, ξn = ξt = 0.8: (a) (Pa), (b) (Pa).

Fig. 32

Distribution of microscopic stresses obtained by applying the displacement condition in the whole domain and subsequent stress relaxation (strategy A), kn = 12.55 GPa, kt/kn = 0.24, αt = αr = 0.4, ξn = ξt = 0.8: (a) (Pa), (b) (Pa).

Fig. 33

Distribution of average stresses obtained by applying the displacement condition in the whole domain (strategy A) before stress relaxation, kn = 12.55 GPa, kt/kn = 0.24, αt = αr = 0.4, ξn = ξt = 0.8: (a) , (b) .

Fig. 34

Distribution of average stresses obtained by applying the displacement condition in the whole domain and subsequent stress relaxation (strategy A), kn = 12.55 GPa, kt/kn = 0.24, αt = αr = 0.4, ξn = ξt = 0.8: (a) , (b) .

Fig. 35

Distribution of average stresses obtained by applying the kinematic conditions on the boundary (strategy B), kn = 12.55 GPa, kt/kn = 0.24, αt = αr = 0.4, ξn = ξt = 0.8: (a) , (b) .

Fig. 36

Histograms of averaged stresses for the arbitrary-shaped sample comparison of distributions of averaged stresses obtained by applying the analytical displacements in the whole domain (strategy A) and on the boundary (strategy B), kn = 12.55 GPa, kt/kn = 0.24, αt = αr = 0.4,ξn = ξt = 0.8: (a) , (b) .

Table 4

Statistics summary of the stress distributions for the arbitrary shaped sample (kt/kn = 0.24, kn = 12.55 GPa, αt = αr = 0.4, ξn = ξt = 0.8, stress values are given in MPa).

Description	initial microstresses displ. domain before relax.		final microstresses displ. domain after relax.		final microstresses displ. boundary		averaged stresses displ. domain before relax. rRVE = 4rmax		averaged stresses displ. domain after relax. rRVE = 4rmax		averaged stresses displ. bound. rRVE = 4rmax
	σ11	σ22	σ11	σ22	σ11	σ22	σ11	σ22	σ11	σ22	σ11	σ22
Mean	−1.93	−3.04	−1.91	−3.02	−1.91	−3.02	−1.75	−2.77	−1.74	−2.74	−1.74	−2.74
Std. deviation	0.135	0.193	0.142	0.183	0.142	0.183	0.029	0.040	0.030	0.038	0.030	0.038
CV	0.070	0.064	0.075	0.061	0.075	0.061	0.017	0.014	0.017	0.014	0.017	0.014
Minimum	−4.68	−7.95	−4.29	−5.636	−4.29	−5.636	−1.88	−2.930	−1.90	−2.86	−1.90	−2.86
Maximum	−0.84	−0.33	0.09	−0.41	0.09	−0.40	−1.59	−2.394	−1.54	−2.34	−1.54	−2.34
σ11/σ22	0.634		0.633		0.633		0.633		0.633		0.633

Summary and conclusions

A methodology for transferring a constant initial stress to an arbitrary shaped DEM domain has been presented. The macroscopic virgin stresses are transferred via the displacements corresponding to the virgin stress/strain field. Introducing a virtual stress free configuration obtained by applying the inverse displacement allowed us to obtain the stressed configuration at initial undeformed geometry.

Two possible strategies of application of kinematic loading conditions have been investigated. In the first strategy, kinematic conditions are applied to the particles in the whole domain and afterwards stress relaxation in the interior of the domain is allowed, oscillations are damped out and quasi-static equilibrium is obtained. In the second strategy, kinematic conditions are applied on the boundary only. The internal particles can accommodate their positions during kinematically controlled loading and during stress relaxation which follows the loading stage. It has been proved that both strategies are equivalent. They produce practically identical final stress states.

Influence of damping value on the results in both strategies and effect of loading velocity in the second strategy have been studied. It has been shown that the value of damping had influence on the effectiveness of energy dissipation but the final stresses were the same for different damping. Low loading velocity in the second strategy allows to realize quasi-static loading conditions. The level of final stresses, however, is the same as for higher loading velocity with subsequent stress relaxation.

Macroscopic stresses have been calculated by two step procedure. In the first step microscopic stresses have been calculated at each particle. A simple correction taking into account reactions has been proposed for stress calculation at constrained particles. Numerical tests have shown that this correction improves micro-stress calculations at the boundary. In the second step macroscopic stresses have been obtained using averaging over representative volume elements. Influence of the RVE size on averaging results has been studied, the RVE radius equal 4rmax,rmax being the maximum particle radius, has been assumed as appropriate for stress averaging.

The methodology of stress transfer was presented for the uniform stress, however it can also be used for the non-uniform stress state, provided that the corresponding strain field is integrable and the displacement field can be recovered up to rigid motion. The presented methodology can be extended to three dimensional problem with initial stress. The methodology can also be adapted for transfer of macroscopic stresses induced by any loading and at any stage of loading. This would require the knowledge of the evolution of displacement field or boundary geometry.

Attention should be paid to the calibration of the microscopic material parameters, in order to obtain the desired macroscopic behavior, which is manifested by specific values for the Young’s modulus and Poisson’s ratio. In the present work, the calibration was based on the dimensionless micro–macro relationships. The relationships depend on particle packing, so the use of packing with different characteristics than that used in the calibration procedure may require some adjustment of microscopic parameters determined from the curves presented in this work. It has been shown that correction of microscopic parameters based on sensitivity analysis leads efficiently to an accurate solution.