Vinicius Ottonio

O Gonçalves

Walter Henrique

S M Talon

Vinicius Kartnaller

Fabricio Venancio

João Cajaiba

Thierry Cabioc'

Jean-Marc Clacens

Frédéric Richard
email: frederic.richard@univ-poitiers.fr

Thierry Cabioc'h

Hydrodeoxygenation of m-cresol as a depolymerized lignin probe molecule: Synergistic effect of NiCo supported alloys

Keywords: m-cresol, hydrodeoxygenation, NiCo Alloy, bimetallic

Introduction

Bio-oil from biomass pyrolysis is becoming attractive to produce non-petroleum based and value-added biofuels and chemicals [START_REF] Huber | Synergies between Bio-and Oil Refineries for the Production of Fuels from Biomass[END_REF][START_REF] Li | Catalytic Transformation of Lignin for the Production of Chemicals and Fuels[END_REF][START_REF] Kan | Lignocellulosic biomass pyrolysis: A review of product properties and effects of pyrolysis parameters[END_REF][START_REF] Dabros | Transportation fuels from biomass fast pyrolysis, catalytic hydrodeoxygenation, and catalytic fast hydropyrolysis[END_REF]. It contains a range of various organic compounds such as phenols, acids, aldehydes, and ketones, depending on the pyrolysis conditions and the biomass source [START_REF] Huber | Synthesis of transportation fuels from biomass: chemistry, catalysts, and engineering[END_REF]. The as-produced bio-oil remains low valuable as a product since it requires a decrease of its oxygen content to be more profitable.

Hydrodeoxygenation (HDO) catalytic process can be used for oxygen removal. In recent years, much efforts were done in order to design a cost-effective catalyst for a HDO efficient process [START_REF] Bu | A review of catalytic hydrodeoxygenation of lignin-derived phenols from biomass pyrolysis[END_REF][START_REF] Arun | Green diesel synthesis by hydrodeoxygenation of bio-based feedstocks: Strategies for catalyst design and development[END_REF][START_REF] Mortensen | A review of catalytic upgrading of bio-oil to engine fuels[END_REF][START_REF] Li | Heterogeneous sulfur-free hydrodeoxygenation catalysts for selectively upgrading the renewable bio-oils to second generation biofuels[END_REF][START_REF] Jin | Catalytic Upgrading of Biomass Model Compounds: Novel Approaches and Lessons Learnt from Traditional Hydrodeoxygenationa Review[END_REF][START_REF] Baloch | Recent advances in production and upgrading of bio-oil from biomass: A critical overview[END_REF][START_REF] Peters | Anisole and Guaiacol Hydrodeoxygenation Reaction Pathways over Selected Catalysts[END_REF].

In this context, the catalyst must provide an acceptable combination of the followed features: stability, sustainability, deoxygenation rate, and selectivity toward more valuable products. Firstly, commercial sulfided solids such as Mo-based sulfides appeared as efficient catalysts for HDO and were extensively studied due to their interesting features as catalysts [START_REF] Furimsky | Catalytic hydrodeoxygenation[END_REF]. However, in the research of more sustainable processes, sulfided catalysts become limited since their use requires an S-containing agent (e.g. H2S or dimethyl disulfide) to be present in the feed to maintain the catalyst activity. Unfortunately, their presence can lead to a possible sulfur contamination of the upgraded bio-oil liquids [START_REF] Ryymin | Insight to sulfur species in the hydrodeoxygenation of aliphatic esters over sulfided NiMo/γ-Al 2 O 3 catalyst[END_REF]. Noble metal catalysts (e.g. containing Pt or Pd) also exhibit attractive catalytic properties such as high activity and stability without the use of unsustainable materials such as sulfur [START_REF] Liu | Hydrodeoxygenation of benzofuran over silica-alumina-supported Pt, Pd, and Pt-Pd catalysts[END_REF][START_REF] Nimmanwudipong | Catalytic conversion of guaiacol catalyzed by platinum supported on alumina: reaction network including hydrodeoxygenation reactions[END_REF]. However, these solids have high cost, meaning that the final product would be too expensive to compete with their petro-based similar compounds.

In the perspective to use an efficient catalyst with an affordable cost, the catalytic properties of nickel-based catalyst materials were recently investigated for the deoxygenation of triglycerides [START_REF] Kordulis | Development of nickel based catalysts for the transformation of natural triglycerides and related compounds into green diesel: a critical review[END_REF]. It is often proposed that active metallic phase containing at least two components exhibits better catalytic characteristics compared to monometallic [START_REF] Alonso | Bimetallic catalysts for upgrading of biomass to fuels and chemicals[END_REF]. It is often reported that the catalytic properties of Ni-Co alloys were more interesting than the monometallics for various reactions such as partial oxidation of CH4 [START_REF] Koh | Hydrogen or synthesis gas production via the partial oxidation of methane over supported nickelcobalt catalysts[END_REF]. CO hydrogenation [START_REF] Ishihara | Effect of alloying on CO hydrogenation activity over SiO2-supported Co-Ni alloy catalysts[END_REF][START_REF] Ishihara | Hydrogenation of carbon monoxide over SiO2-supported Fe-Co, Co-Ni and Ni-Fe bimetallic catalysts[END_REF], steam reforming of various substrates (CH4 [START_REF] Hu | Investigation of steam reforming of acetic acid to hydrogen over Ni-Co metal catalyst[END_REF][START_REF] Zhang | Development of stable bimetallic catalysts for carbon dioxide reforming of methane[END_REF][START_REF] Zhang | Effects of metal content on activity and stability of Ni-Co bimetallic catalysts for CO 2 reforming of CH 4[END_REF][START_REF] Sajjadi | Hydrogen production via CO2reforming of methane over Cu and Co doped Ni/Al2O3 nanocatalyst: impregnation versus sol-gel method and effect of process conditions and promoter[END_REF][START_REF] San-José-Alonso | Co and bimetallic Ni-Co catalysts for the dry reforming of methane[END_REF][START_REF] Das | Facile synthesis of highly disperse Ni-Co nanoparticles over mesoporous silica for enhanced methane dry reforming[END_REF], biomass tar [START_REF] Wang | Catalytic performance and characterization of Ni-Co catalysts for the steam reforming of biomass tar to synthesis gas[END_REF], ethanol [START_REF] Koh | Hydrogen or synthesis gas production via the partial oxidation of methane over supported nickelcobalt catalysts[END_REF][START_REF] Andonova | Structure and redox properties of Co promoted Ni/Al 2 O 3 catalysts for oxidative steam reforming of ethanol[END_REF][START_REF] Chen | Support and alloy effects on activity and product selectivity for ethanol steam reforming over supported nickel cobalt catalysts[END_REF], acetic acid [START_REF] Hu | Investigation of steam reforming of acetic acid to hydrogen over Ni-Co metal catalyst[END_REF] or glycerin [START_REF] Luo | Hydrogen generation from liquid reforming of glycerin over Ni-Co bimetallic catalyst[END_REF]). The promoting effects due to the presence of alloys are mostly attributed to a better dispersion of metallic phase [START_REF] Zhang | Development of stable bimetallic catalysts for carbon dioxide reforming of methane[END_REF][START_REF] Xu | Biogas reforming for hydrogen production over nickel and cobalt bimetallic catalysts[END_REF], to electronic effects [START_REF] Ishihara | Effect of alloying on CO hydrogenation activity over SiO2-supported Co-Ni alloy catalysts[END_REF], or to the enhancement of their reducibility properties [START_REF] Wang | Catalytic performance and characterization of Ni-Co catalysts for the steam reforming of biomass tar to synthesis gas[END_REF][START_REF] Andonova | Structure and redox properties of Co promoted Ni/Al 2 O 3 catalysts for oxidative steam reforming of ethanol[END_REF]. Such behavior was also observed for Ni based catalysts. The association of Ni to different metals such as Cu [START_REF] Khromova | Anisole hydrodeoxygenation over Ni-Cu bimetallic catalysts: The effect of Ni/Cu ratio on selectivity[END_REF][START_REF] Ardiyanti | Catalytic hydrotreatment of fast-pyrolysis oil using non-sulfided bimetallic Ni-Cu catalysts on a δ-Al 2 O 3 support[END_REF], Re [START_REF] Feng | Aqueous-phase hydrodeoxygenation of 4propylphenol as a lignin model to n-propylbenzene over Re-Ni/ZrO2 catalysts[END_REF][START_REF] Yang | Geometric and electronic effects of bimetallic Ni-Re catalysts for selective deoxygenation of m-cresol to toluene[END_REF], Pd [START_REF] Lai | Hydrodeoxygenation of guaiacol over Ni@Pd and Ni@Pt bimetallic overlayer catalysts[END_REF], Pt [START_REF] Lai | Hydrodeoxygenation of guaiacol over Ni@Pd and Ni@Pt bimetallic overlayer catalysts[END_REF][START_REF] Do | Bimetallic effects in the hydrodeoxygenation of meta-cresol on γ-Al 2 O 3 supported Pt-Ni and Pt-Co catalysts[END_REF], Fe [START_REF] Leng | NiFe/γ-Al 2 O 3: a universal catalyst for the hydrodeoxygenation of bio-oil and its model compounds[END_REF][START_REF] Nie | Selective conversion of m-cresol to toluene over bimetallic Ni-Fe catalysts[END_REF][START_REF] Shafaghat | Catalytic hydrodeoxygenation of simulated phenolic bio-oil to cycloalkanes and aromatic hydrocarbons over bifunctional metal/acid catalysts of Ni/HBeta, Fe/HBeta and NiFe/HBeta[END_REF][START_REF] Liu | Hydrodeoxygenation of m-cresol over bimetallic NiFe alloys: Kinetics and thermodynamics insight into reaction mechanism[END_REF][START_REF] Han | The synergistic effect between Ni sites and Ni-Fe alloy sites on hydrodeoxygenation of lignin-derived phenols[END_REF] or Co [START_REF] Huynh | Hydrodeoxygenation of Phenol as a Model Compound for Bio-oil on Non-noble Bimetallic Nickel-based Catalysts[END_REF][START_REF] Huynh | Hydrodeoxygenation of Bio-Oil on Bimetallic Catalysts: From Model Compound to Real Feed[END_REF][START_REF] Zhou | Catalytic conversion of guaiacol to alcohols for bio-oil upgrading[END_REF][START_REF] Zhou | Water-Assisted Selective Hydrodeoxygenation of Guaiacol to Cyclohexanol over Supported Ni and Co Bimetallic Catalysts[END_REF][START_REF] Tran | Vapor-phase hydrodeoxygenation of guaiacol on Al-MCM-41 supported Ni and Co catalysts[END_REF][START_REF] Liu | Significant Promotion of Surface Oxygen Vacancies on Bimetallic CoNi Nanocatalysts for Hydrodeoxygenation of Biomass-derived Vanillin to Produce Methylcyclohexanol[END_REF][START_REF] Hari | Production of diesel fuel by the hydrotreatment of jatropha oil derived fatty acid methyl esters over γ-Al2O3 and SiO2 supported NiCo bimetallic catalysts[END_REF] leads to beneficial changes in their HDO catalytic properties. The presence of Ni-Cu alloys in the corresponding bimetallic catalysts favored the deoxygenation rate of anisole into methylcyclohexane under high pressure (6 MPa) at 280 °C. This promotion was linked to an increase in the nickel content [START_REF] Khromova | Anisole hydrodeoxygenation over Ni-Cu bimetallic catalysts: The effect of Ni/Cu ratio on selectivity[END_REF]. In the same way, Yang et al. [START_REF] Yang | Geometric and electronic effects of bimetallic Ni-Re catalysts for selective deoxygenation of m-cresol to toluene[END_REF] showed a beneficial effect of the presence of Re in Ni-Re catalysts for the HDO of m-cresol under atmospheric pressure at 300°C. This effect was related to an increase of the Ni dispersion due to the presence of Re allowing the formation of active Ni-Re alloys. In the same way, Han et al. [START_REF] Han | The synergistic effect between Ni sites and Ni-Fe alloy sites on hydrodeoxygenation of lignin-derived phenols[END_REF] reported that the bimetallic Ni-Fe catalyst was more active than the monometallic Ni one in HDO of phenol due to the presence of Ni-Fe alloy. In this last case, the bimetallic catalyst presented the highest HDO activity which was assigned to an increase of the adsorption strength of cyclohexanol on Ni-Fe alloy particles. This alcohol is an oxygenated product produced during the HDO of phenol, which can be transformed into cyclohexane by hydrogenolysis. It was also reported in literature that the presence of iron can modify the product distribution observed on Ni-based catalysts during HDO of m-cresol (300 °C, atmospheric pressure) [START_REF] Nie | Selective conversion of m-cresol to toluene over bimetallic Ni-Fe catalysts[END_REF]. Indeed, methylcyclohexanone was the main product over Ni/SiO2 whereas toluene was predominant over Ni-Fe/SiO2. According to the study, this change of selectivity in the presence of iron could suppress the hydrogenation of the phenol aromatic ring. The authors proposed that iron could cause a repulsion with the  electron system of the phenol aromatic ring thereby inhibiting the hydrogenating properties of Ni-Fe particles. Using an acidic support (HZSM5 zeolite), Huynh et al. [START_REF] Huynh | Hydrodeoxygenation of Phenol as a Model Compound for Bio-oil on Non-noble Bimetallic Nickel-based Catalysts[END_REF] demonstrated a beneficial effect of Co when it was associated to Ni for the HDO of phenol (5 MPa, 250 °C). In another case, studying guaiacol as model molecule at 400 °C under atmospheric pressure, Tran et al. [START_REF] Tran | Vapor-phase hydrodeoxygenation of guaiacol on Al-MCM-41 supported Ni and Co catalysts[END_REF] observed that the presence of Co enhanced the HDO activity of Ni based catalysts by favoring the dispersion of nickel metal oxides. Nevertheless, in both cases, the promoting effect of Co was difficult to estimate due to the very high level of conversion of phenol (close to 98 %). Recently, Liu et al. [START_REF] Liu | Significant Promotion of Surface Oxygen Vacancies on Bimetallic CoNi Nanocatalysts for Hydrodeoxygenation of Biomass-derived Vanillin to Produce Methylcyclohexanol[END_REF] reported that bimetallic NiCo alloys supported on Al2O3 were efficient catalysts for the HDO of vanillin, methylcyclohexanol being selectively produced at 2 MPa and 200 °C. Interestingly, the Co/Ni atomic ratio of 2 allowed to obtain a very high HDO performance since the yield in methylcycloghexanol was close to 100%. Such promoting effect was explained by the presence of electron-rich Co 0 species in CoNi nanoparticles which favored the activation of oxygen-containing groups.

