Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media

Accepted Manuscript

Discrete Element Method to simulate the elastic behavior of 3D heterogeneous continuous media W. Leclerc PII: DOI: Reference:

S0020-7683(17)30226-3 10.1016/j.ijsolstr.2017.05.018 SAS 9580

To appear in:

International Journal of Solids and Structures

Received date: Revised date: Accepted date:

21 December 2016 2 April 2017 14 May 2017

Please cite this article as: W. Leclerc, Discrete Element Method to simulate the elastic behavior of 3D heterogeneous continuous media, International Journal of Solids and Structures (2017), doi: 10.1016/j.ijsolstr.2017.05.018

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Discrete Element Method to simulate the elastic behavior of 3D heterogeneous continuous media W. Leclerca de Picardie Jules Verne, LTI, EA 3899, F-02100 Saint-Quentin, France

CR IP T

a Universit´ e

Abstract

M

AN US

The paper aims to investigate the suitability of a discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media. In the present work, the cohesion is modeled using beam elements introduced between each pair of particles in contact within the granular packing describing the investigated domain. Such an approach is an interesting alternative to classical continuous methods which enables to take into account more naturally discontinuities and damaging effects under various solicitations. Besides, the development of coupling strategies and efficient codes open new prospects due to ongoing reduction of calculation cost which was again recently the major drawback of discrete approaches. In the present work, we focus our studies on four points. We first take benefit of the new MULTICOR3D code developed in our laboratory to better define the representative volume element in terms of number of particles in the context of a homogeneous domain. Then, we investigate the definition of the stress tensor in the framework of simple heterogeneous configurations, namely a single cylindrical problem and a spherical inclusion one. The third objective is to highlight the suitability of the discrete element method to predict effective elastic coefficients for the same study cases and a large range of contrast of properties. Finally, we apply the cohesive discrete element approach to the case of an interpenetrating phase composite. In this last case, the influence of interfacial debonding effects are also considered and comparisons are done with experimental measures. Whatever the configuration of the heterogeneous medium, results in terms of stress fields and elastic coefficients are always in quite good agreement with the predictions given by other numerical approaches such as the finite element method and the fast Fourier based technique.

PT

1. Introduction

ED

Keywords: Discrete element method, Composite materials, Stress field, Elastic properties, Numerical modeling

AC

CE

The prediction of the mechanical behavior of composite materials remains a difficult challenge due to various phenomena arising in such a medium. Thus, the complexity of their microstructure, their ability to resist to time degradation, the presence of local defects as void areas and interfacial debonding effects are some factors among other illustrating the magnitude of the task. As a result of these issues, continuous approaches such as the Finite Element Method (FEM) are somewhat limited, especially when thin discontinuities arise and multi-scale effects have to be taken in account. The Discrete Element Method (DEM) is a more flexible numerical tool which more naturally treats damaging and contact effects whatever the heterogeneity of the medium. The major drawback of such an approach is its computational cost which remains prohibitive in comparison with continuous methods such as the FEM or the Fast Fourier Transform (FFT) based approach. However, a great effort has been done in the last few years to develop coupling methods enabling to better take advantage of its benefits without the limitations related to the calculation cost. Thus, among other works, Ben Dhia [1] developed the Arlequin method to consider the intermixing of different mechanical approaches. Jebahi et al. [2] introduced a coupling method between the DEM and the Constrained Natural Element Method (CNEM) and applied it to several dynamic tests. Haddad et al. [3] used a numerical model based ∗ corresponding

author. Tel.: +33 323 503 697. Email address: [email protected] (W. Leclerc)

Preprint submitted to International Journal of Solids and Structures

May 22, 2017

ACCEPTED MANUSCRIPT

on a DEM–FEM approach to determine the thermomechanical behavior of the contact interface between two bodies subjected to a tribological loading. The common idea is to consider a classical approach to model the material in continuous parts and introduce the DEM where damaging and discontinuities occur. Besides, the recent development of new codes such as GranOO [4] and the 3D version of MULTICOR [5], MULTICOR3D which is used to handle up to about 16 million particles in the present contribution, pave the way of efficient Discrete Element (DE) calculations.

M

AN US

CR IP T

Initially, the DEM was mainly used to better model physical phenomena in the context of geomechanics. Thus, Cundall [6] introduced the DEM for rock mechanics applications, and Cleary and Campbell [7, 8] simulated landslides. Other works aimed at simulating granular flows as the silo filling treated by Moakher et al. [9] and its discharge studied by Nicot et al. [10]. In this latter work, the definition of the stress tensor according to the Love-Weber formulation [11, 12] as well as the impact of additional inertial effects were considered. DEM was also used to describe tribological problems such as the dry friction in wheel-rail contact [13]. In this context, the main motivation is related to the modeling of a thin layer called third-body which is mainly composed of debris detached from surfaces due to abrasion forces and is submitted to complex multiphysics phenomena. Numerous works discussed the material damaging and its simulation using a DE approach [14, 15] as well. The most important issue was to define cohesion laws to precisely model the continuity of the material. Thus, pioneer works of Schlangen [16] described the cohesion using beam elements introduced between each pair of particles in contact within a granular packing fitting the computational domain. He proved that beam elements provide more realistic crack patterns than spring ones under various mechanical loadings. More recently, Andr´e et al. [17] considered a hybrid particulate-lattice model using the same paradigm. They performed a calibration process to relate microscopic properties associated to beam elements to expected macroscopic ones and exhibited that DEM is suitable to quantitatively yield the linear elastic behavior of homogeneous continuous media. However, it was verified that such a result is only true when a set of intrinsec geometrical parameters related to granular packings such as the compacity and the coordination number are perfectly monitored. In fact, such an approach is greatly sensitive to the arrangement of granular packings as established in the works of Kumar et al. [18]. Recently, some works exhibited the suitability of this approach to simulate the mechanical behavior of heterogeneous media. Thus, Maheo et al. [19, 20] investigated the case of a 3D unidirectional fiber composites and Haddad et al. [21, 22] explored the notion of Representative Volume Element (RVE) and the definition of the stress tensor in heterogeneous media in the context of a 2D random particulate composite.

AC

CE

PT

ED

In the present work, we consider the simulation of 3D heterogeneous media using a beam element based on EulerBernoulli theory in which torsion, tensile and bending terms are taken into account. This model distinguish itself from the beam element used by Andr´e et al. [17] by the presence of coupling terms between bending and tangential terms. Our first objective is to precisely define the minimum number of particles within the granular packing in the context of a 3D homogeneous medium. To authors’ knowledge, even though this is of crucial importance to ensure that no discretization effect could affect our results no accurate result is currently available in the literature. As a result, a 3D version of MULTICOR code, MULTICORD3D based on a C++ programming was developed to handle several millions of particles and enable more in-depth investigations. The second objective is to explore the definition of the stress tensor in the framework of 3D heterogeneous media, namely single fiber and spherical inclusion models. For that purpose, due to local fluctuations related to a DE approach, a mesoscale is introduced using a tessellation of polyhedra called specimens each of them containing several particles. The influence of the size of specimens is mainly discussed and comparisons are given with stress fields and values determined using the FEM. The third purpose consists in exhibiting the ability of the DEM to yield a suitable elastic response in the present context. We consider the same models as in the case of the definition of the stress field and comparisons are carried out in terms of elastic coefficients with several numerical approaches such as the FEM and the FFT-based technique as well as analytical models. Finally, we aim to apply the investigated DE approach to the case of a real particulate composite the microstructure of which is complex and interpenetrated and makes difficult the generation of a 3D mesh in the framework of Finite Element (FE) simulations. Studies are carried out in terms of stress field and effective Young’s modulus and compared to FE predictions obtained using a structured mesh composed of regular hexahedra, and experimental measures. Debonding effects are also introduced at the interface between the matrix and the inclusions in order to better appreciate the influence of such a defect on the elastic response of the composite. The paper is organized as follows. Section 2 describes the numerical modeling of a homogeneous medium using the DEM and the investigation leading to the required minimum number of particles. The third and fourth sections are 2

ACCEPTED MANUSCRIPT

respectively dedicated to the case of a heterogeneous medium with a cylindrical and a spherical inclusion. In each study, comparisons are performed with several numerical and analytical models in terms of stress fields and effective elastic properties. Finally, section 5 is dedicated to an application to the case of a complex particulate alumina/Al composite. Comparisons are done with FE predictions and experimental measures. The influence of interfacial debonding effects is also investigated. 2. Numerical modeling of a homogeneous medium

CR IP T

