Abstract

This paper deals with the characterization and the numerical modelling of the collapse of composite hollow spherical structures developed to absorb energy during high velocity impacts. The structure is composed of hollow spheres ( mm) made of epoxy resin and mineral powder. First of all, quasi-static and dynamic ( mm·min−1 to  m·s−1) compression tests are conducted at room temperature on a single sphere to study energy dissipation mechanisms. Fracture of the material appears to be predominant. A numerical model based on the discrete element method is investigated to simulate the single sphere crushing. The stress-strain-time relationship of the material based on the Ree-Eyring law is numerically implemented. The DEM modelling takes naturally into account the dynamic fracture and the crack path computed is close to the one observed experimentally in uniaxial compression. Eventually, high velocity impacts ( m·s−1) of a hollow sphere on a rigid surface are conducted with an air cannon. The numerical results are in good agreement with the experimental data and demonstrate the ability of the present model to correctly describe the mechanical behavior of brittle materials at high strain rate.

1. Introduction

Hollow sphere structure (HSS) belongs to cellular solids that have been studied recently for their multiple properties [1].Compared to foam structure, HSS has both closed and open porosity. Its applications can be of various types including energy absorber [2], acoustic damping [3], and thermal insulation [4]. In our case, HSS aims to absorb soft impacts energy on an airliner cockpit. Most of impacts are due to bird-strikes during take-off or landing at high velocity ( m·s−1). Protective structures have to dissipate the kinetic energy of the projectile by all possible means: friction [5], plastic deformation [6], or fracture [7]. Many materials and structures have been tested to absorb energy of external impacts such as aluminium honeycombs [8], polymer [9], and metallic [10] foams and composite materials [7]. CELPACT [11] and MANSART [12] projects have recently focused on cellular structures designed for energy absorption. HSS is now investigated through the SAMBA (Shock Absorber Material for Bird-shield Application) project because of its promises in terms of specific energy dissipated (J·kg−1) during impact [13]. Hollow spheres can also easily fill a sandwich structure to be placed in front of the plane. The commonly used bird-shield structure measures approximately one square meter and ten centimeters thick which represents 10–12 kg.

Hollow spheres used in this study are made of commercial epoxy resin, a thermosetting polymer subjected to crack propagation [14]. In fact, rapid crack propagation (RCP) is observed to be predominant at macroscopic scale to dissipate the available energy stored in the structure [15]. As it is commonly known, most thermosetting polymers have a strain rate and temperature dependence on their mechanical properties (Young’s modulus, yield strain) [16, 17]. An influence of the compression velocity can be predicted.

In the interest of studying the fracture on polymer hollow spheres, a numerical model is investigated. It is based on the discrete element method, often used for rocks or soils modelling [18]. It reveals to be an interesting way to model the mechanical behavior of brittle materials. It allows large displacements and strain and the fracture of the material. It offers an alternative to the FEM by allowing a natural propagation of the crack [19]. Indeed the DEM does not require an adaptive meshing refinement at the crack tip contrary to some finite element methods (FEM) [20, 21]. The use of “classical” continuum models such as the FEM was well adopted to simulate the crushing of elastoplastic hollow spheres [2225]. Despite several methods (elements erosion, cohesive zones, extended finite element method, etc.) developed to represent the fracture of a material [2628], a discrete model is chosen. Continuous material behavior can now be simulated with the DEM by using cohesive beams linking each discrete element [29]. The crack closure can be taken into account by activating and computing contacts between discrete elements.

Before optimizing an assembly of HSS (which corresponds to the macroscopic scale) as shown in Figure 1(a) to dissipate more energy when crushed, the present paper focuses on the investigation of one sphere crushing (mesoscopic scale). The sphere behavior under impact, which remains important to study, is often quasi-invisible at macro-scale. It leads nevertheless to a better understanding of energy dissipation mechanisms. That is why, an experimental and numerical approach is proposed to take into account the fragmentation on a single composite hollow sphere. Experimental part will present the structure response as (force displacement curve) as well as the characterization of the constitutive material for different strain rates and temperatures. A fracture criterion based on the principal stress and changing according to the Ree-Eyring formalism [30] is implemented in the DEM. Eventually, high velocity impacts are conducted with an air cannon. Hollow spheres are fired on a rigid surface, and multicracks are observed and compared to the numerical results.