In the present work, the investigation of nickel and cobalt alloys supported on γ-Al2O3 for the hydrodeoxygenation of m-cresol as bio-oil model compound under high hydrogen pressure was studied for the first time. The interest of studying such materials is based on the non-use of deleterious materials and on their potential effectiveness as catalysts for HDO. The relationship between the solids' properties and their performances as catalyst is discussed further herein.

Experimental

Catalyst preparation

Ni-Co catalysts were prepared by incipient wetness impregnation method of γalumina. Different molar ratios between cobalt and nickel were used in the synthesis using the correct amounts of Ni(NO3)2.6H2O and Co(NO3)2.6H2O as precursors (purchased from Sigma-Aldrich). In order to evaluate the bimetallic catalyst system, Ni/Al2O3 and Co/Al2O3 were also synthesized. The total metal loading (Co + Ni) was set to 10 wt.%. After impregnation, the catalyst was left for maturation during 5 days at 25 ºC. This aging time was already reported to ensure homogeneity between different elements in the case of coimpregnation procedure [START_REF] Ungureanu | Synthesis of highly thermostable copper-nickel nanoparticles confined in the channels of ordered mesoporous SBA-15 silica[END_REF].

The impregnated solid was then calcinated at 500 ºC (1.5 ºC min -1 ) for 5 h in a static oven. Further, the obtained catalysts were pelleted, crushed and sieved (250-315 µm).

Catalyst characterization

The Ni and Co contents of the catalysts were determined by atomic absorption using an AA200 Perkin Elmer equipment. Before analysis, the catalysts were solubilized in an acidic aqueous solution containing both hydrochloric and nitric acids. N2 adsorption-desorption isotherms were measured on a Micromeritics ASAP 2000 analyzer at -196 ºC. Prior to N2 adsorption, the powder samples were degassed under secondary vacuum at 200 ºC overnight. The specific surface area was calculated using the Brunauer-Emmett-Teller (BET) method. The total pore volume was calculated from the adsorbed volume of nitrogen at P/P0 equal to 0.99. H2-TPR measurements were performed for 100 mg of each catalyst in Autochem 2920.

The solid was first pretreated under He flow (20 mL min -1 ) at 200 ºC (10 ºC min -1 ) for 30 min.

Temperature was then adjusted to 50 ºC and the reduction measurements using TCD detector started under 10% H2 in Ar (20 mL min -1 ) until 1000 ºC (5 ºC min -1 ) and kept in isotherm for 90 min.

X-ray diffraction (XRD) analysis were performed on a Bruker AXS D8 Advance diffractometer to determine the crystal structure of nanoparticles. The X-ray diffractometer was operated at 40 kV and 35 mA with a Co-K ( = 1.7889 Å) X-ray source. Each analysis was performed for 2 values from 20 to 100°, with steps of 0.05° and fixed acquisition time of 6s per step. Rietveld refinements of the diffractograms were performed using the Materials Analysis Using Diffraction (MAUD) software [START_REF] Lutterotti | MAUD: a friendly Java program for material analysis using diffraction[END_REF]. Prior to analysis, the Ni-Co solids were reduced in situ at 500 ºC (5 ºC min -1 ) under 4 MPa of H2 during 3 h by a flow of hydrogen (4.7 NL h -1 ). After this procedure, samples were cooled down under hydrogen until room temperature and passivated under atmospheric pressure for 3 h using 5% O2/He (80 mL min - 1 ).

CO chemisorption was determined by pulse method at 30 ºC. The catalyst was reduced in-situ at 500 ºC (3 h) at 4 MPa of H2 as described above for X-ray analysis. The samples were transferred under argon atmosphere to an equipment for CO chemisorption measurements.

First, the solid was pretreated with helium (30 mL min -1 ) at 200 ºC for 30 minutes.

Temperature was then set to 30 °C and CO pulses were injected until saturation. The measurements were detected by TCD.

The acidity of the support and catalysts was determined by NH3-TPD. Experiments were performed on a Micromeritics Autochem 2910 apparatus. All samples were pretreated with pure He at 200°C during 1 h. The temperature was then lowered to 100 °C. The solid was saturated with anhydrous NH3 (10 vol% in He). Thereafter, it was purged under He, then the temperature was raised to 600 °C (heating rate: 5 °C min -1 ). The amount of NH3 desorbed was quantified by TCD.

Carbon content was determined using an elementary analyzer (NA2100 Protein Thermoquest, CE instruments). The catalysts were exposed to 1020 ºC under oxygen atmosphere. The C content was quantified by the analysis of CO2 provided by the combustion of the solid.

Catalytic experiments

Ni-Co catalysts were tested in a continuous fixed-bed reactor for the deoxygenation of m-cresol at 340 °C under 4 MPa of total pressure. Prior to the reaction, each catalyst was reduced in-situ at 500 ºC (5 ºC min -1 ) under 4 MPa of pure hydrogen (4.7 L/h) for 3 h.

Temperature was then set to the reaction temperature (340 ºC). The reactor was fed by a liquid model feed composed of m-cresol (4.65 wt.%) and n-heptane as an internal standard (2.52 wt.%), in dodecane (solvent). The used molar ratio of hydrogen to m-cresol was 486 NL/L and the H2 partial pressure was 3.24 MPa. All chemicals used in this study were purchased from Sigma-Aldrich, except for pure H2 (Air Liquide).

During the reaction, the exit line of the reactor was refrigerated by a Minichiller-Huber at 10 ºC and the condensed samples were analyzed every hour using a Varian 430 chromatograph with a DB1 capillary column (length: 30 m, inside diameter: 0.320 mm, film thickness: 5 µm) and a FID detector.

The conversion and selectivity were calculated based on the following Eqs. ( 1) and ( 2), standard deviation of all values were estimated to 5%: Experiments were carried out during 25 h on stream. Different conversion levels were obtained using several space times ( in g h mol -1 ), changing the reactant flow (Eq. ( 3)).

  0 0 % 100 CRE CRE CRE CRE CC X C   (1) 0 (%) 100 i i CRE CRE C S CC   (2)
CRE W F   (3)
where W is the weight of catalyst and CRE F the molar flow of m-cresol.

The reactant and product concentration can be used to determine kinetic parameters of transformation of m-cresol. Assuming this reaction as 1st-order (Eq. ( 4)):

 

1 ln 1 TOT CRE kX   ( 4
)
where kTOT is the reaction rate constant (mmol g -1 h -1 ).

The turnover frequency (TOF, in h -1 ) was calculated by Eq. ( 5):

TOT k TOF M  ( 5
)
where M is the CO uptake (in mmol g -1 ) obtained from chemisorption experiments.

3 Results and discussion

Catalyst characterization

Ni and Co contents in the prepared catalysts were very close to the nominal values as

shown in Table 1. The N2 adsorption and desorption isotherms for the support and the catalysts reduced in situ at 4 MPa and 340 ºC are presented in Fig. 1. All catalysts possessed a mesoporous structure very similar to their support presenting a type IV isotherm. The fresh catalysts in oxide form also presented the same kind of isotherm (not shown here). Table 1 summarizes the textural properties of the support, the as prepared and reduced catalysts. For supported materials, raw values were recalculated to remove the contribution of oxide phases weight on textural properties (values in brackets in Table 1) allowing to access information on stability of the support during the impregnation-reduction steps. It can be noted that the values of surface area and pore volume of all solids slightly decreased after introduction of oxide phases. This result might be attributed to two factors, the formation of Ni and Co oxides blocking the support pores and the formation of a solid structure in the form of spinel species [START_REF] Andonova | Structure and redox properties of Co promoted Ni/Al 2 O 3 catalysts for oxidative steam reforming of ethanol[END_REF]. The reduction process led to a slight decrease of both SBET and pore volume (less than 3%), showing that the reduction step preserves the textural properties of the catalysts.