The present section describes the methodology of simulation of a homogeneous medium using the Cohesive Discrete Element Method (CDEM). In this approach, a continuous medium is modeled by a granular packing composed of spheres and the cohesion is introduced between each pair of particles in contact using beam elements. The relation between macroscopic and microscopic scales is obtained by a calibration process the charts of which directly provide the parameters associated to the beam element to meet expected macroscopic mechanical properties.

(1)

PT

li j ≤ Ri + R j + εD

ED

M

AN US

2.1. ECD Under several assumptions related to the arrangement of particles, the granular packing can be described as an Equivalent Continuous Domain (ECD) in the sense that this is representative enough of the continuous medium. First, an ECD is characterized by a set of geometric parameters which directly affect the relation between the microscopic properties at the scale of the beam element and the expected macroscopic behavior. In order to avoid such sensitivity effects and better control the mechanical properties associated to the ECD, a natural choice is to meet the assumptions of the ”Random Close Packing” (RCP) [23]. For this purpose, granular packings are generated using the efficient Lubachevsky-Stillinger Algorithm (LSA) [24] which enables an accurate control of intrinsec parameters such as the compacity, the polydispersity and the coordination number. Thus, the compacity in terms of volume fraction is set to 0.64 and a small polydispersity is introduced in order to avoid undesirable directional effects. Thus, in the present work, the particle’s radius follows a Gaussian distribution law and the dispersion is described by a Coefficient of Variation (CoV), e.g. the ratio between the standard deviation and the average radius, set to 0.3. Besides, the coordination number Z which is the average number of bonds per particle is also controlled. For that purpose, we consider a sensitivity parameter ε enabling to increase the range of interaction between particles and adjust Z. Let li j be the center to center interparticle distance between two particles i and j in contact with radii Ri and R j , D a characteristic dimension of the whole packing, for example the side of a cubic sample, we assume that bonds arise when the following relation is verified :

AC

CE

ε is adjusted so that Z is set to 6.2 which is in accordance with results given by Andr´e et al. [17]. Please notice that ε depends on the density of the packing. Thus, ε is 1.7e−4 and 5.0−5 for granular packings constituted of respectively 10,000 and 700,000 particles. Besides, a Delaunay triangulation based on the granular medium could have been an alternative method to densify the system and increase the range of interaction. However, such an approach leads to a coordination number higher than 6.2 which does not respond to our expectations. Second, the macroscopic properties of an ECD are very sensitive to the isotropy of the granular system. Typically, this is verified using statistical tools as the polar plots or the 2-point probability function. In the present work, we consider the distribution of contact angles from which we quantify a geometric anisotropy [17]. The main idea is to extract the contact angles associated to each bond in each direction of the Cartesian coordinate system. Then, a comparison is done in terms of probability distribution with the uniform law using polar plots and a parameter describing the geometric anisotropy which can be seen as a quadratic error : v u t !2 Nc 1 X 1 k P − (2) Ek = Nc i=1 i Nc

where Pki is the probability of occurence of a given class i relative to kth direction (x, y and z directions respectively orthogonal to yz, xz or xy planes), and Nc is the number of class. Figures 1a, b and c illustrate the distribution of 3

ACCEPTED MANUSCRIPT

M

AN US

CR IP T

contact angles for granular packings respectively composed of 10,000; 100,000 and 1,000,000 particles and an Nc parameter set to 32. A comparison is done with the uniform law which is represented by a light blue circle. One can notice that a granular packing composed of 10,000 particles is already sufficient to obtain suitable results without peak or irregularity. The influence of the number of particles on the geometric anisotropy is illustrated on Figure 2. Results exhibit that the denser the packing is the lower the quadratic error is which illustrates the convergence to a theoretical perfectly isotropic system. Thus, the quadratic error comes from 1.3e−3 for 10,000 particles to 3.8−4 for 1,000,000 particles whatever the investigated direction. From this study, we conclude that a system composed of 10,000 particles is already sufficient to obtain a geometric isotropy which is in good agreement with the investigations of Andr´e et al. [17]. A denser system only enables to more accurately converge to the uniform law. However, one has to keep in mind that a small defect in the arrangement of particles could be sufficient to drastically affect the relation between microscopic and macroscopic mechanical properties. We refer to the works of Kumar et al. [18] about this critical point. Besides, an apparent geometric isotropy is necessary but could be insufficient to ensure the mechanical isotropy. Thus, the issue of the influence of the density of granular packings on the mechanical properties is a critical point which will be subsequently studied.

AC

CE

PT

ED

Figure 1: Distribution of contact angles for a granular packing composed of (a) 10,000 particles, (b) 100,000 particles and (c) 1,000,000 particles

Figure 2: Influence of the number of particles on the geometric anisotropy

2.2. Cohesive beam model The cohesion between two spherical particles in contact is modeled using a hybrid particulate-lattice derived from the works of Andr´e et al. [17]. This choice is motivated by the ability of the cohesive beam element to model an effective 4

Figure 3: Cohesive beam model

CR IP T

ACCEPTED MANUSCRIPT

aµ = rµ

AN US

continuous material and to produce more realistic crack patterns than the classical spring model [16]. In this approach, the cohesive contact is described by a set of geometric and physical parameters related to a beam element, namely the length of the beam Lµ , the microscopic Young’s modulus Eµ , the microscopic shear modulus Gµ , the cross-section Aµ and the quadratic moment Iµ (Figure 3). The cross-section of the beam element is a disk of radius aµ which depends on the the radii of both particles i and j in contact Ri and R j and a non-dimensional parameter rµ ∈]0, 1] as follows : Ri + R j 2

(3)

2 r2µ  Ri + R j 4 4 a4µ r4µ  Iµ = π = π Ri + R j 4 64

ED

Aµ = πa2µ = π

M

The cross-section Aµ and the quadratic moment Iµ can be directly connected to rµ parameter as follows :

(4) (5)

AC

CE

PT

Thus, the mechanical behavior of the beam is only described by Eµ , Gµ , rµ and the size of particles. However, this latter is connected to the number of particles and has no effect on the macroscopic mechanical behavior for random granular packings dense enough to meet the assumptions of the RCP. From a practical standpoint, the cohesion is modeled by the six-component vector of generalized forces acting as attractive internal forces. A normal component is considered to counteract the relative displacement in traction between both particles in contact, tangential components are introduced to resist to the tangential relative displacement and moment components enable to prevent bending and torsion effects. The expression of internal cohesive forces comes from the classical Euler-Bernoulli beam theory in which coupling terms between bending and tangential effects are considered. The corresponding system for two particles i and j in contact reads :

5

ACCEPTED MANUSCRIPT

  Fnj→i         Ftj→i         Fbj→i      =    Mnj→i         j→i  Mt       Mbj→i

Kn

0

0

0

0

0

0

0

0

Kt

0

0

0

0

0

Kt

0

Kt Lµ 2

Kt Lµ 2

0

Kt L µ 2

Kt Lµ 2 0

0

0

0

0

Sn

0

0

0

0

0

0

Kt Lµ 2

0

Kt L2µ

Kt L2µ

3

6

0

0

0

Kt Lµ 2

0

0

0

0

Kt L2µ

Kt L2µ

3

6

                                    

 uin − unj    uit − utj    uib − ubj   j  i θn − θn    θti   j  θt   θbi   θbj

(6)

CR IP T

                   

Kn = Kt =

AN US

where superscripts n, t and b refer to the main directions of the local Cartesian coordinates system related to the i, j contact between i and j particles. uni, j and ut,b are respectively the normal and tangential displacements associated to i, j i, j i and j particles. θn and θt,b are the components of rotation respectively in torsion and bending. Kn and Kt are the classical normal and tangential stiffnesses respectively given by: Eµ Aµ Lµ 12Eµ Iµ L3µ

(7) (8)

Sn =

M

Sn is a parameter associated to the torsion motion which reads : 2Gµ Iµ Lµ

(9)

Ii θ¨i = Miext +

PT

j

ED

The linear system of equations exhibit the relation between microscopic and macroscopic scales used to determine internal cohesive forces. Translational and rotational equations of motion for a given particle i are written as follows: X mi u¨ i = Fiext + F j→i (10) X

(11)

CE

j

M j→i

AC

