Periodic cells for large-scale problem initialization

In geotechnical applications the success of the discrete element method (DEM) in simulating fundamental aspects of soil behaviour has increased the interest in applications for direct simulation of engineering scale boundary value problems (BVP’s). The main problem is that the method remains relatively expensive in terms of computational cost. A non-negligible part of that cost is related to specimen creation and initialization. As the response of soil is strongly dependant on its initial state (stress and porosity), attaining a specified initial state is a crucial part of a DEM model. Different procedures for controlled sample generation are available. However, applying the existing REV-oriented initialization procedures to such models is inefficient in terms of computational cost and challenging in terms of sample homogeneity. In this work a simple but efficient procedure to initialize large-scale DEM models is presented. Periodic cells are first generated with a sufficient number of particles matching a desired particle size distribution (PSD). The cells are then equilibrated at low-level isotropic stress at target porosity. Once the cell is in equilibrium, it is replicated in space in order to fill the model domain. After the domain is thus filled a few mechanical cycles are needed to re-equilibrate the large domain. The result is a large, homogeneous sample, equilibrated under prescribed stress at the desired porosity. The method is applicable to both isotropic and anisotropic initial stress states, with stress magnitude varying in space.


Introduction
The discrete element method (DEM) was first proposed to study fundamental aspects of soil response [1] and that is still the main field of application of this numerical method in soil mechanics [2].Typically, small cubic or cylindrical samples are generated to obtain a representative elementary volume (REV), which is later tested to explore different links between micromechanical features and mesomechanical response.Element tests may be performed using boundaries that mimic conditions typical of physical soil tests (e.g.rigid walls, membranes) or else using periodic boundaries (e.g.[5]).Periodic boundaries have the advantage of eliminating boundaryeffects thus offering direct access to material responses.Soil mechanical responses are strongly dependent on state.At the mesosocopic level this is typically described by stress and porosity.Ensemble properties such as particle size distribution (PSD) are also significant.Hence attaining a specified initial state of stress and porosity for given PSD, is a crucial step of DEM modelling in soil mechanics.Different procedures are available for initialization.Those used more often include the Radius Expansion Method [1], the Fixed Point Method [3], the Isotropic Compression Method [4], the Modified Isotropic-compression Method [5] and the multi-layer Under Compaction Method [6].These procedures have been tried and tested with some success to create REV's suitable for element testing, but pose important problems when the method is used to address engineering scale geotechnical problems.Indeed, when large models -large in terms of model dimension to particle size-are assembled the computational cost of system initialization quickly becomes prohibitive.Moreover, specimen homogeneity is difficult to attain [7,8].Because periodic boundaries do not appear naturally in the definition of most engineering-scale BVP the possibilities they offer for model building have not been explored before.
In the following the cell repetition method (CRM), a simple but efficient procedure to initialize large-scale DEM models, is presented.In the CRM periodic cells are first generated matching a desired PSD.The cells are then equilibrated at low-level isotropic stress at the target porosity.Once the cell is in equilibrium, it is replicated in space in order to fill the model domain.For uniform stress states only a few mechanical cycles are needed to reequilibrate the large domain after the domain is thus filled.By means of contact force scaling, the procedure is then extended for cases when the stress state varies in magnitude bot not in orientation (i.e. the principal stress directions are uniform within the domain).This allows to initialize a DEM model under prescribed k0 conditions, applying a gravitational field and simultaneously controlling the horizontal stress (hence k0 where k0= Vch/Vcv).As far as we know this is a novel procedure, with important implications for the simulation of BVP.All the simulations described in this paper were performed using the PFC3D software [9].

Periodic cells and REV
We use the problem of REV generation as an example of the potential of the method for DEM initialization under homogenous stress fields.