2. Materials and Methods

2.1. Hollow Spheres

Hollow spheres are manufactured and patented [32] by ATECA SAS. The constitutive material of the hollow spheres is composed of an epoxy resin with mineral powder to mainly prevent the spheres from sticking together during the manufacture process. It can also increase the structural rigidity and the material yield stress. Powder aggregates are in high proportion in the matrix. The grain size varies from 10 μm to 100 μm as it can be observed with the help of microscope observations; see Figure 1(c). The presence of this reinforcement increases the density of the material from 1200 kg·m−3 (resin only) to 1880 kg·m−3. Hollow spheres are manufactured by coating the resin on expanded polystyrene foam. They are then heated to polymerize the composite layer and to melt the polystyrene. Therefore, the spheres can be considered as hollow, the polystyrene remaining being negligible compared to the sphere’s shell. The hollow sphere diameter depends on the size of the polystyrene ball. It varies within  mm. The thickness ratio is defined by , where is the shell thickness and the average radius (), in Figure 1(b).

Tested hollow spheres present a diameter of  mm, a radius ratio , and a mass  g. For clarity, only results obtained on one geometry are presented although experiments with different radius ratios have been conducted.

2.2. Experimental Procedure

Uniaxial compressive tests on hollow spheres were conducted on two machines: one classical compression machine (Zwick Roell Z250) for quasi-static tests and an original fly wheel (Figure 2) for dynamic tests [31]. The rotational movement of the wheel is transformed into a translational movement. The velocity of the moving plates can be considered as constant due to the high moment of inertia of the wheel. In real conditions of a bird strike, velocity up to 175 m·s−1 can be reported. The maximum macroscopic strain rate of the sandwich structure is then close to approximately 2000 s−1. The macroscopic strain rate of 50 s−1 achievable with the fly wheel is supposed to be a good indication to show the influence of a dynamic solicitation. Conventional compression Hopkinson’s bars have not been able to be directly used because of the high heterogeneity of the HSS. High velocity impacts are conducted with a 35 mm in diameter and 1 m length air cannon. An air booster is used to pressurize an air tank at 60 bars. Hollow spheres are placed at the entry of the cannon and are fired by releasing the pressure with a solenoid valve. The impact surface is made of steel and is considered perfectly rigid. Impact velocities up to 130 m·s−1 were measured by image analysis. Force sensors of 10 kN were used in quasi-static and dynamic compression tests. Force as a function of time for the high velocity impacts was measured with a 150 kN force sensor (Kistler 9377c). A high speed camera (Photron SA-5) is used to capture the hollow sphere crushing and the propagation of cracks. 75000 frames per second with a resolution of pixels is set for uniaxial compression tests. Due to the very short time of the high velocity impacts, the camera is set to 150000 frames per second with a resolution of pixels.

To characterize the constitutive material of the hollow sphere, compression tests were conducted on a compression testing machine (Zwick Roell Z250) equipped with a climatic chamber. High temperature is reached by thermal resistance and low temperature is reached by adding liquid nitrogen. A control loop keeps the temperature constant during the tests. A 250 kN force sensor is used. Hopkinson bars (25 mm in diameter, aluminium) allow conducting high strain rate characterization [33, 34] at room temperature.

2.3. The Discrete Element Method

The DEM was originally developed for granular assemblies [35] but its applications have been expanded these last years to other areas such as tribology [36], fluid mechanics [37], and continuum mechanics [38]. The GranOO workbench [39, 40], developed in the laboratory, is naturally chosen. GranOO is based on an explicit time integration scheme (velocity Verlet scheme) [41] which permits a reasonable computational cost for high velocity simulations. Discrete element spatial position and velocity are calculated using the acceleration determined by the second Newton’s law of motion. The mass of each element is found so that it represents the density of the material modelled. Linear positions and orientations are, respectively, associated with resultant forces and torques. Quaternions are used for elements rotations instead of Euler angles for computing efficiency.