where mi is the elementary mass of the particle i and Ii is the quadratic moment of inertia of the particle i. F j→i and M j→i are respectively the force and the moment of interaction of the particle j on the particle i. Fiext and Miext are respectively the external force and moment acting on particle i. Please notice that a Rayleigh damping is added to internal forces in order to reduce dynamic effects and improve the convergence over time. Such a system is solved using an explicit time integration with a formulation based on the velocity Verlet scheme. This choice is mainly motivated by the simplicity of the scheme and the suitability of an explicit scheme to be handled for massive DEM simulations. 2.3. Determination of macroscopic elastic properties We now aim to determine the elastic properties of the ECD. For that purpose, we consider a cubic pattern also called RVE which describes the homogeneous domain under several assumptions of size and isotropy related to the concept of the RCP. The length of the RVE L is set to 1 centimeter but this has no influence on the macroscopic mechanical behavior. We applied a constant velocity v to a given face of RVE and symmetric boundary conditions are set up in order to improve the accuracy of results in comparison to a clamped configuration (Figure 4). Please notice that the 6

ACCEPTED MANUSCRIPT

CR IP T

Figure 4: Cubic representative pattern and applied boundary conditions

AN US

velocity is directly and suddenly imposed but a faster convergence could be obtained by accelerating the boundary particles to the desired velocity from zero. Tangential εT and longitudinal εL strains as well as the reaction force F in the direction of the imposed velocity are extracted from the numerical simulation in order to evaluate the macroscopic Young’s modulus E M and the Poisson’s ratio ν M . Please notice that F is obtained by summing the contributions of all particles belonging to the opposite face of this submitted to the imposed velocity. Thus, let S be the area of the cross-section of the pattern, E M and ν M read as follows : F SεL εT = − εL

EM = νM

(12) (13)

PT

ED

M

The reaction force and consequently the macroscopic elastic properties are sensitive to dynamic effects which are only reduced by the Rayleigh damping. In fact, results only converge for a very high number of time steps which affects the efficiency of the DEM. One solution consists in averaging the results over time in order to accelerate the process for a given accuracy. In the present work, values are only averaged over all previous time steps but a moving average could also be envisaged. In a first step, we aim to estimate the convergence of extracted properties as a function of real time and the number of time steps for two given configurations. The first material is an alumina the microscopic parameters of which are Eµ =2536 GPa and rµ =0.545 which correspond to a macroscopic Young’s modulus of 350 GPa and a Poisson’s ratio of 0.24. The second one is a silica glass the microscopic parameters of which are Eµ =220.2 GPa and rµ =0.722 which correspond to a macroscopic Young’s modulus of 65 GPa and a Poisson’s ratio of 0.2. DE c calculations are performed using an adjustable time step which is determined as a function of critical time steps ∆tcrit associated to each beam element as follows : c ∆t = Dt min(∆tcrit ) c∈ζ

(14)

AC

CE

where ζ designates the set of contacts within the granular packing and Dt is a security ratio which is set to 0.5 in both cases. Thus, ∆t=2.15e−5 s in the case of the alumina and ∆t=5.21e−5 in the case of the silica glass. For information purposes, the granular packing used for DE calculations is composed of 700,000 particles which is a priori sufficient to meet the assumptions of the RCP and ensure the isotropy of the medium. This point will be later discussed. Figures 3a and b illustrate the evolution of both macroscopic Young’s modulus and Poisson’s ratio as a function of the real time for both materials. Green dashed lines are the average properties over time and the asymptotic result is represented by a horizontal black solid line. Results show slow convergences with oscillations decreasing in amplitude over time. Besides, the numerical outputs seem to be very similar in both configurations even if macroscopic Young’s moduli are different. To make comparisons, we consider a convergence criterion based on a relative error of 1e−3 with respect to the asymptotic value. Thus, more than 50,000 time steps are required to reach the convergence whatever the property and the configuration when results are not averaged over time. Conversely, all average coefficients converge for less than 20,000 time steps, namely respectively 17,600 and 5,200 time steps in the case of the Young’s modulus for alumina and silica glass samples. In a second step, we aim to determine the required minimum number of particles to prevent all effects related to the density of the granular packing. For this purpose, the macroscopic elastic properties are estimated for both materials and several densities from 10,000 to 2,000,000 particles. Figures 6a and b illustrate 7

ACCEPTED MANUSCRIPT

M

AN US

CR IP T

the results obtained for the same convergence criterion based on a relative error of 1e−3 as previously seen. Please notice that each data dot is determined using a set of 15 granular packings and error bars represent a prediction interval of 95%. Numerical outputs exhibit a slow convergence whatever the property and the sample. Besides, from a qualitative standpoint, the material configuration has a very low impact on the convergence. From this scope of results, we estimate a required minimum number of particles close to 700,000 DE. This value is much higher than the prediction given by Andr´e et al. [17] which is close to 10,000 DE. However, we hypothesize three explanations. First, the numerical test used by the previous authors is quite different since they used a cylindrical pattern when we consider a cubic one. Second, coupling terms between bending and tangential terms are not considered in Andr´e’s model. As a result, the relation between the microscopic scale and the macroscopic one could be affected. Finally, the previous authors could have underestimated the required density since their convergence criterion is only based on results obtained using granular packings composed of less than 20,000 particles. Nevertheless, the number of 700,000 DE is in quite good adequation with an estimate of 7,000 DE obtained for a 2D configuration in a previous work [21] since 7,0003/2 ≈600,000 DE.

AC

CE

PT

ED

Figure 5: Evolution of macroscopic elastic properties as a function of the real time : (a) case of the alumina sample and b) case of the silica glass sample

Figure 6: Influence of the number of particles on the macroscopic elastic properties (a) case of the alumina sample and b) case of the silica glass sample

2.4. Calibration process As one can observe in previous examples, the microscopic Young’s modulus Eµ does not correspond to its macroscopic counterpart E M and the Poisson’s ratio ν M is a priori unknown. Consequently, a calibration process has to be set up 8

ACCEPTED MANUSCRIPT

M

AN US

CR IP T

in order to identify the local parameters leading to expected macroscopic properties E M and ν M . A large range of configurations are tested with Eµ ∈[2GPa,1000GPa] and rµ ∈[0.1,0.9] and tensile tests with symmetric boundary conditions are performed in order to estimate the macroscopic properties. For information purposes, DE calculations are carried out using a cubic pattern composed of 700,000 particles according to the outcome of the previous study. In this context, the local shear modulus Gµ has no impact on the macroscopic behavior and is consequently not considered in the present process. Figure 7a illustrates the influence of rµ parameter on the non-dimensional Young’s modulus Ea which is the ratio between the macroscopic Young’s modulus E M and the microscopic one Eµ . One can notice that the numerical results do not depend on Eµ but only on rµ . This result is in good agreement with analytical formulations based on Voigt and best-fit hypotheses for spring elements [18, 25]. Figure 7b shows the evolution of the Poisson’s ratio ν M as a function of rµ parameter. As previously seen, Eµ has no impact on the numerical outputs, what is also in accordance with the formulations of Voigt and best-fit hypotheses. Shapes and trends of different curves are very similar to those obtained for cylinders in 2D [21]. From a practical standpoint, polynomial regressions are drawn which provide to a user the microscopic parameters Eµ and rµ leading to the expected macroscopic mechanical behavior.

ED

Figure 7: Influence of rµ parameter on (a) the non-dimensional Young’s modulus Ea and b) the Poisson’s ratio ν M

AC

CE

PT

2.5. Stress field

Figure 8: Introduction of a mesoscopic scale using a tessellation of polyhedral specimens

The stress field within the ECD is obtained using the definition of the stress tensor given by Love-Weber [11, 12]. In this formulation, the stress tensor is determined at a mescoscopic scale using a tessellation of polyhedra called specimens under the assumption that granular and continuous media produce identical internal and external forces [10]. Thus, the average stress associated to each specimen of volume V and supposedly composed of a sufficient

9

ACCEPTED MANUSCRIPT

number of particles is written using tensorial notations as follows: σiVj =

1 X i j rα dα V i

(15)

α∈KV

AN US

CR IP T

where rα is the contact reaction force and dα is the branch vector joining the center of the contacting grains. KVi is the set of contacts in the specimen. To take into account the average of oscillation effects due to the harmonic behavior of the ECD, we extend this definition to a space-time analysis. This consists in considering an expression of the mean stress tensor averaged in space and time [26]. The accuracy and the continuity of the stress field are sensitive to the heterogeneity of the domain, the chosen polyhedron, the average number of particles by specimen and the corresponding standard deviation. Typically, a random network of tetrahedral specimens is more flexible and consequently more practical to approximate complex discontinuities. Conversely, a structured grid of hexahedra is easier to handle, more adapted for homogeneous media, especially for cubic domains, and avoids some variabilities in terms of size and number of particles by specimen. To authors’ knowledge, no consensus exists about the issue of the number of grains within a specimen but a recent work [21] exhibits that a minimum number of 30 particles is necessary to ensure suitable stress fields in 2D heterogeneous media. Such a definition is indeed difficult to set up due to multiple factors affecting the results, among them, the heterogeneity of the medium, the applied solicitations, the dimension of the problem, the model of cohesive link and the expected accuracy.

