Modelling polymeric deformable granular materials-from experimental data to numerical models at the grain scale

Polymeric deformable granular materials are widely used in industry and the understanding and the modelling of their shaping process is a point of interest. This kind of materials often presents a viscoelasticplastic behaviour and the present study promotes a joint approach between numerical simulations and experiments in order to derive the behaviour law of such granular material. The experiment is conducted on a polystyrene powder on which a confining pressure of 7MPa and an axial pressure reaching 30MPa are applied. Between different steps of the in-situ test, the sample is scanned in an X-rays microtomograph in order to know the structure of the material depending on the density. From the tomographic images and by using specific algorithms to improve the images quality, grains are automatically identified, separated and a finite element mesh is generated. The long-term objective of this study is to derive a representative sample directly from the experiments in order to run numerical simulations using a viscoelactic or viscoelastic-plastic constitutive law and compare numerical and experimental results at the particle scale.


Introduction
The use of polymeric granular materials having a deformable behaviour is of common use in the industry.The behaviour of such materials, when mechanically loaded under relatively large stresses, and at a given temperature and moisture content, is ruled by grain kinematics but also, in a very significant manner, by the deformability of the grains.
In the present study, an analysis of the mechanical behaviour under mechanical loading is developed, using in complementary ways, X-rays tomography and in situ triaxial compression test.The aim is to be able to observe each grain of the experimental sample at each step of the experiment.The approach developed here consists in generating digitized grains from the images obtained through tomography, meshing these grains and integrating these meshed geometries into a finite element model.After applying consistent boundary conditions and mechanical loading, a grain-by-grain comparison of indicators such as contact areas and their orientations should allow an accurate analysis of mechanical behaviour of an average set of grains by using a 3D multi-particle finite element method (MP-FEM).That method has been already used on deformable materials [1][2][3].Some studies also took into account the same kind of material, but they used only 2D models [4,5].

Material
The deformable granular material used in that work is a polymeric powder composed of polystyrene grains supplied by Goodfellow (UK).The size of the grains is between 150μm and 500μm and their polyhedral shape is very random due to the freeze-fracturing process (see paragraph 2.2).

Compression
The experiment performed in this study was a triaxial compression test.Three tests were performed at different confining pressures: 1, 2 and 7MPa.The results presented here were obtained by applying a constant confining pressure of 7MPa while the axial stress reached 30MPa.The tested powder (approximately 2.5g) was poured into a flexible membrane made of latex, having a thickness of 400μm in average and a diameter of approximately 11mm.The length of the sample before loading was approximately 22mm.The confining pressure was applied using a hydraulic pump and held to its set point with an accuracy of ±0.05MPa.The axial pressure was applied by a displacement-controlled piston and was measured by a 5kN force sensor.The compressive load was applied in subsequent 700μm displacement steps at a velocity of 2μm/s.Figure 1 shows the typical curves of the experiment.X-ray tomography scans are performed at each step,

Microtomography
The scans were performed on the 3SR-lab microtomograph (manufactured by RX Solutions, with a Hammatsu L12161-07 source, and a Varian flat panel detector).The pixel size was set to 9μm and the beam was obtained by applying a voltage of 80kV and a current of 125μA.3D images were captured using 1440 projections per scan and averaged with 8 images per projection.Figure 2 shows the quality of a non-processed image.

Density and mechanical analysis
The deviatoric stress is plotted on figure 1 as a function of the axial strain.The figure also shows the almost constant increase of the mean relative density of the sample under compression and, thus, the decrease of the volumetric strain.The measurement of the relative density was made by thresholding the 3D images and calculating the ratio of the volume of voxels corresponding to the solid phase and the total volume.This process was done on subvolumes, the size of which was increased until the measured density stabilised so as to define the size of a representative element.Then, the determination of the relative density was calculated on all the subvolumes of this size within the sample.The volumetric strain was calculated as: Where V 0 t , V 0 g and D 0 are, respectively, the total volume, the volume of grains and the density before loading; V t , V g and D the same quantities at a given compression step.It is considered, here, that V g equals V 0 g .It is observable that, when using a confining pressure of 7MPa, the deformable behaviour of the grains allows the sample to reach a relative density of 77% with an axial strain of 25%.No dilatancy phenomenon is observable, which is explained by the ability of the particles to adapt to the load by changing their shape.

Image processing
As shown in figure 2, the quality of the reconstructed images was not suitable for direct use.Indeed, the contrast between the grains and the porosity was poor and a significant noise was visible.Figure 3 shows the histogram (number of voxels as a function of the grey scale value) of the original volume.The presence of two peaks was due to the bimodal character of the sample (solid material and voids).The width of the peaks was mainly due to the experimental noise (visible on the corresponding picture).To improve the contrast between the two phases, the noise was reduced by applying a blurring filter (3D total variance).The effect of the filter was to achieve a better separation between the peaks and to improve their magnitude as seen on the lower histogram.As a result, one clearly observes the decrease of noise inside the grains (dark dots) and outside (bright dots) on the lower image.The signal being much closer to a binary image, the determination of the threshold in the processing algorithm was made easier and more accurate.Note that if the blurring filter efficiently reduced the noise, it also blurred the edges of the grains.