Discrete elements can interact only by contact or can be connected by springs or 3D beams. As established in [29], DEM models using cohesive beams to link discrete elements are appropriate to model continuous materials. The idea is to build a hollow sphere made of thousands of discrete elements; see Figure 3. Appendix A gives the methodology to find the microscopic parameters (i.e., physical parameters of cohesive beams) of the model.

The aptitude of the DEM to simulate the fracture of materials is possible by breaking the beams according to a failure criterion. The simplest fracture model is based on the Rankine stress criterion applied for each beam:where is the maximal equivalent Rankine stress, the maximum normal stress, and the maximal shear stress. But, as it was shown before in [42], this criterion failed to simulate accurately the failure of brittle material. A novel failure criterion was recently developed [42], no longer on each beam but on each discrete element. It consists in an equivalent Cauchy stress tensor (similar to what can be obtained with the finite element method) expressed at the center of each discrete element. The stress tensor is calculated taking into account the forces generated by all the beams connected to the element under consideration. It is computed as follows:where(i) is the tensor product between two vectors;(ii) is the equivalent Cauchy stress tensor of the discrete element ;(iii) is the volume of the discrete element ;(iv) is the interaction force between and adjacent DE;(v) is the relative position vector of the center of the two adjacent discrete elements and in a global basis .

The equivalent Cauchy stress tensor can be diagonalized in order to find the principal stresses in the principal basis :

A failure criterion based on the Cauchy stress tensor can be defined, the maximum principal stress criterion, computed here:

In order to take into account the time-temperature superposition principle, the strain rate has to be evaluated on each element. The strain calculation in DEM is not as direct as it can be for continuous modelling (e.g., finite element). The Lagrangian formulation for small deformations is indeed not well adapted in discrete domains. Specific models have been implemented to overpass this problem [43]. Most of these techniques rely on an equivalent discrete-continuum approach. Delaunay triangulation [44, 45] or Dirichlet tessellation [46] has made it possible to compute a local strain in discrete arrangements. The construction of a domain around the elements is nevertheless necessary which increases the computational time. A simpler method is proposed in this paper; it relies on the Cauchy stress tensor. Even if the stress tensor and the strain tensor in discrete domains show a strong duality, it is not complete since expression for both tensors does not belong to the same domain. The strain tensor takes into account the space between elements whereas the stress tensor is calculated for one element. Since only an equivalent strain rate is sought, the true strain is not essential. Similar reasoning was done for the DEM modelling of concrete [47] where an interaction strain between elements is calculated to implement a strain rate dependent model. In that case, it is possible to compute the equivalent strain in each principal direction:

The strain rate is then calculated knowing the strain on each element at each iteration:with the time step.

The failure criterion becomes

The temperature is considered constant in the simulation; no heat conduction is performed.

Once the element considered reaches the stress limit, the element breaks. All its beams are disabled and they do not interact anymore in the simulation. The broken element remains in the domain giving the right microscopic brittle behavior so that the energy conservation is verified. A brittle calibration shall be carried out in the same manner as the elastic calibration to find the limit stresses values.

3. Numerical Calibration

Elastic and fracture parameters are used in the DEM model to reproduce the elastobrittle behavior of the material. These parameters are numerically calibrated according to the properties of the material: Young’s modulus, Poisson ratio, and the yield stress.

3.1. Experimental Compression Tests for Numerical Calibration

In order to calibrate the numerical model, preliminary quasi-static and dynamic compressive tests on the material’s hollow sphere are conducted. Compression cylindrical specimens were provided to determine directly the mechanical properties such as Young’s modulus and the failure force of this type of structure. These specimens were moulded and they underwent the same baking cycle as hollow spheres.