AC

CE

PT

ED

M

We aim to determine the variability of the stress field within a homogeneous medium. For that purpose, we consider a structured tessellation of regular hexahedra and we limit our investigations to the case of the silica glass sample. The cubic pattern is submitted to the same boundary conditions as those seen in subsection 2.3 and a velocity of 0.00001m/s is imposed in the x direction. σ xx stress field is determined for a strain set to 3.75e−3 using a full range of discretizations from 343 to 21,952 hexahedral specimens and the same granular packing composed of 700,000 particles. We precise that the number of particles by specimen increases from about 30 to about 2040 DE when the number of specimens decreases from 21,952 to 343 hexahedra. Figure 9a shows the variability of the σ xx stress field in the cutting plane y=L/2 for a tessellation of the cubic pattern composed of 3375 specimens. Results exhibit neither privileged directions nor hot spots so that we can assimilate the variability of the stress field to a random noise. Minimum and maximum values which are 191.7 MPa and 344.4 MPa respectively are not located in the investigated cutting plane but close to the loaded face. Besides, extremum values are not symmetrically dispersed compared to the average value. This denotes that stress values could be lognormally distributed but this point should be more in-depth investigated. Figure 9b illustrates the influence of the number of DE by specimen on extremum and average σ xx stress values. One can notice that the greater the number of DE is the lower the variability is since extremum values tend to close in the average value which is 243.3 MPa in the present context. Besides, prediction intervals at a 5% level are also depicted and their magnitude decreases as a function of the number of particles by hexahedral specimen. Figure 9c shows the influence of the number of particles on the CoV which is obtained as the ratio between the standard deviation and the mean value of the stress. This enables to better exhibit the variability of the stress field and provides some criteria to choose a suitable tessellation. As expected, the greater the density is, the lower the variability is. A first indicator level at 10% is reached for about 80 particles per specimen and a second one, at 5%, is reached for about 1200 particles. Thus, a minimum of about 80 particles per specimen is required to ensure a statistical variability less than 10%. However, this result will have to be tempered in the context of heterogeneous stress fields for which the size of specimens also plays an important role due to local stress gradients. 3. Case of a heterogeneous medium with a cylindrical inclusion The present section is dedicated to a heterogeneous medium composed of a single fiber drowned in a matrix. A cylindrical model is considered and investigations are performed in three steps. First, we evaluate the suitable density of particles which is affected by the heterogeneity of the microstructure. Then, we determine stress fields for a ceramic/metal material. Finally, elastic coefficients are estimated for a full range of configurations and compared to numerical and analytical models.

10

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

3.1. Volume conservation

ED

Figure 9: (a) Variability of the σ xx stress field for a tessellation composed of 3375 specimens (Cutting plane y=L/2), (b) Influence of the number of DE by specimen on σ xx stress values and (c) Influence of the number of DE by specimen on the CoV.

AC

CE

PT

We consider a cubic pattern of length L=1cm composed of a single cylinder the axis of which is parallel to the y direction. The radius of the cylinder is equal to one quarter of the length of the cube so that the volume fraction of the cylinder φc is theoretically equal to 0.196. Figure 10a illustrates a granular packing composed of 100,000 particles in which DE belonging to the inclusion are depicted in blue. Please notice that a DE is associated to a given phase according to the position of its center. Besides, φc is estimated by summing the volume of each particle belonging to the inclusion and dividing the result by the contributions of all particles. Thus, spurious irregularities arise at the interface between the matrix and the inclusion, and the number of particles within the cylinder has to be large enough to reduce these effects and meet the theoretical φc . Figure 10b shows the evolution of φc as a function of the total number of particles within the cubic pattern. One can notice that the volume fraction slowly converges to the theoretical value. Strictly speaking, this latter is never respected and we introduce an error threshold of 1% with respect to the theoretical value to estimate the suitable number of particles. Under this assumption, a minimum total number of 1,000,000 DE has to be respected which is a little bit more important that the density estimated for a homogeneous case. This represents about 200,000 particles associated to the sole cylindrical inclusion. Consequently, according to this finding, all calculations are performed in the present section using granular packings composed of 1,000,000 particles. 3.2. Stress fields Stress fields are now determined using the Love-Weber formulation discussed in subsection 2.5. We consider the f 72.5 GPa and Poisson’s ratio 0.33. The corresponding case of a metallic fiber of macroscopic Young’s modulus E M f f microscopic Young’s modulus Eµ is 12,020 GPa and the rµ parameter is 0.169. The matrix phase is composed of 11

CR IP T

ACCEPTED MANUSCRIPT

Figure 10: (a) Geometric model and granular packing, and (b) Volume fraction of the cylinder φc as a function of the number of particles

AC

CE

PT

ED

M

AN US

ceramic material of macroscopic Young’s modulus EmM 350 GPa and Poisson’s ratio 0.24. The corresponding microm scopic Young’s modulus Em µ is 2,536 GPa and the rµ parameter is 0.545. Besides, the microscopic Young’s modulus EΓµ associated to beam elements joining two particles belonging to two different phases throughout the interface Γ beq tween the matrix and the fiber is given by the geometric average of microscopic Young’s moduli so that EΓµ = Eµf Em µ. In a first study, we consider symmetric boundary conditions as already seen in the case of a homogeneous material and we apply an imposed velocity of 1e−5 m/s in the x direction. DE calculations are carried out using a time step of 1.815e−5 s and the stress state is determined at a mesoscale introduced using a tessellation of pentahedral specimens for a strain equal to 3.63e−4 corresponding to 20,000 time steps. In a DE approach, microscopic scale effects lead to local variabilities in the stress field. That is why, the number of specimens has to be carefully chosen in order to sufficiently smooth the results and avoid a loss of informations due to an excessive smoothing where local stress gradients could arise. Results given in subsection 2.5 in the context of a homogeneous stress field provide a minimum bound which is chosen in terms of CoV. Our choice is to consider a criterion of 10% which leads to an average density of about 80 particles per specimen. An upper bound is defined according to a second criterion based on the dimensions of specimens. We define the maximal length of each specimen as the maximum distance between two points of the polyhedron. An overall average is performed and provides a parameter which we designate as the average length of specimens. Our choice is to limit it to a value of about 0.17L which, roughly speaking, signifies that we have about 6 specimens in the width of the cubic pattern. In the present work, reference results are provided by the FEM and help us to determine the suitable specimen within these bounds. However, one has to keep in mind that the choice of the specimen strongly depends on the nature of the stress field and the presence of local stress gradients both existence and location of which are a priori unknown. Thus, a possible method could be to consider a tessellation in the middle of the interval defined by lower and upper bounds and a posteriori proceed to a magnification/refinement step according to a criterion based on stress gradient. Table 1 shows the influence of the number of pentahedral specimens on the extremum and average σ xx stress values. The greater the number of specimens is the greater the maximun value is. Thus, this increases from 139.5 to 236 MPa when the number of specimens increases from 516 to 13,908. For comparison purposes, FE estimates are performed using a structured mesh composed of 5.8 million hexahedra and the same material configuration and boundary conditions. The estimated FE maximum stress value is 216.3 MPa and the minimum one is 37.3 MPa. Hence, we notice that DE results obtained for 3984 pentahedral specimens are in quite good agreement with those given by FE calculations, namely a maximum σ xx value equal to 193.3 MPa and a minimum σ xx value equal to 35.0 MPa. Besides, the average stress value which is obtained using a volume weighting is 99.8 MPa is in perfect agreement with the value extracted from the reaction force. Figure 11a and b illustrate the σ xx stress field in the cutting plane y=L/2 for DE and FE simulations respectively. For information purposes, the DE stress field is estimated using a tessellation of 3984 pentahedra according to the previous discussion. From a qualitative standpoint, stress fields are very similar with a maximum value located close to A position of coordinates (L/2,3L/4,L/2) and a minimum one located close to C position of coordinates (3L/4,L/2,L/2). From a quantitative 12

ACCEPTED MANUSCRIPT

Number of pentahedral specimens Average number of DE per specimen Average length of specimens Maximum σ xx value (MPa) Minimum σ xx value (MPa) Average σ xx value (MPa)