However, the surface areas of both Co/Al2O3 and Ni/Al2O3 were similar, close to 195 m 2 g -1 .

Fig. 2 shows the TPR profiles of Ni/Al2O3, Co/Al2O3 and NiCo/Al2O3 with different molar ratios. The H2 uptakes obtained from TPR measurements as well as the theoretical value assuming fully reduction of the metals can be found in Table 2. From these values, Ni and Co species were completely reduced to metallic species. The order of H2 consumption followed during the reduction process. Indeed, it was proposed that the reducibility of Ni and Co species was simultaneously improved in NiCo samples [START_REF] Xu | Biogas reforming for hydrogen production over nickel and cobalt bimetallic catalysts[END_REF][START_REF] Liu | Significant Promotion of Surface Oxygen Vacancies on Bimetallic CoNi Nanocatalysts for Hydrodeoxygenation of Biomass-derived Vanillin to Produce Methylcyclohexanol[END_REF]. In addition, the reduction of spinel NiAl2O4 species occurs at much higher temperature with the increase of cobalt content as already reported [START_REF] Andonova | Structure and redox properties of Co promoted Ni/Al 2 O 3 catalysts for oxidative steam reforming of ethanol[END_REF]. In order to verify that all metallic phases were reduced under the reduction step carried out prior the catalytic test (500 ºC, 4 MPa), additional TPR experiments were carried out for each pre-treated sample. The reactor was open in a glovebox filled with Argon and the reduced catalysts were transferred under Ar atmosphere to perform TPR. Over all reduced samples, no consumption of hydrogen was detected indicating that all oxide phases were totally reduced at 500 °C under 4 MPa of hydrogen pressure. angle gradually with increasing Co amount. Fig. 5 shows a linear relationship between the lattice parameter (calculated by Rietveld refinements) and the composition of metallic phase estimated by the Ni/(Ni+Co) atomic ratio. This result is in agreement with the Vegard's law [START_REF] Cullity | Elements of x-ray diffraction[END_REF] and confirms the presence of Ni-Co alloy in bimetallic samples. In addition, the size of metallic nanoparticles seems to increase with the introduction of cobalt, in agreement with results reported by Wang et al. [START_REF] Wang | Catalytic performance and characterization of Ni-Co catalysts for the steam reforming of biomass tar to synthesis gas[END_REF].

As expected, the acidity of the support and catalysts, determined by NH3-TPD, was close to 1.0 µmol g -1 (Table 2), irrespective the presence and the nature of the metallic phase.

The value of CO uptake for Ni/Al2O3 is about twice higher than the one of Co/Al2O3 (Table 2), in agreement with results obtained by Takanabe et al. [START_REF] Takanabe | Titania-supported cobalt and nickel bimetallic catalysts for carbon dioxide reforming of methane[END_REF]. Interestingly, the bimetallic solid with a Ni/Co molar ratio close to 3 (NiCo(3:1)/Al2O3), shows the highest CO uptake. This alloy composition seems to favour the amount of sites allowing to adsorb CO.

HDO of m-cresol

The HDO of m-cresol was carried out at 4 MPa and 340 ºC in order to evaluate the effect of alloying nickel and cobalt. For each catalyst, experiments were performed for 25 h on stream and space time was changed every 5 h in order to study the evolution of the reaction products with the conversion of m-cresol (Fig. 6). In a run with high value of space time using the bare support as catalyst, no conversion of the reactant was observed. All catalysts were stable under our experimental conditions during the whole set of experiments, despite the presence of carbon over spent catalysts (Table 3). Fig. 7 shows a linear correlation between -ln (1-XCRE) and space time (τ) indicating that the reaction can be modeled by a first-order kinetic, as already reported for several types of catalysts [START_REF] Gonçalves | Hydrodeoxygenation of m-cresol over nickel and nickel phosphide based catalysts. Influence of the nature of the active phase and the support[END_REF][START_REF] Gonçalves | Kinetics of the hydrodeoxygenation of cresol isomers over Ni2P/SiO2: Proposals of nature of deoxygenation active sites based on an experimental study[END_REF]. Thus, the overall rate constant (kTOT) determined for each catalyst (Table 4) showed that NiCo (3:1)/Al2O3 was the most active of the studied series. It was about 1.2 times more active than Ni/Al2O3 and 9 times more active than Co/Al2O3, showing a promoting effect of the presence of Ni-Co alloy for HDO of m-cresol. In addition, the Ni-Co alloys composition had an effect on their HDO catalytic properties: in our studied samples, the higher the nickel content is, the higher the HDO activity is. Indeed, the order of activity among Ni-Co catalysts followed NiCo (3:1)/Al2O3 > NiCo(1:1)/Al2O3 > NiCo(1:3)/Al2O3 . Interestingly, this trend was the same that the one obtained from the CO uptake (Table 2), as shown in Fig. 8. The evolution of both kTOT and CO uptakes as a function of the active phase composition (determined by the Ni/(Ni + Co) ratio) followed the same trend. These results seem to indicate that the sites allowing to adsorb CO are the active sites in HDO of m-cresol. Consequently, such technique seems to be well appropriate to measure the amount of HDO active sites.

Recently, it was proposed that the active sites present on CoNi nanoparticles could be surface electron-rich Co 0 species, which could interact with oxygen atom and favor the deoxygenation rate [START_REF] Liu | Significant Promotion of Surface Oxygen Vacancies on Bimetallic CoNi Nanocatalysts for Hydrodeoxygenation of Biomass-derived Vanillin to Produce Methylcyclohexanol[END_REF]. Due to the lower amount of active sites estimated by CO adsorption for Ni/Al2O3 (35 mol g -1 ) compared to NiCo(3:1)/Al2O3 (46 mol g -1 ), the TOF value of Ni/Al2O3 was about 1.1 times higher than the one of NiCo(3:1)/Al2O3 (Table 4), even if the total activity (estimated by kTOT) of NiCo(3:1)/Al2O3 was 1.2 times higher than Ni/Al2O3. In addition, the TOF values reported in Table 4 indicates that the activity of Ni/Al2O3 was higher than the one of Co/Al2O3, in agreement with results obtained for HDO of phenol [START_REF] Huynh | Hydrodeoxygenation of Phenol as a Model Compound for Bio-oil on Non-noble Bimetallic Nickel-based Catalysts[END_REF].

In Under the same experimental conditions (340 °C, 3 MPa), the HDO of m-cresol was already reported over various supported catalytic systems, such as molybdenum sulfide [START_REF] Gonçalves | Hydrodeoxygenation of Cresols Over Mo/Al2O3 and CoMo/Al2O3 Sulfided Catalysts[END_REF],

molybdenum oxide [START_REF] Gonçalves | Effect of the support on the hydrodeoxygenation of m-cresol over molybdenum oxide based catalysts[END_REF] and nickel phosphides [START_REF] Gonçalves | Hydrodeoxygenation of m-cresol over nickel and nickel phosphide based catalysts. Influence of the nature of the active phase and the support[END_REF][START_REF] Gonçalves | Kinetics of the hydrodeoxygenation of cresol isomers over Ni2P/SiO2: Proposals of nature of deoxygenation active sites based on an experimental study[END_REF][START_REF] Gonçalves | Effect of P/Ni ratio on the performance of nickel phosphide phases supported on zirconia for the hydrodeoxygenation of m-cresol[END_REF]. Only CoMoS/Al2O3 was slightly more active than the NiCo(3:1)/Al2O3 (about 1.1 times more active). Nevertheless, the selectivity of this sulfide phase was different since toluene was the main product over CoMo/Al2O3.

In addition, except to Co/Al2O3, the particle size increases linearly with the amount of Co (Fig. 5). This trend was rather different that the one observed for the HDO activity, indicating that the particle size of the metallic nanoparticles (monometallic or alloy) was

probably not a key parameter to explain their catalytic properties, as previously suggested [START_REF] Huynh | Hydrodeoxygenation of Phenol as a Model Compound for Bio-oil on Non-noble Bimetallic Nickel-based Catalysts[END_REF][START_REF] Takanabe | Titania-supported cobalt and nickel bimetallic catalysts for carbon dioxide reforming of methane[END_REF]. °C. Under these experimental conditions, m-cresol was mainly deoxygenated into methylcyclohexane ( 56 mol%) and methylcyclohexenes ( 33 mol%). Toluene was also detected but in low quantities (less than 2 mol%). An oxygenated compound, methylcyclohexanone, was also observed in the reactional mixture; its amount being close to 10 mol%. As indicated in Table 4, the presence of Co favored the total deoxygenation of mcresol since the HDO values were higher than 90 mol% for the three bimetallic catalysts and for Co/Al2O3, whereas deoxygenation selectivity was lower for Ni/Al2O3 (equal to 84 mol.%)

showing that the presence of Co 0 species favored the activation of m-cresol. Such behavior was already reported by Liu et al. [START_REF] Liu | Significant Promotion of Surface Oxygen Vacancies on Bimetallic CoNi Nanocatalysts for Hydrodeoxygenation of Biomass-derived Vanillin to Produce Methylcyclohexanol[END_REF] and attributed to the presence of surface oxygenvacancy which decorated bimetallic CoNi nanoparticles. In addition, the proportion of methylcyclohexane was the highest on NiCo(3:1)/Al2O3 (56 mol%) and the lowest on Co/Al2O3 (41 mol.%). Our results are in accordance with results reported by Huynh et al. [START_REF] Huynh | Hydrodeoxygenation of Phenol as a Model Compound for Bio-oil on Non-noble Bimetallic Nickel-based Catalysts[END_REF] using phenol as probe molecule, cyclohexane being the main product (250°C, 5 MPa).

However, over Co/Al2O3, methylcyclohexenes become the major products, showing the lower hydrogenating properties of Co compared to Ni. In the same way, Co seems more selective into toluene than Ni. According to Table 4, the degree of deoxygenation (HDO %) is favored by increasing the Co content. Zhou et al. [START_REF] Zhou | Catalytic conversion of guaiacol to alcohols for bio-oil upgrading[END_REF] observed that the selectivity into cyclohexanol in the HDO of guaiacol was improved using nickel catalyst promoted by cobalt; this catalyst being supported on carbon nanotubes; meanwhile pure nickel mostly produced catechol, a product obtained without oxygen removal. Similarly, for the hydrodeoxygenation of guaiacol under 1 MPa and at 300 ºC, Mochizuki et al. [START_REF] Mochizuki | Deoxygenation of guaiacol and woody tar over reduced catalysts[END_REF] showed that the degree of deoxygenation over Co/SiO2 was higher, leading mainly to aromatic and naphthenic compounds, than Ni/SiO2, which gave cyclohexanol as main product.

Fig. 10 shows the behavior of the selectivity from the transformation of m-cresol as function of space times. At low conversion levels (low contact time values) alkene isomers (methylcyclohexenes) were preferentially formed, their selectivity decreased by increasing the conversion of the phenolic reactant. Regarding the evolution of selectivity with m-cresol conversion, 3-methylcyclohexanone followed the same trend than methylcyclohexenes. At higher space times, methylcyclohexane selectivity increases until it becomes the major product. The selectivity into toluene was not modified in the range of contact time studied for each catalyst.