To calculate Young’s modulus strain gauges were placed on cylindrical specimens of 30 mm in diameter and 20 mm in height to measure the true strain. The geometry was chosen to ensure a good measurement of the strain according to the strain gauges used. Young’s modulus is only evaluated in quasi-static test (the crosshead velocity was set at 5 mm·min−1) and at room temperature (°C). A value of  GPa was found which is three times more than the pure epoxy resin. As expected, the reinforcements stiffen significantly the material while reducing the deformation at failure. The composite material showed a brittle behavior when compressed as expected for thermosetting polymer matrices at ambient temperature. After the elastic phase, samples present fractures (surfaces oriented at 90° to 45° as a function of the compression direction is observed, in Figure 6(a)). The elastic curve goes up to the failure of the material and no plasticity seems to occur. The fracture surfaces do not exhibit permanent deformation.

Compressive tests in temperature were conducted on cylindrical specimen of 12.7 mm in diameter and 25.4 mm in height (ASTM D695). Temperatures range from −40°C to +80°C by step of 20°C. The crosshead velocity was set to 5 mm·min−1 ( s−1) and 500 mm·min−1 ( s−1).

High strain rate tests ( s−1) were conducted on Hopkinson’s bars at room temperature. Cylindrical specimens of 10 mm in diameter and 5 mm in height were used.

Results in terms of failure stress as a function of temperature and for different strain rate are presented in Figure 4. The material shows a high sensibility to temperature: failure stress ranges from 50 MPa to more than 200 MPa between −40°C and 80°C for a strain rate of  s−1. The strain rate effect is also important; the failure stress in dynamic  s−1 is more than twice as high compared to quasi-static compression.

In order to predict the yield stress of polymers, physical or phenomenological laws have been studied, especially for thermosetting resin [16, 48, 49]. The original Eyring law considers only the process of pure flow and neglects nonviscous effect (isotherm adiabatic heating or high elasticity). Validated for moderate strain rates [50], this theory is not suitable for very high strain rates. The generalized theory of Ree-Eyring [30] is more adapted and was successfully applied for epoxy resins [51]:with(i) the yield stress,(ii) the absolute temperature,(iii) the strain rate,(iv) the Boltzmann constant,(v) the universal gas constant,(vi) and the activation volumes for each process of and ,(vii) and the elementary shear,(viii) and the rate constants,(ix), , and the activation energies for each process of and .

Ree-Eyring constants are identified by numerical optimization. The Downhill Simplex Method [52] is used in Python to find the constants that will fit to the experimental data. Constants initial values of a similar material [51] are set before optimization. New values found are close to the initial ones and are presented in Table 1. The failure stress as a function of strain rate and temperature according to the Ree-Eyring law with the constants identified is plotted in Figure 5.

3.2. Brittle Calibration

Prior to the brittle calibration, the elastic calibration is conducted according to the methodology presented in Appendix B. This process remains the same for any type of elastic material and was successfully achieved before. Microscopic values found, corresponding to the material studied here, are shown in Appendix B.

Brittle calibration is then carried out. The brittle parameter (stress limit) is scanned over a wide range of values in order to find the right microscopic value leading to the correct macroscopic value. The geometry tested was chosen to be the same as in the experiment (a cylinder with  m,  m is chosen). Firstly, the strain rate effect is not taken into account as only the correspondence between macroscopic and microscopic yield stress values is sought. Therefore, force displacement curves can be directly compared between experimental and numerical results. The sample is compressed between two plates. A very low loading velocity is applied to be in quasi-static condition (kinetic energy is negligible). The top and bottom faces of cylinders have been numerically smoothed so that numerical roughness (which is much larger than the cylindrical specimen at macroscopic scale) does not interfere with the simulation. To permit a crack closure and the sliding between cracked faces, the contact is activated between each element which are linked by a broken beam. The stiffness contact, computed by the penalty method, is equal to the rigidity of the beam in compression:with the beam section, Young’s modulus of the beam, the radius of the beam, the mean radius of the two elements connected, and the radius ratio.

