Karthik Ram Ramakrishnan
email: karthik.ramakrishnan@eng.ox.ac.uk

Stéphane Corn
email: stephane.corn@mines-ales.fr

Nicolas Le Moigne

Patrick Ienny

Pierre Slangen

Experimental assessment of low velocity impact damage in flax fabrics reinforced biocomposites by coupled high-speed imaging and DIC analysis

Keywords: A, Natural Fibers A, Polymer-matrix composites (PMCs) B, Impact behaviour C, Damage mechanics

Structural composite materials with good impact resistance are required in several industrial sectors to design parts submitted to crashes or falling objects. This work investigates the impact behaviour of flax fabrics-epoxy and commingled flax-polypropylene (PP) fibres composite plates manufactured using compression moulding process. An instrumented drop tower was used to conduct low velocity impact tests at several energies and highspeed imaging coupled with DIC analysis was used to assess the damage evolution. A crack-tracking algorithm was implemented in ImageJ to identify the initiation and propagation of cracks. A comparison of the forcedisplacement histories and high-speed camera results for both composites were used to identify the different modes of damage and the critical energy for complete penetration. The proposed method is able to depict the successive events that lead to the breakage of biocomposites and it will be a valuable tool to investigate the impact behaviour of other materials as well.

Introduction

Biocomposites with natural fibre reinforcements are considered a potential replacement to traditional structural composites due to their low density, high specific properties, low cost of raw material, relative abundance, and positive environmental profile [START_REF] Shah | Can flax replace E-glass in structural composites? A small wind turbine blade case study[END_REF]. Bledzki and Gassan [START_REF] Bledzki | Composites reinforced with cellulose based fibre[END_REF], Faruk et al. [START_REF] Faruk | Biocomposites reinforced with natural fibers: 2000-2010[END_REF] and Pickering et al. [START_REF] Pickering | A review of recent developments in natural fibre composites and their mechanical performance[END_REF] have published extensive review articles on the mechanical behavior of polymer composites reinforced with natural fibres. In particular, it has been reported that composites reinforced with plant fibres from flax plant (Linum usitatissimum) in thermoplastic and thermoset matrices have exhibited good specific mechanical properties comparable to those of glass fibre composites [START_REF] Yan | Flax fibre and its composites -A review[END_REF][START_REF] Pil | Why are designers fascinated by flax and hemp fibre composites? Compos[END_REF]. However, a lack of data on durability properties of these composites for loading conditions such as fatigue and impact has limited the use of natural fibres for high performance applications [START_REF] Bensadoun | Low velocity impact properties of flax composites[END_REF]. The biocomposite structure may encounter impact loads during their service life due to dropped tools, collisions with loading and unloading fixtures, ballast or cargo and the internal damage produced by such impacts may severely reduce their residual mechanical properties. The difference in impact resistance between biocomposites and traditional composites was investigated by Bos et al. [START_REF] Bos | Mechanical properties of short-flax-fibre reinforced compounds[END_REF] using unnotched Charpy tests and it was reported that the impact strength of short flax fibre compounds were much lower than that of glass fibre. In contrast, Puech et al. [START_REF] Puech | Investigating the impact behaviour of short hemp fibres reinforced polypropylene biocomposites through high speed imaging and finite element modelling[END_REF] showed that PP-hemp biocomposites exhibit higher absorbed energies (up to 40%) than PP-glass composites due to higher strain at break. Recently, several researchers have studied the impact resistance of biocomposites with non-woven short flax fibres in thermoplastic matrix [START_REF] Hoang | Non-woven flax fiber reinforced polypropylene : mechanical properties in static and low velocity impact behaviours[END_REF] and thermoset matrix [START_REF] Habibi | Influence of low-velocity impact on residual tensile properties of nonwoven flax/epoxy composite[END_REF]. However, the use of natural fibres in structural applications requires higher reinforcement efficiency that can only be achieved with twistless long continuous fibres. The literature on impact resistance of composite laminates with long continuous natural fibres is still limited. There have been some studies on woven hemp fibres [START_REF] Scarponi | Low-velocity impact behaviour of hemp fibre reinforced bio-based epoxy laminates[END_REF], jute fibres [START_REF] Dhakal | Influence of temperature and impact velocity on the impact response of jute/UP composites[END_REF], and unidirectional and cross-ply flax fabric composites [START_REF] Liang | Impact behaviour of flax/epoxy composite plates[END_REF][START_REF] Ravandi | Low velocity impact performance of stitched flax/epoxy composite laminates[END_REF][START_REF] Sy | Damage evolution in unidirectional and cross-ply flax/epoxy laminates subjected to low velocity impact loading[END_REF]. Liang et al. [START_REF] Liang | Impact behaviour of flax/epoxy composite plates[END_REF] found that the fracture mechanism in flax/epoxy composites included delaminations that occurred at low energy level, followed by the development of intra-laminar transverse cracks and macrocracks due to fibre failure. Ravandi et al. [START_REF] Ravandi | Low velocity impact performance of stitched flax/epoxy composite laminates[END_REF] investigated the low velocity impact response of continuous flax/epoxy composite laminates and found that through-thickness stitching had an adverse effect on the impact resistance as it facilitated the propagation of in-plane cracks. Sy et al. [START_REF] Sy | Damage evolution in unidirectional and cross-ply flax/epoxy laminates subjected to low velocity impact loading[END_REF] compared the impact performance of unidirectional and cross-ply flax/epoxy laminates and found that the cross-ply laminate had higher penetration threshold energy and impact toughness compared to the unidirectional laminate. One of the limitations of these studies is that there is very limited literature on the comparison of the impact resistance of long continuous natural fibre reinforced composites with thermoset and thermoplastic matrix. In this paper, a comprehensive study on the impact behaviour of woven flax fibre-reinforced composites with epoxy and polypropylene matrix is conducted.

The experimental characterization of impact behaviour of composites is typically conducted using a drop tower. The force-displacement curve obtained from the impact experiment is used to identify the evolution of damage as the creation of damage or cracks results in the sudden loss of stiffness of the specimen and a corresponding drop in the force history. Non-destructive inspection methods such as visual inspection and ultrasonic C-scan are used for the assessment of damage but it is typically conducted after the impact event. For instance, Scarponi et al. [START_REF] Scarponi | Low-velocity impact behaviour of hemp fibre reinforced bio-based epoxy laminates[END_REF] used Ultrasonic C-scan to evaluate the damage area in lowvelocity impact experiments of plain weave hemp fabric reinforced bio-based epoxy composites. Dhakal et al. [START_REF] Dhakal | Influence of temperature and impact velocity on the impact response of jute/UP composites[END_REF] used acoustic emission technique to characterize the failure modes of woven jute fibre reinforced unsaturated polyester composite. Santulli [START_REF] Santulli | Falling weight impact damage characterisation on flax / epoxy laminates[END_REF] developed a method to assess the impact damage in flax-epoxy composites using a combination of hysteresis cycles of the force-deflection curves and optical microscopy. However, the force-displacement data is rarely correlated to the actual creation and propagation of damage. Puech et al. [START_REF] Puech | Investigating the impact behaviour of short hemp fibres reinforced polypropylene biocomposites through high speed imaging and finite element modelling[END_REF] conducted such study using high-speed imaging of the impact response of short hemp fibres reinforced PP. It was found from the video tracking analysis that for a given cumulated crack length, PP-hemp composite absorbed much more energy than PP-glass, which was assumed to be related to differences in failure mechanisms [START_REF] Puech | Investigating the impact behaviour of short hemp fibres reinforced polypropylene biocomposites through high speed imaging and finite element modelling[END_REF]. Flores et al. [START_REF] Flores | Highspeed 3D digital image correlation of low-velocity impacts on composite plates[END_REF] have shown that high-speed imaging can be used to record strains and deformation as well as crack initiation and propagation during low velocity impact of carbon fibre composites. It is clear that high-speed imaging of the dynamic impact event should provide important information about the impact damage of composites.

Digital Image Correlation (DIC) method is a powerful and flexible optical metrology technique based on digital image processing and numerical computing for the full-field measurements of displacement or strain [START_REF] Pan | Digital image correlation for surface deformation measurement: historical developments, recent advances and future goals[END_REF]. Typically, image analysis techniques such as DIC were limited to mostly quasi-static and low strain rate applications for the measurement of surface deformations and has been widely used for the identification of materials constitutive parameters from mechanical tests such as tensile tests. Analysis of stresses in real time during a transient failure event, such as dynamic crack initiation and growth, was challenging due to a combination of spatial and temporal resolution demands such as synchronization difficulties and poor image resolution at high rates [START_REF] Kirugulige | Measurement of transient deformations using digital image correlation method and high-speed photography: application to dynamic fracture[END_REF]. Recent improvements in high speed camera technology have allowed the use of DIC of high-speed images in dynamic event testing [START_REF] Tiwari | Application of digital image correlation in impact testing[END_REF][START_REF] Reu | The application of high-speed digital image correlation[END_REF][START_REF] Pan | Full-field transient 3D deformation measurement of 3D braided composite panels during ballistic impact using singlecamera high-speed stereo-digital image correlation[END_REF][START_REF] Cuynet | Damage characterisation of flax fibre fabric reinforced epoxy composites during low velocity impacts using high-speed imaging and Stereo Image Correlation[END_REF]. Reu and Miller [START_REF] Reu | The application of high-speed digital image correlation[END_REF] provided an overview of the applications of high-speed DIC used for many experiments including fracture and impact testing. Recently, Pan et al. [START_REF] Pan | Full-field transient 3D deformation measurement of 3D braided composite panels during ballistic impact using singlecamera high-speed stereo-digital image correlation[END_REF] and Cuynet et al. [START_REF] Cuynet | Damage characterisation of flax fibre fabric reinforced epoxy composites during low velocity impacts using high-speed imaging and Stereo Image Correlation[END_REF] have also demonstrated the application of high speed imaging for full field measurement in impact testing of composites using a single camera stereo-DIC system and two cameras system, respectively.

The aim of this work is to study the impact behaviour of flax fabrics reinforced biocomposites and to assess the evolution of damages using high-speed imaging coupled with DIC analysis. The high-speed imaging enables the user to better depict the successive events occurring upon the crack propagation. A method to correlate the force-displacement measurements to the cracks propagation during impact test is proposed.

Materials and methods

Manufacturing of composite specimen