Segmentation
The aim of the segmentation was to identify and separate the grains.Grains with realistic shapes and an accurate contact detection can be achieved by this means.For this paper, the segmentation was first tested on the two grains presented on figures 3 and 4. The segmentation algorithm was based on the diffusion of markers.The script uses the random walker function from  the python module of scikit-image [6].The first task consisted in selecting what is solid material and what is void.An automatic threshold was chosen for that purpose.Two methods were tested: Kittler [7] and Otsu [8] methods.Even if both of them gave realistic and adequate binary images, they led to slightly different values of threshold.To take advantages of both Otsu and Kittler methods, an averaged threshold was defined (equation 2): where T , T Otsu , T Kittler are respectively the obtained threshold value, the one from Otsu and the one from Kittler.The variable α is a parameter chosen between 0 and 1 in order to get images as realistic as possible on a visual aspect.In most cases α was approximately 0.4.
Once the images were binarised, the residual noise had to be removed.To do that, opening (erosion and dilation) and closing (dilation and erosion) processes were performed to remove all the black or white dots that were isolated.It was then necessary to divide each image into three categories: -the "sure background", obtained by dilation of the solid material; -the "unknown region", which is not the "sure foreground" or the "sure background".
Every marker in the sure foreground was labelled and the random walker was applied in a specific order: the labels were growing in the unknown region depending on the neighbourhood and, then, the background propagated toward the grains.The propagation rate depends on the blurred original images and a given "penalization coefficient".Some inconsistencies appeared near the contact areas depending on the order of the labelled markers and the distance between the markers and the contact areas: the boundary of one grain could be taken as the boundary of another grain.To avoid that, an opening process was made on the isolated grains, so that the very small defects around the grain were removed.The surfaces of the new grains were then smoother and the volume could be considered unchanged.
Contacts were characterized as the borders between labels after propagation of the markers.

Toward a mechanical discrete analysis
As mentioned in paragraph 3.2, it was possible to obtain the volume of the grains from the tomography images.The aim is to introduce a certain number of grains (approximately 5×5×5) into a finite element code.For the purpose of using the numerical grains, it was first needed to mesh them.
The mesh generation was performed thanks to the free matlab/octave-based toolbox iso2Mesh [9], partially based on the software project CGAL 3.5 3D mesher [10].The grains were meshed by using the so-called "cgalmesh" method, which allowed to mesh the surface, to rearrange the elements of the surface and to mesh the volume.This method also gives the opportunity to manage the size of the triangular surface faces by choosing the radius of the Delaunay spheres.The maximum volume of the tetrahedral elements can also be chosen in order to manage the volumetric density of elements.Figure 6 shows the mesh generation of the two grains studied in the previous paragraphs when using a radius of the Delaunay spheres of 13.5μm and a maximum volume of 729μm 3 .Considering these parameters, slightly more than 30500 nodes and 13500 tetrahedral quadratic elements were sufficient to accurately mesh the two grains.
The digitized and meshed neighbouring grains will be introduced into a MP-FEM code where a viscoelastic or viscoelastic-plastic constitutive law will be considered for the material behaviour and a Coulomb friction model will be used.Periodic boundary conditions or flexible membranes will be chosen as boundary conditions.Contact areas and their orientations should be efficient indicators to compare the real and the numerical granular material at a same density.The analysis of the displacement or strain fields will also be taken into account.

Conclusion
This study presented a methodology designed for modelling the mechanical behaviour of a granular material using X-ray tomography as an experimental characterization tool at the grain scale.Indeed, from the 3D images of the grain assembly structure, it is possible to detect and isolate the grains and to generate a 3D mesh for each grain.
A region of meshed particles can then be considered and introduced into a finite element code to simulate its mechanical behaviour.
To this end, we used already existing tools or methods, in order to remove as much as possible the experimental noise without damaging the borders of the grains and to perform an accurate detection of the contact zones between the grains.The next step will consist in studying the mechanical behaviour of the granular material by using, in a MP-FEM code, a certain number of meshed grains and applying a load comparable to the one obtained with the triaxial compression test.The constitutive law will be based on viscoelastic or viscoelastic-plastic models.The comparison between the real region (observed in microtomography) and the numerical one (resulting from simulation) should allow to understand the mechanical behaviour of the granular material.

Figure 1 .
Figure 1.Deviatoric stress, volumetric strain and mean density (with 3D density maps) measured during the triaxial compression test as functions of the axial strain

Figure 3 .
Figure 3. 3D blurring filter effect on images and histograms

Figure 4 .
Figure 4. Evolution of two grains during compression

Figure 5 .
Figure 5. Identification and segmentation of two grains

Figure 6 .
Figure 6.Mesh generation of the two grains with iso2mesh