Friction coefficients from 0.25 to 0.5 are encountered in literature for polymers; see [53] for pure epoxy in contact with steel. The failure force value varies from 1% with a fluctuation of the friction coefficient value between 0.3 and 0.5. The authors decide therefore to set a friction coefficient of 0.3 for both contacts between each discrete element (cracked faces) and between discrete elements and the plates. Static and dynamic friction coefficients are supposed to be the same as their values are close. The result with the principal stress criteria is plotted in Figure 7. The macroscopic yield stress is found to be . For example, at °C and  s−1,  MPa and  MPa. Qualitative results can be seen in Figure 6 where elements colour shows their horizontal displacement after failure ( mm). Failure crack paths are well represented by the DEM model: crack surfaces are oriented at 90° to 45°. The force decreases faster after the main failure which is due to the broken elements that are ejected from the crack surfaces. Indeed, the third body created by the fracture is thicker (minimum the size of the discrete elements) and more mobile in the simulation. This induces a less important stress take-up.

Numerical compression tests have been conducted at different velocities to validate the strain rate dependent model (based on the general theory of Ree-Eyring). Cylindrical geometries (12.7 mm in diameter and 25.4 mm in height) with 2000 discrete elements where used. This small number of elements does not give accurate results but allows decreasing the computational time to run low velocities compression tests. Results for both the strain rate dependent and strain rate independent models are presented in Figure 8. Below 10 m·s−1, the failure force is quasi-constant for the strain rate independent model but increases accordingly to the Ree-Eyring law for the strain rate dependent model.

Above 10 m·s−1, the specimen is not in an equilibrium state, inertia effects are not negligible anymore and the failure stress increased for both models. The equilibrium time can be calculated with this equation:with Young’s modulus, the length of the specimen, and the density of the material. A value of μs is found in this case. Considering that the failure is reached at a strain of 1.5%, the velocity threshold occurs at  m·s−1 which is consistent with the results found for both models.

4. Results

4.1. Compressive Tests on Hollow Spheres

Quasi-static (Figure 9) and dynamic compression tests on hollow spheres have been conducted and compared in Figure 10. The dispersion is plotted as the area between the maximum and minimum force curve during compression. As expected, the structure of a hollow sphere presents strain rate dependence. The increase in failure force is more than 50% higher in dynamic compression ( m·s−1). A significant dispersion is observed both in quasi-static and dynamic compression. The standard deviation (SD) for the failure force is 110 N in quasi-static and 160 N in dynamic compression. And for the displacement at failure, it is, respectively, 0.24 mm and 0.20 mm. This dispersion may have two origins: dispersion associated with the accuracy of force and displacement sensors and the dispersion due to geometrical defaults of hollow spheres induced by the manufacturing process. The two first causes can be quantified. In quasi-static compression, the compressive speed is the same and less than accurate. In dynamic compression, the measurement is estimated to be 5% accurate. Geometrical defaults of hollow spheres are important. Indeed, a mass dispersion of about % can be measured and hollow spheres sphericity is not perfect. The failure force is therefore highly dependent of the orientation of the sphere between plates.

The failure process of the sphere remains the same in quasi-static or dynamic compression. It can be described in three stages:(1)Elastic phase: at the beginning of the compression ( of relative longitudinal displacement, where is the displacement of the moving plate and the diameter of the sphere), a linear elastic behavior can be considered as no drop of force occurs.(2)Failure: a great drop of force is observed for which corresponds to rapid crack propagation. The circumferential stress in the hollow sphere’s shell is mainly positive (tensile stress) far from the point of load application when the structure is under uniaxial compression [54]. The opening mode I of the crack is then expected and is confirmed here. Before failure, the structure reaches its maximum force which is on average 1037 160 N in dynamic and 752 110 N in quasi-static compression. The crack propagates vertically (between the two plates) on a half perimeter, in Figures 9 and 11. The curves plotted in the paper only represent the first two phases. The displacement of the plates was stopped after the main failure.(3)Multiple cracks: once the hollow sphere is vertically cracked, the resultant force is close to zero. Then walls of the sphere are bending and buckling causing the proliferation of cracks oriented in every direction.