Biotex Flax-PP supplied by Composites Evolution, UK was chosen for the thermoplastic composite. Flax/thermoplastic fabric consists of twistless flax fibres and polypropylene (PP) fibres mixed in the ratio of 40/60 (v/v). The commingled fibres are then woven in a balanced 2x2 twill architecture to form fabrics with areal density of 400 g/m 2 [START_REF] Ramakrishnan | Optimized manufacturing of thermoplastic biocomposites by fast induction-heated compression moulding: Influence of processing parameters on microstructure development and mechanical behaviour[END_REF]. For the thermoset composite, dry flax woven fabrics (Biotex Flax, 400 g/ m 2 ) also supplied by Composites Evolution and a two-part epoxy resin (1800 epoxy resin and 1805 aliphatic amine hardener) supplied by Resoltech Advanced Technology Resins, France were chosen. The resin and the hardener are mixed in the ratio of 100:17 by mass. 10 plies of the flax-PP fabric were used for the thermoplastic composite (flax-PP) while 5 plies of dry flax fabric were used for the thermoset composite (flaxepoxy), so that both biocomposites plates have comparable areal density of flax fibres and thicknesses.

The first step in the manufacturing of the flax-PP biocomposite is the

Experimental setup for low velocity impact testing

Low velocity impact testing with velocities below 10 m/s is typically accomplished using a drop tower. Applications for flax composites in the transport industry are expected to involve this velocity regime. The tests were conducted based on standard test method for drop weight impact testing of composites (ASTM D7136) using a CEAST 9350 instrumented drop tower (Instron, UK) shown in Fig. 1. A hemispherical impactor of diameter 20 mm is secured to a carriage assembly weighing 3.1 kg. A composite target plate is clamped in a fixture with 40 mm internal diameter. The carriage is released from an assigned height and falls along guideposts and impacts the target.

A Phantom V2512 high-speed camera (Vision Research, Inc., USA), coupled with long distance macro-zoom lens, was used for slow motion observation of the bottom surface of the composite plates during the impact event. The scheme shown in Fig. 1 b present the camera view and the side view of the experiment. The target plate is placed on a sample holder at the top of a U shaped column and the camera is set up facing and centered on the column. The sample holder has a circular hole to let the impactor go through the sample after rupture and this allows uninterrupted view of the bottom surface of the target plate. A mirror setup placed at an angle below the composite plate allowed the camera to safely capture the impact event and also to shine about 15,000 lm from superLEDs cold light device (VERITAS Constellation 120, USA). The height of the U shaped column is fixed and is defined by K while the height of the center of the camera is given by H, and M is the mirror height placed inside the U shape column. The mirror height M is kept fixed, as it also enables 15000 lx superlum array to shine on the object for short exposure time. The camera is tilted to the correct angle. This angle is reached when the maximum contrast is reached on the whole image (planar sharp focus). Some accurate adjustments of pan, tilt, rotation and slide can be necessary to optimize contrast. As H > M, then the mirror angle is adjusted to reach H at the working distance defined by the macrolens zoom (taking into account the reflection path). Optical calibration is performed with optical caliper reference. The scale factor or magnification, which depends on the focal length and the distance between the camera and object, was fixed to 69 μm per pixel. Since the high speed camera coupled with the macro zoom lens has a mass of about 8 kg, a Manfrotto heavy weight tripod 475B (Manfrotto Via Valsugana, 100-36022 Cassola (VI) Italy) was used to ensure stable set-up. This coupling is done by using friction control Pan-Tilt-Rotation Manfrotto tripod head, able to carry the camera weight. The quality of the images has to be perfect to ensure DIC results and so the set-up is improved by adding fine coarse translation Manfrotto 454 micropositioning sliding plate between the 3D head and the camera. The images were acquired at a frame rate of 39,000 fps and exposure time of 25 μs with a resolution of 786 × 786 pixels. Photoelectronic switches were used to measure the initial velocity (and rebound velocity) of the impactor and also acted as the trigger for the image acquisition in rolling buffer FIFO memory mode.

Image sharpness of the rear surface of the target is determined by maximizing the contrast of the live images while carefully sliding the optical head. In order to create a high contrast zone for the purpose of skeletonisation of the images to monitor the propagation of cracks during the impact of the composite plates, their bottom surface was painted with a thin layer of white powder. Additionally, a black paint was sputtered over the layer of white paint on the bottom surface of one set of the flax composite plates to give it a granular random fly-speckled pattern as shown in the (Fig. 2), which is required for the "pattern matching" based DIC measurements.

The drop height was varied from 0.1 m to 0.3, 0.35, 0.4, and 0.5 m to give a range of kinetic impact energies of 3 J to 15 J. The initial energies were chosen to represent a wide range of damage intensities. The upper bound energy was chosen in order to get experimental results up to the full perforation of the plate. Three specimens were tested at each impact condition to ensure the repeatability of the experiments. The energy absorbed by the composite plate was calculated from the area of the hysteresis loop by integrating the force-displacement curve. The total energy is defined as the energy at peak displacement (i.e. without considering the rebound part of the curve). Recovered energy can then be calculated as the difference between the total energy and absorbed energy. It should be noted that not all of the initial kinetic energy of the impactor is available to the plate and some of it is lost in other dissipative losses. Post-impact assessment of the surface deformations and failure mechanisms in the impacted samples were evaluated using an optical 3D measurement system, InfiniteFocus G5 by Alicona (Alicona Imaging Gmbh, Austria). Damage assessment was also conducted by observing cross-sections of the flax-PP and flax-epoxy composite sample using optical microscope. The impacted sample was embedded in coloured resin to preserve the deformed shape and the damage and then the cut section was mounted in a clear epoxy resin. The mounted sample was polished with successively finer abrasive papers and finished with diamond polishing paste.

Digital image correlation of high-speed camera images with CinEMA

Digital image correlation (DIC) analysis was conducted on the images obtained from the high-speed camera. DIC is based on the principle of comparing small zones (or subsets) of the digital images of the specimen surface coated with a random pattern. The image correlation is performed by matching subsets from the reference or undeformed sample and subsequent deformed states to compute full-field displacements and strains. In the 2D DIC of the x-y plane corresponding to the rear surface of the plate, the in-plane displacement fields designated by Ux(x,y) and Uy(x,y) correspond to the displacement that maximizes a similarity criterion between reference and deformed subsets. The image post-processing is based on a three-step DIC approach developed inhouse called CinEMA. The definition of the parameters used for the DIC and the details of the algorithm is described in [START_REF] Ivanov | Strain Localisation in iPP/MWCNT Nanocomposites Using Digital Image Correlation: Strain Localisation in iPP/MWCNT Nanocomposites[END_REF][START_REF] Christmann | Mechanical behaviour at large strain of polycarbonate nanocomposites during uniaxial tensile test[END_REF]. In the first step, the rigid body displacement is calculated by a 2D cross correlation of the selected subimages. In the second step, the DIC processing leads to an incremental displacement evaluation. In-plane displacement fields Ux(x,y) and Uy(x,y) were evaluated at each point of the virtual mesh in the region of interest (ROI) of the sample including the crack zone. The peak of the correlation function was detected at a subpixel accuracy by interpolation and this process is repeated for the entire image to get fullfield in-plane displacements. In the third step, the calculation of the displacement gradient around the considered point allows the evaluation of the in-plane Green-Lagrange strain tensor E.

Crack-tracking algorithm of high-speed camera images

High-speed video frames of the bottom surface of a flax-PP plate for a 6 J impact is shown in Fig. 3. The initiation and propagation of macrocracks at different moments of the impact are observed. It can be seen that at time of 0.77 ms after impact, there is localized damage beneath the point of impact. Radial cracks begin to propagate mainly along the two fibre yarns directions, i.e. 0 • and 90 • , corresponding to x and y axes.

The DIC analysis is not suitable when discontinuities such as large macroscopic cracks occur, as the measured displacements are no longer associated with the deformation of the material. In order to overcome this limitation, a crack-tracking method is proposed. A macro was written in ImageJ software [START_REF] Ferreira | ImageJ User Guide[END_REF] to detect cracks and measure the cumulative length from each image. The crack-tracking algorithm involves several steps as illustrated in Fig. 4. In the first step, an undamaged reference image is converted to 8-bit and a combination of median and Gaussian blur filters are applied to the image. Next, each image from the high-speed camera was treated with similar filters and subtracted from the reference image to isolate the cracks. The image is converted to mask and binarised and the cracks are clearly visible as black in a white background. A skeletonisation algorithm is then used after different steps of dilation and erosion to highlight the cracks and to trace the crack length with 1-pixel wide lines. The total length of the macro-cracks is measured at each time interval during impact and the evolution of damage observed in the plate was correlated to its force-displacement response. A composite image with the cracks overlaid on the original image is then created to verify that the tracked cracks correspond to the actual damage.

Experimental results and discussion

Drop tower impact testing

The sensor instrumentation in the drop tower allows the measurement of force and displacement history of the impactor during impact. during a 3 J impact. The force curve has an asymmetric bell-shape corresponding to the loading and unloading of the plate. The force response starts with the linear increase up to 1000 N which is followed by a nonlinear material behaviour up to peak force of approximately 1550 N. After the peak is reached, there is an unloading phase due to rebound or fracture of the composite. Another significant feature observed from the curves is that the tests are highly repeatable with the three samples showing nearly identical time-histories of force. Fig. 5 b shows the history of the mechanical energy calculated from integrating the force--displacement data. It can be seen that the kinetic energy of the impactor is transferred to the plate during the loading phase of the impact. After the peak is reached, there is a rebound phase where the energy is returned to the impactor. In an ideal elastic impact case, there is a symmetric unloading and all the energy is returned to the impactor. However, it can be seen that there is some energy dissipation and not all the energy transferred to the plate is returned to the impactor. This energy gap results from plastic deformation, viscous dissipation and various damage mechanisms in the composite.

Fig. 6 shows the typical force-displacement curves for impact at different initial impact energies for both the flax-PP and flax-epoxy plates. The impact energy varies from 3 J corresponding to low damage up to full penetration impact at 15 J. A representative sample is shown for each energy as the earlier figure showed that the response is highly repeatable. The initial linear part of the curve is identical for all the loading conditions while the figure shows that the residual displacement after impact increased with impact energy. For impact at 15 J, there is complete perforation of the plate and there is virtually no rebound. The tail in the force response for the perforated sample of flax-epoxy is a Fig. 9. Visual and profilometer image of post-impact deformation in flax-PP composite after 6 J and 12 J impact. Fig. 10. Visual and profilometer image of post-impact deformation in flax-epoxy composite after 6 J and 12 J impact. consequence of the residual contact force at the end of impact due to friction between the impactor and the hole in the plate.

A comparison of the typical force-displacement data for a flax-PP and a flax-epoxy composite plates for 6 J and 12 J impact energies is shown in Fig. 7. It can be seen that the two materials exhibited very similar response, suggesting a fibre-dominated behavior of the composites. However, at 6 J, the initial stiffness of the flax-epoxy plate is marginally higher than that of the flax-PP composite due to higher stiffness of the epoxy matrix compared to PP. Besides, it is also evident from the curves that the flax-PP composite has higher maximum displacement for both load cases, which is due to the ductile behaviour of the PP matrix when compared to epoxy.