REV Identification: traditional approach
The REV size is selected to be as small as possible to ease computation, but big enough so that the observed mechanical response is stable with further size increases.This criterion is usually applied by trial and error comparing performance of different-sized cells.Small (with N = 1118 particles), medium (N = 8884) and large periodic cells (N = 29885), were prepared from a PSD with mean radius of 0.5 mm and standard deviation of 0.2.Using the Isotropic Compression Method [4], by setting the contact friction to zero, a dense state (porosity of 0.365) at 100 kPa of isotropic pressure was obtained (Fig. 1).Triaxial compression tests were then performed on the three cells using a linear contact model with normal contact stiffness of 1000 MPa, kn/ks ratio of 2.5 and contact friction coefficient (P) of 0.25.The three samples are first loaded isotropically to 1 MPa and then sheared until the axial strain reaches a value of 30%.The results are shown in Fig. 2. All samples tested showed a similar response and the large fluctuations clearly visible for the small sample reduce as the cell size (particle number) is increased.Depending on the application, either the medium or the large models would be selected as REV.What is here most relevant is to note that the equilibrated initial state for this test was obtained after 35, 85 and 240 minutes for the small medium and large samples respectively.

REV from periodic cell repetition
Now that the REV minimum dimension has been identified, to exploit the idea of the CRM, the equilibrated small sample (which was discarded as a REV) is replicated 8 and 27 times in space in order to fill a cubic cell with twice and three times its initial size respectively.These two new samples are now tested under the same triaxial conditions as the precedent ones.The results are shown in Fig. 3.It appears that the new samples, created from replicas of the precedent small cell are both behaving as REV's.The time needed to generate these samples in an equilibrated condition from replicas of the small cell resulted to be 5 and 10 minutes respectively.If we add the 35 minutes needed to generate the small cell, it appears that we have obtained medium and large samples in 40 and 45 minutes, or, respectively, 47% and 19% of the time spent with the standard procedure.Fig. 2. Deviatoric stress (q) axial strain (Ha) response of the small, medium and large DEM models under triaxial compression at cell pressure of 1 MPa.Fig. 3. Stress-strain triaxial compression response of the medium (8k particles), large (30k particles), 8 times replicated small and 27 times replicated small periodic cell samples.

Periodic cells and BVP
To illustrate the computational efficiency of the CRM and a possible application of this approach, the CRM is here used to initialize a large scale 2D domain, of base B=DL and height H=EL, under a 10g gravitational field.The intent is to use the DEM as a virtual centrifuge to simulate spudcan penetration test.Spudcan penetration is a challenging problem in offshore geotechnical engineering, for which centrifuge model testing has been regularly employed [12] and that is currently simulated through a variety of continuum-based methods [13].The approach followed to initialize this BVP for DEM based simulation is schematized in Fig. 4 and can be summarized as follows: x Generate a LxL 2D periodic cell and equilibrate a sample under a stress state such that i) the effective vertical stress corresponds to the vertical stress the ground should have at a depth of H and ii) the horizontal stress corresponds to the target k0.x Replicate the cell D times in the horizontal direction and E times in the vertical direction.x For each replica added, loop through all the contacts and scale the contact force, CF by a ratio z/H where z is the distance from the ground level to the position of the contact.x Install gravity and fix a thin layer of particles in contact with the horizontal boundaries.x Cycle to equilibrate.

Conclusions
In this paper a novel approach for fast initialization of DEM models has been presented.The cell repetition method (CRM) builds upon a simple idea: using equilibrated periodic DEM models to fill larger spaces.This idea is complemented by force scaling to initialize anisotropic stress fields of magnitude variable in space.Although the anisotropic case has been illustrated with a 2D example, extension to 3D appears straightforward.The computational efficiency, control of initial conditions and homogeneity of the generated specimen give the CRM characteristics that make it very attractive for simulations requiring large scale DEM models.

Fig. 4 .
Fig. 4. a) 2D cell and b) model geometry definition, c) 2511 disk periodic cell and d) corresponding CF for the equilibrated sample.

Fig. 5 Fig. 5 .
Fig.5illustrates some aspects of the equilibrated models, namely the distribution of vertical particle stresses and the contact force distribution.The time required to prepare the model was of roughly 30 minutes.The geometrical dimension of the problem and the number of particles used to generate the model are summarized in Error!Reference source not found.. To compare the computational advantage of using the CRM, a new sample, having the same number of disks was generated at a higher porosity and equilibrated under gravitational filed.The time required to equilibrate the sample under 1g resulted to be around 11000 minutes (just more than 1 week).It is clear that for such application the CRM is

Fig. 6 .
Fig. 6. model showing a,c) the disks coloured by particle displacement and b,d) contact forces scaled by CF magnitude after 0.1m and 0.9m of spudcan penetration in 10g gravitational field respectively.