4.2. Numerical Modelling of Hollow Spheres

Once the hollow sphere was created (see Appendix C for more details), micro-parameters found in the elastic and brittle calibration were used (Table 2). The velocity of the moving plate was set to 1.3 m·s−1, close to the average impact velocity for experimental tests. The principal stress criterion with the strain rate dependence was employed and compared to the experimental results. The force displacement curve of the dynamic experimental data and the numerical result is plotted in Figure 12. The force at failure is 1010 N, close to the experimental results. The displacement at failure is 1.29 mm, 6% more, and the energy at failure is 0.71 J, 9% less than in the experiments.

In addition to quantitative results, the numerical modelling permits to obtain excellent qualitative results. The crack propagates in the same way: the crack is frank, rectilinear, and directed vertically between the two plates (Figure 13). The mean crack tip velocity computed is approximately 200 m·s−1, close to the experimental results ( m·s−1). The position of the crack tip as a function of time was measured directly on experimental and numerical images.

Microcracks are also observed at the top and bottom of the hollow sphere; see Figure 13. It shows that, before the main crack, hollow spheres are subjected to permanent deformation and damage. The numerical model takes into account microcracks and allows quantifying them, one can notice the number of broken elements associated with a small drop of force which is not observable with experimental data, in Figure 14. It is now possible to analyse the position of the microcracks from inside the hollow sphere. Cracks are concentrated at the bottom and top center from inside the hollow sphere but are more spread at the outside face, as one can see the experimental data. Indeed, hollow sphere’s shells are bending and microcracks are observed where tensile stress occurs.

4.3. Validation for High Strain Rate

Eventually, high velocity tests are conducted with an air cannon. Qualitative results of both the strain rate dependent and strain rate independent model compared to the experimental results are presented in Figure 15. In the experiment, 33 μs after impact, several cracks (8 to 10) are observed near the point of impact. Cracks propagate as far as 2/3 of the hollow sphere (66 μs). These scraps of spheres then bend and break (99 μs). The remaining dome crashes in turn. The strain rate independent model fails to reproduce the crushing of the hollow sphere at high strain rate: the hollow sphere breaks instantly when it impacts. On the contrary, the strain rate dependent model takes into account the different phases of crushing.

The impact force has been experimentally measured on a reduced number of tests. Indeed, an impact location outside the center of the force sensor caused nonaccurate results that cannot be exploited. Figure 16 presents the impact force (in red) of a hollow sphere impacting a rigid surface at 70 m·s−1.

The DEM model is then used to reproduce the impact at the same velocity. The noisy computed force is due to the progressive impact of all the discrete elements. For clarity of results, a moving average technique is conducted for the numerical data. Raw results (in transparency) and smooth results are presented in Figure 16. The maximal force is close for both the strain rate dependent model and the experiment (respectively, 3350 N and 3696 N for the smooth results). The decrease of force after the peak is similar until 150 μs after impact. The strain rate independent model has a lower maximal force (2444 N, compared to the experimental results).

The numerical results are in good agreement with the experimental data and justify the use of a strain rate dependent law.

5. Discussion and Conclusions

The fracture behavior of composite hollow sphere under quasi-static and dynamic compression was investigated using an experimental and a numerical approach. Experiments have been conducted at room temperature and at different compressive velocities (5 mm·min−1 to 130 m·s−1). A numerical model based on the discrete element method which takes into account the strain rate dependence is proposed.

Compression tests on composite hollow spheres have shown that the compression velocity influences the response of the material/structure. As expected for polymer materials [51, 55], the more the compressive velocity, the more the force at failure (+50%) and the energy at failure (+30%). The material is subjected to dynamic fracture to dissipate the free available energy stored in the structure: a principal crack propagates at high speed ( m·s−1) and initiates near the sphere summits.