The absorbed energy calculated from the integration of force--displacement curves are shown in Fig. 8. It can be seen that the absorbed energy increases almost linearly with increasing initial impact energy. However, the recovered energy starts to decrease if perforation failure is initiated, reaching zero for complete penetration of the composite sample. We can estimate the critical energy for penetration failure based on the recovered energy as approximately 9 J for the flax-PP and flax-epoxy composites. A comparison of the energy absorbed by the flax-PP and flax-epoxy composite for different initial impact energies shows that the energy absorbed by the two materials are very similar, with the flax-PP composite absorbing marginally more energy especially for the 3 and 6 J impacts. At higher energies, the two composites show a similar damage response with fibre failure and plate penetration.

Visual assessment of impact damage

A comparison of the visual inspection of the bottom surface of the flax-PP composite samples impacted at energies of 6 J and 12 J is shown in Fig. 9. The profile and topography of the indentations produced by the impact were analyzed by measuring the depth and width of the surface deformations from the Alicona profilometer. Since the polypropylene matrix is a ductile material, extensive deformation is observed in the flax-PP composite as evidenced by the bulged surface and fibre pull-out in the profilometer images. The peak residual deformation of 1.8 and 7.5 mm are observed for the plates impacted at 6 J and 12 J respectively. At lower energies, there is matrix cracking around the point of impact and as the energy of impact increases, fibre fracture forms a crack that is aligned with the fibre directions of 0 and 90 • . There is extensive fibre failure of the composite for initial energy of 12 J for the flax-PP composite and it can be seen that the crack nearly reaches the edge of the Similar analysis of the flax-epoxy composite is shown in Fig. 10 and the profilometer data for the 6 J impact shows that due to the brittle behaviour of the epoxy matrix, the permanent deformation is not as pronounced as in the thermoplastic composite. For the 6 J impact, the peak residual displacement is only 1.1 mm. The fibre fracture at the centre of the impact zone is similar to the flax-PP composite. It can been seen that there is complete penetration of the flax-epoxy composite for the 12 J impact.

The depth profiles of the rear of the impacted plate obtained from the profilometer is compared for the 6 J and 12 J impacts of flax-PP and flaxepoxy composites in Fig. 11. Since the deformation of the plate is symmetric, the depth is plot along the central line of only one axis. It can be seen that the depth is highest at the center of the impacted area (20 mm), with 0 and 40 mm representing the clamped edge. The overall depth profile is similar for the flax composites with PP and epoxy matrix but it can be seen that the epoxy composite has less permanent deformation than the ductile PP matrix for the 6 J impact while at the higher energy there is complete perforation in both plates seen as discontinuities in the depth profile.

Optical microscope images at different magnifications of the crosssection of the flax-PP composite impacted at 12 J are shown in Fig. 12. It can be seen that there is a large impact crater and fibre fracture within the composite. Higher magnification images show the presence of matrix damage in addition to the fibre failure with evident matrix cracks through the flax fibre bundles. It is known that the limited wettability of natural fibres towards apolar and viscous matrices such as PP leads to poor impregnation and the formation of micro and macropores. Ramakrishnan et al. [START_REF] Ramakrishnan | Optimized manufacturing of thermoplastic biocomposites by fast induction-heated compression moulding: Influence of processing parameters on microstructure development and mechanical behaviour[END_REF] calculated the porosity volume fraction from X-ray computed tomography analysis and reported that the pressure and temperature of the manufacturing process have a strong influence on the formation of pores. It has also been reported [START_REF] Moigne | Surfaces and Interfaces in Natural Fibre Reinforced Composites Fundamentals , Modifications and Characterization[END_REF] that surface treatments can be successfully used to increase the interactions between natural fibres and polymer matrices, and significantly enhance the interfacial adhesion leading to better stiffness and possibly modifying impact strength properties.

In the case of the cross-section of the flax-epoxy composites in Fig. 13, it can be seen in the higher magnification image that transverse matrix cracks (highlighted by arrow) lead to delamination type failure between the plies of the composite, which was not visible in the case of flax-PP composite. This is consistent with the complex fracture mechanisms observed in the micrographs of impacted flax epoxy plates [START_REF] Liang | Impact behaviour of flax/epoxy composite plates[END_REF], where a combination of delamination and transverse cracking were reported near the specimen rear face. This highlights a more cohesive behavior for flax-epoxy composite with crack propagation perpendicular to the impact direction, i.e. delamination. It was also reported that the damage in the rear face of flax/epoxy laminates is principally fibre breakage with limited delamination suggesting a fibre-dominant rather than matrix-controlled damage [START_REF] Sy | Damage evolution in unidirectional and cross-ply flax/epoxy laminates subjected to low velocity impact loading[END_REF].

Results of digital image correlation analysis

The displacement contours and strain distributions obtained from the Digital Image Correlation for the flax -PP composites impacted at 6 J are shown in Fig. 14. The image is at time of 0.54 ms from the beginning of the impact (21st image), which corresponds to the initiation of the cracking as clearly identified by the DIC program. It can be seen that the Ux and Uy displacement fields are qualitatively similar showing that the 2D-DIC is able to compute the full-field displacement of the plate. The strain contours show the initiation of cracks as strain localization and it Fig. 13. Optical microscope images of the cross-section of flax-epoxy composite plate in the damage zone after 12 J impact. can be seen that at the point of initiation, the cracks are not symmetric in the × and y directions.

Fig. 15 shows the strain history from the center of flax-PP composite for 6 J impact case. It can be seen that the strain values in the xxdirection and yy-direction closely follow one another, supporting the biaxial loading response of the composite plate. The response diverges from each other at the initiation of the cracks, as the cracks are not formed in the x and ydirections at the same moment. This is supported in the strain contours of Fig. 14. However, the response stabilizes after the crack propagates to a certain length.

A similar DIC analysis was conducted on flax-epoxy composite for impact of 6 J and results are shown in Fig. 16. The major strain contour is for image 14 after the first contact and just before the beginning of cracks. It can be seen that the strain response of the plate is symmetric at this time with a sharp peak at the centre of the plate and gradually decreasing to the edges. It should be noted again that the subset matching algorithm of the DIC method is not capable of calculating strain histories when there are discontinuities due to the cracks and fibre failures. Therefore the DIC method is applicable only for lower energy impacts before the displacement discontinuities induced by the creation of macro-cracks on the surface.

The advantage of the DIC analysis is that it is possible to obtain fullfield displacement measurement of the impacted plate which is useful when there is strain localization. In heterogeneous materials such as composites, it is important to have the local deformation history in addition to the global response. This is a promising methodology for studying the damage evolution in biocomposites. However, there are still some significant challenges regarding the use of high-speed camera and DIC for impact damage assessment. One of the main problems is that the out-of-plane displacement of the specimen during impact leads to a change in the image-to-object distance [START_REF] Pan | Two-dimensional digital image correlation for inplane displacement and strain measurement: A review[END_REF]. It is possible that the magnification of the images will change considerably due to the out-ofplane displacement, which will lead to errors in the calculation of inplane displacements. Even though out-of-plane displacement during the impact can cause a strain error in the 2D DIC measurements. A stereo DIC based on the principle of binocular stereovision can help overcome this problem. Pankow et al. [START_REF] Pankow | Three-dimensional digital image correlation technique using single high-speed camera for measuring large out-of-plane displacements at high framing rates[END_REF] have shown that large out-of-plane displacements can be captured by a 3D DIC system. It would be interesting to adapt this stereo-correlation system for impact assessment in the future. Nevertheless, Reu and Miller [START_REF] Reu | The application of high-speed digital image correlation[END_REF] pointed out that despite its drawbacks, two-dimensional information is more than adequate in a number of applications.

Crack-tracking during impact

Fig. 17 shows the evolution of cumulative crack length correlated to the force-time history obtained from the sensors. The force response exhibits an initial linear region followed by softening phase where the stiffness changes and it can be observed that the initiation of crack corresponds to the stiffness change at critical force of 1100 N. The crack length continues to increase as the crack propagates along the principal directions and reaches a maximum at the end of the loading phase. It should be noted that this is the total crack length and not the length of each radial arm. The kinematics of the crack propagation can be observed with this crack tracking method which is a significant difference to published literature such as Sy et al. [START_REF] Sy | Damage evolution in unidirectional and cross-ply flax/epoxy laminates subjected to low velocity impact loading[END_REF] who only report the post-impact total crack length.

A comparison of the crack evolution curves for flax-PP and flaxepoxy composites impacted at different energies is shown in Fig. 18. It can be seen that the crack length is initiated approximately 0.6 ms after the start of impact and increases to a maximum length which depends on the impact energy and also the materials. Indeed, there are some significant differences between the evolution of crack lengths for the flax composite with thermoplastic and thermoset matrix. The flax-PP composite reaches a maximum crack length of 29 mm, 47 mm and 51 mm for 3, 6 and 9 J impacts, respectively. On the other hand, the flax-epoxy composite reaches maximum crack length of 38, 40 and 61 mm for the same initial impact energies. It is noticeable that cracks propagation is much faster in the brittle epoxy matrix compared to the more ductile PP matrix, maximum crack length being thus reached at shorter time.

The total crack length for the different energy impacts is plotted against the absorbed energy in Fig. 19. This shows that for a given crack length, the flax-PP and flax-epoxy composites absorb similar amount of energy, indicating that despite the differences in failure mechanisms for the two composites, there is no significant effect on the quantity of energy absorbed per length of cracks by changing the matrix from epoxy to PP (with very different mechanical properties). This proves the potential of commingled thermoplastic biocomposite with shorter manufacturing cycle and lower density to have similar impact performance as the flaxepoxy composite. The ductility of the PP matrix has the potential to further improve the impact performance if the fibre-matrix interface can be improved by optimized manufacturing conditions and fibre treatments.

Conclusions

The aim of the paper was to compare the impact behaviour of biocomposites reinforced with flax fibres in thermoplastic and thermoset matrices. A compression moulding process was used to manufacture flax-PP and flax-epoxy composite plates with similar fibre volume fraction and areal density. A drop tower system was used to conduct impact tests at multiple energies, and force and displacement histories were measured during impact. Additionally, a high-speed image acquisition system was used to observe the rear surface of the target to evaluate the damage kinematics in the composite. The different modes of damage and the critical energy required for complete penetration were identified for the composite. The key findings of the paper are:

1. Comparison of the impact response shows that the performance of the flax fibre reinforced composite was similar for the thermoset and thermoplastic matrix. The flax-PP composite exhibits more ductile behavior and absorbs slightly more energy than the thermoset composite, while a more cohesive behavior was observed for flaxepoxy composite with the appearance of delamination type failure. 2. The optical instrumentation implemented in this study shows that the initiation of damage cannot be clearly identified simply from the force-displacement curve (as it is usually performed). Therefore, the assessment of damage requires the support of a local analysis, which was implemented in this work. 3. A method was proposed to use high-speed imaging for in-situ observation of the creation and propagation of damages in natural fibre-reinforced composites. A macro was implemented in ImageJ software so as to measure cracks lengths from images obtained from Phantom high-speed camera. The cumulative crack length was evaluated at each time interval and the evolution of damage observed in the plate during impact was correlated to the force--displacement response. The maximum crack length for flax-PP composite was 29 mm, 47 mm and 51 mm for 3, 6 and 9 J impacts, respectively, while the flax-epoxy composite reached maximum crack length of 38, 40 and 61 mm for the same initial impact energies. 4. In literature, the DIC technique has been used previously to measure the mechanical properties in quasi-static strain rates but their use for high-speed image processing has been limited. In our study, we applied DIC technique during impact on samples coated with white paint and with speckled patterns. The displacement contours and strain distributions were calculated using in-house program CinEMA.

From the strain contours, it was possible to identify strain localization and the initiation of the cracking. 5. One of the main outcomes of this work is that, in addition to the postmortem damage assessment of the impacted composites, it is possible to follow the damage kinematics over time with the high-speed camera.

The proposed method better depicts the successive events that lead to the breakage of composites and can be a valuable tool to study the impact behaviour of not only natural fibre-reinforced composites but other materials as well. Future study will utilize high-speed imaging with stereo-correlation to take into account out-of-plane displacement. The development of a Finite Element model of the impact behaviour of flax biocomposites and a comparison of the simulated and experimental results is also under study.

CRediT authorship contribution statement

Fig. 1 .

 1 Fig. 1. (a) Drop tower setup used for low velocity impact testing monitored with high speed camera, (b) scheme of the experimental set-up.

Fig. 2 .

 2 Fig. 2. Typical view of the high-speed camera of the target plate with (a) white powder and (b) fly-speckled pattern.

Fig. 3 .

 3 Fig. 3. Typical high-speed video frames showing initiation and propagation of macro-cracks in flax-PP during 6 J impact.

Fig. 5 a

 5 shows the typical force history for three flax-PP composite plates

Fig. 4 .

 4 Fig. 4. Procedure for crack-tracking from the high-speed camera images.

Fig. 5 .

 5 Fig. 5. Typical force-time and energy-time history for 3 J impact of flax-PP composite.

Fig. 6 .

 6 Fig. 6. Comparison of the typical force-displacement curves for different initial impact energies (a) flax-PP composite and (b) flax-Epoxy composite.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Comparison of typical Forcedisplacement curves for (a) 6 J and (b) 12 J impact.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Depth profile for flax-PP and flax-epoxy composites after 6 J and 12 J impact.

Fig. 14 .

 14 Fig. 14. (a) and (b) : Contours of Ux and Uy displacements; (c) and (d) : contours of xx and yy strains components; for flax-PP composite. (time = 0.54 ms; 21st image).

Fig. 15 .

 15 Fig. 15. Lagrangian strains calculated by DIC method for 6 J impact on flax-PP composite, with cropped views of the plate bottom surface.

Fig. 16 .

 16 Fig. 16. Strain contours in the xx and yy directions for flax-epoxy composite (6 J impact).

Fig. 17 .

 17 Fig. 17. (a) Forcedisplacement curves for 6 J impact on Flax-PP composite, along with total crack length; (b) Real time high-speed image of the plate at maximum crack length; (c) Skeletonised image of the crack.

Fig. 18 .

 18 Fig. 18. Comparison of crack lengths for 3, 6 and 9 J impact of (a) flax-PP and (b) flax-epoxy composite.

Fig. 19 .

 19 Fig. 19. Comparison of total crack lengths vs. absorbed energy for flax-PP and flax-epoxy composites.

  Ram Ramakrishnan: Conceptualization, Methodology, Investigation. Stéphane Corn: Conceptualization, Methodology, Investigation, Supervision. Nicolas Le Moigne: Conceptualization, Methodology, Investigation, Supervision. Patrick Ienny: Conceptualization, Methodology, Investigation. Pierre Slangen: Conceptualization, Methodology, Investigation.

Acknowledgements

This project was funded by the European Union's Horizon 2020 Programme (HERMES project, grant agreement n • 636520). KRR also acknowledges Mr. Jere Ahkiola for his assistance with profilometer and optical microscope imaging.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Hao Xiong

Zhen-Yu Yin
email: zhenyu.yin@polyu.edu.hk

François Nicot

Antoine Wautier

• Miot

Marie • Félix

Guillaume Veylon

Pierre Philippe

Félix Darve

A novel multi-scale large deformation approach for modelling of granular collapse

Keywords: Multiscale approach, Granular material, Micromechanics, Meso-scale, Smoothed particle hydrodynamics (SPH), Granular collapse

Collapse of granular material is usually accompanied by long run-out granular flows in natural hazards, e.g., rock/debris flow and snow avalanches. This paper presents a novel multi-scale approach for modelling granular column collapse with large deformation. This approach employs the smoothed particle hydrodynamics (SPH) method to solve large deformation boundary value problems while using a micromechanical model to derive the non-linear material response required by the SPH method. After examining the effect of initial cell size, the proposed approach is subsequently applied to simulate the flow of granular column in a rectangular channel at a low water content by varying the initial aspect ratio. The numerical results show good agreement with various experimental observations on both collapse process and final deposit morphology. Furthermore, the mesoscale behavior is also captured owing to the advantages of the micromechanical model. Finally, it was demonstrated that the novel multi-scale approach is helpful in improving the understanding of granular collapse and should be an effective computational tool for the analysis of real-scale granular flow.

Introduction

The transformation of granular material from a stationary state into a dynamic flow state can be a prototype found in many natural disasters and industrial processes. The understanding of granular flow mechanisms plays a crucial role in minimizing natural hazards (e.g., landslides, snow avalanches, and soil liquefaction) and in optimizing industrial processes (e.g., mineral processing, ceramic, food processing, pharmaceutical manufacture, and civil engineering applications). Many theories and assumptions have been proposed in the attempt to explain the apparent high mobility of granular material, including the bubbles at the bottom of granular flow [START_REF] Shreve | Leakage and fluidization in air-layer lubricated avalanches[END_REF][START_REF] Han | Functional Catastrophe Analysis of Progressive Failures for Deep Tunnel Roof Considering Variable Dilatancy Angle and Detaching Velocity[END_REF], sand fluidization [START_REF] Hungr | Entrainment of debris in rock avalanches: an analysis of a long run-out mechanism[END_REF][START_REF] Jianbo | Experimental investigation on granular flow past baffle piles and numerical simulation using a µ(I)-rheology-based approach[END_REF], grain segregation-induced friction decrease [START_REF] Phillips | Enhanced mobility of granular mixtures of fine and coarse particles[END_REF]. However, these theories and assumptions have not been fully confirmed by experiments, and there are still many mechanisms of granular flow that are still unclear.

It has been recognized that debris flows and granular flows may behave similarly: for instance, they can sustain shear stresses with very slow deformation due to lasting, frictional grain contacts, and they can flow rapidly, sustaining inelastic grain collisions [START_REF] Iverson | The physics of debris flows[END_REF][START_REF] Neto | Simulation of debris flow on an instrumented test slope using an updated Lagrangian continuum particle method[END_REF]. Thus, research has mainly focused on small scale laboratory experiments and numerical simulations of granular materials [START_REF] Tan | New simple method for calculating impact force on flexible barrier considering partial muddy debris flow passing through[END_REF][START_REF] Tan | Experimental study on impact and deposition behaviours of multiple surges of channelized debris flow on a flexible barrier[END_REF][START_REF] Chen | Influence of matric suction on nonlinear time-dependent compression behavior of a granular fill material[END_REF][START_REF] Chen | Crushing and flooding effects on one-dimensional time-dependent behaviors of a granular soil[END_REF][START_REF] Lajeunesse | Granular slumping on a horizontal surface[END_REF][START_REF] Crosta | Numerical modeling of 2-d granular step collapse on erodible and nonerodible surface[END_REF][START_REF] Lube | Collapses of two-dimensional granular columns[END_REF][START_REF] Lube | Axisymmetric collapses of granular columns[END_REF][START_REF] Rondon | Granular collapse in a fluid: role of the initial volume fraction[END_REF][START_REF] Fern | Granular column collapse of wet sand[END_REF][START_REF] Zhao | Mechanical behavior of intact completely decomposed granite soils along multi-stage loading-unloading path[END_REF][START_REF] Brezzi | Collapse of granular-cohesive soil mixtures on a horizontal plane[END_REF][START_REF] Kermani | Simulation of quasi-static axisymmetric collapse of granular columns using smoothed particle hydrodynamics and discrete element methods[END_REF]. Among them, two classic experiments are widely performed: rectangular channel flowing tests and column flowing tests. In the former, granular collapse is obtained by putting the granular material in a rectangular channel and quickly removing a vertical side boundary, whereas, in the latter, granular material is in a hollow cylinder pipe. This study only discusses the former.

There are many possible solutions to simulate the granular column collapse by using numerical methods, e.g., the finite element method (FEM) using the arbitrary Lagrangian-Eulerian (ALE) technique [START_REF] Crosta | Numerical modeling of 2-d granular step collapse on erodible and nonerodible surface[END_REF][START_REF] Wu | Numerical modelling of granular column collapse using coupled eulerian-lagrangian technique with critical state soil model[END_REF], the discrete element method (DEM) [START_REF] Utili | 3d dem investigation of granular column collapse: evaluation of debris motion and its destructive power[END_REF][START_REF] Li | Modelling the shearing behaviour of joints using an improved shear box genesis approach in particle flow code (2d) and its validation[END_REF], the particle finite element method (PFEM) [START_REF] Zhang | Numerical simulation of a flow-like landslide using the particle finite element method[END_REF], the smoothing particle finite element method (SPFEM) [START_REF] Jin | An edge-based strain smoothing particle finite element method for large deformation problems in geotechnical engineering[END_REF], the material point method (MPM) (So lowski and Sloan, 2015), and the smooth particle hydrodynamics (SPH) method [START_REF] Peng | A sph approach for large deformation analysis with hypoplastic constitutive model[END_REF][START_REF] Neto | Continuum hydrodynamics of dry granular flows employing multiplicative elastoplasticity[END_REF][START_REF] Wang | Smooth particle hydrodynamics studies of wet granular column collapses[END_REF][START_REF] Peng | Loquat: an open-source gpu-accelerated sph solver for geotechnical modeling[END_REF][START_REF] Bui | Lagrangian meshfree particles method (sph) for large deformation and failure flows of geomaterial using elastic-plastic soil constitutive model[END_REF][START_REF] Yin | Novel sph simsand-based approach for modeling of granular collapse[END_REF]Jin et al., 2019b,a). Furthermore, some multi-scale approaches have also been employed to simulate the granular column collapse problem, e.g., a hierarchical coupling scheme to integrate MPM with DEM packings for multi-scale modelling of large deformation in granular materials [START_REF] Liang | Multiscale modeling of large deformation in geomechanics[END_REF]. These methods have contributed to the investigation of granular flows, but problems remain. On one side, although ALE technology can be used to solve large deformation problems in continuum methods, it does not really simulate the flow of granular materials because it lacks a unified formulation for going from solid to fluid in continuum mechanics. On the other side, the discretization methods require lots of calculations and are not suitable for large-scale boundary value problems. However, micromechanical models try to give the best compromise. The micromechanical models have no explicit constitutive behavior but the implicit description by statistical description of grain interactions, e.g. [START_REF] Zhao | Integrating a micromechanical model for multiscale analyses[END_REF][START_REF] Zhu | Three-dimensional constitutive relations for granular materials based on the dilatant double shearing mechanism and the concept of fabric[END_REF][START_REF] Xiong | A three-dimensional micromechanically based model[END_REF][START_REF] Nicot | The H-microdirectional model: accounting for a mesoscopic scale[END_REF][START_REF] Zhu | A micromechanics-based elastoplastic damage model for granular materials at low confining pressure[END_REF][START_REF] La Ragione | A simplified model for inelastic behavior of an idealized granular material[END_REF][START_REF] Zhao | An evolution law for fabric anisotropy and its application in micromechanical modelling of granular materials[END_REF].