516 1938 0.156L 139.5 54.4 97.2

1072 933 0.125L 162.4 37.8 98.9

2060 485 0.101L 187.2 34.1 101.7

3984 251 0.079L 193.3 35.0 99.9

13908 72 0.054L 236.0 30.1 99.8

Table 1: Influence of the number of specimens on extremum and average σ xx stress values

AN US

CR IP T

standpoint, stress values extracted at A, B, C and D positions for the DE simulation are in quite good agreement with those given by the FE approach. Thus, relative errors with respect to FE estimations are close to 10%. In a second

Figure 11: σ xx stress field (a) Case of DE simulations and (b) Case of a FE simulations (Cutting plane y=L/2)

AC

CE

PT

ED

M

approach, we consider an imposed velocity of 1e−5 m/s in the y direction instead of the x direction. The material configuration and the boundary conditions are not modified so that the time step is still equal to 1.815e−5 s and the stress state is still determined for a strain of 3.63e−4 . Table 2 shows the influence of the number of pentahedral specimens on extremum and average σyy stress values. Please notice that we consider the same tessellations as those used for the previous study under the same hypotheses of size and density of particles. The minimum σyy stress value decreases from 40.6 MPa to 16.6 MPa when the maximum value increases from 148.3 MPa to 278.7 MPa and the number of pentahedral specimens increases from 516 to 13,908. FE simulations are performed to meet the present assumptions using the same mesh of hexahedral elements as previously used. We obtain a minimum stress value of 23.8 MPa and a maximum one of 131 MPa which are quite close to the values given by a tessellation of 1072 elements, namely 27.5 MPa and 158.6 MPa respectively. Furthermore, the average stress value which is 113.2 MPa for this discretization is in quite good agreement with the value extracted from the reaction force, namely 114.0 MPa. Figures 12a and b shows the σyy stress field in the cutting plane x=L/2 for DE and FE simulations respectively. One can notice that in the case of the FE solution the stress field is quasi-homogeneous within each phase. Thus, the σyy stress value is 113.1 MPa in the ceramic part and 25.2 MPa in the metallic one. The σyy stress field given by the DE approach is estimated using a tessellation of 1072 pentahedra according to the previous findings. From a qualitative standpoint, both stress fields are very similar except for small fluctuations which are still noticeable in the case of DE predictions. Comparisons are done at A’(L/2,0,3L/4), B’(L/2,0,L/2), C’(L/2,L,L/2) and D’(L/2,L,3L/4) positions. Thus, σyy stress values at A’and D’ positions are 116.9 MPa and 114.2 MPa respectively and are in good agreement with the FE prediction of 113.1 MPa. Small fluctuations are more noticeable in the case of predictions given at B’ and C’ positions which are 28.5 MPa and 29.0 MPa respectively when the FE prediction is 25.2 MPa. Globally, relative errors with respect to FE estimations are always less than 15% which highlights the ability of the DE approach to yield a suitable stress field in the present configuration. 3.3. Effective elastic properties Effective Young’s and shear moduli are estimated using the DE approach and comparisons are performed using several numerical and analytical models. DE predictions are carried out according to the methodology described in section 2 for a homogeneous material and using a granular packing composed of 1,000,000 particles. In the present configuration, the medium follows a linear elastic transverse isotropic behavior. In other words, the elastic coefficients depend 13

ACCEPTED MANUSCRIPT

Number of pentahedral specimens Average number of DE per specimen Average length of specimens Maximum σyy value (MPa) Minimum σyy value (MPa) Average σyy value (MPa)

516 1938 0.156L 148.3 40.6 115.2

1072 933 0.125L 158.6 27.5 113.2

2060 485 0.101L 191.5 24.7 113.8

3984 251 0.079L 193.3 20.9 113.9

13908 72 0.054L 278.7 16.6 114.0

AN US

CR IP T

Table 2: Influence of the number of specimens on extremum and average σyy stress values

Figure 12: σyy stress field (a) Case of DE simulations and (b) Case of a FE simulations (Cutting plane x=L/2)

PT

ED

M

on the direction, transversely or longitudinally to the cylinder. Thus, the longitudinal Young’s modulus EL refers to the Young’s modulus in the y direction and the transverse Young’s modulus ET refers to the Young’s modulus in the z or x direction since both directions are equivalent. Practically, from a numerical standpoint, ET is obtained by averaging the results obtained for both directions. The longitudinal shear modulus GL refers to the shear modulus in the z direction which is independent of Young’s moduli. We remind that a transversely isotropic material is characterized by five independent coefficients. However, for a sake of conciseness, we only consider three of them and longitudinal and transverse Poisson’s ratios are not studied. DE results are compared to predictions given by two main numerical approaches, the FFT-based method and the FEM. The FFT-based method consists in solving the Lippmann-Schwinger’s equation in Fourier space using an iterative algorithm [27, 28]. In the present work, calculations are performed using the Eyre-Milton scheme [29] and a voxelized map of the cubic pattern consisted of 16,7 million voxels. Two FE approaches are also considered, namely the double-scale homogenization method (2SFEM) [30] which is based on periodic boundary conditions and the classical FEM for which boundary conditions are the same as those considered in the DE approach. For information purposes, all FE calculations are carried out using a structured mesh composed of 34 million tetrahedral elements. Finally, DE predictions are also compared to classical Voigt and Reuss bounds in the case of Young’s moduli.

AC

CE

Effective young’s moduli EL and ET are estimated in the three directions of the cubic pattern using symmetric boundary conditions and an imposed velocity of 0.00001 m/s as described in Figure 4. The longitudinal shear modulus GL is determined using a simple shear test and antisymmetric boundary conditions. Properties are evaluated as a function f of the contrast of properties cr =E M /EmM . Please notice that two configurations are considered. The first one is the case of a fiber less stiff than the matrix so that cr is lower than 1. The second one is the case of a fiber stiffer than the matrix so that cr is higher than 1. Figures 13, 14 and 15 illustrate the influence of the contrast of properties cr on non-dimensional Young’s and shear moduli which are obtained by dividing the effective coefficients by their counterparts associated to the matrix phase, namely EmM and GmM respectively. Results exhibit a good agreement between all numerical approaches and theoretical estimates whatever cr . In the context of EL , one can notice that the results are in perfect agreement with the Voigt estimate. Besides, in the case of ET , all numerical predictions are welllocated between Voigt and Reuss bounds. These findings highlight the consistency and the accuracy of predictions. For example, in the context of a contrast of 100 (cr =1/100 or cr =100), relative differences with respect to the value given by the FEM are 3.70%, 1.24% and 0.45% respectively for EL , ET and GL when cr < 1. Relative differences are 0.40%, 4.17% and 3.09% respectively for EL , ET and GL when cr > 1. Globally, all relative differences are less than 5% whatever the contrast and the numerical approach. These results show the ability of the DE approach to predict 14

ACCEPTED MANUSCRIPT

CR IP T

elastic coefficients in the present configuration. The next section is dedicated to the case of a spherical inclusion in order to investigate the suitability of the DEM to consider other microstructures.

PT

ED

M

AN US

Figure 13: Case of a cylindrical inclusion: non-dimensional longitudinal Young’s modulus as a function of the contrast of properties

CE

Figure 14: Case of a cylindrical inclusion: non-dimensional transverse Young’s modulus as a function of the contrast of properties

4. Case of a heterogeneous medium with a spherical inclusion

AC

The case of a heterogeneous medium composed of a single spherical inclusion embedded a matrix is now studied. Investigations are carried out according to the same steps as previously done for a cylindrical inclusion. Thus, the issue of volume conservation is first studied, then stress fields are determined and compared to FE results. Finally, both Young’s and shear modulus are estimated and compared to predictions given by several numerical and analytical models for a large range of contrast. 4.1. Volume conservation The heterogeneous medium is modeled using a cubic pattern of length L=1cm composed of a centered spherical inclusion. The radius of the sphere is equal to L/3 so that the volume fraction of the sphere φs is theoretically equal to 0.155. Figure 11a shows a section view of a granular packing composed of 100,000 particles in which DE belonging to the inclusion are depicted in blue. As previously done for a cylindrical inclusion, the numerical φs is estimated by summing the volumes of each sphere associated to the spherical inclusion and dividing the result by the sum of 15

CR IP T

ACCEPTED MANUSCRIPT

Figure 15: Case of a cylindrical inclusion: non-dimensional longitudinal shear modulus as a function of the contrast of properties

AC

CE

PT

ED

M

AN US