The numerical model using DEM has been developed to simulate the dynamic fracture and to consider the time-temperature superposition principle of the polymer. The constants of the Ree-Eyring law are identified according to compression tests at different strain rates ( s−1 to  s−1) and temperatures (−40°C to +80°C). The strain rate is successfully computed in DEM through the stress tensor. By taking into account the strain rate dependence on the yield stress of the material, the model allows simulating accurately the fracture of the material for dynamic compression, both on a cylindrical specimen and on one hollow sphere. The maximal force numerically measured before the principal crack propagates in the hollow is close to the experimental data (difference of 3%).

The model highlights the phenomena of microcracks coalescence prior to the main failure. It makes it possible to quantify the energy dissipated by fracture before and after the failure of the hollow sphere. Moreover, the DEM shows all its potential for high velocity impacts ( m·s−1 to 130 m·s−1) where multicracks, large displacements, and contacts contribute majorly to the results. It is shown that the use of a strain rate dependent law is essential to properly simulate the ruin of the structure. This study highlights a DEM approach to model a dynamic fracture mechanism with taking into account the time dependency of polymer materials.

Eventually, the objective is to develop a macroscopic model of the hollow spheres crushing. Thousands of hollow spheres will interact between each other. Each sphere will be modelled by one discrete element and the interaction between each element will contain the constitutive law identified with the model presented in this paper. A promising multiscale model will eventually take advantage of the DEM to simulate fracture on the microscopic scale (one to five hollow spheres) and to describe the mechanical behavior of the whole structure on another higher scale. Optimization could be performed at the macroscopic scale by finding the right microscopic parameters to dissipate more energy during an impact.

Appendix

A. Discrete to Continuous Model

The beams have geometrical and material parameters which will be referenced as microscopic parameters. Beam’s length is equal to the distance between two elements linked. The radius ratio (where is the mean radius of the elements linked) is used to calculate the radius of the beam. This dimensionless parameter allows the model to be independent of the number of discrete elements in the simulation. Material parameters include Young’s modulus and Poisson ratio of the beam. Beams are computed using the Euler-Bernoulli beams [56] with all 3D solicitations (traction-compression, torsion, and bending) considered. A macroscopic elastic behavior emerges from these microscopic parameters and permits to describe the material studied.

The microscopic Young modulus and the radius ratio have to be found for each material elastic behavior. The Poisson ratio value of the beam has almost no influence on the macroscopic rigidity [29]. Indeed, torsion deformation is negligible compared to traction, compression, and bending deformations. Its value will be set to . A parametric analysis is conducted to find the microscopic parameters. This elastic calibration phase requires to(i)find the microscopic beam’s radius ratio giving the macroscopic Poisson’s ratio,(ii)find the microscopic Young’s modulus giving the macroscopic Young’s modulus.

The model can work properly in the elastic region of the material as long as the following parameters are checked [29]: the numerical assembly of elements has to be isotropic with no preferable orientation. The cardinal number (average number of bonds per element) has to be close to . The volume fraction () has to be constant along the geometry and roughly equal to for a good representation of the material continuity. For simple geometries, the convergence is obtained with a minimum of elements. All these parameters depend on the creation of the geometry which is performed before the calibration. A special module has been developed to create spherical or cylindrical random compact assemblies of discrete elements in three steps.(i)Stage I: random free filling with discrete elements in the volume considered (box, sphere, etc.) until no more element can be added (without interpenetration with another element).(ii)Stage II: new discrete elements which are inserted one by one. At each iteration a DEM calculus is used to compute the reorganisation of all the discrete elements due to the possible interpenetration with the freshly inserted new discrete element. When the kinetic energy of the discrete domain lowers under a small given threshold, a new discrete element can then be randomly inserted into the discrete domain. Stage II stops when the cardinal number reaches the setting value (the default value is ).(iii)Stage III: relaxing the boundaries of the domain by lowering the stiffness of the walls to decrease the average interpenetration of the elements (smaller than 10−5%).