In this paper, a novel multi-scale framework is proposed. Under this framework, the SPH method is used to solve the large deformation boundary value problem while the 3D-H model is used to derive the nonlinear constitutive material response of granular material required by SPH for each point. The 3D-H model is also extended to the condition of low water content by considering the capillary forces in water bridges between connected grains. Firstly, five simulations with different initial cell size are compared. Then, the proposed framework is subsequently applied to simulate the rectangular channel flowing tests at a low water content by varying the initial aspect ratio of the granular column. By comparing with various experimental observations, it verified that the proposed framework is able to predict the collapse process and the final deposit morphology. Furthermore, the mesoscopic analysis is performed owing to the advantages of the 3D-H model.

2 Multi-scale framework for granular flows

Fundamentals of the SPH method

The SPH method was adopted in this study to model large deformation of granular collapse. The SPH method was initially proposed by [START_REF] Gingold | Smoothed particle hydrodynamics: theory and application to non-spherical stars[END_REF][START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF] in astrophysics. It was then introduced to other fields owing its further developments. The derivation of the SPH method can be divided into two stages. The first stage implements a kernel approximation, where all the physical variables Smoothing function, W Integration domain i j Fig. 1: Particle approximations using particles within the integration domain of the smoothing function W for particle i.

are approximated using an integral transformation of discrete quantities involving the convolution of the physical variable with a chosen kernel function.

The kernel approximation of a function f (x) at position x = (x, y, z) in physical domain Ω is expressed as

f (x) = Ω f (x )W (x -x , h)dx (1)
where W is the kernel function and h is the smoothing length.

The kernel function W has to satisfy certain requirement such as normalization condition, delta function property and compact support condition [START_REF] Monaghan | An introduction to sph[END_REF][START_REF] Peng | Loquat: an open-source gpu-accelerated sph solver for geotechnical modeling[END_REF]. The cubic spine kernel is employed in this study. The kernel approximation of the gradient of the function ∇f (x) can also be obtained in the same manner to give

∇f (x) = - Ω f (x )∇ x W (x -x , h)dx (2)
where ∇ x indicates that the derivatives are evaluated at x . If we assume that the continuous field f (x) comes from a piecewise constant field f , then the kernel representation for a function and its gradient as exhibited above can be approximated by a Riemann summation over the particles (i.e. the subdomains Ω j where f takes constant value f (x j ), x j being the center of Ω i ) used in the simulation, thus

f (x) = N j=1 f (x j )W (x i -x j , h) m j ρ j (3) ∇f (x) = - N j=1 f (x j )∇ j W (x i -x j , h) m j ρ j (4)
where ρ j and m j are the density and mass of particle j (m j /ρ j = |Ω j |), respectively, and N is the number of particles. Note that ∇ i W (x i -x j , h) = -∇ j W (x i -x j , h) due to the symmetry of the kernel function.

The second stage of the SPH method involves solving the equations governing granular impact; namely, the mass and momentum balance equations. These equations can be discretized with the kernel approximation as follows:

             dρ i dt = N j=1 m j v α i -v α j • ∂W (x i -x j , h) ∂x α i dv α i dt = N j=1 m j σ αβ i + σ αβ j ρ i ρ j + Π ij δ αβ ∂W (x i -x j , h) ∂x β i (5)
Fig. 2: General homogenization scheme of 3D-H model [START_REF] Cambou | Homogenization for granular materials[END_REF].

where ρ i , v i , p i , σ i , m i are the density, velocity, pressure, stress and mass of particle i; d(•)/dt denotes the time derivative of a physical quantity. The term Π ij in the momentum equation is the Monaghantype artificial viscosity. This term is introduced to smooth potential unphysical oscillations, to prevent unphysical particle-particle penetrations, and to stabilize the numerical solutions. The Monaghan-type artificial viscosity is incorporated into the momentum equation as follows [START_REF] Monaghan | An introduction to sph[END_REF]:

Π ij =      -ac ij φ ij + βφ 2 ij ρ ij v ij • x ij < 0 0 v ij • x ij ≥ 0 (6)
where

φ ij = (h ij v ij • x ij )/(|x ij | 2 + ϕ 2 ), c ij = 1 2 (c i + c j ), ρ ij = 1 2 (ρ i + ρ j ), h ij = 1 2 (h i + h j ), v ij = v i -v j , and x ij = x i -x j .
Here c i is the speed of sound associated with the particle i; α and β are constant coefficients that have values of about 1.0; and, φ = 0.1h ij is applied in order to avoid singularity when two particles come too close.

Multi-scale constitutive model: 3D-H model

In Equation 5, the stress tensor has to be specified. This stress is linked to the strain field through a constitutive equation. In the present work, we chose to account for this constitutive equation in an implicit way with use of a micromechanical model. The 3D-H model [START_REF] Xiong | A three-dimensional micromechanically based model[END_REF] was previously proposed by the authors and has been successfully implemented within a finite element code to solve geotechnical engineering problem (Xiong et al., 2019a,b). Based on the kinematic hypothesis, the 3D-H model enables the derivation of the macroscopic stress tensor from the macroscopic strain tensor according to the following steps (Figure 2): (1) Kinematic localisation: The meso-structure (shown in Figure 4) is a connection between the macroand meso-scale. The dimension of the meso-structure can be characterised by the vector: L = [l 1 , l 2 , l 3 ] T , wherein l 1 , l 2 , l 3 represent the lengths along directions n, t, w, respectively (see Figure 5a and Figure 6a). Thus the kinematic localisation assumption gives:

     δl 1 l 1 = t n • δε • n = P • δε • P -1 11 δl 2 l 2 = t t • δε • t = P • δε • P -1 22 δl 3 l 3 = t w • δε • w = P • δε • P -1 33 (7)
where δε is the incremental macroscopic strain tensor, and P is the rotation matrix from global frame (x 1 , x 2 , x 3 ) to local frame (n, t, w) (see Figure 3).

The kinematic localisation defined by Equation 7 is a homogenisation process. It is analogous to the usual Voigt approximation in the field of continuous media. The reader should note that the localisation process goes from the macro-to the meso-scale, not to the micro-scale, as assumed in the Voigt approximation.

F a 1 F a 1 F 2 F 2 F 2 F 2 G 2 G 2 G 2 G 2 1 4 9 8 10 7 F b 1 F b 1 F 3 F 3 F 3 F 3 G 3 G 3 G 3 G 3
Fig. 4: The 3D mesostructure and its decomposition procedure in the 3D-H model.

Moreover, it has been widely used in granular materials, such as in [START_REF] Cambou | Homogenization for granular materials[END_REF][START_REF] Nicot | The H-microdirectional model: accounting for a mesoscopic scale[END_REF][START_REF] Nicot | A multi-scale approach to granular materials[END_REF].

In the 3D-H model, the hexagon deformation described by the vector δL = [δl 1 , δl 2 , δl 3 ] T is derived from the macroscopic strain tensor. Then, the term δL is used as the known variable to compute the relative displacement at each contact and to compute contact forces.

(2) Meso-structure behaviour: The meso-structure (Figure 4) can be decomposed into two independent hexagon patterns: Hexagon A (Figure 5) and Hexagon B (Figure 6), both being similar. The geometrical configuration and external forces applied to the meso-structure are symmetrical; thus for each hexagon, only two grains need to be analysed. For Hexagon A, as shown in Figure 5, only grains 1 and 2 are analysed.

The contact between grains 1 and 2 is denoted by contact 1, whereas the contact between grains 2 and 3 is denoted by contact 2. Then, the kinematic relations read (for Hexagon A):

   δu 1 n = δd 1 δu 1 t = d 1 δα 1 δu 2 n = δd 2 (8)
where u i n and u i t represent the normal and tangential relative displacements at contact i (For Hexagon A, i = 1 or 2, for Hexagon B, i = 3 or 4), respectively. As depicted in Figure 5a, the geometrical description for Hexagon A gives:

l 1 = d 2 + 2d 1 cos α 1 l 2 = 2d 1 sin α 1 (9)
The force balance of grain 1 along direction n and of grain 2 along directions w and n, together with the moment balance of grain 2, reads:

F a 1 =2(N 1 cos α 1 + T 1 sin α 1 ) (10a) F 2 =N 1 sin α 1 -T 1 cos α 1 (10b) N 2 =N 1 cos α 1 + T 1 sin α 1 + G 2 (10c) G 2 =T 2 (10d)
where N i and T i represent the normal and tangential contact forces of contact i, respectively. The elasticperfect plastic inter-particle contact law reads in an incremental form, for a given contact i:

   δN i = k n δu i n δT i = min T i + k t δu i t , tan ϕ g (N i + δN i ) × T i + k t δu i t T i + k t δu i t -T i (11)
After simplifying (see more details in Appendix A), the contact law (for Hexagon A) can be rewritten as follows:

   δN 1 = -k n δd 1 δN 2 = -k n δd 2 δT 1 = B 1 δα 1 -A 1 δd 1 + C 1 (12)
Term C 1 differs from zero only during a transition from an elastic to a plastic regime. For very small strain increments, as considered throughout this paper, term C 1 can therefore be neglected.

To obtain the incremental evolution of the external forces, δd 1 , δd 2 and δα 1 need to be expressed as a function of the mesoscopic strain. Three equations are therefore required. Compatibilities (Equations 9) provide two relations. The third is the balance equation of grain 2 along direction n (Equation 10c). By taking the inter-particle contact law (Equations 12) into account, we reach the following algebraic system expressing the incremental changes in δd 1 , δd 2 and δα 1 with respect to the incremental changes in δl 1 and δl 2 :

  2 cos α 1 1 -2d 1 sin α 1 2 sin α 1 0 2d 1 cos α 1 cos α 1 + A 1 k n (sin α 1 + 1) -1 F 2 -B 1 (sin α 1 +1) k n     δd 1 δd 2 δα 1   =   δl 1 δl 2 0   ( 13
)
where A 1 , B 1 are given in A.