contributions of all particles. Thus, the numerical φs strongly depends on the density of the granular packing and the total number of particles has to be chosen in order to ensure that numerical and theoretical volume fractions are in good agreement. Figure 16b depicts the influence of the total number of DE on φs . Results exhibit two main conclusions. First, φc slowly converges to the theoretical value. Second, this latter is never reached even for a granular packing composed of 2,000,000 particles. That is why, as done in the case of a cylindrical inclusion, we consider an error threshold of 1% with respect to the theoretical value to determine the suitable number of particles. Thus, a minimum number of particles of about 2,000,000 DE has to be used to ensure the volume conservation. This value is higher than those obtained in the context of homogeneous medium and a cylindrical inclusion, namely 700,000 DE and 1,000,000 DE respectively. For information purposes, the number of particles within the spherical inclusion is close to 310,000 DE. From now on, in the present section and except otherwise stated, all calculations are performed using granular packings composed of 2,000,000 particles.

Figure 16: (a) Geometric model and granular packing, and (b) Volume fraction of the sphere φs as a function of the number of particles

4.2. Stress fields As previously done in the case of a cylindrical inclusion, stress fields are estimated using the Love-Weber formulation discussed in subsection 2.5. Again, we consider the context of the ceramic/metal material the elastic coefficients of which are given in subsection 3.2. Two configurations are investigated. In both cases symmetric boundary conditions are set up and an imposed velocity of 1e−5 m/s is applied. In the first study, the velocity is applied in the x direction 16

ACCEPTED MANUSCRIPT

AN US

CR IP T

whereas this is applied in the y direction in the second problem. All DE calculations are performed using a time step of 1.815e−5 s and stress fields are again evaluated with the help of a tessellation, this time composed of tetrahedral specimens, for a strain equal to 3.63e−4 corresponding to 20,000 time steps. As previously done in the context of the cylindrical inclusion, a set of tessellations is defined so that the average number of particles per specimen is more than 80 particles per specimen and the average length of specimens is less than 0.17L. The choice of the suitable tessellation is limited to a comparison with the FEM but another approach based on the detection of local stress gradients could be used in a general context without reference. Table 3 illustrates the influence of the number of tetrahedral specimens on extremum and average stress values. We precise that σ xx stress values are determined in the first configuration and σyy stress values are determined in the second one. Besides, for comparison purposes, FE calculations are carried out using a structured mesh composed of 5.8 million hexahedra and the same material configuration and boundary conditions. In the present context, stress fields given by FEM are exactly the same whatever the direction of solicitation x or y since both geometry and mesh are perfectly symmetric. Consequently, from now on, we refer to only one FE calculation which gives a maximum value of 172.7 MPa and a minimum one of 27.6 MPa. From a DE standpoint, in the first configuration, the maximum σ xx stress value increases from 154.1 MPa to 355.8 MPa while the minimum σ xx stress value decreases from 37.8 MPa to 1.8 MPa and the number of tetrahedra increases from 1038 to 25152. In the second case, the maximum σ xx stress value increases from 163.1 MPa to 349.5 MPa when the minimum σ xx stress value decreases from 36.6 MPa to 10.4 MPa. Globally, both configurations lead to very close extremum values what respects the material and geometric symmetry. From this scope of results, we deduce that the tessellation composed of 8978 tetrahedra provides the most suitable values which are in quite good adequation with the FE results. We also notice that the average stress values determined in both cases, namely 108.3 MPa and 113.1 MPa respectively for a velocity imposed in x and y direction, do not perfectly match. This highlights the variability effects related to a DE approach which do not exist in a FEM.

ED

M

Number of tetrahedral specimens Average number of DE per specimen Average length of specimens Maximum σ xx value (MPa) Minimum σ xx value (MPa) Average σ xx value (MPa) Maximum σyy value (MPa) Minimum σyy value (MPa) Average σyy value (MPa)

1038 1927 0.182L 154.1 37.8 111.6 163.1 36.6 110.8

3144 1276 0.141L 228.3 30.7 105.3 227.9 32.6 116.5

8978 223 0.087L 202.2 27.1 108.3 214.6 25.9 113.1

25152 80 0.071L 355.8 1.8 108.4 349.5 10.4 113.1

Table 3: Influence of the number of specimens on extremum and average stress values

AC

CE

PT

Figures 17a and b show σ xx and σyy stress fields obtained by DEM using a tessellation of 8978 tetrahedra in cutting planes y=L/2 and x=L/2 respectively. From a qualitative standpoint, results are visually very similar and close to the FE reference given in Figure 17c. Thus, whatever the configuration and the numerical approach, extremum values are located at top and bottom positions of the inclusion while minimum ones are located at both edges of the sphere. From a quantitative point of view, indicators are given at A(L/2,L/2,5L/6), B(L/2,L/2,L/2), C(5L/6,L/2,L/2), D(5/6,L/2,5L/6), E(L/2,5L/6,L/2) and F(L/2,5L/6,5L/6) positions. At A one, σ xx and σyy stress values are respectively 164.3 MPa and 166.4 MPa and in quite good agreement with the FE prediction which is 167.6 MPa. At C and E positions, stress values given by the DE approach are respectively 38.8 MPa and 31.2 MPa and quite well conform the FE estimation of 32.2 MPa. All relative errors with respect to FE predictions are always less than 20% which exhibits the good adequation between FE and DE approaches. We can conclude that the DE approach is able to estimate suitable stress values whatever the geometry of the inclusion, cylinder or sphere. 4.3. Effective elastic properties We aim to evaluate effective Young’s and shear moduli using the DE approach for a large range of contrasts. For that purpose, we consider the methodology described in section 2 and we use a granular packing composed of 2,000,000 particles according to the discussion done in subsection 4.1. In the present context, the heterogeneous medium follows a linear elastic isotropic behavior. In other words, all directions are equivalent and the model only depends on two independent parameters. For a sake of clarity, elastic coefficients are first estimated in the three directions of the Cartesian coordinate system then the results are averaged. Thus, we only consider two mechanical properties: the 17

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 17: (a) σ xx stress field with DE simulations (Cutting plane y=L/2), (b) σyy stress field with DE simulations (Cutting plane x=L/2) and (c) σ xx/yy stress field with FE simulations (Cutting plane x/y=L/2)

AC

CE

PT

ED

M

effective Young’s modulus E and the effective shear modulus G. DE predictions are compared to several numerical approaches, namely the FFT-based method, the double-scale homogenization method (2SFEM) and the classical FEM based on the same boundary conditions and solicitation as those considered in the DEM. Please notice that FFT calculations are carried out using a regular grid composed of 16,7 million voxels and all FE predictions are obtained using a structured mesh composed of 34 million tetrahedral elements. Besides, DE predictions are also compared to an analytical estimate: the Mori-Tanaka model [31]. We precise that all calculations are performed for a constant Poisson’s ratio set to 0.24 which corresponds to a rµ parameter equal to 0.545. Figures 18a and b show the evolution of the non-dimensional Young’s modulus as a function of the contrast of properties. One can notice that DE predictions are in quite good agreement with the other numerical approaches. Thus, relative differences with respect to the prediction given by the FEM are 2.4% and 2.8% for cr =1/100 and cr =100 respectively. Figures 17a and b illustrate the case of the non-dimensional shear modulus. In this present context, relative differences are more important, namely 6.5% and 7.6% for cr =1/100 and cr =100 respectively. However, all DE predictions conform well the estimates given by FFT and FE approaches with a relative difference always less than 8% whatever the numerical reference. That is why, we can conclude that whatever the geometry of the inclusion, cylinder or sphere, and the contrast of properties the CDEM is well-adapted to determine the elastic coefficients of heterogeneous media. The next section is dedicated to an application to a more complex microstructure. 5. Application to the case of a particulate alumina/Al composite The DE approach is applied to the case of a particulate alumina/Al composite. In a first step, the complex interpenetrated microstructure of this material is modeled using the cherry-pit model. Comparisons with FE simulations are performed in terms of stress field and effective Young’s modulus. Then, debonding effects are taken into account in the DEM at the interface between metallic and ceramic phases and results finally compared to a set of experimental measurements.

18

CR IP T

ACCEPTED MANUSCRIPT

ED

M

AN US

Figure 18: Case of a spherical inclusion: non-dimensional Young’s modulus as a function of the contrast of properties

AC

CE

PT

Figure 19: Case of a spherical inclusion: non-dimensional shear modulus as a function of the contrast of properties

Figure 20: Example of granular packing composed of 500 particles with interconnection areas and definition of the cherry-pit model