From previous results over several active phases (CoMoS, Ni, Ni2P, MoOx) [START_REF] Gonçalves | Hydrodeoxygenation of m-cresol over nickel and nickel phosphide based catalysts. Influence of the nature of the active phase and the support[END_REF][START_REF] Gonçalves | Kinetics of the hydrodeoxygenation of cresol isomers over Ni2P/SiO2: Proposals of nature of deoxygenation active sites based on an experimental study[END_REF][START_REF] Gonçalves | Hydrodeoxygenation of Cresols Over Mo/Al2O3 and CoMo/Al2O3 Sulfided Catalysts[END_REF][START_REF] Gonçalves | Effect of the support on the hydrodeoxygenation of m-cresol over molybdenum oxide based catalysts[END_REF][START_REF] Gonçalves | Effect of P/Ni ratio on the performance of nickel phosphide phases supported on zirconia for the hydrodeoxygenation of m-cresol[END_REF],

it can be admitted that the transformation of phenolic compounds proceeds mainly via two parallel routes. Products of m-cresol ring hydrogenation are favored over NiCo catalysts, and are formed by the so-called HYDdrogenation route, which proceeds from the adsorption of m-cresol leading first to 3-methylcyclohexanone. The hydrogenation of its C=O bond leads to 3-methylcyclohexanol. This alcohol was never observed in these experimental conditions owing to its high reactivity through dehydration reaction. Such results were also described in the literature when acid supports like Al2O3 were used for HDO of alkylphenols [START_REF] Gonçalves | Effect of the support on the hydrodeoxygenation of m-cresol over molybdenum oxide based catalysts[END_REF][START_REF] Romero | Hydrodeoxygenation of 2-ethylphenol as a model compound of bio-crude over sulfided Mo-based catalysts: Promoting effect and reaction mechanism[END_REF][START_REF] Bouvier | Effect of H2S and CO on the transformation of 2-ethylphenol as a model compound of bio-crude over sulfided Mobased catalysts: propositions of promoted active sites for deoxygenation pathways based on an experimental study[END_REF].

The acidity of the catalytic support alone, determined by NH3-TPD, was close to 1.04 mol g -1 , which was practically unmodified in the presence of reduced Ni and/or Co nanoparticles (Table 2). Methylcylohexene isomers (MCHes) were hence produced from dehydration of 3methylcyclohexanol, which may occur on alumina acid sites. Interestingly, 1methylcyclohexene was always the main alkene isomer instead of the expected isomers, being about 63 mol% of all alkene isomers (Table 4). This value was close to the percentage calculated from thermodynamics (65 mol%), suggesting isomerization until thermodynamic equilibrium [START_REF] Peereboom | Experimental and calculated thermodynamic data for the isomeric methylcyclohexenes and methylenecyclohexane[END_REF]. This observation was also reported for catalysts exhibiting acidic properties [START_REF] Gonçalves | Hydrodeoxygenation of Cresols Over Mo/Al2O3 and CoMo/Al2O3 Sulfided Catalysts[END_REF][START_REF] Gonçalves | Effect of the support on the hydrodeoxygenation of m-cresol over molybdenum oxide based catalysts[END_REF].

By the HYD route, methylcyclohexane is then obtained via the alkene hydrogenation.

The other pathway leads to the direct scission of C-O bond forming toluene being called direct deoxygenation (DDO route). The proposed reaction pathways for the transformation of mcresol is shown in Scheme 1.

The addition of nickel in the Co/Al2O3 strongly enhances the hydrogenating properties of the solid, which can be seen by the ratio of methylcyclohexane to methylcyclohexenes that increases with Ni content (Table 4). This result is expected since nickel is known for its ability to hydrogenate even hard refractory aromatics such as benzene [START_REF] Mokrane | Benzene hydrogenation over aluminasupported nickel nanoparticles prepared by polyol method[END_REF]. From the aromatic selectivity order, NiCo catalyst tends to break the C-O bond leading to toluene by increasing the Co content, even if it is in a small extent. Indeed, the synergistic effect between cobalt and nickel forming an alloy increases the metal dispersion and leads to a stronger catalytic site compared to the physical mixture. Ni/(Co+Ni)

  concentrations of m-cresol in the feed and in the collected liquid sample; i C is the molar concentration of a given product.

  Co/Al2O3 > NiCo (1:3)/Al2O3 > NiCo (1:1)/Al2O3 > NiCo (3:1)/Al2O3 > Ni/Al2O3, which is the expected order since NiO has the lowest hydrogen consumption and the spinel Co3O4 the highest one. In the case of Co/Al2O3, the first reduction peak is attributed to the reduction of Co3O4 into CoO (360 ºC) and the second one to the reduction of CoO to metallic Co (574 ºC)[START_REF] Ji | Comparative study on the formation and reduction of bulk and Al2O3-supported cobalt oxides by H2-TPR technique[END_REF][START_REF] Jacobs | Fischer-Tropsch synthesis: support, loading, and promoter effects on the reducibility of cobalt catalysts[END_REF]. For Ni/Al2O3, a wide reduction peak is observed at 665 ºC for the reduction of NiO small crystallites into Ni(0). A small shoulder at 824 ºC is due to the reduction of spinel NiAl2O4[START_REF] Andonova | Structure and redox properties of Co promoted Ni/Al 2 O 3 catalysts for oxidative steam reforming of ethanol[END_REF][START_REF] Zhang | Effect of preparation method on structural characteristics and propane steam reforming performance of Ni-Al 2 O 3 catalysts[END_REF][START_REF] Zhang | Selective catalytic reduction of nitric oxide with propane over Ni-Al 2 O 3: effect of Ni loading[END_REF]. It can be seen that the addition of nickel seems to favor the reduction of Co3O4 into CoO. Indeed, the peak observed at 360 °C without Ni shifts to lower temperature in the presence of Ni, which appeared at 287 °C for NiCo(3:1)/Al2O3. For the three solids having both Ni and Co, the reduction peaks of NiO and CoO species appearing between 400 and 800 °C became much closer; a single peak at around 550 °C is noted for the sample with the Ni to Co ratio equal to 3:1. This result indicates the formation of a bimetallic alloy Ni-Co

Fig. 3

 3 Fig. 3 displays the XRD diffraction patterns of Ni, Co and CoNi catalysts after being

  order to highlight the promoting effect in HDO of the presence of Ni-Co alloy instead of neighboring Ni and Co metallic phases, the catalytic properties of NiCo(3:1)/Al2O3 were compared to the ones of a mechanical mixture composed by 75 wt.% of Ni/Al2O3 and 25 wt.% of Co/Al2O3. As highlighted in Fig. 7, the former was about 1.7 times more active than the latter, showing the advantage of the presence of Ni-Co alloys to favor the HDO catalytic properties of bimetallic Ni-Co solids, especially when NiCo(3:1) phase was formed.

Fig. 9

 9 Fig. 9 displays the influence of time on stream on both conversion and selectivity for the transformation of m-cresol (τ = 5 g h mol -1 ) over NiCo (3:1)/Al2O3 under 4 MPa and 340

4 Conclusion

 4 In this work, the hydrodeoxygenation of m-cresol was carried out at 340 °C under 4 MPa over monometallic Ni and Co nanoparticles and various bimetallic Ni-Co systems, all metallic phases were supported on Al2O3. It was demonstrated that NiCo alloys were formed in bimetallic systems. The catalyst having a Ni/Co atomic ratio close to 3 (NiCo(3:1)/Al2O3) was the most active one. A synergistic effect between Ni and Co due to the presence of NiCo alloy was highlighted. Nickel enhanced the hydrogenating properties of the catalyst whereas the addition of cobalt favored the deoxygenation rate of m-cresol due to the presence of CoNi nanoparticles which favored the formation of oxygen vacancy acting as deoxygenation active site. All catalysts remained stable during experimental conditions despite carbon deposition probably due to the high hydrogen pressure.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. N2 adsorption and desorption isotherms of the bare support and the reduced catalysts