Differentiating Equation 10a, Equation 10b and combining with Equations 12 gives:

δF a 1 = -k n cos α 1 δu 1 n + k t sin α 1 δu 1 t -F 2 δα 1 δF 2 = -k t cos α 1 δu 1 t -k n sin α 1 δu 1 n + F a 1 δα 1 (14)
Thus, combining Equations 8, Equations 13 and Equation 14, the incremental constitutive relation for Hexagon A can be expressed as follows: where:

1 |D| a K a 11 K a 12 K a 21 K a 22 δl 1 δl 2 = δF a 1 δF 2 (15) d 1 d 2 d 1 d 2 d 1 d 1 α 1 l 2 l 1 (a) Dimension of Hexagon A 1 4 5 3 6 2 F a 1 F a 1 F 2 F 2 F 2 F 2 G 2 G 2 G 2 G 2 (b) External forces 2 N 2 N 1 T 1 G 2 F 2 1 F a 1 N 1 T 1 N 1 T 1 (c)
                                           K a 11 = 2 F 2 sin α 1 -k n d 1 cos 2 α 1 -k t d 1 sin 2 α 1 K a 12 = (k t d 1 sin α 1 -F 2 ) A 1 k n sin α 1 + A 1 k n + 3 cos α 1 -cos α 1 (B 1 sin α 1 + B 1 -F 2 + 2k n d 1 sin α 1 ) K a 21 = 2 (k t -k n ) d 1 sin α 1 cos α 1 -2F a 1 sin α 1 K a 22 = (F a 1 -k t d 1 cos α 1 ) A 1 k n sin α 1 + A 1 k n + 3 cos α 1 -sin α 1 (B 1 sin α 1 + B 1 -F 2 + 2k n d 1 sin α 1 ) |D| a = 2 k n (B 1 sin α 1 + A 1 d 1 cos α 1 )(sin α 1 + 1) - 2 k n (F 2 sin α 1 + k n d 1 cos 2 α 1 + 2k n d 1 ) (16)
Similarly, the incremental constitutive relation for Hexagon B can also be obtained. Consequently, superimposing Hexagon A and Hexagon B, the total incremental force along direction n is δF 1 = δF a 1 +δF b 1 . The incremental constitutive relation of the 3D meso-structure is finally obtained.

(3) Stress averaging: Averaging the mesoscopic stress σ taking place within all the meso-structures in the specimen of volume V can be performed as follows:

σ = 1 V ω(θ, ϕ, ψ)P -1 V σ(n, t, w)P sin ϕdθdϕdψ ( 17
)
where σ is the macro-stress tensor operating on the specimen scale, V is the volume of meso-structure along the local frame (n, t, w), V is total volume of the specimen, which is computed from V as follows:

V = ω(θ, ϕ, ψ) V sin ϕdθdϕdψ (18)
ω(θ, ϕ, ψ) is a probability density function, which fulfills the following relation: For an isotropic specimen, ω 0 = 1/4π 2 . V σ with respect to the local frame can be computed from the local variables (Figure 5 and Figure 6) using the Love-Weber formula [START_REF] De Saxcé | About the numerical simulation of the dynamics of granular media and the definition of the mean stress tensor[END_REF][START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF][START_REF] Christoffersen | A micromechanical description of granular material behavior[END_REF][START_REF] Mehrabadi | On statistical description of stress and fabric in granular materials[END_REF]:

2π 0 π 0 2π ω(θ, ϕ, ψ) sin ϕdθdϕdψ = 1 (19) d 3 d 4 d 3 d 4 d 3 d 3 α 2 l 3 l 1 (a) Dimension of Hexagon B 1 4 9 8 10 7 F b 1 F b 1 F 3 F 3 F 3 F 3 G 3 G 3 G 3 G 3 (b) External forces 7 N 4 N 3 T 3 G 3 F 3 1 F b 1 N 3 T 3 N 3 T 3 (c)
       V σ 11 = 4N 1 d 1 cos 2 α 1 + 4T 1 d 1 cos α 1 sin α 1 + 2N 2 d 2 +4N 3 d 3 cos 2 α 2 + 4T 3 d 3 cos α 2 sin α 2 + 2N 4 d 4 V σ 22 = 4N 1 d 1 sin 2 α 1 -4T 1 d 1 cos α 1 sin α 1 V σ 33 = 4N 3 d 3 sin 2 α 2 -4T 3 d 3 cos α 2 sin α 2 (20)
Notably, the local void ratio is related to the opening angle, which is not related to local anisotropy. The opening angle α 1(2) is a geometrical parameter (Figure 5 and Figure 6). The opening angle, together with the components l 1 , l 2 , l 3 , determine the initial shape of the hexagons as well as the local void ratio of the meso-structure. For a virgin specimen, the initial opening angle is denoted as α 0 . Then, α 1 = α 2 = α 0 . The initial void ratio e 0 (n, t, w) of each meso-structure belonging to the local frame (n, t, w) can be estimated using the initial opening angle α 0 as follows:

e 0 = - 4 π cos 3 α 0 - 6 π cos 2 α 0 + 4 π cos α 0 + 6 π -1 (21)
In summary, the total number of input parameters of the 3D-H model is four (in case of initial microstructure isotropy). Three material parameters plus one geometry parameter. All parameters can be obtained by conventional drained triaxial test. The validations of this constitutive model have been conducted in our previous works [START_REF] Xiong | A three-dimensional micromechanically based model[END_REF](Xiong et al., , 2019a[START_REF] Xiong | Programming a micro-mechanical model of granular materials in julia[END_REF], in which an excellent agreement between the 3D-H model and experimental data on dry sand has been achieved; therefore, these will not be repeated here.

Extension of 3D-H model at a low water content

At low water contents, water in a granular material forms individual capillary bridges between pairs of grains. These capillary bridges modify the constitutive behaviour of granular materials by introducing some cohesion. In this section, we show how the 3D-H model used so far can be extended to take into account capillary forces directly at contact scale for low water contents. The macroscopic effect of the capillary cohesion thus introduced is then investigated on the same numerical set up of Figure 13.

From a mechanical point of view, a capillary bridge between two grains can be taken into account in the form of an attractive force between the two connected grains. This capillary force can be computed with use of a surface energy minimization method. Thus, the total surface energy E s of two grains connected by a volume of water is given by

E s = γ lg A lg + 2 i=1 γ sl i A sl i + 2 i=1 γ sg i A sg i (22)
where γ lg = 0.073 N.m -1 is the surface tension of air/water interface at 20 • C given in [START_REF] Molenkamp | Interactions between two rough spheres, water bridge and water vapour[END_REF], A lg is the area of the liquid/gas interface, γ sl i and A sl i are respectively the surface tension and the area of the interface between grain i ∈ {1, 2} and water, and γ sg i and A sg i are respectively the surface tension and the area of the interface between grain i and air.

The surface tensions of the interfaces are related to the wetting angle θ by the Young-Dupré equation Yuan and Lee (2013):

γ sl i -γ sg i = -cos θ i γ lg (23)
Thus the surface energy can be expressed by:

E s = γ lg A lg -γ lg 2 i=1 cos θ i A sl i + C ( 24
)
where C = 2 i=1 γ sg i A s i is a constant, with A s i = A sg i + A sl i the total area of grain i. The geometries of the liquid/gas and the solid/liquid interfaces are determined numerically with the surface energy minimization method implemented in Surface Evolver [START_REF] Brakke | The surface evolver[END_REF]. The numerical procedure is briefly detailed in the following. First, a rough mesh of the interfaces is created. With use of a gradient descent method, the nodes of the mesh are moved in order to minimize the total energy while keeping constant the volume of water. Then, the mesh is refined and the gradient descent procedure repeated until a precise geometry of the capillary bridge is obtained.

Eventually, the capillary force acting between the two grains is obtained from the virtual work principle. As presented in Figure 7, an incremental displacement δd is imposed to one of the grain in the direction of the bridge, and the surface energy of this new configuration E s (δd) is calculated and compared with the original one. In a closed, static system without gravity, there are neither kinetic nor potential and internal energy. The incremental work δW of the capillary force F c is thus related to the surface energy variation as:

δW = F c .δd = δE s ( 25
)
This means that the capillary force F c can be expressed as a function of the variation of surface energy E s (δd) -E s (0) and the incremental displacement δd:

F c = E s (δd) -E s (0) δd (26)
By repeating the whole procedure for varying liquid volume, the dimensionless capillary force

F * c = F c γ sl r
and the half-filling angle β (as defined in Figure 7) can be computed as functions of the dimensionless liquid volume V * = V /R 3 where R is the grain radius. Note that, in the present work, the two grains are of equal radius and in contact. The reason for this restriction lies in the fact that the purpose of these simulations is to provide an empirical expression of the capillary force acting in between the contacting grains of the ten grain meso-structure of the H-model (Figure 10). Figure 8 shows that the capillary force in a liquid bridge can be approximated fairly well with a 4th degree polynomial while Figure 9 shows that the half filling angle can be approximated really accurately with a power law:

F * c = 0.701V * 4 -3.167V * 3 + 5.150V * 2 -4.290V * + 5.380 β = 68.0V * 0.293 (27)
For low water contents, water in a meso-structure forms independent capillary bridges between grains two by two, as presented in Figure 10.

We suppose that the water content w c provided at macroscale is the same in all the meso-structures (whatever their local orientations). This hypothesis can be seen as a localisation of the water content, as for the strain localisation. Under this hypothesis, the dimensionless water volume V * w in a given meso-domain is: Fig. 8: Capillary force in a liquid bridge between two grains in contact as a function of the water volume calculated with the surface energy minimization method (solid line) and approached by a 4th degree polynomial (dashed line).

V * w = w c V * s γ s γ w (28)
Fig. 9: Half-filling angle in a liquid bridge between two grains in contact as a function of the water volume calculated with the surface energy minimization method (solid line) and approached by a power law (dashed line).

Fig. 10: Individual capillary bridges between grains in the H-model meso-structure. Vertical bridges are coloured in blue while inclined bridges are coloured in green with V s being the volume of the 10 grains contained into the meso-domain V * s = V s r 3 = 40 3 π, and γ s γ w = 2.4 the ratio of the solid and liquid specific weight.

The total volume of water V * w lying inside the dodecahedron come from the contribution of the four vertical bridges and eight inclined bridges visible in Figure 10. To keep the symmetry of the meso-structure, we assume that all vertical bridges have the same volume V * v and all inclined bridges have the same volume V * i . As a result, the total volume of water splits as:

V * w = 8V * i + 4V * v (29)
In the above equation, we have one degree of freedom to split the volume V * w between V * i and V * v . We suppose that the distribution of water volume between the different capillary bridges is the distribution which minimizes the total surface energy in the meso-structure. However, this distribution has to insure that the bridges are not in contact, which would involve the merging of capillary bridges into more complex coalesced bridges (this particular point will be deeply investigated in a forthcoming paper). To carry on the computation, it is hypothesized that grain interpenetration is neglected so that the two hexagonal patterns of a cell are identical. Thus, it is assumed that:

   α = α 1 + α 2 2 d 1 = d 3 = d 2 = d 4 = 2r (30)
The dimensionless total surface energy in a meso-structure E meso * s can therefore be written as:

E meso * s = 8E bridge * s (V * i ) + 4E bridge * s (V * v ) = 8E bridge * s ( V * w -4V * v 8 ) + 4E bridge * s (V * v ) (31)
The derivative of the total surface energy with respect to the volume of one vertical bridge is:

dE meso * s dV * v = -4 dE bridge * s dV * v ( V * w -V * v 4 ) + 4 dE bridge * s dV * v (V * v ) (32)
Figure 11 shows the surface energy in a bridge. The minimum of E meso * s , given by

dE meso * s dV * v = 0 corre- sponds to the case where V * i = V * v = V * w 12 .
To insure that the bridges do not touch each other, the filling angles of the vertical bridges β v and of the inclined bridges β i , as defined in Figure 12, should fulfil the following conditions:

     cos β i > cos β i,coal = 1 + cos 2 α 2 β v > β v,coal = π -α -β i (33)
where β i,coal is the limit half-filling angle of the inclined bridges corresponding to the contact of the inclined bridge and β v,coal is the limit half-filling angle of the vertical bridges corresponding to the contact between inclined and vertical bridges. The conditions of contact defined by Equation 33 are checked for V * i = V * w /12. In case of contact, V i is decreased to the maximal value which insure the independence of the capillary bridges, in order to obtain the geometrically compatible configuration with the minimal surface energy.

At this point, volume in each bridge is known, and the capillary force at each contact of the cell can be deduced from the Equation 27. These capillary forces, noted N cap 17, thus, using Love-Weber formula:

V σ cap 11 = 4(N 1 + N cap 1 )d 1 cos 2 α 1 + 4T 1 d 1 cos α 1 sin α 1 + 2(N 2 + N cap 2 )d 2 +4(N 3 + N cap 3 )d 3 cos 2 α 2 + 4T 3 d 3 cos α 2 sin α 2 + 2(N 4 + N cap 4 )d 4 V σ cap 22 = 4(N 1 + N cap 1 )d 1 sin 2 α 1 -4T 1 d 1 cos α 1 sin α 1 V σ cap 33 = 4(N 3 + N cap 3 )d 3 sin 2 α 2 -4T 3 d 3 cos α 2 sin α 2 (34) σ = 1 V ω(θ, ϕ, ψ)P -1 V σcap (n, t, w)P sin ϕdϕdθdψ (35)
When considering a large deformation problem such as sand column collapse, an invariant stress tensor with respect to rigid-body rotation must be applied. Accordingly, the Green-Naghdi rate which gives an objective measure of the stress rate was adopted in this study [START_REF] Green | A general theory of an elastic-plastic continuum[END_REF][START_REF] Neto | Continuum hydrodynamics of dry granular flows employing multiplicative elastoplasticity[END_REF][START_REF] Neto | Simulation of debris flow on an instrumented test slope using an updated Lagrangian continuum particle method[END_REF][START_REF] Liang | Multiscale modeling of large deformation in geomechanics[END_REF][START_REF] Jin | An edge-based strain smoothing particle finite element method for large deformation problems in geotechnical engineering[END_REF][START_REF] Systèmes | Abaqus 6.14 analysis user manual[END_REF]: where Ω = Ṙ • R T , R is the rigid body rotation in the polar decomposition of the deformation gradient F . The differences between Green-Naghdi rate and Jaumann rate are significant only if finite rotation of a material point is accompanied by finite shear.

σ G = σ -Ω • σ + σ • Ω T (36)

Implementation of 3D-H model in SPH

The implementation of the 3D-H model in the SPH codes constitutes a complete multi-scale procedure.

The SPH method is an effective method to solve the large deformation boundary value problem. SPH searches whether particles i and j interact with each other. Once determined that the two particles interact, the statistical representative volume element (RVE) of the 3D-H model embedded in the SPH particles begins to work. As illustrated before, the 3D-H model localizes the incremental macro strain tensor by using Equation 7. Then, the meso-structures are decomposed into two perpendicular hexagons, in which the constitutive behaviour is able to be solved (Equation 13). Hereafter, the incremental macro stress tensor of the RVE is obtained by stress averaging (Equation 17). Finally, the position and velocity of SPH particles are updated. The entire calculation process of the proposed multi-scale approach is shown in Figure 13.

3 Multi-scale modelling of granular collapse

Model setup

The numerical example chosen for the demonstration of the proposed multi-scale approach is a classical dynamic problem: the two-dimensional granular column collapsed to frictional ground under gravity. The geometry and the boundary conditions of the granular column are shown in Figure 14. L i and H i denote the initial length and initial height, respectively. The initial aspect ratio is defined as a = H i /L i . After granular columns collapse to finial deposit profile, L f and H f denote the run-out distance and deposit height, respectively. Except for the ground, other boundary walls are frictionless walls. The contact between SPH particles and the frictional wall is described by the classical Coulomb friction law with a friction coefficient of 0.3. Dimension and discretization parameters for all numerical simulations are summarized in Table 2. The parameters used in the extended 3D-H model is reported in Table 1. The process of numerical simulation is divided into two stages: (1) SPH particles are generated and are reached the self-weight balance;

(2) move up the right wall with a speed of 1 m/s. It is worth noting that the initial stress of SPH particles is set to balance with the self-weight. Due to the particularity of the micro-mechanical model, all RVEs need to be confined up to the initial stress to obtain the state variables at the microscopic scale, just like DEM. During the calculation, six representative SPH particles (A-F marked by red cross in Figure 14a) are selected for meso-scale analysis.

Effect of initial cell size

Negligence of initial cell size can sometimes be an embarrassment in SPH, especially in large deformation problems and dynamic problems. Thus the effect of initial cell size of the proposed multi-scale approach is first examined. Five tests with different numbers of SPH particles are performed. Figure 15 shows a comparison of the final deposition profile between different initial cell size tests. The results show that there are slight differences that are negligible.

Performance of the proposed approach

The problem of granular column collapse has been extensively studied in the literature. For granular flows, the assessment of the final run-out distance is crucially significant, as it determines the extent of the regions affected by the avalanches or landslides. In this section, our simulations have been compared with previous numerical results from [START_REF] Staron | Study of the collapse of granular columns using two-dimensional discrete-grain simulation[END_REF][START_REF] Crosta | Numerical modeling of 2-d granular step collapse on erodible and nonerodible surface[END_REF][START_REF] Utili | 3d dem investigation of granular column collapse: evaluation of debris motion and its destructive power[END_REF] and experimental observations from [START_REF] Lajeunesse | Granular slumping on a horizontal surface[END_REF][START_REF] Lube | Collapses of two-dimensional granular columns[END_REF] in the literature. Following the dimensional analysis by [START_REF] Utili | 3d dem investigation of granular column collapse: evaluation of debris motion and its destructive power[END_REF], the final normalized deposit height is defined as [H] = H f /L i while the final normalized run-out distance is defined as [L] = (L f -L i )/L i . Figure 16 shows the final normalized deposit height [H] and the final normalized run-out distance [L] versus different initial aspect ratio. As

,QLWLDODVSHFWUDWLRa shown in Figure 16a, it can be observed that the final normalized run-out distance predicted by this study shows good agreement with experimental observations and numerical results from a qualitative viewpoint, including [START_REF] Lube | Collapses of two-dimensional granular columns[END_REF] from experimental tests performed under plane strain conditions, the 3D DEM numerical results from [START_REF] Utili | 3d dem investigation of granular column collapse: evaluation of debris motion and its destructive power[END_REF], and the 2D FEM numerical results from [START_REF] Crosta | Numerical modeling of 2-d granular step collapse on erodible and nonerodible surface[END_REF]. The difference between experimental observations from [START_REF] Lajeunesse | Granular slumping on a horizontal surface[END_REF] and the results obtained by this study probably stems from the usage of glass beans. Since particle angularity tends to reduce run-out, the unrealistically large run-out distance is observed by [START_REF] Lajeunesse | Granular slumping on a horizontal surface[END_REF]. Similarly, unrealistically large run-out distance is also observed by [START_REF] Staron | Study of the collapse of granular columns using two-dimensional discrete-grain simulation[END_REF], due to the fact that 2D DEM simulations are employed. Meanwhile, the final normalized deposit height shown in Figure 16b has similar observations as the final normalized run-out distance. The predicted results are consistent with the results obtained by [START_REF] Lube | Collapses of two-dimensional granular columns[END_REF][START_REF] Utili | 3d dem investigation of granular column collapse: evaluation of debris motion and its destructive power[END_REF][START_REF] Crosta | Numerical modeling of 2-d granular step collapse on erodible and nonerodible surface[END_REF], but lower than [START_REF] Lajeunesse | Granular slumping on a horizontal surface[END_REF][START_REF] Staron | Study of the collapse of granular columns using two-dimensional discrete-grain simulation[END_REF]. It can be concluded that the proposed approach is able to predict the final deposit morphology of collapse of granular column. The prediction of the final deposit profile by the proposed approach was discussed previously, below we will focus on the collapse process. In Figure 17, the evolutions of normalized granular spreading length versus normalized time is plotted for various initial aspect ratio. Following [START_REF] Lube | Collapses of two-dimensional granular columns[END_REF][START_REF] Crosta | Numerical modeling of 2-d granular step collapse on erodible and nonerodible surface[END_REF], time has been normalized by [t] = t/ H i /g, which can be thought of as the time taken by a single particle in free fall to travel from the centre of the column to the base. As shown in the figure, the development of the curve of all the results can be roughly divided into four stages: (1) 0

[L] = (L f -L i )/L i /DMHXQHVVHHWDO /XEHHWDO 6WDURQ +LQFK &URVWDHWDO 8WLOLHWDO 7KLVVWXG\ (a) Final normalized run-out distance ,QLWLDODVSHFWUDWLRa [H] = H f /L i /DMHXQHVVHHWDO /XEHHWDO 6WDURQ +LQFK &URVWDHWDO 8WLOLHWDO 7KLVVWXG\ (b) Final normalized deposit height
t H i /g [L] = (L f -L i )L
[t] < 1, the onset of collapse with an initial transient acceleration; (2) 1 [t] < 2, the development of collapse with a constant velocity flow; (3) 2 [t] < 3 stabilizing of granular flow with a gradual deceleration; (4) [t] ≥ 3 final deposit of granular flow. By comparing the numerical results obtained by this study with the results from literatures, it can be found that the predicted results by the proposed approach located in a reasonable range. This verifies the correctness of the proposed multi-scale approach from another point of view.

Certainly, the verification of Figure 17 is only from the perspective of the spreading length of granular column collapse. For the sake of concision, the test a=1.0 is selected for a detailed comparison of the collapse process. Figure 18 shows the comparison of failure progress of granular collapse between numerical results obtained by proposed approach and experimental observations from [START_REF] Nguyen | A new sph-based approach to simulation of granular flows using viscous damping and stress regularisation[END_REF] for initial aspect ratio a=1.0 for both dry and unsaturated materials. The results of dry material are shown on the left column while the unsaturated case with water content w c = 10% are shown on the right column. It is worth noting that the velocity fields obtained by numerical calculations are all extracted in strict accordance with the time points recorded in the experiment. In general, the simulated predictions are basically consistent with the observed morphological evolution. Since the experimental results are based on dry conditions, the numerical results of dry case is closer to the experimental results. It also implies that the proposed approach can not only predict the final deposit state, but also capture the whole process of granular column collapse problem. The comparison between the dry and unsaturated conditions shows that the capillary forces delay the collapse and reduces the run-out distance of the granular column. This observation is qualitatively in agreement with experimental results obtained by [START_REF] Gabrieli | Discrete particle simulations and experiments on the collapse of wet granular columns[END_REF] and reproduced numerically by [START_REF] Wang | Smooth particle hydrodynamics studies of wet granular column collapses[END_REF]. Quantitative comparison are out of the scope of the present paper but will be performed in forthcoming works.

Meso-scale analysis

Taking advantage of the 3D-H model, a meso-scale analysis can be performed. The key microscopic behaviours hiding behind the macroscopic scale are helpful to understand and interpret the micro-mechanisms governing the overall response. Six representative SPH particles (A-F, marked as red cross in Figure 14a) are selected. Figure 19

σ I n (ϕ) = ω σ n dθdψ σ I t (ϕ) = ω σ t dθdψ (37)
where σ n and σ t are normal and tangential stress on the meso-scale, thus σ I n (ϕ) and σ I t (ϕ) can be considered as a function of the last Euler angle ϕ, it should be noted that ϕ = 0 • coincides with the positive z-axis direction in the global coordinate system.

The first observation is that all mesoscopic variables are isotropically distributed when t=0 s, as mentioned earlier. The initial magnitude of the mesoscopic variables at each selected point is different, depending on the initial stress at which the point is located. As the collapse progressed, the distribution of σ I n (ϕ) and σ I t (ϕ) gradually became anisotropic, indicating that force fabric gradually evolved into anisotropy. This is a natural evolution of the 3D-H model, with no parameter control. In response to the large deformation dynamics problem, such as granular column collapse, the proposed approach is dedicated to simulating macroscopic stress-strain response based on the evolution of microscopic and mesoscopic information. Naturally, In granular materials, the change in contact fabric is due to the disconnection of existing contacts and the re-contact of new contacts. In the current version of the 3D-H model, when a meso-structure fails, it is not lost, but stored in the system with a 'no-contribution' label. If the global deformation develops and makes this opening contact re-contact, the failed meso-structure is reactivated. However, we do mention that in the current version of the model, no new local meso-structure can appear during a loading program. The meso-scale analysis helps to compare the observations by micro-CT or results by DEM simulations and further improve the model, especially in an actual boundary value problem.

σI n (ϕ)3D σI t (ϕ)3D α 2 • (a) Point A σI n (ϕ)3D σI t (ϕ)3D α 2 • (b) Point B σI n (ϕ)3D σI t (ϕ)3D α 2 • (c) Point C
σI n (ϕ)3D σI t (ϕ)3D α 2 • (a) Point D σI n (ϕ)3D σI t (ϕ)3D α 2 • (b) Point E σI n (ϕ)3D σI t (ϕ)3D α 2 • (c) Point F

Conclusions

This paper presents a new multi-scale coupling approach for simulation of granular column collapse without using the traditional phenomenological constitutive models. The SPH method is employed to solve large deformation boundary value problems while the 3D-H model is used to derive the necessary constitutive response for the representative volume element embedded in each material point of the SPH. A preliminary attempt of the 3D-H model for a low water content condition has been made by taking into account the capillary force between connected grains. The main results are summarized as below:

1. The impact of capillary forces for small water contents is noticeable and impacts both the dynamic of the collapse and the run-out distance of the granular column. Nevertheless, the impact for small water contents is shown to be quantitatively limited. The maximum water content that can be presently reached is limited by the pendular regime approximation. Capillary forces calculation in 3D-H model for higher water contents, for which pendular bridges merge, is currently being investigated. 2. With regard to the large deformation problem, the multi-scale approach can not only capture the final deposition profiles of the granular material after large deformation, but also capture the microstructure failure process from small deformation to large deformation. 3. The final deposit profiles of granular column predicted by the proposed approach quantitatively match the DEM analysis by [START_REF] Utili | 3d dem investigation of granular column collapse: evaluation of debris motion and its destructive power[END_REF] and experimental results by [START_REF] Lube | Collapses of two-dimensional granular columns[END_REF]. The evolution of final normalized run-out distance [L] and deposit height [H] versus initial aspect ratio α are consistent with the results by [START_REF] Lajeunesse | Granular slumping on a horizontal surface[END_REF][START_REF] Staron | Study of the collapse of granular columns using two-dimensional discrete-grain simulation[END_REF][START_REF] Crosta | Numerical modeling of 2-d granular step collapse on erodible and nonerodible surface[END_REF]. 4. Different four stages have been observed in the evolution of normalized run-out distance [L] versus normalized time [t]. The whole process of collapse is consistent with the experimental observations by [START_REF] Nguyen | A new sph-based approach to simulation of granular flows using viscous damping and stress regularisation[END_REF]. 5. Owing to the advantages of the proposed multi-scale approach, the meso-scale analysis is performed within six material points in SPH. It shows how the mesoscopic variables evolve during such a dynamic large deformation problem. It also helps to the further improvements of the micromechanical model by comparing with the observations by micro-CT or numerical results by DEM.

Finally, the coupling multi-scale approach provides an effective computational tool for the analysis of granular flows.

Acknowledgment

The financial support provided by the GRF project (Grant No. 15209119) and the RIF project (Grant No. PolyU R5037-18F) from the Research Grants Council (RGC) of Hong Kong are gratefully acknowledged. We also gratefully acknowledge the CNRS International Research Network GeoMech for having offered the opportunity to make this project possible through a long term collaboration of all the authors (http://gdrmege.univ-lr.fr/).

A Contact law

This elastic-perfect plastic model includes a Mohr-Coulomb criterion and can be expressed under the following incremental formalism:

δN i = knδu i n δT i = min T i + ktδu i t , tan ϕg N i + knδu i n × T i +k t δu i t T i +k t δu i t -T i (38)
where: i = 1, 2, 3, 4 denotes the identifier of contact number. According to Equations 8, Equations 38 can be rewritten as follows:

   δN i = -knδd i δT i = ktd i δα j elastic regime δT i = tan ϕg (N i -knδd i ) ξ i -T i plastic regime (39) where: ξ i is the sign of T i + ktd i δα j ; j = 1 when i = 1, 2; j = 2 when i = 3, 4; plastic regime is reached when ktd i δα j + T i tan ϕg (N i -knδd i ), otherwise it is in elastic regime. To facilitate the derivation, I p i and I e i are introduced as indicator functions of the contact state, expressed as follow:

I p i =
1 in plastic regime 0 in elastic regime ; I e i = 1 -I p i (40)

Thus, the constitutive relations can be expressed as:

δN i = -knδd i δT i = B i δα j -A i δd i + C i (41)
where:

   A i = I p i knξ i tan ϕg B i = I e i ktd i C i = I p i (ξ i tan ϕgN i -T i )

Fig. 3 :

 3 Fig. 3: Global and local coordinate system transformation by employing Euler angles in 3D conditions.

Fig. 5 :

 5 Fig. 5: Mechanical description of Hexagon A.

  Fig. 6: Mechanical description of Hexagon B.

Fig. 7 :

 7 Fig. 7: Calculation of capillary force in a liquid bridge based on the Virtual Work Principle.

Fig. 11 :

 11 Fig. 11: Surface energy in a liquid bridge between two grains in contact, as a function of the water volume.

Fig. 12 :

 12 Fig. 12: Definition of the vertical filling angle β v and the inclinated filling angle β i .

  bridges are added to the contact forces. At mesoscale, V σ are replaced by V σ cap in Equation 20 and Equation

Fig. 13 :Fig. 15 :

 1315 Fig. 13: Flowchart of the proposed multi-scale approach.

Fig. 16 :

 16 Fig. 16: Comparison between the numerical results (w c = 0%) of this study and results from literatures.

Fig. 17 :

 17 Fig. 17: Normalized granular spreading length versus normalized time of experimental and numerical (w c = 0%) results.

Fig. 18 :

 18 Fig. 18: Comparison of failure progress of granular collapse between simulations and experiments for initial aspect ratio (a = 1): left column shows the results of dry material while right column shows the results of unsaturated material with water content w c = 10%.

Fig. 19 :

 19 Fig.19: Meso-scale variable distributions of 3D-H model inside SPH particles A, B and C at three states. The selected SPH particle positions are marked as red cross in Figure14a.

Fig. 20 :

 20 Fig.20: Meso-scale variable distributions of 3D-H model inside SPH particles D, E and F at three states. The selected SPH particle positions are marked as red cross in Figure14a.

Table 1 :

 1 Parameters selected in calibration and prediction phases

			kn/rg(Pa)	kt/kn	e 0	ϕg( • )	wc
			1.90 × 10 9	0.6	0.53	25	10%
	Test	Aspect ratio Initial length Initial height particle number Initial cell size Water content
		a	L i (mm)	H i (mm)		-	(mm)	wc
	cell500	1.0	100	100		500	10	0%
	cell1183	1.0	100	100		1183	7.5	0%
	cell2048	1.0	100	100		2048	6.3	0%
	cell4000	1.0	100	100		4000	5	0%
	cell8125	1.0	100	100		8125	4	0%
	a0.5	0.5	100	50		1024	6.3	0%
	a1.0	1.0	100	100		2048	6.3	0%
	a1.5	1.5	100	150		3072	6.3	0%
	a3.0	3.0	100	300		6144	6.3	0%
	a7.0	7.0	100	700		14336	6.3	0%
	a9.0	9.0	100	900		18432	6.3	0%
	dry	1.0	100	100		8125	6.3	0%
	unsaturated	1.0	100	100		8125	6.3	10%

Table 2 :

 2 Dimension and discretization parameters for all numerical simulations