5.1. Context We consider the context of an Interpenetrating Phase Composite (IPC) the microstructure of which is characterized by the interpenetration of two phases. Such a material has potentially superior properties in comparison with classical 19

ACCEPTED MANUSCRIPT

AN US

CR IP T

materials with embedded inclusions due to the continuity and the percolation of each phase. In the present work, the mechanical response of an IPC composed of a ceramic matrix and a metallic network of interconnected spheres is investigated. The material is elaborated using a specific manufacturing process based on a scaffold of PMMA particles. In a first step, an aqueous solution of acetone (RPE 99.8% Carlo Erba) is used to dissolve the surface of PMMA particles and force the coalescence. Such a process leads to the generation of circular interconnections the size of which is measured and controled by the apparent shrinkage of the material. Then the polymer scaffold is R infiltrated by an alumina slurry (CT1200SG, Almatis ). This step is followed by a debinding treatment (1◦ C/min, ◦ 30h at 200 C), the elimination of PMMA particles and the sintering of the alumina (5◦ C/min, 1h at 1670◦ C). Finally, the obtained macroporous alumina foam is infiltrated by a melt aluminum alloy by capillarity through open pores. Microstructural characterizations conducted using X-Ray computed tomography (OIS Engineering Limited, HMX 225) exhibit the homogeneity of the dispersion of infiltrated pores and local void areas due to closed pores. More details about the manufacturing process are provided in previous works [32, 33]. The modeling of IPC is still quite difficult due to their complex microstructure which is difficult to reproduce and requires a high number of nodes in the context of a FE simulation. In the present case, the particulate structure of the metallic phase enables us to model the microstructure of the IPC using the cherry-pit model [34]. In this approach, each spherical particle is partially penetrable and interconnection sizes are controlled by a λ parameter that lies in the interval [0,1]. Thus, each particle has two different phases, namely a hard core of diameter d=λD encompassed by a fully penetrable concentric shell of diameter D. In the present work, the cherry-pit model is based on a granular packing of hard particles in point contact generated using the LSA. The soft phase is modeled by enlarging the hard core according to λ parameter which is chosen so as to meet experimental measurements of interconnection sizes. Indeed, the dimensionless contact radius R∗c which is obtained by dividing the average contact radius of circular contact areas Rc by the average radius of particles R=D/2, is related to λ parameter as follows: q (16) λ = 1 − R∗c 2

AC

CE

5.2. Stress field

PT

ED

M

In a previous work [33], thermomechanical simulations were performed for a full range of λ parameter using FE and FFT approaches and compared to experimental measurements. In these simulations, the heterogeneous material was supposed perfect in the sense that no void areas or interfacial effects were taken into account. Results exhibited that experimental measurements are much lower than numerical predictions with relative differences coming from 20 to 55% for the effective Young’s modulus. Other findings showed that local void areas have a very little impact on mechanical coefficients since these are limited to a low volume ratio less than 5% with respect to the volume of the sample. As a result, we suspect interfacial effects to affect the mechanical behavior of the alumina/Al composite. We now aim to simulate the mechanical response of the heterogeneous medium using the DE approach and estimate the influence of debonding effects.

Figure 21: σ xx stress field (a) Case of DE simulations and (b) Case of a FE simulations (Cutting plane y=L/10)

20

ACCEPTED MANUSCRIPT

ED

M

AN US

CR IP T

A periodic cubic pattern of the microstructure of the particulate alumina/Al composite is modeled. In a first step, the network of heterogeneities is generated using an initial granular packing composed of 300 monodisperse spherical inclusions in point contact. Then, each one is enlarged to take into account interconnection areas so that the final volume fraction of inclusions is 68.02% which corresponds to a λ parameter equal to 0.975. For information purposes, the number of inclusions is chosen according to the findings of a previous work [33] in order to ensure the representativity of the medium and avoid directional effects. Besides, the length of the cubic pattern L is set to 2 centimeters. Before starting calculations, we aim at assessing the suitable number of particles to ensure the volume conservation within the cubic pattern. In a first approach, an assessment can be given by the result of 310,000 particles per inclusion obtained in the context of a single spherical inclusion. This leads to a total number of about 90 million particles for the set of 300 inclusions and a total number of about 140 million particles. However, it turns out that volumetric compensation effects arise between inclusions with excessive or insufficient quantities of particles and the suitable number is in fact probably much lower. That is why, we consider the same methodology as in previous studies dedicated to single inclusion problems. We investigate several granular packings composed of from about 1.6 million particles to about 125.9 million particles and, in each case, we directly approximate the volume fraction of inclusions by summing the contributions of each particle. Results exhibit that a granular packing composed of about 15.8 million particles leads to an estimated volume fraction of inclusions of 68.00% which is dense enough to respect a volumetric error threshold of 1%. Besides, this value is in good agreement with the 2D estimate of 65,000 DE given by Haddad et al. [21] since 65,0003/2 ≈16,700,000 DE. We consider the following configuration in order to determine and compare the stress field estimated using the Love-Weber formulation with the results given by the FEM. Young’s moduli and Poisson’s ratios are those given in subsection 3.2 for a ceramic/metal material. Symmetric boundary conditions are set up and an imposed velocity of 1e−5 m/s is applied in the x direction. DE calculations are performed using a time step of 1.44e−5 s. σ xx stress field is evaluated with the help of a tessellation of 241,980 tetrahedral elements for a strain equal to 3.42e−4 corresponding to 47,500 time steps. Comparisons are done with FE calculations based on a structured mesh composed of 5.8 million regular hexahedra and the same material configuration and boundary conditions. Figure 20a and b exhibit σ xx stress fields in the cutting plane y=L/10 given by DE and FE approaches respectively. From a qualitative point of view, both stress fields are very similar with hot spots located between the metallic inclusions. From a quantitative standpoint, comparisons are done at A(0.21L,0.07L), B(0.49L,0.63L) and C(0.75L,0.83L) positions. Thus, σ xx stress values at A and C positions are respectively 29.3 MPa and 149.4MPa for DE calculations versus 32.6 MPa and 137.4 MPa for DE ones. It shows a quite good agreement between both approaches. 5.3. Effective Young’s modulus and debonding effects

AC

CE

PT

We now aim to estimate the effective Young’s modulus E of the alumina/Al particulate composite and perform some comparisons with FE predictions and experimental measurements. Practically, we consider the same calculation as previously seen except that results are averaged over time as discussed in subsection 2.3. For a sake of clarity, we suppose that the convergence is reached when the relative error with respect to a previous value determined 2,500 time steps earlier is less than 0.001 two times in a row. This leads us to the total of 47,500 time steps. Experimental values are determined by Resonant Frequency by Damping Analysis (RFDA) using a set of 6 samples of the particulate composite elaborated according to the experimental process described in 5.1. For each sample, the dimensionless contact radius R∗c is related to the shrinkage after the sintering step and the volume fraction of the metallic phase is estimated by hydrostatic weighing. Thus, the average R∗c parameter is 0.2±0.03 and the average volume fraction is 0.682±0.05. These values are in good agreement with the chosen numerical model for which λ is set to 0.975 which corresponds to R∗c =0.218 and the volume fraction is also 0.682. In previous investigations, the microscopic Young’s modulus EΓµ associated to beam elements joining two particles belonging to two different phases throughout their interface Γ was given by the geometric average of local Young’s moduli. In the present study, we introduce an interfacial coefficient CΓ ∈ [0, 1] in order to model debonding effects as follows : q f (17) EΓµ = CΓ Em µ Eµ

f where we remind that Em µ and Eµ are respectively the microscopic Young’s modulus of the matrix and the microscopic Young’s modulus of the fiber. Thus, a CΓ parameter equal to 1 leads to the same interfacial microscopic Young’s modulus as considered in previous studies and a CΓ parameter less than 1 leads to a loss of cohesion at the interface.

21

ACCEPTED MANUSCRIPT

Please notice that a CΓ equal to 0 corresponds to a perfectly debonded interface without any cohesion link. Table 4 illustrates the effective Young’s modulus predicted using DE calculations for CΓ =1, 0.5, 0.2 and 0.1 respectively, the estimate given by the FE approach without debonding and the average of experimental measures. FE and DE Effective Young’s modulus E (GPa)

DE (CΓ =1) 137.0

DE (CΓ =0.5) 131.2

DE (CΓ =0.2) 125.0

DE (CΓ =0.1) 119.5

FE 129.5

Exp. measures 97.4±8.0

Table 4: Comparison of several numerically predicted E values with experimental measures