Fig. 4 :Fig. 5 :Fig. 6 .Fig. 7 .

 4567 Fig. 4: XRD diagram of Co/Al2O3 (a) and NiCo(3:1)/Al2O3 (b) samples (experimental data

Fig. 8 .

 8 Fig. 8. Influence of the Ni/(Ni+Co) ratio of Al2O3 supported catalysts on the kTOT values (■,

Fig. 9 .

 9 Fig. 9. Conversion and selectivity as function of time on stream for NiCo (3:1)/Al2O3 at τ = 5 g h mol -1 at 4 MPa and 340 ºC for HDO of m-cresol. Conversion ( ); methylcylohexanone

Fig. 10 .

 10 Fig. 10. Conversion and selectivity as function of space time for NiCo/Al2O3 catalysts at 4

Fig. 1 . 3 Fig. 2 .Fig. 3 .Fig. 4 :Fig. 5 .Fig. 6 .Fig. 7 .(Fig. 8 .

 132345678 Fig. 1. N2 adsorption and desorption isotherms of the bare support and the reduced catalysts

Table 1 .

 1 Textural properties

		Content
					Oxide samples			Reduced samples
		(wt.%)
				SBET	Pore volume	Average pore	SBET	Pore volume	Average pore
	Solid	Ni	Co
				(m 2 g -1 )	(cm³ g -1 )	diameter (nm)	(m 2 g -1 )	(cm³ g -1 )	diameter (nm)
	Al2O3	-	-	251	0.77	12	-	-	-
	Ni/Al2O3	9.7	-	197 (225) a	0.54 (0.62) a	11	193 (220) a	0.53 (0.60) a	10
	NiCo(3:1)/Al2O3	7.4	2.7	196 (225) a	0.52 (0.60) a	10	190 (218) a	0.51 (0.58) a	10
	NiCo(1:1)/Al2O3	5.1	5.3	191 (220) a	0.51 (0.59) a	10	187 (215) a	0.51 (0.59) a	10
	NiCo(1:3)/Al2O3	7	2.7	186 (212) a	0.49 (0.56) a	10	183 (209) a	0.49 (0.56) a	10
	Co/Al2O3	-	10.1 194 (223) a	0.52 (0.60) a	11	192 (220) a	0.51 (0.58) a	11

a: in brackets are values determined after correction due to the contribution of the weight gain consecutive to the introduction of oxide phases

Table 2 .

 2 Properties of the catalysts

		CO uptake	H2 theoretical	H2-TPR	NH3-TPD	a (Å) b	D c
	Catalyst				uptake
		(µmol g -1 )	uptake a (µmol g -1 )	uptake (µmol g -1 )	(µmol g -1 )		(nm) c
	Al2O3	0	0	0	1.04	-	-
	Ni/Al2O3	35	1.7	1.8	1.06	3.5347	3.5
	NiCo(3:1)/Al2O3	46	1.8	2.0	nd	3.5374	4.7
	NiCo(1:1)/Al2O3	30	2.0	2.2	1.05	3.5412	6.4
	NiCo(1:3)/Al2O3	22	2.1	2.3	nd	3.5448	8
	Co/Al2O3	18	2.3	2.4	1.03	3.548	6.5
	a: calculated from NiO + H2 → Ni 0 + H2O and Co3O4 + 4H2 → 3Co 0 + 4H2O
	b: lattice parameter determined from XRD using Rietveld refinements.

Table 3 .

 3 Carbon content in the spent catalysts after 25 h on stream at 4 MPa and 340 ºC.

	Catalyst	C (wt %)
	Ni/Al2O3	4.7
	NiCo(3:1)/Al2O3	5.8
	NiCo(1:1)/Al2O3	4.7
	NiCo(1:3)/Al2O3	5.0
	Co/Al2O3	5.6

Table 4 .

 4 Product distribution of m-cresol at 340 °C under 4 MPa of total pressure over NiCo/Al2O3.

	Catalyst	τ	XCRE	kTOT b	TOF		Selectivity (mol %)		1-MCHe/	MCH/	HDO
		(g h mol -1 )	(%)	(mmol g -1 h -1 )	(h -1 )	Tol 3-MCHnone MCHes	MCH	MCHes c	MCHes d	(mol %) e
	Ni/Al2O3	5.0	19.2	42.6	1217	1.4	15.3	34.1	49.0	0.63	1.4	84
	NiCo(3:1)/Al2O3	5.0	22.3	50.4	1096	1.5	10.1	32.7	55.7	0.63	1.7	90
	NiCo(1:1)/Al2O3	9.0	19.9	27.8	927	1.7	9.2	34.8	54.3	0.64	1.6	91
	NiCo(1:3)/Al2O3	20.0	22.3	13.2	600	2.0	7.9	39.9	50.2	0.63	1.3	92
	Co/Al2O3	41.0	20.8	5.7	317	6.1	3.9	49.0	41.0	0.64	0.8	96
	Ni/Al2O3 +	8.0	20.7	29.0	nd	1.6	12.6	34.5	51.3	0.643	1.5	87
	Co/Al2O3 a
	a: Physical mixture of 0.75 Ni/Al2O3 and 0.25 Co/Al2O3
	b: standart deviation of rate constants close to 5%
	c 1-MCHe/MCHEs: ratio of 1-methylcyclohexene on all methylcyclohexene isomers
	d: MCH/MCHes: ratio of methylcyclohexane on all methylcyclohexene isomers
	e: HDO: percentage of deoxygenated products

Acknowledgments

This research was founded by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), the Institute of Material and Environment Chemistry of Poitiers (IC2MP) and the University of Poitiers. Vinicius O.O. Gonçalves is grateful to the Brazilian program Science without Borders and to Université de Poitiers for the Invited Professor Grant. This research was supported by the CAPES-COFECUB program Ph-S_Sv 912/18. Walter Henrique S.M. Talon thanks the Scientific Initiation Program from Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ). Jean-Marc Clacens and Frédéric Richard acknowledge financial support from the European Union (ERDF) and "Région Nouvelle Aquitaine".


Jérôme Duriez
email: jerome.duriez@inrae.fr

Cédric Galusinski

A Level Set-Discrete Element Method in YADE for numerical, micro-scale, geomechanics with refined grain shapes

Keywords: Discrete Element Method (DEM), Level Set-DEM (LS-DEM), particle's shape, Fast Marching Method (FMM)

A C++-Python package is proposed for 3D mechanical simulations of granular geomaterials, seen as a collection of particles being in contact interaction one with another while showing complex grain shapes. Following the socalled Level Set-Discrete Element Method (LS-DEM), the simulation workflow stems from a discrete field for the signed distance function to every particle, with its zero-level set corresponding to a particle's surface. A Fast Marching Method is proposed to construct such a distance field for a wide class of surfaces. In connection with dedicated contact algorithms and Paraview visualization procedures, this shape description eventually extends the YADE platform for discrete simulations. Its versatility is illustrated on superquadric particles i.e. superellipsoids. On computational aspects, memory requirements possibly exceed one megabyte (MB) per particle when using a double numeric precision, and time costs, though also significant, appear to ˚Corresponding

Introduction

Geomaterials very often show a discrete nature which controls their solidlike strains or fluid-like strain rates while being under stress, e.g. granular soils. A proper description of that mechanical behavior is of interest to countless geo-engineering problems, e.g. the safe design of large rockfill dams [START_REF] Deluzarche | Discrete numerical modelling of rockfill dams[END_REF], possibly rising in the order of hundreds of meters after piling up decimetric pieces of rock, or the forecast of snow mechanical stability and avalanches [START_REF] Hagenmuller | Microstructure-based modeling of snow mechanics: a discrete element approach[END_REF]. Unlike the equivalent continuum descriptions classically used in engineering practice, numerical modelling approaches based on the Discrete Element Method (DEM, [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] duly respect this granular nature by describing the time evolution of a discrete set of particles, the so-called Discrete Elements (DE), in mechanical interaction. While DEM approaches often serve for qualitative studies in discrete geomechanics (e.g. [START_REF] Guo | The signature of shear-induced anisotropy in granular media[END_REF][START_REF] Duriez | Revisiting the existence of an effective stress for wet granular soils with micromechanics[END_REF], they also are more and more often deployed for quantitative modelling (e.g. Aboul [START_REF] Aboul Hosn | Discrete numerical modeling of loose soil with spherical particles and interparticle rolling friction[END_REF], possibly in a multiscale framework where a DEM description of a Representative Elementary Volume eventually substitutes constitutive relations in a FEM-like model [START_REF] Miehe | Homogenization and two-scale simulations of granular materials for different microstructural constraints[END_REF][START_REF] Guo | A coupled FEM/DEM approach for hierarchical multiscale modelling of granular media[END_REF].

On that quantitative point of view, the predictive abilities of DEM may appear as variable depending on the loading conditions [START_REF] Aboul Hosn | Discrete numerical modeling of loose soil with spherical particles and interparticle rolling friction[END_REF]. As a matter of fact, they certainly often suffer from a spherical shape assumption (adopted e.g. by [START_REF] Duriez | A discrete modeling-based constitutive relation for infilled rock joints[END_REF][START_REF] Guo | The signature of shear-induced anisotropy in granular media[END_REF][START_REF] Duriez | Revisiting the existence of an effective stress for wet granular soils with micromechanics[END_REF][START_REF] Aboul Hosn | Discrete numerical modeling of loose soil with spherical particles and interparticle rolling friction[END_REF], since such spheres constitute a very strong simplification of material particles while particle's shape has a known influence on the macroscopic behavior (e.g. [START_REF] Cho | Particle shape effects on packing density, stiffness, and strength: Natural and crushed sands[END_REF]. Therefore, various DEM strategies towards a better shape description have been introduced, such as the use of rigid aggregates of spheres, so-called clumps, that should mimic real shapes (e.g. [START_REF] Garcia | A clustered overlapping sphere algorithm to represent real particles in discrete element modelling[END_REF][START_REF] Mede | A medial axis based method for irregular grain shape representation in DEM simulations[END_REF]; or the direct consideration of polyhedra (e.g. [START_REF] Eliáš | Simulation of railway ballast using crushable polyhedral particles[END_REF][START_REF] Gladkyy | DEM simulation of polyhedral particle cracking using a combined Mohr-Coulomb-Weibull failure criterion[END_REF].

Clumps offer the advantage to accommodate straightforward and computationally cheap contact algorithms designed for spheres, but still present some unrealistic local roundness. On the other hand, polyhedra resort to more complex algorithms, which still remain restricted, most often, to convex surfaces [START_REF] Dubois | Numerical modeling of granular media composed of polyhedral particles[END_REF]. One can also note the potential particles (PP) or potential blocks (PB) approaches by [START_REF] Houlsby | Potential particles: a method for modelling non-circular particles in DEM[END_REF]; [START_REF] Boon | A new algorithm for contact detection between convex polygonal and polyhedral particles in the discrete element method[END_REF][START_REF] Boon | A new contact detection algorithm for three-dimensional non-spherical particles[END_REF], which both describe particles' surfaces resorting to the zero-level of a so-called potential. Each scalar potential is given in a set of closed-form expressions with a variable number of shape parameters, leading to rounded (PP-case) or angular (PB-case) surfaces that are necessarily convex. Then, the so-called Level Set Discrete Element Method (LS-DEM) has been recently proposed by [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF], in 3D, as another DEM extension towards realistic shapes. In LS-DEM, and with a limited similarity to PP and PB approaches, every DE's surface is implicitly described as the zero-level set of the specific signed distance function to that surface. Contributing to its generality, no closed-form equation or convexity assumptions are required in LS-DEM since Level Set and Fast Marching Methods [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF][START_REF] Sethian | Level set methods and fast marching methods[END_REF] are available to construct distance fields for arbitrary, possibly concave, surfaces and a wide class of scientific applications (e.g. [START_REF] Yang | Assessing and visualizing uncertainty of 3D geological surfaces using level sets with stochastic motion[END_REF]. [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF][START_REF] Kawamoto | All you need is shape: Predicting shear banding in sand with LS-DEM[END_REF] actually illustrated the capabilities of the LS-DEM to describe real soil grains, shapes being acquired through X-ray computed tomography, as well as its promising features to reproduce observed behaviors both qualitatively and quantitatively.

The present contribution then proposes an independent and original implementation of LS-DEM into the existing YADE open-source platform ( Šmilauer et al., 2015), which is often used for geo-mechanical simulations (e.g. [START_REF] Duriez | A discrete modeling-based constitutive relation for infilled rock joints[END_REF][START_REF] Boon | A new contact detection algorithm for three-dimensional non-spherical particles[END_REF][START_REF] Duriez | Revisiting the existence of an effective stress for wet granular soils with micromechanics[END_REF][START_REF] Aboul Hosn | Discrete numerical modeling of loose soil with spherical particles and interparticle rolling friction[END_REF][START_REF] Pirnia | ICY: An interface between COMSOL multiphysics and discrete element code YADE for the modelling of porous media[END_REF]. Example usages are furthermore provided for complex, superquadric shapes, alongside discussing computational costs in comparison with the polyhedral shape description.

Section 2 first recalls Level Set and Fast Marching Methods serving to establish distance fields for arbitrary surfaces. Then, Section 3 describes how LS-DEM uses the particles' distance fields for DEM simulations of granular soils, usually following here the initial guidelines of [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF] or [START_REF] Duriez | Precision and computational costs of Level Set-Discrete Element Method (LS-DEM) with respect to DEM[END_REF]. An original LS-DEM code is proposed accordingly and summarized in Section 4. Section 5 and 6 present a direct application to non-spherical, superquadric, shapes with illustrative simulations and a computational comparison with the use of convex polyhedra.

Level Set and Fast Marching methods

Level set formalism

Level set approaches [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Sethian | Level set methods and fast marching methods[END_REF] see interfaces Sptq as the zero-level set of a function φ t p x, tq being defined from R d ˆR into R, with R d covering the whole space of a d dimensionnality.

Evolving contours (resp. surfaces) can then be described for d " 2 (resp.

d " 3). While the interfaces evolve, propagating with a normal velocity v " F p x, tq n (where n stands for the outwards normal), all level sets evolve with an extended velocity parallel to the gradient of the level set function φ t p x, tq. Since φ t p x, tq is constant along Sptq (equal to 0 @t), the nullity of the material derivative along the interface front leads to the following level set equation:

Bφ t Bt `F || ∇φ t || " 0 (1)
Eq. ( 1) conforms the formalism of Hamilton-Jacobi partial derivative equations [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF], with the Hamiltonian H t , as a function of the spatial derivative(s) of φ t , being equal to:

H t pφ t x q " F |φ t x | for d " 1 (2) H t pφ t x , φ t y q " F b φ t x 2 `φt y 2 for d " 2 (3)
H t pφ t x , φ t y , φ t z q " F b φ t x 2 `φt y 2 `φt z 2 for d " 3 (4)
where f x , f y , f z stand for the spatial derivatives of any scalar function f with respect to x, y, z.

The signed (shortest) distance to Sptq is a typical choice for the function φ t , with the convention of a negative, resp. positive, distance when being inside, resp. outside, of Sptq. Doing so, and for a constant and uniform speed F p x, tq " F , one can relate φ t to T p xq, the arrival time of S at x:

φ t p x, tq " F pT p xq ´tq (5)
Inserting Eq. ( 5) into the level set equation ( 1), one easily re-obtains the so-called Eikonal equation:

F || ∇T || " 1 (6)
With respect to the use of φ t and the level set Eq. ( 1), the consideration of T and the Eikonal Eq. ( 6) forms another description of evolving interfaces, adapted to the case of a constant and uniform sign for the normal velocity.

Doing so, the current interface Sptq is the t-level set of T and no time variable enters the partial differential equation ( 6). That stationary perspective can be finally complemented by the consideration of φp xq, the distance to Spt " 0q:

φp xq " φ t p x, 0q " F T p xq (7)
with the following form for the Eikonal equation:

|| ∇φ|| " 1 ô Hpφ x , ..q " 1 (8)
Similar to Eqs. ( 1)-( 4), Eq. ( 8) can be cast in the form of a Hamiltonian Hpφ x , ..q " 1 with:

Hpφ x q " |φ x | for d " 1 (9) Hpφ x , φ y q " b φ x 2 `φy 2 for d " 2 (10) Hpφ x , φ y , φ z q " b φ x 2 `φy 2 `φz 2 for d " 3 (11)

A Fast Marching Method for the stationary perspective

Looking for the distance field φ to a given, constant, surface S, the Eikonal equation ( 8) can be efficiently solved using a so-called Fast Marching Method (FMM, [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF][START_REF] Sethian | Level set methods and fast marching methods[END_REF]. Space being discretized on a grid, the Eikonal equation makes the φ-value at some gridpoint x i being directly dependent upon surrounding φ-values at adjacent gridpoints, as can be seen from finite difference expressions for the spatial derivatives in Eqs. ( 9) to (11). Accounting for the monotonous nature of φ, which strictly increases (in absolute value) when x goes away of S, the FMM eventually gives the full discrete field φp x i q starting from an initial set of gridpoints being along, or close to, the surface and serving as boundary conditions. In more details, the FMM recursively applies Eq. ( 8), in the form of Eqs. ( 9) or (10) or (11) depending on gradient's dimensionnality, and adopting gradient expressions decentred to low and known φ-values. Recursive applications actually go in a downwind direction away from the surface, until the whole spatial grid has been handled. The key point of the FMM is to go through the grid points in the right order, following at each step the minimal value of distance.

The FMM for instance directly applies to any surface S showing a scalar inside/outside function f p xq, being positive (resp. negative) for x located outside (resp. inside) the surface and null along the surface. In such a case, boundary conditions gridpoints are easily identified as all gridpoints being outside of the surface and having a grid neighbor inside and they can be assigned the following φ-value:

φp xq " f p xq || ∇f p xq|| for x close to S (12)
By construction, Eq. ( 12) is a first order approximation to φ, obviously obeying φ " 0 along S and also verifying the Eikonal equation ( 8) close to S, provided that ∇p1{|| ∇f ||q is finite. This constitutes the initialization of the distance function on the grid points close to the interface, before applying the recursive operations of the FMM.

Figure 1 illustrates the distance output of such a FMM procedure, herein implemented in a DistFMM C++ class presented in Section 4.1, when applied to the following "flake-like" inside/outside function:

f p xq " r ´rR `∆R sinp5θq sinp4ϕqs (13)
In Eq. ( 13), pr, θ, ϕq refer to spherical coordinates with θ P r0; πs measured from z axis and ϕ P r0; 2πs measured in ( x, y) plane.

LS-DEM formulation

For any DEM mechanical simulation to progress in time, it is first necessary to describe the shapes of the bodies, i.e. DEs, and detect their pos- 13), with pR; ∆Rq " p3; 1.5q, after executing a FMM on a 0.1-spaced isotropic grid sible contact interactions with neighbors. Then, contact-scale constitutive relationships express forces and torques so that rigid motion equations can finally be integrated. The following sections detail these three steps.

LS-DEM shape description

Following [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF][START_REF] Kawamoto | All you need is shape: Predicting shear banding in sand with LS-DEM[END_REF], a discrete signed distance field on a body-centered regular grid, possibly obtained from the previous FMM, is the first LS-DEM ingredient. That (grid ; distance field) pair is independently defined for every DE in a local coordinate system and first serves for defining the DE inertial quantities (mass m and inertia matrix I " I αβ , α, β P tx, y, zu) summing contributions from grid voxels v making up the body's volume V as per the following discrete-form equations:

m " ρ ÿ vPV V v " ρ N vox V v (14)
x "

1 N vox ÿ vPV x v (15)
I xx " ρ ÿ vPV " py v ´yq 2 `pz v ´zq 2 ‰ V v ( 16
)
I yy " ρ ÿ vPV " px v ´xq 2 `pz v ´zq 2 ‰ V v ( 17
)
I zz " ρ ÿ vPV " px v ´xq 2 `py v ´yq 2 ‰ V v ( 18
)
I xy " ´ρ ÿ vPV px v ´xq ˆpy v ´yqV v ( 19
)
I xz " ´ρ ÿ vPV px v ´xq ˆpz v ´zqV v (20)
I yz " ´ρ ÿ vPV py v ´yq ˆpz v ´zqV v (21)
In the above equations, ρ is the material mass density, x v " px v , y v , z v q the middle point of a voxel and x " px, y, zq the body's center of mass. Eqs.

(15),( 19)-( 21) serve for verification purposes since the body-attached local frame is expected to be inertial and x, I xy , I xz , I yz to be nil. By a simple convention, a grid voxel v of volume V v is herein said to be part of the body's volume V when its lowest corner is inside the surface, showing a zero or negative distance value. While smoother choices have been proposed by [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF][START_REF] Kawamoto | All you need is shape: Predicting shear banding in sand with LS-DEM[END_REF], it will be verified in Section 5.2 that the present choice does not inhibit precision for grids being fine enough, i.e. showing a spacing g grid at least ten times smaller than a grain's characteristic size l grain .

For the purpose of LS-DEM contact algorithms that will be described in Section 3.2 below, a second LS-DEM ingredient adds to the distance field, in the form of a set of N n boundary nodes tN i , i P r0; N n ´1su discretizing each body's surface S. Generally speaking, boundary nodes should count in the order of thousands and their positions are defined at the intersection of S i.e.

φp xq " 0 and N n half-lines i.e. rays λ v, with v a direction and λ a positive abscissa, that stem from the center of mass. Due to the adopted tri-linear interpolation of the discrete distance field within the grid extents, φp x " λ vq is a cubic polynomial in λ whose coefficients depend upon grid distance values and ray tracing boundary nodes corresponds to solve its positive roots (see e.g. [START_REF] Lin | An efficient volume-rendering algorithm with an analytic approach[END_REF]. Since rays should provide an appropriate discretization of spherical angles pθ, ϕq, the corresponding directions pθ, ϕq are chosen to follow a spiral path instead of a simple rectangular discretization of r0; πs ˆr0; 2πs in order to avoid a possible (shape-dependent) concentration of nodes at the poles θ " 0rπs. More details about the spiral path or the choice of boundary nodes number are given by [START_REF] Duriez | Precision and computational costs of Level Set-Discrete Element Method (LS-DEM) with respect to DEM[END_REF] and in the next sections.

Contact law

From the distance fields and the set of boundary nodes, contact detection between two bodies 1 and 2 relies on a master-slave algorithm whereby nodes N 1 i of body 1 are tested in the distance field φ 2 of body 2 (see also Figure 9). In order to increase precision, the body 1 is chosen as the smallest one in volume, which enables one to explore distance fields with the greatest surface density in nodes. Contact is detected as soon as one node N 1 i verifies φ 2 pN 1 i q ď 0. LS-DEM belongs to the wide class of "soft" DEM whereby small overlaps, φ 2 pN 1 i q ă 0, are possible: these overlaps would in reality materialize through slight changes in shape which are neglected in soft DEM approaches. After identifying the set of contacting boundary nodes, a unique contact point is herein chosen from the node N c showing the greatest interpenetration depth u n , which also gives the contact normal as the local gradient of φ 1 :

u n " ´minpφ 2 p Ý Ý Ñ ON i q, Ý Ý Ñ ON i P S 1 q " ´φ2 p ÝÝÑ ON c q ě 0 (22) n " ∇φ 1 p ÝÝÑ ON c q (23)
That final consideration of a unique contacting point, also adopted by [START_REF] Li | Capturing the inter-particle force distribution in granular material using LS-DEM[END_REF], currently restricts the proposed LS-DEM implementation to convex shapes. For a pair of contacting bodies with concave shapes, multiple contact points would occur but these could be easily detected with the same master-slave algorithm. As such, [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF][START_REF] Kawamoto | All you need is shape: Predicting shear banding in sand with LS-DEM[END_REF]) also addressed concave shapes by defining a mechanical interaction at each contacting boundary node. While being more general, this choice nevertheless poses the risk to make the macroscopic behavior, e.g. the bulk stiffness, to directly depend upon the chosen number of boundary nodes in case a physically unique contact area would involve more than one boundary node.

A classical contact law for cohesionless materials finally expresses the interaction force after decomposing the latter in a normal, along n, and a tangential component. A repulsive normal force F n first arises due to the interpenetration depth u n , as per a linear elastic model with a k n stiffness:

F n " k n u n n (24)
Along the tangential direction, a linear elastic-plastic relationship governs the shear force variations. Denoting k t the shear stiffness and µ the friction coefficient, the following Eqs. ( 25)-( 26) describe the variations of the shear force F t , starting from 0:

d F t " d ˜|| F t || F t || F t || ¸" || F t ||d ˜ F t || F t || ¸`dp|| F t ||q F t || F t || (25) dp|| F t ||q F t || F t || " k t d u t enforcing || F t || ď µ|| F n || (26)
while updates in the shear force direction, F t {|| F t ||, are applied in order to follow changes in the tangent plane's orientation, e.g. a change in contact normal ( Šmilauer et al., 2015).

Motion integration

As for general DEM, the translation and rotation of each DE in space, under resultant force F and torque Γ (computed at the center of mass), finally follow Newton-Euler equations for rigid bodies with v and ω the linear and angular velocities:

m d v dt " p1 ˘Dq F (27) I d ω dt ` ω ^I ω " p1 ˘Dq Γ (28)
The above Newton-Euler equations are classically damped using a numerical coefficient D, which modifies the resultant force and torque so that kinetic energy always decreases (or is led to increase by a smaller extent) as soon as F ‰ 0. Eq. ( 28), relating the variation in ω with Γ, is expressed in local axes where the inertia tensor I is constant. Denoting Rptq the rotation matrix passing from local axes to global ones, and Ω the antisymmetric matrix such that Ω x " ω ^ x, @ x, Eq. ( 28) is finally supplemented with:

dR dt " R Ω (29)
These equations ( 27)-( 29) are then integrated over the time steps through an explicit algorithm common to any non-spherical shape in YADE ( Šmilauer et al., 2015).

Proposed implementation

Source code

The present C++ and Python implementation inserts LS-DEM into the 2020.01a version, i.e. the git commit 9964f53, of the YADE platform ( Šmilauer et al., 2015). Figure 3 Boost.Math library is adopted for this purpose, being preferred over canonical formulae for numerical stability. Moreover, the distance cubic polynomial Set of Discrete Elements is first turned dimensionless (with respect to the grid spacing) for the relative magnitude of its coefficients to be unaffected by the unit system, insuring a constant behavior of the root finding algorithm whatever the user's choice in this aspect. It is recalled the distance field, its grid and the boundary nodes all refer to a reference local (inertial) frame for each particle.

• LevelSet • DistFMM • RegularGrid • VTKRecorder Contact treatment • Bo1 LevelSet Aabb • Ig2 LevelSet LevelSet ScGeom

Explicit motion integration

The distance field at the heart of the shape descriptor can be directly In the end, the set of kinematic variables for a LS-DEM interaction is equivalent in nature to those used for spheres in general DEM and it can be stored in the existing ScGeom class. Constitutive properties k n , k t and µ also correspond to the pre-existing FrictPhys contact model in YADE. This enables LS-DEM simulations to adopt a pre-existing contact law, namely

Law2 ScGeom FrictPhys CundallStrack. Motion integration as described in previous Section 3.3 has also been readily available from the NewtonIntegrator class.

Code usage

The (modified or classical) YADE platform starts in the form of a Python3 interactive interface, invoked from install/bin/yadelevelSet in the proposed installation procedure (see the "Computer code availability" section).

Instead of an interactive session, scripts prepared beforehand can be as The latter enables users to directly insert any distance field they would have otherwise acquired, for instance from computed tomography [START_REF] Vlahinić | Towards a more accurate characterization of granular media: extracting quantitative descriptors from tomographic images[END_REF]. In all cases, grid spacing g grid is input through a spacing attribute while a nNodes attribute of levelSetBody() controls the boundary nodes number N n .

Documentation can be obtained for any class or attribute in the usual interactive Python manner, typing e.g. LevelSet? or levelSetBody?. An HTML version of the documentation can also be built executing make doc from the compilation folder.

Code validation

The implementation is first validated for what concerns the FMM in DistFMM class. Applying the procedure on a sphere of radius R with a known distance field φ th p xq " r ´R, numerical precision can be quantified, looking e.g. at the average relative error on all gridpoints (excluding those with φ th " 0) or at the relative error at the center, as follows:

err avg " average ˆ"ˇˇˇˇφ p x i q ´φth p x i q φ th p x i q ˇˇˇ, x i | φ th p x i q ‰ 0 *˙( 30) err ctr " minpφq `R R (31)
Considering Eq. ( 30), another average error is also analyzed for the more complex flake-like shape previously presented in Figure 1. While no exact distance field is known for such a surface, one could attempt a reconstruction of its inside/outside function f , Eq. ( 13), solving another variant of the Eikonal equation, with a non-unit speed, i.e.:

|| ∇φ|| " || ∇f || (32)
By initializing the FMM, close to S, with values of f : φp x i q " f p x i q and solving for Eq. ( 32), one should indeed await φ th " f as an exact solution.

For these two examples of a FMM application, Figure 4 illustrates how the FMM results approach their respective φ th with a decreasing grid spacing g grid i.e. an increasing grid resolution r g " 2R{g grid . While the precision is somewhat worse for the flake-like surface, in line with an increasing complexity of the problem, it always linearly scales with the grid resolution, in accordance with the first order expression of ∇φ in the numerical method. GHz. Denoting N gp the total number of gridpoints, with N gp " Opr g 3 q, a Opr g 6 q " OpN gp 2 q complexity appears, in accordance with classical Level Set Methods. [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF] actually proposed a lighter complexity for the FMM, through adopting a heap sort when searching the minimum φ-value for propagating the distance field. For the purpose of LS-DEM, the FMM will apply only once per DE, at the very beginning of a simulation and the present time cost in the order of a second for few tens of grid voxel per particle length is actually acceptable, considering the final time cost of a complete LS-DEM simulation. Figure 5 finally illustrates that the FMM computation of distance for the more complex, non-spherical, flake-like shape logically shows the same time costs. nodes are appropriately chosen [START_REF] Duriez | Level set representation on octree for granular material with arbitrary grain shape[END_REF][START_REF] Duriez | Precision and computational costs of Level Set-Discrete Element Method (LS-DEM) with respect to DEM[END_REF].

A direct application of LS-DEM to superquadric shapes

The versatility of LS-DEM to address complex shapes is now illustrated on superellipsoids, also known as superquadric ellipsoids. These surfaces are first presented from an analytical point of view before that their LS-DEM description is introduced with its corresponding precision and eventually compared with the possible use of convex polyhedra.

Superellipsoids surfaces

Superellipsoids [START_REF] Barr | Superquadrics and angle-preserving transformations[END_REF][START_REF] Barr | Rigid physically based superquadrics[END_REF] form a versatile class of surfaces which can be used as more complex shape models of granular soils (see e.g. [START_REF] Wang | Superellipsoid-based study on reproducing 3D particle geometry from 2D projections[END_REF]. They generalize ellipsoids through two additional exponents ǫ e and ǫ n that enter their surface equation together with three different radii r x , r y , r z . In a local frame, the surface equation namely reads:

f px, y, zq " ˜ˇˇˇx r x ˇˇˇ2 ǫe `ˇˇˇy r y ˇˇˇ2 ǫe ¸ǫe ǫn `ˇˇˇz r z ˇˇˇ2 ǫn ´1 " 0 (33)
Figure 6 illustrates five different superellipsoids, with their corresponding shape parameters presented in Table 1. Table 2 also details their volume and inertia properties, as obtained from closed form expressions given by [START_REF] Barr | Rigid physically based superquadrics[END_REF]. One can here observe how the ǫ n exponent modifies the z-variation of cross-sections in px, yq planes. For instance, adopting ǫ n Ñ 0 induces fairly constant cross-sections and a wider distribution of matter for extreme values along the "north-south" axis z, see Shapes A or C. On the other hand, the ǫ n " 1 case corresponds to a rounded variation of these cross sections when progressing along z (Shape B). Some singularity, i.e. a sharpness at z " 0, would appear for ǫ n ě 2, alongside concavity in a plane tangent to z for ǫ n ą 2. For a given ǫ n , ǫ e controls the contour's roundness in the px, yq plane of these cross-sections. While ǫ e " 1 corresponds to perfectly round (circles or ellipses) contours, decreasing ǫ e towards 0 induce edges that tend to align with the x and y axes, see Shape A vs C. Alternate edges and sharpnesses would be obtained in the px, yq plane at ǫ e " 2, just before concavity in that plane, for ǫ e ą 2.

LS-DEM description of superellipsoids

Previous DEM descriptions of superellipsoids have already been proposed

by [START_REF] Podlozhnyuk | Efficient implementation of superquadric particles in Discrete Element Method within an open-source framework[END_REF] or [START_REF] Weinhart | Fast, flexible particle simulations -an introduction to MercuryDPM[END_REF], for instance. In those [START_REF] Barr | Rigid physically based superquadrics[END_REF] studies, contact detection involves a minimization procedure that endows the shape equation ( 33) with an approximated distance nature, following the potential approach by [START_REF] Houlsby | Potential particles: a method for modelling non-circular particles in DEM[END_REF]. Such a minimization is then performed by an iterative numerical method, at each DEM iteration. On the other hand, the generic workflow of LS-DEM is herein proposed to directly apply to superquadrics, considering true distance quantities and avoiding the need for an iterative procedure, outside the consideration of boundary nodes.

The LS-DEM description of a superellipsoid particle nevertheless logically shows a finite precision, with for instance the inertial quantities depending on the chosen resolution for the grid carrying φ, as per the above Section 3.1.

Quantifying now the grid resolution as r g " 2 minpr x , r y , r z q{g grid , Figure 7 then compares the obtained LS-DEM volume with the expected volume presented in Table 2. It shows that using at least ten grid cells per particle's length leads to satisfactory results with an error on the volume being smaller than few %. A similar precision is achieved for inertia components, as shown in Figure 15 in the Appendix. While this analysis is merely geometric in nature, a direct connection between errors in describing particles' volumes and bias in mechanical results was proposed by [START_REF] Mede | A medial axis based method for irregular grain shape representation in DEM simulations[END_REF] when using clumps.

The influence of grid spacing on inertial quantities directly relates to the voxellised nature of the present description of particle's volume, in connection with the sign of discrete φ-values φp x i q. Section 4.3 previously illustrated how the grid spacing also affects the precision in the actual values of those, after solving through a FMM the Eikonal equation. A last impact of grid spacing onto the LS-DEM precision exists through the tri-linear interpolation used to evaluate distance at any location other than a gridpoint, such as a boundary node for the purpose of contact detection. From the present and past results (Duriez and Galusinski, 2020; Duriez and Bonelli, 2021), using r g in the order of few tens (10 to 50) appears to be an adequate compromise between precision and computational (memory) costs, on all aspects.

Time costs in comparison with convex polyhedra

LS-DEM time costs are now briefly illustrated in comparison with the use of convex polyhedra as initially implemented in the YADE platform by [START_REF] Eliáš | Simulation of railway ballast using crushable polyhedral particles[END_REF]. Describing such Polyhedra shapes in YADE relies on the CGAL library, used here in its 4.11 version [START_REF] Kettner | 3D polyhedral surface[END_REF]. That external library determines for instance a possible overlapping volume between two convex polyhedra for the purpose of contact treatment.

Such shapes may actually also apply to the present five superellipsoids, after locating the polyhedra's vertices along the superquadric surface. Previously determined LS-DEM boundary nodes (with r g " 50) can be used for such a purpose. These vertices, through their connecting edges and plane portions (facets) making the polyhedra's surface, govern the precision in describing a superellipsoid shape even though, by the present construction, the obtained particles volumes are always smaller than the exact volumes of the considered (convex) superellipsoids. Figure 8 illustrates how the number of vertices controls the obtained volume and the necessity to use hundreds of polyhedra vertices in order to limit the error on the volume below few %. to a greater precision. As for the use of convex polyhedra, the corresponding time cost appears as proportional to N v 1.7 , then close to OpN v 2 q. With N v being checked to be itself proportional to the number of edges, N e , or planar facets, N f , making up each polyhedral surface, this OpN v 2 q " OpN e 2 q time complexity is actually consistent with the consideration of all possible edge pairs adopted by [START_REF] Eliáš | Simulation of railway ballast using crushable polyhedral particles[END_REF] for that contact algorithm. Mostly, the polyhedral time cost is several orders of magnitude higher than its LS-DEM counterpart for N n " N v . In spite of the incomplete equivalence between N n and N v , these important differences in time costs clearly suggest LS-DEM might be lighter to use in terms of time, especially if a high fidelity is desired at the particle scale since this would here require hundreds of polyhedra vertices (Figure 8).

Results

After executing 56 000 DEM iterations over 1.4 s of model time, a final equilibrium state can be observed in Figure 12, for what concerns the average coordination number z c or the vertical load exerted on the ground wall F , compared in a F rel ratio with the expected weight F th that corresponds to the theoretical solid volumes of all particles (Table 2):

F rel " F ρ|| g|| n DE ř i"1 V th piq (34)
In this illustrative simulation, most dissipation of the initial gravitational energy is artificial, coming from the numerical damping mentioned in the above Section 3.3 used with D " 0.3. here appears as optimal since finer particle descriptions eventually lead to the same results though with higher computational costs, as discussed in the following. 4. It is to note though that those memory requirements could be reduced in the future, adopting octree structures to carry the distance field instead of regular grids [START_REF] Duriez | Level set representation on octree for granular material with arbitrary grain shape[END_REF]. Table 4: Total (whole simulation) RAM usage for the discharge simulations of Figure 12 lighter LS-DEM costs with respect to polyhedra, even though those costs would logically be even more reduced with ideal spherical shapes [START_REF] Duriez | Precision and computational costs of Level Set-Discrete Element Method (LS-DEM) with respect to DEM[END_REF]. Time costs can anyway be significantly decreased using simple OpenMP parallel computing in a shared memory paradigm. Doing so, loops over interactions (for contact treatment in InteractionLoop) or bodies (for motion integration in NewtonIntegrator) are split into different OpenMP threads which are simultaneously executed by different CPU cores. With respect to the sequential case, additional supervisory operations become necessary in order to avoid simultaneous access to the same variable in memory from different threads. Nevertheless, OpenMP execution of the present discharge simulation, using the optimal choices r g " 20 and

N n " 2000, appears as very beneficial, with a significant, linear and nearly optimal, speedup as depicted in Figure 14. Speedup is here measured repeating 3 times each parallel execution as well as the sequential one, on a server machine with two Intel Xeon Platinum 8270, 2.7 GHz, processors with 26 cores and 36 MB of cache memory each, i.e. a total of 52 cores and 104 threads, together with 1.5 TB 2.9 GHz RAM. From all these simulations, 9

parallel / sequential timing ratios are computed and depicted in Figure 14 through their average and standard deviation. cmake,g++,boost,Qt,freeglut3,libQGLViewer,eigen,gdb,sqlite3,Loki,

  Figure1: Level Set description of a flake-like surface defined by Eq. (13), with pR; ∆Rq "

Figure 2 :

 2 Figure 2: Shape description in LS-DEM: a regular grid t x G u carrying the distance field φ, together with boundary nodes tN i u (2D view for clarity)

  illustrates the LS-DEM workflow exposed in the previous section together with the most noticeable new (or modified) C++ classes responsible for execution. Looking from the lsYade root folder of the proposed source code, the files pkg/dem/LevelSet.*pp introduce the new shape descriptor LevelSet. That class includes the discrete distance field as a LevelSet.distField attribute. The regular grid carrying the distance field is LevelSet.lsGrid, which is an instance of the RegularGrid class. Boundary nodes are stored in Lev-elSet.boundNodes and computed (once, at the beginning of a simulation) solving for cubic roots during the ray tracing procedure mentioned in the above Section 3.1. A Newton-Raphson algorithm proposed by the external

Figure 3 :

 3 Figure 3: New or modified C++ classes for a LS-DEM workflow in YADE

  passed from the user (see next section) or also obtained from the Fast Marching Method proposed in DistFMM class. In the latter case, DistFMM.phi() is to execute once the grid and the boundary conditions are defined as grid and phiIni class attributes, respectively. Various predefined functions are proposed to build appropriate boundary conditions expected in phiIni. In addition to a distIniSE function intended to compute the distance to superquadric shapes detailed in Section 5.1, a versatile PhiIniPy function may be based upon any user-defined Python function that discriminates between the inside and the outside of a surface and outputs boundary condition values for the FMM, such as shown in Figure 1. Visualization of LevelSet-shaped bodies relies on vtk exports of the discrete distance field for each DE in current configuration, thanks to a modified version of files pkg/dem/VTKRecorder.*pp. Actual display is typically done from Paraview software (Ayachit, 2019), using its Python interface and a provided pvVisu function defined in examples/levelSet/pvVisu.py. For the purpose of alternate vizualisation methods at the user's discretion, a Lev-elSet.marchingCubes method is also available from YADE interface and gives a triangulated description of a particle's surface as per the Marching Cubes algorithm (Lorensen and Cline, 1987). The files pkg/dem/LevelSetInteraction.*pp finally implement the contact algorithms described in previous Section 3.2. The Bo1 LevelSet Aabb is first responsible to compute an axis-aligned bounding box (Aabb) used in YADE for a first, crude and fast, detection of possible contacts. At the beginning of a simulation, that class first loops over the whole distance field to compute the 8 corners of the corresponding Aabb in local axes, stored in LevelSet.corners. Then, the current Aabb in model (global) axes is easily determined following rigid transformations. In case of overlap between Aabb, precise contact detection subsequently resorts to the other class Ig2 -LevelSet LevelSet ScGeom, that implements the master-slave contact detection based on boundary nodes, identifying the contact point, if any, and the associated kinematic variables (normal vector, interpenetration depth) between two LevelSet-shaped bodies. Similar classes enable contact interaction between a LevelSet-shaped body and existing Wall or Box shapes, often adopted in YADE to simulate rigid boundaries.

  well passed as argument and launched in the same manner than classical Python scripts. Examples of YADE scripts using the new LS-DEM features can be found in the source code at lsYade/examples/levelSet/*.py (levelSetBody.py in particular) and also lsYade/scripts/checks-andtests/checks/checkLSdem.py. The latter actually serves as a new regression test into the YADE platform (Haustein et al., 2017), to insure stability of the LS-DEM features in the future. These examples illustrate the definition of LevelSet bodies through a new levelSetBody() YADE function. That function proposes level set descriptions of pre-defined analytical shapes (from boxes and spheres to superellipsoids, see next Section 5), together with the possibility of a direct assignment of the regular grid with its distance field.

  Figure 4: Influence of the grid resolution on the FMM precision, with reference to spherical or flake-like (Figure 1) surfaces

Figure 5 :

 5 Figure 5: Time cost of the FMM according to the number of gridpoints per space axis 3 a N gp , with N gp the total number

Figure 6 :

 6 Figure 6: Five possible superquadric shapes

Figure 7 :

 7 Figure 7: LS-DEM precision in describing superellipsoid volumes, with reference to Figure 6

Figure 8 :

 8 Figure 8: Precision in describing superellipsoid volumes when using convex polyhedra, with reference to Figure 6

Figure 9 :Figure 10 :

 910 Figure 9: Two contacting superellipsoids described using LS-DEM (left, with 51 boundary nodes) or convex polyhedra (right, with 107 vertices per body)

  of LS-DEM is proposed in examples/levelSet/discharge.py as the discharge under gravity ( g " ´g z, with g " 9.8m/s 2 ) and into a rigid container (L x ˆLy ˆLz " 0.25 2 ˆ8 m 3 ) of n DE " 1000 superellipsoids with equal proportions of the previous five shapes A to E (Figure11). Similar dynamic simulations could serve to study rock falls and slides up to an obstacle, or the angle of repose of granular geomaterials conveyed in industrial processes.

Figure 11 :

 11 Figure 11: Views of the initial cloud of superellipsoids (left: in whole, right: close-up), with lateral and ground walls of the container not shown

Figure 12

 12 Figure 12 also illustrates the possible influence of LS-DEM discretization parameters N n and r g . Using just N n " 50 boundary nodes, together with r g " 20, for instance prevents stabilization because contacts are hardly detected and too easily lost. On the other hand, choosing pr g " 20; N n " 2000q

Figure 12 :

 12 Figure 12: Dynamics of the discharge illustration (with only a fraction of datapoints for the r g " 20; N n " 50 case on (b), for readability)

Figure 13 :

 13 Figure 13: LS-DEM memory requirements per discrete element definition. Each data point is obtained from a different discharge simulation. The planar fit is colored according to memory (also on the z-axis) and is obtained after bilinear regression, following the expression a ˆNgp `b ˆNn with a " 0.0841 ˆ10 ´4 MB and b " 0.2637 ˆ10 ´4 MB

Figure 14 :

 14 Figure 14: OpenMP scalability of the LS-DEM discharge simulation for r g " 20 and N n " 2000. Sequential time cost is 28696 s ˘106 s (« 8 h) from average and standard deviation on 3 runs, being reduced to 392 s ˘5 s (« 6.5 min) using 100 CPU cores

Figure 15 :

 15 Figure 15: LS-DEM precision in describing inertia components for the five superellipsoids of Figure 6

Table 1 :

 1 Shape parameters of the five superellipsoids shown in Figure6

	Shape Half-extents (cm) Curvature exponents
		r x	r y	r z	ǫ e	ǫ n
	A	0.58 1	0.83 0.1	0.5
	B	0.42 1	0.83 0.1	1
	C	same as Shape B	1	0.5
	D	0.5 0.7	1	1.4	1.2
	E	0.4	1	0.8	0.4	1.6
	Shape Volume (cm 3 ) Inertia components (cm 5 )
			V th	I th xx {ρ	I th yy {ρ	I th zz {ρ
	A		3.353	1.649 0.9751	1.358
	B		1.852	0.7456 0.3417 0.5770
	C		1.914	0.7996 0.4389 0.5153
	D		1.093	0.2773 0.2350 0.1304
	E		1.086	0.1283 0.3184 0.2625
	Table 2: Geometric properties of the considered superellipsoid shapes. Inertia
	components are obtained following

Table 3

 3

	lists the simulation's parameters, with

Table 3 :

 3 LS-DEM parameters for the discharge simulation

Acknowledgements

The authors acknowledge financial support of the French Sud region to the LS-ENROC project; Stéphane Bonelli (INRAE, Aix Marseille Univ, RE-

method, be it in terms of memory or execution time. Time costs are nevertheless beneficial with respect to a polyhedral description of complex shapes, as already available in YADE, and they can be furthermore reduced through OpenMP parallel computing with a significant speed-up. As for the memory requirements, these could also decrease in the future using a more appropriate data structure than the current regular grid [START_REF] Duriez | Level set representation on octree for granular material with arbitrary grain shape[END_REF].

Perspectives lie in user-friendly LS-DEM simulations in YADE for multiscale investigations in granular mechanics. A particular multiscale avenue is formed by the hierarchical modelling approaches where the DEM serves as an alternative to phenomenological (e.g. elasto-pastic) stress-strain constitutive relations in structure-scale FEM simulations (e.g. [START_REF] Guo | A coupled FEM/DEM approach for hierarchical multiscale modelling of granular media[END_REF].

Appendix

Confirming the analysis made on volumes in Section 5.2 (Figure 7), Figure 15 illustrates how LS-DEM achieves to describe inertia components of superellipsoids with a very good precision, provided the grid resolution is fine enough i.e. includes more than 10 grid voxels per particle length.

Conflict of Interest

We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.