Since the filling is a random process, all the compact domains created are different for the same geometry. Discrete elements are contained in the same geometry but their arrangement is statistically independent.

B. Elastic Calibration

The first step in the elastic calibration is to ensure that the discrete meshing is fine enough to retrieve the macroscopic elastic properties obtained experimentally. The numerical convergence is assumed when the macroscopic parameters such as Young’s modulus and the Poisson ratio converge to one value, Figure 18. To achieve this convergence, several discrete cylindrical domains from 500 to 35000 elements were created; see Figure 17. Each of them was generated five times in order to use the average for each discrete domain refinement and to quantify the dispersion.

A numerical compressive test on each specimen is then performed. A load force is imposed on each side of the cylinder (on each center of elements belonging, respectively, to the top and bottom set of the cylinder), progressively applied and stabilized. Only elastic properties of the model are investigated so plates are not simulated and also frictions between plates and the cylinder are not considered. The rate of loading forces is chosen to ensure a quasi-static solicitation and therefore to have a negligible kinetic energy compared to the total energy. The macroscopic Young’s modulus and the Poisson ratio are estimated knowing the force applied at the top and bottom and the change in length and radius of the cylinder. According to the charts in (a) and (b) in Figure 18, a discrete element number greater than 15000 is found to get converged: Young’s modulus and Poisson ratio values are within 2.5%. Microscopic parameters are scanned, interpolated, and finally selected according to Section 2.3; see the charts in (c) and (d) in Figure 18 and Table 2.

C. Creation of Discrete Element Hollow Spheres

In order to create representative hollow sphere with discrete elements, a convergence analysis of the elastic phase was performed. A minimum number of discrete elements through the thickness are required to simulate continuous material depending on the geometry. Indeed, for thin structures like hollow spheres, a great number of elements are needed to keep enough elements in the thickness. As it can be observed in Figure 19, the computing time to create this geometry (with a thickness ratio ) varies from 4 s for 1 element in the thickness to 2600 s for 4.2 elements in the thickness. For the geometry of hollow sphere studied (), the computing time is much more important. With the available laboratory computing resource at this time, it is not reasonable to have more than 2 or 3 elements in the thickness. The cardinal number is approximately equal to 5.7 (the optimal value is 6.2 for isotropic representation) for the hollow sphere structure and is hardly influenced by the number of discrete elements through the thickness. A simulation of a thin structure where the analytical result (based on the plate theory) is known has to be conducted for different DEM refinement to validate the use of a limited number of elements.

The simulation of a circular plate composed of discrete elements has been chosen by convenience. The contour of the plate is clamped and the structure is subjected to a normal force at its center; see Figure 20. The deflection of the center of the plate obtained analytically is given in where is the radius of the circular plate, its thickness, and Young’s modulus (N/m2). DEM circular plates were built with a mean diameter of discrete elements ranging from one to one-fifth of the plate thickness. Five geometries for each discretization have been created to account for the random packing. It has been shown that a number of discrete elements in the thickness of 2 allow an accuracy of % (the accuracy increases with the number of elements to reach % for 5 elements in the thickness) on the analytic . It is supposed to be a good compromise for calculation cost. In fact, hollow sphere may be so thin () that it becomes impossible to create compact DEM geometries with more than two elements in the thickness wall. A ratio of two elements in the thickness has been retained to model the hollow sphere.

Disclosure

Part of this work was presented at the 5th International Carbon Composites Conference (2016) in Arcachon, France [57].

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was carried out in the frame of the FUI project SAMBA funded by BPI, Midi-Pyrénées and Aquitaine regions, and approved by the Aerospace Valley. The SAMBA project (for Shock Absorber Material for Bird-shield Application), backed by STELIA Design and Research Division, concerns the development of new generation bird impact shields for commercial and business aircraft. The authors gratefully acknowledge the members of the project: STELIA (AirbusGroup), CEDREM, ESTEVE, NIMITECH, HUTCHINSON, Institut Clément Ader, and especially V. Fascio from ATECA SAS for the provision of hollow spheres and for her contribution.