CR IP T

predictions without debonding (CΓ =1) are in quite good adequation since the relative error with respect to FE estimate is 5.8%. Besides, the lower CΓ is the lower the effective Young’s modulus E is what is in good agreement with a loss of cohesion. However, whatever the value of the interfacial coefficient, estimates remain much higher than the average of experimental measures. As a result, we can conclude that the sole interfacial debonding does not explain the important gap between experimental and numerical values. We hypothesize that this result is more probably due to combined effects of several phenomena as the interfacial debonding, the local damaging and the presence of closed pores. Moreover, we notice that the definition of EΓµ sensitively affects the prediction of elastic coefficients. Thus, the definition of EΓµ could be discussed and a reverse average could also be envisaged.

AN US

Conclusions

AC

CE

PT

ED

M

The present work was dedicated to the investigation of the suitability of a DEM to simulate the elastic behavior of heterogeneous continuous media. First, the DE approach was studied in the framework of a homogeneous medium. Results showed that a minimum of 700,000 particles have to be handled in order to prevent discretization effects. They also exhibited that the stress field which is based on the formulation of Love-Weber, has to be estimated using a tessellation of polyhedra the size of which has to be carefully chosen to smooth the local variability related to the multiscale behavior of the DEM. Second, the case of a heterogeneous medium composed of a single cylindrical inclusion embedded in a matrix was studied. In the case of a ceramic/metal material, stress fields were in quite good agreement with those given by the FEM. Then, effective elastic coefficients were determined for a large range of contrast of properties cr ∈ [1/100, 100]. Results exhibited a quite good agreement with several numerical models with relative differences less than 5%. In a third step, we similarly explored the case of a heterogeneous medium composed of a single spherical inclusion embedded in a matrix. In this context, whatever cr , DE predictions well conform the estimates obtained using the FEM and the FFT-based approach with relative differences less than 8%. Finally, the DEM was applied to the framework of an alumina/Al composite the microstructure of which is interpenetrated and consequently difficult to handle in a FEM. Numerical outputs were compared in terms of stress field and effective Young’s modulus to FE estimates and experimental measures determined by RFDA. Both numerical approaches are in quite good adequation with a relative difference close to 5% in the case of the effective Young’s modulus, but experimental values turned out to be much lower. As a result, interfacial debonding effects were introduced in the DEM using a coefficient affecting the stiffness of the interface. Results exhibited a better but unsatisfactory agreement which led us to assume that combined effects should be considered with local damaging and the presence of void areas. We are currently working on the simulation of damaging effects in 3D heterogeneous media which has already been set up in 2D in a previous work. In a next future, we expect to introduce a thermomechanical coupling using a dilatation model and develop the model for non-linear mechanical behaviors. Conflicts of interest None declared. Bibliography References [1] Ben Dhia H, Rateau G. The Arlequin method as a flexible engineering design tool. Int J Numer Methods Eng 2005;62:1442-1462.

22

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

[2] Jebahi M, Andr´e A, Dau F, Charles JL, Iordanoff I. Simulation of Vickers indentation of silica glass. J Noncryst Solids 2013;378:15-24. [3] Haddad H, Guessasma M, Fortin J. A DEM–FEM coupling based approach simulating thermomechanical behaviour of frictional bodies with interface layer. Int J Solids Struct 2016;81:203-218. [4] Andr´e D, Charles JL, Iordanoff I, N´eauport J. The GranOO workbench, a new tool for developing discrete element simulations, and its application to tribological problems. Adv Eng Softw 2014;74:40-48. [5] Sanni I, Bellenger E, Fortin J, Coorevits P. A reliable algorithm to solve 3D frictional multi-contact problems : application to granular media. J Comput Appl Math 2010;234:1161-1171. [6] Cundall PA, Strack ODL. Discrete numerical model for granular assemblies. Geotech 1979;29:47-65. [7] Cleary PW, Campbell CS. Self-lubrication for long run-out landslides: examination by computer simulation. J Geophys Res 1993;98(B12):21911-24. [8] Campbell CS, Cleary PW, Hopkins MA. Large-scale landslide simulations: global deformation, velocities, and basal friction. J Geophys Res 1995;100(B5):8267-83. [9] Moakher M, Shinbrot T, Muzzio FJ. Experimentally validated computations of flow, mixing and segregation of non-cohesive grains in 3D tumbling blenders. Powder Technol 2000;109:58-71. [10] Nicot F, Hadda N, Guessasma M, Fortin J, Millet O. On the definition of the stress tensor in granular media. Int J Solids Struct 2013;50:25082517. [11] Love AEH. A treatise of mathematical theory of elasticity. Cambridge University Presse, Cambridge, 1927. [12] Weber J. Recherches concernant les contraintes intergranulaires dans les milieux pulv´erulents, application a` la rh´eologie de ces milieux. Bulletin de Liaison des Ponts et Chauss´ees, 1966. [13] Bucher F, Dmitriev AI, Ertz M, Knothe K, Popov VL, Psakhie SG, Shilko EV. Multiscale simulation of dry friction in wheel/rail contact. Wear 2006;261:874-884. [14] Hentz S, Donz´e FV, Daudeville L. Discrete element modelling of concrete submitted to dynamic loading at high strain rates. Comput Struct 2004;82(29-30):2509-2524. [15] Tan Y, Yang D, Sheng Y. Discrete element modelling of fracture and damage in the machining of polycristalline sic. J Eur Ceram Soc 2009;29(6):1029-1037. [16] Schlangen E, Garboczi EJ. Fracture simulations of concrete using lattice models : Computational aspects. Eng Fract Mech 1997;57(2):319332. [17] Andr´e D, Iordanoff I, Charles JC, Neauport J. Discrete element method to simulate continuous material by using the cohesive beam model. Comput Methods Appl Mech Eng 2012;213:113-128. [18] Kumar R, Rommel S, Jauffr`es D, Lhuissier P, Martin C. Effect of packing characteristics on the discrete element simulation of elasticity and buckling. Int J Mech Sci 2016;110:14-21. [19] Maheo L, Dau F, Andr´e D, Charles JL, Iordanoff I. A promising way to model cracks in composite using Discrete Element Method. Compos Part B 2015;71:193-202. [20] Le Bd, Dau F, Charles JL, Iordanoff I. Modeling damages and cracks growth in composite with a 3D discrete element method. Compos Part B 2016;91:615-630. [21] Haddad H, Leclerc W, Guessasma M. Application of DEM to predict the elastic behavior of particulate composite materials. Granul Matter 2015;17:459-473. [22] Leclerc W, Haddad H, Guessasma M. On the suitability of a Discrete Element Method to simulate cracks initiation and propagation in heterogeneous media. Int J Solids Struct 2017;108:98-104. [23] Donev A, Cisse I, Sachs D, Variano EA, Stillinger FH, Connelly R, et al. Improving the density of jammed disordered packings using ellipsoids. Sci 2004;303:990-993. [24] Lubachevsky BD, Stillinger FH. Geometric properties of random disk packings. J Stat Phys 1990;60:561-583. [25] Liao CL, Chang TP, Young DH, Chang CS. Stress-strain relationship for granular materials based on the hypothesis of best fit. Int J Solids Struct 1997;34:4087-4100. [26] De Saxc´e G, Fortin J, Millet O. About the numerical simulation of the dynamics of granular media and the definition of the mean stress tensor. Mech Mater 2004;36:1175-1184. [27] Moulinec H, Suquet P. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Comput Methods Appl Mech Eng 1998;157:69-94. [28] Michel J, Moulinec H, Suquet P. Effective properties of composite materials with periodic microstructure : a computational approach. Comput Methods Appl Mech Eng 2001;172:109-143. [29] Eyre D, Milton G. A fast numerical scheme for computing the response of composites using grid refinement. Eur Phys J Appl Phys 1999;6:41-47. [30] Sanchez-Palencia E. Non-Homogeneous Media and Vibration Theory, Springer, Berlin, 1980. [31] Mori A, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall 1973;2:571-574. [32] Ferguen N, Cogn´e C, Bellenger E, Guessasma M, P´elegris C. A numerical model for predicting effective thermal conductivities of alumina/Al composites. J Compos Mater 2013;47:3311-3321. [33] Leclerc W, Ferguen N, P´elegris C, Haddad H, Bellenger E, Guessasma M. A numerical investigation of effective thermoelastic properties of interconnected alumina/Al composites using FFT and FE approaches. Mech Mater 2016;92:42-57. [34] Torquato S. Random heterogeneous materials. New-York: Springer; 2002.

23