R Troncy

Gilles Bonnet
email: bgbonnet@univ-lr.fr

F Pedraza
email: cfpedraza@univ-lr.fr

Microstructural characterization of NiAl-Al2O3 composite materials obtained by in situ aluminothermic reduction of NiO for potential coating applications

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Nowadays, the demand for materials with high performances at elevated temperature is in constant increase. Nickel based superalloys possess excellent mechanical properties at high temperature due to their high melting point. However, their resistance to oxidation by hot gases and to so-called hot corrosion appears limited [START_REF] Goward | Current research on the surface protection of superalloys for gas turbine engines[END_REF][START_REF] Stringer | High-temperature corrosion of superalloys[END_REF]. To overcome these disadvantages, the use of intermetallic compounds was developed. NiAl is the most attractive due to capacity to develop a protective aluminum oxide layer through oxidation [START_REF] Brumm | The oxidation behaviour of NiAl -I. Phase transformations in the alumina scale during oxidation of NiAl and NiAl-Cr alloys[END_REF][START_REF] Choi | Transmission electron microscopy observations on the phase composition and microstructure of the oxidation scale grown on as-polished and yttriumimplanted β-NiAl[END_REF]. The coating acts as an aluminum reservoir whereby aluminum is consumed by oxidation at the surface of the coating to form an alumina protective layer [START_REF] Goward | Current research on the surface protection of superalloys for gas turbine engines[END_REF][START_REF] Stringer | High-temperature corrosion of superalloys[END_REF][START_REF] Nicholls | Designing oxidation-resistant coatings[END_REF][START_REF] Goward | Mechanisms of formation of diffusion aluminide coatings on nickelbase superalloys[END_REF]. Further, Al diffuses into the substrate, which lowers the overall Al content in the coating. Below a critical concentration, other non-protective oxides such as NiAl2O4 and NiO may form resulting in the rapid degradation of the material. For instance, during the oxidation of Ni(Pt)Al systems at high temperature (100 h at 1050°C), 2 to 8% of Al are consumed by the formation of the aluminum oxide while 20 to 35% of Al are "lost" by diffusion to the superalloy substrate [START_REF] Vialas | Effect of Pt and Al content on the long-term high temperature oxidation behavior and interdiffusion of Pt-modified aluminide coating deposited of Ni-base superalloys[END_REF]. Diffusion barriers can be employed to prevent diffusion of aluminum in the substrate. These diffusion barriers can be made of refractory elements or rare-earth ones whose effectiveness to slow down the diffusion of aluminum has been demonstrated [START_REF] Matsumura | Oxidation behavior of a Re-base diffusion-barrier/NiAl coating on Nb-5Mo-15W at high temperatures[END_REF][START_REF] Wang | Diffusion barrier behaviors of (Ru, Ni) Al/NiAl coatings on Ni-based superalloy substrate[END_REF][START_REF] Narita | Diffusion barrier coating system concept for high temperature applications[END_REF]. However, such diffusion barriers mechanically weaken the coated components and increase the coating prices.

Our approach to maintain a sufficient amount of Al is to introduce aluminum micro-containers into the coating so that Al is released once a concentration gradient is established between the reservoir and the Al-impoverished coating matrix. The micro-containers are small volumes (tens of µm) of an Al-rich NixAly intermetallics compound entrapped by aluminum oxide. The micro-containers shall thus have an Al-rich intermetallic core (Al3Ni or Al3Ni2) encapsulated in an alumina shell. The intermetallic core can be synthesized by Self Propagating High-Temperature Synthesis (SHS). SHS is a complex thermally activated 2 physicochemical process that takes place at elevated temperatures [START_REF] German | Sintering-Theory and Practice[END_REF]. A great number of factors and process conditions can exert a significant influence on the results, in particular temperature, time, size of the particles, heating and cooling rates, pressure and atmosphere [START_REF] Nieroda | Development of the method for the preparation of Mg2Si by SPS technique[END_REF][START_REF] Chmielewski | Processing, microstructure and mechanical properties of Al2O3-Cr nanocomposites[END_REF]. The relationship between these different parameters and the reaction process between the three elements (Ni, Al and O) was reported by Levin & al. and results in multiphase (metal, liquid and oxide) components as shown in Fig. 1 [START_REF] Levin | Phase diagrams for ceramists[END_REF]. [START_REF] Levin | Phase diagrams for ceramists[END_REF].

The alumina shell can be produced in situ by an aluminothermic process following the thermite reaction (eq. 1) [START_REF] When | Characterization of thermochemical properties of Al nanoparticles and NiO nanowire composites[END_REF]: The redox reaction between NiO and Al results in the release of a large amount of energy. The products of such a reaction are solid, liquid and sometimes gaseous, inducing an intense sputtering and destruction of the samples [START_REF] Vrel | Reactions kinetics and phase segregation in the 3NiO + 2Al Al2O3 Thermite system[END_REF][START_REF] Vrel | SHS reactions in the NiO-Al system: Influence of stoichiometry[END_REF]. However, the exothermic effect can be inhibited by the addition of an inert diluent, i.e. of a heat absorber. Vrel et al. studied the effects of the addition of alumina and of the limitation of one of the two reagents on the aluminothermic reaction [START_REF] Vrel | Reactions kinetics and phase segregation in the 3NiO + 2Al Al2O3 Thermite system[END_REF][START_REF] Vrel | SHS reactions in the NiO-Al system: Influence of stoichiometry[END_REF]. In addition, Biswas et al. [START_REF] Biswas | A study of self-propagating hightemperature synthesis of NiAl in thermal explosion mode[END_REF] reported that the heating rate and the particle size could be modified to limit the amount of energy released and produce different intermetallic compounds. For instance, low heating rates and a large size of particles can cancel the explosion of the powders.

Based on the above, the Al-Ni SHS and the aluminothermic reactions will be combined in this work to synthesize NiXAly-Al2O3 phases that could serve as self-healing coatings. Different parameters including particle size, heating rates, amount of nickel oxide and contact area between particles will be studied.

Experimental procedure

Commercially available Ni microparticles underwent a preoxidation treatment at 700°C in air for different times to form NiO oxide shells over the Ni powders. This oxidation temperature of 700°C was chosen following the works of Peraldi et al. [START_REF] Peraldi | Correlations between growth kinetics and microstructure for scales formed by high-temperature oxidation of pure nickel, I. Morphologies and Microstructures[END_REF][START_REF] Peraldi | Correlations between growth kinetics and microstructure for scales formed by high-temperature oxidation of pure nickel, II. Growth Kinetics[END_REF], who described the microstructure of the oxide scales grown on pure Ni as function of temperature and time. The porous NiO formed at 700°C for short times is expected to facilitate the Al flow during the SHS reactions. The particles size of the powders and the time of preoxidation are given in Table 1. [START_REF] Vrel | Reactions kinetics and phase segregation in the 3NiO + 2Al Al2O3 Thermite system[END_REF][START_REF] Vrel | SHS reactions in the NiO-Al system: Influence of stoichiometry[END_REF]. Moreover, it allows to keep a metal core surrounded by a porous NiO to form the desired composite compounds [START_REF] Biswas | A study of self-propagating hightemperature synthesis of NiAl in thermal explosion mode[END_REF][START_REF] Peraldi | Correlations between growth kinetics and microstructure for scales formed by high-temperature oxidation of pure nickel, I. Morphologies and Microstructures[END_REF][START_REF] Peraldi | Correlations between growth kinetics and microstructure for scales formed by high-temperature oxidation of pure nickel, II. Growth Kinetics[END_REF].

A powder mixture containing the preoxidized nickel particles and aluminum particles was blended in the stoichiometric atomic ratio to form NiAl according to the Al-Ni binary diagram (Fig. 2.a). In addition, the preoxidized nickel particles were deposited over the Al microparticles (Fig. 2.b) to simulate the interfacial reactions of a potential coating. Both configurations (mixed and stacked powders) were placed in alumina crucibles for the DSC experiments (SETARAM Labsys Evo 1600 thermal analysis). The heating rates were fixed at 2 and 25°C/min with a maximum temperature of 1300°C; the cooling rate was fixed at 50°C/min (to keep the microstructures and crystal phases obtained at high temperature). Upon annealing, the reaction chamber was ventilated with a flow of Ar (20 mL/min) to limit oxidation phenomena and prevent secondary reactions in the system, although some air can be trapped between the powders when poured into the crucible. For each test, energy changes and reaction temperatures were determined using the Calisto software coupled to the DSC technique. To identify all the transformations, additional heat treatments were performed in the same conditions but the tests were interrupted at 750 and 1000°C. After the heat treatment, the crucibles were prepared for microscopy observation of the cross-sections. They were thus mounted in a phenolic resin, ground till SiC P4000 and finally polished with a 1 µm diamond suspension (Struers). The characterization of the cross-sections was performed with a SEM (FEI QUANTA 200F Field Emission Gun) coupled to an EDAX detector for chemical analysis at 20 kV in low vacuum (0.9 mbar). Raman micro-spectrometry (Jobin Yvon LabRam HR800, λ=632.82 nm) was also carried out to identify the different oxide phases.

Results

Effect of the preoxidation time of the Ni powder

Fig. 3 shows the cross-section of the preoxidized nickel particles. The Raman spectra of Fig. 4 clearly demonstrate the unique formation of NiO on all particles at different oxidation times at 700°C [START_REF] Mironova-Ulmane | Raman scattering in nanosized nickel oxide NiO[END_REF]. The smaller nickel particles (Ni-2) present thicker oxide scales on their surface although the oxidation time is shorter (15 and 40 min) than for the coarse particles (1 and 4 h). Some of the grain boundaries of the coarse particles (Ni-1) are also oxidized. At increased magnifications (see inset in Fig. 3), the surface of Ni-1 exhibits a duplex scale (platelet and compact) after 1 h but the oxide is thicker and porous after 4 h of oxidation. It thus appears that the oxides get thicker with increasing time and number of cycles of oxidation. However, the thickness of the oxide layer varies among the particles (0.3 to 3 µm for 1 h and 1 to 5 µm for 4×1 h). Additional characterization of the coverage with NiO as a function of particle size was performed with "Image J" software (Fig. 5). It can be observed that the smaller the particles, the greater the oxidation. Also, the dependence of oxidation with the particle size decreases with increasing the number of oxidation cycles. Moreover, sintering of the particles caused by oxide growth is observed for all the oxidation times and particle sizes, in particular for the smallest ones (Ni-2). This indicates that the crushing between oxidation cycles was effective in avoiding sintering. In addition, the large particles seem to break along the oxidized grain boundaries and thus form smaller particles upon the repetitive crushing.

Ni-1.a (coarse particle preoxidized for 1 h) Ni-1.b (coarse particle preoxidized 4×1 h) Ni-2.a (small particle preoxidized for 3×5 min) Ni-2.b (small particle preoxidized for 4×10 min)

Reactivity of Al-Niox microparticles

Differential scanning calorimetry

The DSC results of the reference powders (Al; metallic Ni-1 and Ni-2; preoxidized Ni-1 and Ni-2) heated at 25°C/min till 1300°C are shown in Fig. 6. The data are summarized in Table 2. Coarse Ni (Ni-1) undergoes an event starting at 346°C that corresponds to the Curie temperature (358°C) [START_REF] Serway | Physics for Scientists and Engineers[END_REF]. In contrast, the same thermal event is more energetic and starts and finishes at a larger temperature interval when the particles are small (Ni-2). For Al, the onset and the maximum melting temperatures also differ (665 and 677°C for, respectively, the 2 and 25°C/min heating rates) but the overall endotherm is relatively similar.

The differences with respect the theoretical melting temperature (660.5°C [START_REF] Okamoto | Al-Cr (Aluminum-Chromium)[END_REF]) can be thus ascribed to the different heating rates [START_REF] Pedraza | Potential thermal barrier coating systems from Al microparticles. Mechanisms of coating formation on pure nickel[END_REF] that change the manner heat is absorbed. In the mixed Al + Ni powders, the first exothermic peak at 620±2°C (label 1 in Fig. 7) is most likely related to a solid-state reaction between Al and NiO [START_REF] Udhayabanu | Mechanical activation of aluminothermic reduction of NiO by high energy ball milling[END_REF][START_REF] Zhu | In-situ processing of NiAl-alumina composites by thermite reaction[END_REF]. This peak is followed by an endothermic one at 671±6°C, which corresponds to the melting point of Al (label 2 in Fig. 7) [START_REF] Okamoto | Al-Cr (Aluminum-Chromium)[END_REF]. For the coarse oxidized particles (Ni-1), the reaction between Al and Ni occurs at 852±14°C (label 3 in Fig. 7) and forms NixAly intermetallic compounds [START_REF] Padmavardhani | Synthesis and microstructural characterization of NiAl-Al2O3 functionally gradient composites[END_REF]. Right after, a last exothermic peak (label 4 in Fig. 7) starting above 1000°C can be related to the stabilization of the NiAl intermetallic compound and to the transformation of metastable γ-Al2O3 into stable α-Al2O3 [START_REF] Udhayabanu | Mechanical activation of aluminothermic reduction of NiO by high energy ball milling[END_REF][START_REF] Padmavardhani | Synthesis and microstructural characterization of NiAl-Al2O3 functionally gradient composites[END_REF]. Interestingly, the maximum temperature for NiAl formation is lower in the least preoxidized coarse particles than the ones oxidized for longer times (4×1 h).

For the small preoxidized particles (Ni-2), the same reactions occur until 750°C (Fig. 7, Table 3). However, a clear separation between NiAl and Al2O3 cannot be made at higher temperatures because the phenomena are now significantly exothermic. Further, the least oxidized particles (3×5 min) do not

Microstructural development

The microstructures and the composition of the mixed powders obtained at 750, 1000 and 1300°C are shown in Fig. 8 and Fig. 9. These temperatures correspond to the main exothermic reactions observed in the DSC curves of Fig. 7 (peaks 1, 3 and 4). Therefore in the rest of the study, we will focus on the coarse particles oxidized 4×1 h (Ni-1.b) because the oxide layers are thicker and more homogeneous (Figs. 3 and5), which should be favorable for the formation of micro-reservoirs. Concerning the small oxidized particles, the particles oxidized 4×10 min (Ni-2.b) will be discarded. Indeed, the oxide layers are conversely too important (Figs. 3 and5), which would consume a large part of aluminum, likely impeding the formation of an Al-rich intermetallic core rich for the micro-reservoirs.

Reactivity of preoxidized Ni over pure Al Differential scanning calorimetry

The DSC results of the different preoxidized Ni particles obtained with two heating ramps (2 and 25°C/min) are shown in Fig. 10 and the thermal data are summarized in Table 4. Reduction of NiO and melting of Al occur at respectively, 620±2°C and 670±2°C. For the two heating rates applied to the coarse particles, and for the heating rate of 25°C/min applied to the small ones, many peaks (label 3 in Fig. 10) appear at temperatures ranging between 750 and 1300°C.

Microstructural development

The microstructure and composition of the layout "preoxidized Ni over Al" heated in Ar following different temperature ramps (2 and 25°C/min) till 1300°C are respectively shown in Figs. 11 and12. The crucible can be divided into three main sections from the top to the bottom. A summarized description of these different sections is given in Table 5.

The top of the crucibles coincides with the area where the preoxidized Ni particles did not react with Al.

No noticeable difference can be observed between the samples (Figs 11.c and 11.d) except for the small oxidized particles for which the recrystallization of the NiO grains occurred forming compact facetted grains (Figs 12.c and 12.d).

The central section of the crucible corresponds to the interface area between preoxidized Ni and Al particles. This zone is composed of different layers of nickel, intermetallic compounds and alumina. The coarse Ni particles form a compact Ni-rich NiAl phase when heated at 25°C/min (label 1, Fig. 11.e).

Underneath, an Al-rich NiAl phase with Al2O3 segregated at the grain boundaries occurs (label 2 on Fig.

11.e). In between, many pores can be observed. The Al content increases further to the bottom (label 3 on Fig. 11e). For a heating ramp of 2°C/min (Fig. 11.f), the NiAl compounds are also observed but the Al2O3 appears more segregated in the Ni-rich NiAl than in the Al-rich counterpart. Further, the interface is not homogeneous and displays cracks rather than pores. Ni/Al2O3/Ni-rich NiAl/Al-rich NiAl are distributed on one side and Ni/Al2O3 Al-rich NiAl on the other but this can be due to the evolution of the solid powders upon melting.

NiAl (label 2 on Fig. 12.e) can be also observed. The interfaces between these clusters are composed by Al2O3.

For the bottom of the crucibles, there is barely any significant difference between the heating ramps of 25 and 2°C/min. In contrast, the compounds formed differ between the coarse and small Ni particles. Indeed, Al-rich NiAl particles (label 1 in Fig. 11.g) surrounded by Al3Ni2 and simultaneous Al2O3 formation (label 2 in Fig. 11.g) are developed from the coarse particles. In contrast, the small Ni particles develop both the Ni-rich and the Al-rich NiAl compounds (labels 1 and 2, respectively in Fig. 12.g). ). This difference can be due to the stresses generated by the crushing of the particles, which in turn foster nickel and oxygen interdiffusion and/or provoke the spallation of the columnar grains [START_REF] Evans | Stress effects in high temperature oxidation of metals[END_REF]. In addition, the increase in oxidation time brings about the intergranular oxidation of the Ni particles. The action of crushing effectively disassembles the grains and causes the reduction of the grain sizes and the creation of new un-oxidized surfaces.

The greater reactivity of the smaller Ni particles (Ni-2) is due to their greater specific surface. Such a greater reactivity is demonstrated through the oxide thickness that ranges between 0.5 and 2 µm after just 15 min (i.e. 3×5 min) (Fig 5 .c) and between 1.5 and 4 µm after 40 min (i.e. 4×10 min) (Fig 5 .d) at 700°C. The oxidation of the small particles is so fast that some of them appear completely oxidized (Fig.

3.d).

The small particles preoxidized for 40 min were thus consequently discarded as they would not allow the formation of Al-rich NixAly cores in the particle by reaction with Al. Similarly, the coarse particles preoxidized for 1 h did not develop sufficiently thick nickel oxide to results in Al2O3 shells upon the aluminothermic reaction. The coarse and small particles preoxidized for, respectively, 4 h (i.e. 4×1 h) and 15 min (i.e. 3×5 min) were then used to study the formation of the NixAly + Al2O3 phases.

Layout 1-Mixture of preoxidized Ni and Al particles

The DSC of the Al and preoxidized Ni powders show two peaks at temperatures corresponding to the melting of Al and to the Curie temperature of Ni. When these two powders are mixed together and heated, various additional phenomena appear. The first exothermic peak at 620±2°C is most likely related to the solid-state aluminothermic reaction given in equation ( 1) [START_REF] Padhan | Mechanical activation on aluminothermic reduction and magnetic properties of NiO powders[END_REF]. The heat released upon this highly exothermic reaction causes the melting of Al and the stop of the solid-state aluminothermic reaction. Whereas by solid state or by molten Al, the whole NiO is consumed by Al to form Al2O3 and Ni (Figs. [START_REF] Chmielewski | Processing, microstructure and mechanical properties of Al2O3-Cr nanocomposites[END_REF].b and 13.c). When the temperature exceeds the melting of Al, the dissolution of Ni can readily occur [START_REF] Bouchaud | Slurry aluminizing mechanisms of Ni-based superalloys incorporating an electrosynthesized ceria diffusion barrier[END_REF][START_REF] Galetz | The role of combustion synthesis in the formation of slurry aluminization[END_REF]. Ni and Al are thus quickly consumed and form Ni3Al (eq. 2) (Fig 13 . c).

3Ni(s) + Al(l) Ni3Al(s) (eq. 2)

With increasing temperature between 740 and 1000°C, further Al inward diffusion occurs and the coarse Ni particles thus develop a concentration gradient between Ni at the core, the intermediate layer of Ni3Al and the external layer of NiAl (Fig. 13.c). The formation of NiAl is accompanied of a strong exothermic peak shown in Fig. 7 following equation (

The residual Ni is fully transformed with the subsequent heating till 1300°C. The final microstructure is made of clusters of NiAl particles with many cracks (Fig. 13.d). In contrast, no alumina resulting from the aluminothermic reaction could be detected under the SEM, probably because it is either too thin and/or it came up to the surface given its much lower density than that of the Ni-Al intermetallic compounds.

For the small particles, the exothermic peak at 620°C associated with the NiO reduction is also observed but the energy released is lower than with the coarse particles. As a result, a large amount of the NiO grown on the small Ni particles has not reacted yet at about 656°C (Fig. 13.f). Between 750 and 1000°C, the exothermal reaction can be associated with the solid-liquid reduction reaction that brings about an exothermic reaction (Fig. 7). At this range of temperature three types of particles are formed: NiAl/Al2O3, Ni/Ni3Al/NiAl and Ni/NiO/Ni multi-layered microstructures (Fig. 13.g). For the first one, the Ni formed the periphery of the surface is the result of the aluminothermic reaction. These features would imply that the reaction is controlled by the solid-state diffusion of species through the products layers because of the incomplete dissolution of NiO. For the second and the third ones, all the NiO is consumed to form Ni and Al2O3. Moreover, for the third case of Ni(NiO)Ni particles, the Al2O3 was expelled.

The full consumption of NiO, Ni and Al at 1050°C results in Ni3Al and Al2O3 and the metastable γ-Al2O3 transforms into the stable α-Al2O3 that releases heat (Fig. 7) [START_REF] Pedraza | Potential thermal barrier coating systems from Al microparticles. Mechanisms of coating formation on pure nickel[END_REF][START_REF] Eisenreich | On the mechanism of low temperature oxidation for aluminum particles down to the nano-scale[END_REF]. Moreover, isolated alumina particles can also be found inside the particles of Ni3Al.

Layout 2-preoxidized Ni powders over the Al particles

The comparison of the results obtained with the mixed particles with that for the layout "preoxidized Ni over Al" clearly points out the tremendous impact of the contact between the reactive particles. Indeed, for the same particle size, the melting enthalpy of Al is much lower with the mixed powders (∼ 65 J g -1 ) than with the double layer configuration (∼ 100 J g -1 ). Therefore, the greater contact area between Al and NiO results in the formation of more Al2O3 absorbing heat and hampering melting of Al. In contrast, when the preoxidized Ni powder is deposited over the Al particles, a front of reactions and of heat propagates from the interface, resulting in different compounds and microstructures.

In this latter configuration (Figs. 14.a and 15.a), the influence of the heating rate appears negligible. Indeed, the DSC curves present the same reaction peaks (within a couple of degrees of difference) and the SEM images exhibit similar zones. However, it has been noted that the extent of the Al-rich areas is greater than that of the Ni-rich area when a slow heating rate (2°C/min) has been applied and vice-versa. This can be related to the fact that diffusion of aluminum is favored below 1000°C [START_REF] Brossard | Modelling of aluminized coating growth on nickel[END_REF] (inward diffusion) while the outward diffusion of Ni is promoted above 1000°C [START_REF] Brossard | Modelling of aluminized coating growth on nickel[END_REF]. The system therefore spends a longer time at low temperature with slow heating rates.

When the coarse preoxidized Ni particles are heated till 550°C, the solid/solid aluminothermic reaction begins [START_REF] Pedraza | Potential thermal barrier coating systems from Al microparticles. Mechanisms of coating formation on pure nickel[END_REF]. The external NiO is reduced by Al, bringing about the formation of alumina and of nickel (Fig. For the small particles the effect of the heating rate is very significant as demonstrated by the different DSC peaks (Fig. 10) and SEM microstructures at the center of the crucibles (Fig. 12). These differences can result from the preoxidation of the powders whereby the oxide growth at 700°C is ruled by the faster diffusion rate of the Ni 2+ ions than that of the O 2-ions. As a consequence, pores and cavities are generated at the oxide/metal interface [START_REF] Jeangros | Oxidation mechanism of nickel particles studied in an environmental transmission electron microscope[END_REF]. Further the particles are subjected to intermediate air quenching and crushing (3x15 min), which induce stress in the oxide and in the metallic Ni particle. These phenomena are enhanced on the small particles vs. the coarse ones because of their greater oxide thickness and less spherical morphology. Overall, the small particles become more reactive than the coarse ones.

When the temperature of 550°C is reached, the solid/solid aluminothermic reaction begins [START_REF] Padmavardhani | Synthesis and microstructural characterization of NiAl-Al2O3 functionally gradient composites[END_REF], NiO and Al react at the interface of the preoxidized Ni particles (Fig. When the temperature reaches 950°C, a third exothermic peak appears when heating at 25°C/min (Fig. 10 peak labeled 3'). This peak results from the transformation of the Ni phase into Ni3Al and NiAl but the reaction does not take place continuously. This reaction can be explained by the difference in the thermal expansion coefficient between Ni and NiO and the potential residual stress and cavities that cause cracking and spallation of the NiO layer. Therefore, Al becomes free to react with Ni and form Ni3Al and NiAl. The NiO expelled continue its reduction in the molten Al to form Al2O3 that gets trapped between the grains of the intermetallic compounds as they grow. This step will be repeated until the intermetallic layer thickens enough to hamper the arrival of Al. The heat produced during this stage and the very low amount of trapped oxide do not prevent the sintering of the particles to each other to form a compact structure. The Ni3Al layer will subsequently enrich with Al to form NiAl (Fig. 15.d).

In contrast, cracking and spallation of NiO do not occur when heating at 2°C/min. Therefore, when the temperature reaches 750°C, the solid/liquid aluminothermic reaction occurs and the remainder of the NiO near the interface is finally converted into an adherent Al2O3. Simultaneously, the NiO "barrier" over the Ni particles being reduced, the Ni core can now react directly with the molten Al (Fig. 15.c). The greater amount of alumina formed at this stage aggregates around the particles and prevent their sintering. As a result, clusters of intermetallic compounds with significant porosity are obtained (Fig. At the surface of these clusters, almost no Al2O3 was observed, which suggests that Al2O3 either dissolves and/or that is transported to the surface of the melt given its lower density than the intermetallic compounds. These results, therefore, suggest that the synthesis of self-regenerating coating could be obtained by e.g. spraying slurries containing mixtures of preoxidized Ni and metal Al or multilayering them onto a substrate. Subsequent heat treatments with relatively fast heating ramps (25°C) till moderate temperatures (∼ 750°C) assisted by the aluminothermic reaction between NiO and Al could consolidate the formation of the micro-containers trapped in the coating.

Fig. 1 :

 1 Fig. 1: Ternary phase diagram of Ni-Al-O at 1638 K (1365°C) [14].

Fig. 2 :

 2 Fig. 2: Schematic drawings representing the two different layouts: (a) mixture of preoxidized Ni and of Al particles; (b) preoxidized Ni over the Al particles.

Fig. 3 :

 3 Fig.3: SEM (BSE mode) cross-section images on the coarse (Ni-1.a & Ni-1.b) and small (Ni-2.a & Ni-2.b) preoxidized Ni particles at 700°C in air for different times.

Fig. 4 :

 4 Fig. 4: Raman spectra of the oxides formed on the coarse (Ni-1) and small (Ni-2) particles after oxidation at 700°C in air fordifferent times[START_REF] Mironova-Ulmane | Raman scattering in nanosized nickel oxide NiO[END_REF].

Fig. 5 :

 5 Fig. 5: Area proportion of NiO (SNiO/SNi ×100) as a function of particle size on the coarse (Ni-1.a & Ni-1.b) and small (Ni-2.a & Ni-2.b) preoxidized Ni particles.

Fig. 6 :

 6 Fig. 6: DSC thermograms for the reference powders (Al, Ni and Ni oxidized) upon heating at 2 and 25°C/min till 1300°C in Ar.

Fig. 7 :

 7 Fig. 7: DSC thermograms of the mixtures of Al with the different Ni particles upon heating 25°C/min till 1300°C in Ar.

  In the coarse Ni particles (Ni-1.a and Ni-1.b), all the NiO formed by preoxidation is consumed by the aluminothermic reaction at 671±6°C. As a result, the particles present a homogeneous composition of about 30Al-70Ni (at %) (Figs. 8.a and 8.b). At 1000°C, quite similar concentration gradients of Al and Ni between the core and the external areas are observed in Ni1.a and Ni1.b (labels 1 to 4 in Figs. 8.c and 8.d). Concentrations of Al and Ni of 21-79; 37-63; 48-52 and 58-42 at % are respectively retrieved from the core to the surface. After the full treatment at 1300°C, all the particles show a fully homogenized NiAl phase with 46Al-54Ni (at %). Moreover, cracks and Kirkendall porosity appear.

Fig. 8 :Fig. 9

 89 Fig. 8: SEM (BSE mode) of the mixture between Al and coarse Ni preoxidized particles (Ni1.a & Ni2.b) after DSC treatments at 25°C/min for different temperatures.

Fig. 10 :

 10 Fig. 10: DSC thermograms of the preoxidized Ni over Al upon heating (2 and 25°C/min) till 1300°C in Ar.

Fig. 12 :

 12 Fig. 12: SEM (BSE mode) cross-sections of the layout 2 with small preoxidized Ni particles (Ni-2.a over Al) after DSC treatments at 25 and 2°C/min. (a) and (b) overview; (c) and (d) top; (e) and (f) center; (g) and h) bottom of the crucible. Heat treatment (crosssection)

1

 1 Effect of the preoxidation timeTwo different Ni powders were oxidized at 700°C under air. The coarse particles (Ni-1) were crushed every hour to avoid the formation of different NiO morphologies and sintering of the oxidized powders. Indeed, porous equiaxial and duplex layer (equiaxial and columnar oxides) forms at the periphery of the Ni coarse particles during oxidation (Fig.3.a). According to the works of Peraldi et al., a porous NiO scale forms at 700°C for short times[START_REF] Peraldi | Correlations between growth kinetics and microstructure for scales formed by high-temperature oxidation of pure nickel, I. Morphologies and Microstructures[END_REF][START_REF] Peraldi | Correlations between growth kinetics and microstructure for scales formed by high-temperature oxidation of pure nickel, II. Growth Kinetics[END_REF]. With increasing the oxidation time to 3 h, a duplex scale grows with the porous and equiaxed layer surmounted by columnar grains but the thickness remains almost constant (about 1 µm)[START_REF] Peraldi | Correlations between growth kinetics and microstructure for scales formed by high-temperature oxidation of pure nickel, I. Morphologies and Microstructures[END_REF][START_REF] Peraldi | Correlations between growth kinetics and microstructure for scales formed by high-temperature oxidation of pure nickel, II. Growth Kinetics[END_REF]. In our case, the formation of the duplex layer appears along the first hour of oxidation with oxide thicknesses ranging from 0.3 to 3 µm (Figs. 3.a and 5.a). The existence of these two types of oxide morphologies within the system may be due to the larger active surface of the Ni powders compared to the average grain size (2 mm) of the Ni studied by Peraldi et al.[START_REF] Peraldi | Correlations between growth kinetics and microstructure for scales formed by high-temperature oxidation of pure nickel, I. Morphologies and Microstructures[END_REF][START_REF] Peraldi | Correlations between growth kinetics and microstructure for scales formed by high-temperature oxidation of pure nickel, II. Growth Kinetics[END_REF]. When the exposure time is extended to 4 h with crushing of the powder every hour, the scale morphology is exclusively porous with thicknesses ranging from 1 to 5 µm (fig3.b and 5.b

Fig. 13 :

 13 Fig. 13: Schematic model of NiO reduction and formation of intermetallic particles for a mixture of (a-d) coarse preoxidized Ni (4h) and Al and (e-h) small preoxidized (15 min) Ni and Al heated at f 25°C/min.

  14.b). The resulting particles close to the interface (between Al and Niox particles) are then exposed to the flux of Al, which in turn allows the formation of intermetallic compounds (Ni3Al and NiAl) successively towards the core (Fig.14.b). This stage progressively continues along the crucible until the onset of Al melting is reached (645±1°C). Part of this molten aluminum may be transported upwards by capillarity between the nickel particles. Therefore, the nickel particles get embedded in the molten aluminum (Fig.14.c). Moreover, the melting of aluminum induces a preferential dissolution of Ni at the grain boundaries and the subsequent formation of cracks, either by shrinkage upon quenching or through the formation of denser phases. Beyond 750°C, aluminum keeps on diffusing through the Ni particles. Therefore, the NixAly intermetallics get progressively richer in Al at the expense of Ni-rich NixAly. In addition, the grains of the newly formed intermetallic phases grow during this stage causing the disappearance of porosity. The molten aluminum between these grains will be gradually "absorbed" by the Ni particles leaving behind alumina trapped between the grains (Fig.14.c).Between 1000 and 1300°C, the microstructure barely evolves while the composition slightly changes due to favored nickel diffusion. This results in Al-rich and Ni-rich NiAl compounds (Fig.14.d). Quenching to room temperature provokes solidification of molten Al and the segregation of Al3Ni2 between the grains of Al-rich NiAl (Fig.14.d).

Fig. 14 :

 14 Fig. 14: Schematic drawings of the aluminothermic reduction and reactions between the coarse preoxidized Ni over Al particles (double layer configuration) heated at 25°C/min.

  15.b). As opposed to the coarse particles, the thickness of the oxide layer is now quite significant and the aluminothermic reaction is thus momentarily superseded by the onset of the Al melting at 648±1°C. Molten Al flows upwardly by capillarity between the Ni preoxidized particles NiO grains causing the entrapment of nickel particles in the molten Al and the enrichment of the NiO in Al (Fig.15.c). Beyond 750°C, the second aluminothermic reaction (liquid/solid) begins. The remainder of the NiO near the interface is thus finally converted into an adherent Al2O3 (Fig.15.c). Simultaneously, the NiO "barrier" over the Ni particles is transformed in Al2O3 creating cracks, the Ni core is directly exposed to the molten Al (Fig.15.c). The greater amount of alumina formed at this stage aggregates around the particles and prevents their sintering. As a result, clusters of intermetallic compounds with a lot of porosity are obtained (Fig.15.c).

  15.c) Deeper in the crucible and irrespective of the heating rate, clusters of Al3Ni2 and Al3Ni are formed upon cooling (Fig. 15.d).

Fig. 15 :

 15 Fig. 15: Schematic drawings of the aluminothermic reduction and reactions between the small preoxidized Ni an Al (double layer configuration) heated at 2 and 25°C/min

Table 1 :

 1 Main features of initial Ni, Al and preoxidized Ni particles.

	The coarse powder (Ni-1) was heated at 700°C for 1h (Ni-1.a) and 4h (Ni-1.b) in air (intermediate

crushing in a mortar every hour). The small-sized powder (Ni-2) underwent the same preparation but the oxidation time was of 15 min (Ni-2.a) and 40 min (Ni-2.b) (intermediate crushing every 5 and 10 min, respectively). The intermediate crushing step was realized to increase the oxide thickness and avoid sintering. Such partial oxidation is assumed to present two simultaneous advantages, i.e. to reduce the energy release during the thermite reaction and to lower the content of one of the two reagents like in the works of Vrel et al.

Table 2 :

 2 Thermal data of the effects observed for the reference powders (Al, Ni and Ni oxidized) upon heating at 2 and 25°C/min till 1300°C in Ar.

Table 3 :

 3 Thermal data and potential transformations in the mixtures of Al with the different Ni particles.

		Tmax (°C)	ΔT(°C)	ΔH (J g -1 )	Tmax (°C)	ΔT(°C)	ΔH (J g -1 )
	Transformation	Al + Ni-1.a (coarse oxidized Ni particles: 1 h)	Al + Ni-1.b (coarse oxidized Ni particles: 4×1 h)
	NiO reduction	620	526-656	-73	618	540-656	-78
	Al melting	671	656-712	63	669	656-711	69
	NiAl formation	838	742-988	-596	867	740-1000	-576
	γ-Al2O3 α-Al2O3	1095	988-1252	-54	1080	993-1142	-52
			Al + Ni-2.a			Al + Ni-2.b
		(small oxidized Ni particles: 3×5 min)	(small oxidized Ni particles: 4×10 min)
	NiO reduction	621	555-655	-66	621	551-655	-72
	Al melting	668	655-724	61	670	655-727	81
	NiAl formation and γ-Al2O3 α-Al2O3	1023	754-1244	-1457	1036	742-1136	-1290
	NiAl formation	---	---	---	1192	1136-1272	-314

Table 4 :

 4 Thermal data and the potential transformations in the preoxidized Ni over Al layout

		Tmax (°C)	ΔT(°C)	ΔH (J g -1 )	Tmax (°C)	ΔT(°C)	ΔH (J g -1 )
	Heating ramp		2°C/min			25°C/min
	Transformation		Ni-1.b / Al			Ni-1.b / Al
	NiO reduction	612	529-649	-50	622	577-648	-9
	Al melting	665	649-671	99	678	648-732	104
	NiAl formation	1086	1056-1124	-85	1275	1123-1292	-52
			Ni2.a / Al			Ni2.a / Al
	NiO reduction	556	511-549	-76	612	560-648	-14
	Al melting	663	649-670	108	675	648-766	95
	Ni3Al and NiAl
	formation	---	---	---	1081	757-1300	-137
	γ-Al2O3 α-Al2O3

Table 5 :

 5 Summary of the different zones observed metallographically

	Heating ramp	25°C/min		2°C/min
	Particles	Ni-1.a / Al	Ni-2.a / Al	Ni-1.b / Al	Ni-2.a / Al
	Zone 1 (top)	Ni particles with NiO	Ni particles with NiO recrystallized	Same as Ni 1.a at 25°C/min	Same as Ni 2.a at 25°C/min
	Zone 2 (center)	-Ni particles (10 at % Al) -NiAl (Ni rich) -NiAl (Al rich) with Al2O3 grain boundaries	-Ni particles with NiO Interface of Al2O3 -Ni3Al -NiAl (Ni rich) with grain boundaries of Al2O3	Same as 25°C/min	-Ni particles with NiO -Interface of Al2O3 -cluster of NiAl particles with grain boundaries of Al2O3
			-cluster of NiAl (Al rich)
	Zone 3 (bottom)	-grains with a NiAl core and Al3Ni2 at the Al3Ni2 edge + Al2O3 in grain	particles with grain boundaries of Al2O3 with grain boundaries of -cluster of Al3Ni2 particles -interface of Al2O3	Same as Ni 1.a zone developed 25°C/min but more	Same as Ni 2.a 25°C/min
			Al2O3


Yue Tong

Wan-Qing Shen

Jian-Fu Shao
email: jian-fu.shao@polytech-lille.fr

An adaptive coupling method of state-based peridynamics theory and finite element method for modeling progressive failure process in cohesive materials

Keywords: State-based peridynamics theory, Finite element method, Adaptive coupling, Crack propagation, Progressive failure, Cohesive materials

Introduction 1

Failure in cohesive geological and cement-based materials is relevant to the transition from diffuse 2 micro-cracking to localized macro-fracturing [START_REF] Bažant | Instability, ductility, and size effect in strain-softening concrete[END_REF][START_REF] Shao | A microcrack-based continuous damage model for brittle geomaterials[END_REF][START_REF] Zhao | Analysis of localized cracking in quasi-brittle materials with a micro-mechanics based friction-damage approach[END_REF].

3

The macroscopic failure is generally a progressive process due to the internal cohesion of those materials 4 [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF][START_REF] Rots | Computational modeling of concrete fracture[END_REF][START_REF] Moës | Extended finite element method for cohesive crack growth[END_REF]. Unlike propagation modeling of a single 5 existing crack, modelling of the progressive failure from the onset, coalescence and localization of multiple 6 cracks is still open issue.

7

Considerable efforts have been provided and different kinds of numerical methods have been developed 8 during recent decades. Not exhaustively but representatively, some of them are here mentioned. The en-9 riched finite element method (EFEM) [START_REF] Oliver | Modelling strong discontinuities in solid mechanics via strain softening constitutive equations, part 1: fundamentales[END_REF] with elementary enrichments and the eXtended finite 10 element method (XFEM) [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Stolarska | Modelling crack growth by level sets in the extended finite element method[END_REF] with node enrichments have been widely used for dealing with cracking problems with strong discontinuities, even with hydromechanical coupling [START_REF] Zeng | Numerical study of hydraulic fracture propagation accounting for rock anisotropy[END_REF][START_REF] Zeng | Study of hydraulic fracturing in an anisotropic poroelastic medium via a hybrid edfm-xfem approach[END_REF]. Some comparative studies of different kinds of numerical methods can be found in [START_REF] Jirásek | Comparative study on finite elements with embedded discontinuities[END_REF][START_REF] Dias-Da Costa | A comparative study on the modelling of discontinuous fracture by means of enriched nodal and element techniques and interface elements[END_REF]. But in those previous studies, the choice of specific enrichment functions and crack propagation criteria was a key issue and not an easy task [START_REF] Wu | Phase field modeling of fracture[END_REF][START_REF] Wu | Computational modeling of localized failure in solids: Xfem vs pf-czm[END_REF].

Further, the ability of those methods in dealing with crack coalescence and branching as well as multiple cracks propagation is still not fully demonstrated. More recently, the phase-field methods have been developed by approximating the sharp crack topology by the regularised one [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]. This class of methods attracted more and more attention due to their efficiency for describing the transition from diffuse damage to localized cracking. However, the development of the real crack is represented by the localized damage band which requires an extra scale length parameter to characterize its width [START_REF] Molnár | 2d and 3d abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF]Gravouil, 2017, Wang et al., 2019).

On the other hand, the peridynamics (PD) theory has been developed as an extension of classical continuum mechanics framework [START_REF] Silling | Reformulation of elasticity theory for discontinuities and long-range forces[END_REF]. Instead of solving partially derivative equations, which is particularly delicate with the presence of crack singularities, integral motion equations are solved in PD theory. Therefore, this theory is particularly efficient for dealing with cracking problems [START_REF] Madenci | Peridynamic theory and its applications[END_REF]. The cracking process of continuum is directly related to the damage of internal material bonds or links. Therefore, the initiation and propagation of cracks can be naturally described without introducing any extra crack propagation criteria or crack-tracking method. In terms of interactions between internal material points, two different formulations of the PD theory are proposed. The bond-based PD theory was first formulated by only considering pair-wise interactions between two neighboring points [START_REF] Silling | Reformulation of elasticity theory for discontinuities and long-range forces[END_REF]. More recently, some improvements have been proposed in the bond-based PD theory by introducing rotation effect in order to avoid the restriction of Poisson's ratio (Zhu andNi, 2017, Ni et al., 2019). In order to describe more complex deformation mechanisms, the state-based PD theory has further been developed [START_REF] Silling | Peridynamic states and constitutive modeling[END_REF]. The motion of each material point depends on the interactions with all other points inside a limited surrounding zone [START_REF] Silling | Peridynamic states and constitutive modeling[END_REF]. The classical concepts of local stress and strain tensors are generalized to the non-local force and deformation states. The constitutive equations are thus described by the relations between such states. In addition, the state-based PD theory is able to conveniently describe distortional and volumetric deformations other than the single volumetric deformation in the bond-based one [START_REF] Silling | Reformulation of elasticity theory for discontinuities and long-range forces[END_REF]Lehoucq, 2010, Madenci and[START_REF] Madenci | Peridynamic theory and its applications[END_REF]. Thus, the state-based PD is suitable for modeling the deformation and failure of geological materials, for which the shearing induced volumetric dilatancy is an important feature. However, owing to the non-local formulations, the numerical methods based on the PD theories are computationally time-consuming. In view of studying large scale boundary values problems, it is generally convenient to combine the PD theories in cracking zones with the finite element method (FEM) in elastic or plastic zones. A number of works have been reported on the combination of bond-based PD theory and FEM [START_REF] Macek | Peridynamics via finite element analysis[END_REF], Shojaei et al., 2016[START_REF] Zaccariotto | Coupling of fem meshes with peridynamic grids[END_REF]. More recently, the state-based PD theory was also been combined with the FEM framework [START_REF] Madenci | A state-based peridynamic analysis in a finite element framework[END_REF][START_REF] Bie | A coupling approach of state-based peridynamics with node-based smoothed finite element method[END_REF].

In terms of description of bond damage, different kinds of criteria have been formulated. For instance, a critical bond stretch criterion for the prototype micro-elastic brittle (PMB) material was first proposed and implemented in the bond-based PD theory in [START_REF] Silling | A meshfree method based on the peridynamic model of solid mechanics[END_REF]. This criterion has been widely applied to simulate brittle fracture problems [START_REF] Javili | Peridynamics review[END_REF][START_REF] Diehl | A review of benchmark experiments for the validation of peridynamics models[END_REF]. Some improvements have been obtained by using a trilinear bond model [START_REF] Yang | Investigation on mode-i crack propagation in concrete using bond-based peridynamics with a new damage model[END_REF]. Similar critical bond stretch criteria have also been introduced into the state-based PD theory for brittle materials (Madenci andOterkus, 2014, Zhang and[START_REF] Zhang | A state-based peridynamic model for quantitative fracture analysis[END_REF]. However, all those brittle models are not able to correctly describe the progressive failure process in cohesive materials such as concrete (Petersson, 1981, Hoover and[START_REF] Hoover | Cohesive crack, size effect, crack band and work-of-fracture models compared to comprehensive concrete fracture tests[END_REF]. In particular, the progressive softening behavior on the post-peak regime is badly reproduced.

More recently, a new bond continuum damage model has been formulated in the framework of the bondbased PD theory [START_REF] Tong | A new bond model in peridynamics theory for progressive failure in cohesive brittle materials[END_REF]. In order to overcome the limitations of the bond-based PD theory, the bond continuum damage model is here extended to the state-based PD theory, in order to better modeling the progressive failure process in cohesive geological materials, especially for the compressive shearing conditions. Further, an adaptive switching strategy is also developed for coupling with the standard finite element method. The capability of the proposed bond damage model as well as the coupling strategy will be verified through three representative laboratory tests. This paper is organized as follows. The general framework of the ordinary state-based PD theory is first presented in section 2. In section 3, a new bond continuum damage model for cohesive materials is introduced into the ordinary state-based PD theory. In section 4, an adaptively strategy for coupling the PD theory and finite element method is established. In section 5, the effectiveness of the coupling method considering the damage or not is verified.

The ordinary state-based PD theory

The state-based PD theory can be seen as a generalized framework for the classical solid mechanics [START_REF] Silling | Peridynamic states and constitutive modeling[END_REF][START_REF] Silling | Peridynamic theory of solid mechanics[END_REF]. The kernel is the use of the state concept for mapping deformation at a material point by including the information of all internal bonds linking to other material points within its horizon zone. As shown in Figure 1, every material point x interacts directly with other material points x within its horizon H x determined by the prescribed horizon radius δ. Due to the external loading, the solid body deforms. Accordingly, the material points x, x at the initial configuration are moved to y, y in the deformed configuration. Now the vector X is defined as the reference vector state mapping the initial position of one bond as:

X x -x = x -x (1)
and Y as the deformation vector state mapping the bond into its deformed image as:

Y[x, t] x -x = y(x , t) -y(x, t) (2)
In these definitions y(x, t) = x + u(x, t), y(x , t) = x + u(x , t) , u(x, t) and u(x , t) are the displacement vectors associated to x and x respectively. Besides, let T define the force density vector state. The force vector states related to x and x are respectively T[x, t] and T[x , t]. If the force state T is aligned with the deformation state Y, the PD theory is named as the ordinary state-based. Otherwise, it is called as the non-ordinary state-based PD theory. In this paper, the ordinary state-based PD theory is applied. In this case, the motion equation of a material point is defined as: as [START_REF] Silling | Peridynamic states and constitutive modeling[END_REF]:

ρ(x) ü(x, t) = H x {T[x, t] x -x -T[x , t] x -x }dV x + b(x, t) (3)
where ρ(x) is the mass density, V x the volume associated with the material point x , and b(x, t) the external body force density. For numerical implementation, the integral formulation (3) is approximation by the discrete for a given PD material point x k :

ρ(x k ) ü(x k , t) = N j=1 (T x j -x k -T x k -x j )V j + b(x k , t) (4)
in which N is the number of material points (x j ) within the horizon of x k , V j is the volume of x j . For the case of static or quasi static problems, (4) is simplified to:

0 = N j=1 (T x j -x k -T x k -x j )V j + b(x k ) (5)
As in the classical continuum mechanics, the force state can be obtained from the derivation of an energy functional W(x) with respect to the deformation state [START_REF] Madenci | Peridynamic theory and its applications[END_REF]. For instance, the force states for the material points x k and x j can be expressed as follows:

T x j -x k = 1 V j ∂W(x k ) ∂(|y j -y k |) y j -y k |y j -y k | (6)
and

T x k -x j = 1 V k ∂W(x j ) ∂(|y k -y j |) y k -y j |y k -y j | (7)
For isotropic and elastic materials, the generalized form of the strain energy density function W(x) at any PD material point x k is defined as:

W(x k ) = aθ(x k ) 2 + b N j=1 ω k j (|y j -y k | -|x j -x k |) 2 V j (8)
The term θ(x k ) is related to volume dilatation and calculated by:

θ(x k ) = d N j=1 ω k j (|y j -y k | -|x j -x k |) y j -y k |y j -y k | x j -x k |x j -x k | V j (9)
In these expressions, a, b and d are three PD parameters which are related to the macroscopic elastic constants of materials and the PD horizon radius. ω k j is a non-dimensional influence function which reflects the degree of non-locality between material points so that the interaction between points decreases as the distance increases [START_REF] Madenci | Peridynamic differential operator and its applications[END_REF].

With the assumption of small deformation, the following simplifying relations can be adopted:

|u(x j ) -u(x k )| << |x j -x k | (10) y j -y k |y j -y k | x j -x k |x j -x k | ≈ 1 (11)
Further, by defining the influence function ω k j as:

ω k j = δ |x j -x k | ( 12
)
the force vector states in ( 6) and ( 7) can be rewtitten as:

T x j -x k = ( 2ad 2 δ 2 |x j -x k | N m=1 |y m -y k | -|x m -x k | |x m -x k | V m + 2bδ |y j -y k | -|x j -x k | |x j -x k | ) x j -x k |x j -x k | (13)
and

T x k -x j = ( 2ad 2 δ 2 |x k -x j | N n=1 |y n -y j | -|x n -x j | |x n -x j | V n + 2bδ |y k -y j | -|x k -x j | |x k -x j | ) x k -x j |x k -x j | (14)
where x m and x n denote the material points within the respective horizon of x k and x j . Now by using the equivalence condition between the classical continuum mechanics and PD theory of the elastic strain energy respectively for the shear strain part W(x) and the volumetric strain part θ(x), the PD parameters a, b and d can be identified [START_REF] Madenci | Peridynamic theory and its applications[END_REF]:

                               a = 1 2 (k - 5 3 µ), b = 15µ 2πδ 5 , d = 9 4πδ 4 , 3D a = 1 2 (k -2µ), b = 6µ πhδ 4 , d = 2 πhδ 3 , 2D a = 0, b = E 2Aδ 3 , d = 1 2Aδ 2 , 1D . ( 15
)
where k, µ and E are the macroscopic bulk modulus, the shear modulus and the Young's modulus, h is the thickness for 2D geometry, A is the cross sectional area for 1D geometry. And k and µ are expressed of Young's modulus E and Poisson's ratio ν by:

k =                                E 3(1 -2ν) , 3D E 2(1 -ν) , plane stress E 2(1 + ν)(1 -2ν)
, plane strain

, µ = E 2(1 + ν) (16)
As in the bond-based PD theory (Foster et al., 2011, Zhang and[START_REF] Zhang | A state-based peridynamic model for quantitative fracture analysis[END_REF], the stretch (deformation)

of the bond linking a pair of two points x k and x j is denoted as s k j and calculated by:

s k j = |y j -y k | -|x j -x k | |x j -x k | (17)
Substituting ( 17) for ( 13) and ( 14), the force states are explicitly expressed as functions of the bond stretches:

T x j -x k = ( 2ad 2 δ 2 |x j -x k | N m=1 s km V m + 2bδs k j ) x j -x k |x j -x k | (18)
and

T x k -x j = ( 2ad 2 δ 2 |x k -x j | N n=1 s jn V n + 2bδs jk ) x k -x j |x k -x j | (19)
The interaction between two materials points can now be defined by a generalized pairwise force density function as follows:

f x j -x k = T x j -x k -T x k -x j (20)
Substituting ( 18) and ( 19) for (20), one gets:

f x j -x k =         4bδs k j + 2ad 2 δ 2 |x j -x k | ( N m=1 s km V m + N n=1 s jn V n )         x j -x k |x j -x k | (21)
Like the that defined in the bond-based PD theory, the generalized pairwise force density function is also a function of stretches of paired points. However, as a fundamental difference, it not only depends on the own stretches of the paired points but also on the stretches associated with the points within their horizons.

Based on the above relation between the pairwise force and stretch, it appears convenient to introduce a local damage criterion as a function of the bond stretch to describe progressive failure of materials. That means that when the stretch of a bond reaches a critical value, the bond is progressively broken and the corresponding pairwise force decreases. The effect of bond damage on the bond force is here described by introducing a history dependent scalar function ζ k j , which represents the status of the invented bond connecting x k and x j . Consequently, the generalized pairwise force density function is updated by taking into account the bond damage:

f x j -x k = ζ k j         4bδs k j + 2ad 2 δ 2 |x j -x k | ( N m=1 ζ km s km V m + N n=1 ζ jn s jn V n )         x j -x k |x j -x k | (22)
In the basic model for prototype micro-elastic brittle (PMB) material [START_REF] Silling | Reformulation of elasticity theory for discontinuities and long-range forces[END_REF]Askari, 2005, Madenci and[START_REF] Madenci | Peridynamic theory and its applications[END_REF], the status function ζ k j is defined as follows:

ζ k j =                  1, s k j < s c 0, otherwise . ( 23
)
where s c is the critical bond stretch. According to this model, the interaction force between the paired points increases linearly with the stretch and drops abruptly as the stretch reaches the critical value, as shown in Figure 2a. The value of s c can be determined by considering that the work required to break all bonds across a newly created crack surface is equivalent to the macroscopic critical energy release rate G c , one gets:

s c =                        G c (3µ + ( 3 4 ) 4 (k -5µ 3 ) δ , 3D G c 6 π µ + 16 9π 2 (k -2µ) δ , 2D . (24)
With a fixed value of Poisson's ratio ν=1/4 for 3D cases or 2D plane strain cases and ν=1/3 for 2D plane stress cases, the terms (k -5µ 3

) and (k -2µ) in ( 24) vanish and one gets the critical bond stretch value used in the bond-based PD theory.

Finally, the macroscopic damage state at any PD point x k is quantified by the scalar variable ϕ(x k ) ∈ [0, 1], which defines the ratio of the number of broken bonds to the total number of bonds:

ϕ(x k ) = 1 - N j=1 ζ k j /N (25)

Bond damage model for cohesive materials

The basic elastic brittle model illustrated in Figure 2a has widely used in modeling brittle failure in elastic materials generally subjected to tensile stresses [START_REF] Javili | Peridynamics review[END_REF][START_REF] Diehl | A review of benchmark experiments for the validation of peridynamics models[END_REF]. However, as mentioned above, the failure of cohesive materials under compressive stresses is a progressive process through the initiation and propagation of cracks in cohesive zones [START_REF] Planas | Cohesive cracks versus nonlocal models: Closing the gap[END_REF][START_REF] Li | Eigenvalue analysis of size effect for cohesive crack model[END_REF][START_REF] Zi | Eigenvalue method for computing size effect of cohesive cracks with residual stress, with application to kink-bands in composites[END_REF][START_REF] Hoover | Cohesive crack, size effect, crack band and work-of-fracture models compared to comprehensive concrete fracture tests[END_REF]. At the macroscopic scale, one obtains a smooth decrease of loading capacity in the post-peak regime [START_REF] Petersson | Crack growth and development of fracture zones in plain concrete and similar materials[END_REF][START_REF] Reinhardt | Tensile tests and failure analysis of concrete[END_REF][START_REF] Bažant | Concrete fracture models: testing and practice[END_REF]. Obviously, the basic elastic brittle model is not able to correctly describe such failure process. Therefore, a new bond damage model is here proposed in the framework of the state-based PD theory. To this end, the bond status function for any paired points x k and x j is first modified as follows:

ζ k j = f x j -x k f max x j -x k (26)
where f max x j -x k denotes the peak value of the bond force density while the bond stretch reaches the elastic limit. Compared with that defined in the bond-based PD theory, the new status function directly reflects the current degradation state of the bond force.

Furthermore, the status function is here expressed as a continuous function of the bond stretch. According to the result obtained in a concrete beam bending test [START_REF] Rots | Computational modeling of concrete fracture[END_REF], the global load-deflection response can be described by using an exponential form of the stress-crack opening relation. Inspired by this results, the exponential law has widely used in modeling tensile softening in concrete materials and structures [START_REF] Bazant | Cohesive crack with rate-dependent opening and viscoelasticity: I. mathematical model and scaling[END_REF][START_REF] Winkler | Application of a constitutive model for concrete to the analysis of a precast segmental tunnel lining[END_REF][START_REF] Grassl | Damage-plastic model for concrete failure[END_REF][START_REF] Unger | Modelling of cohesive crack growth in concrete structures with the extended finite element method[END_REF][START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF][START_REF] Le | Localised failure mechanism as the basis for constitutive modelling of geomaterials[END_REF]. Based on those previous studies, an exponential damage evolution function is here introduced to describe the local mechanical behavior of bonds for cohesive materials. The new status function is then defined as:

ζ k j =                    1, s ≤ s 0 e -k 1 s -s 0 s 0 + k 2 s -s 0 s , s > s 0 . ( 27
)
where s The bond stretch elastic limit s 0 , the bond force reduction parameter k 1 and the residual bond force parameter k 2 are now determined to completely calibrate the new bond damage model for the state-based PD theory. Similarly to the bond-based theory, the value of s 0 is also related to the macroscopic critical fracture energy G c . As depicted in Figure 3, the formation of a new crack surface requires the breakage of all bonds between the material points x k + (above the crack surface) and x j -(below the crack surface). Let K -denote the number of material points within the horizon of x k + below the crack surface and intersecting with the crack surface. And let J + represent the number of material points within the horizon of x j -above the crack surface and intersecting with the crack surface. By making use of ( 22) and ( 27), the critical energy density needed to eliminate the interaction force between x k + and x j -is given by: The last term in the bracket of (28) represents the horizon volume of the bond and can be related to the horizon radius by:

w = s 0 f ξds = 4bδξs 2 0 1 2 + 1 k 1 (1 -e -k 1 s-s 0 s 0 ) + k 2 ( s -s 0 s 0 -ln s s 0 ) + 2ad 2 δ 2 s 2 0        1 2 + 1 2k 1 (1 -e -2k 1 s-s 0 s 0 ) + 2k 2 k 1 (1 -e -k 1 s-s 0 s 0 ) + k 2 2 ( s 2 -s 2 0 ss 0 -2ln s s 0 )                K - m=1 V m + J + n=1 V n         (28) in which ξ = |x k + -x j -|.
        K - m=1 V m + J + n=1 V n         =                  πδ 3 2 , 3D case 4hδ 2 3 , 2D case . ( 29
)
Subsequently, the sum of energy used to break all the bonds crossing the unit area of a new crack surface is assumed to be equal to macroscopic fracture energy. Thus one gets for 3D cases:

G c = δ 0 2π 0 δ z cos -1 z ξ 0 wξ 2 sinφdφdξdθdz = 2πbδ 6 5 s 2 0 1 + 2 k 1 (1 -e -k 1 s-s 0 s 0 ) + 2k 2 ( s -s 0 s 0 -ln s s 0 ) + π 2 ad 2 δ 9 8 s 2 0        1 + 1 k 1 (1 -e -2k 1 s-s 0 s 0 ) + 4k 2 k 1 (1 -e -k 1 s-s 0 s 0 ) + 2k 2 2 ( s 2 -s 2 0 ss 0 -2ln s s 0 )        (30)
and for 2D cases:

G c = 2h δ 0 δ z cos -1 z ξ 0 wξdφdξdz = bhδ 5 s 2 0 1 + 2 k 1 (1 -e -k 1 s-s 0 s 0 ) + 2k 2 ( s -s 0 s 0 -ln s s 0 ) + 8ad 2 h 2 δ 7 9 s 2 0        1 + 1 k 1 (1 -e -2k 1 s-s 0 s 0 ) + 4k 2 k 1 (1 -e -k 1 s-s 0 s 0 ) + 2k 2 2 ( s 2 -s 2 0 ss 0 -2ln s s 0 )        (31)
According to ( 30) and ( 31), the bond stretch elastic limit s 0 cannot be directly determined from G c due to the unknown values of k 1 and k 2 . Thus, some simplifications are here made. It is first assumed that the energy related to the residual bond force is a very small part in the total energy. This leads to set k 2 to 0. It is further considered that the bond failure occurs when s is significantly larger than s 0 . Thus the values of

(1 -e -k 1 s-s 0 s 0 ) and (1 -e -2k 1
s-s 0 s 0 ) are approximatively equal to 1. Based on these simplifications, the value of s 0 can be calculated by:

s 0 =                          G c 2πbδ 6 5 (1 + 2 k 1 ) + π 2 ad 2 δ 9 8 (1 + 1 k 1 ) , 3D G c bhδ 5 (1 + 2 k 1 ) + 8ad 2 h 2 δ 7 9 (1 + 1 k 1 ) , 2D . ( 32
)
The above relations can be further simplified. To this end, the coefficients g 1 and g 2 are introduced as follows:

g 1 =                  2πbδ 6 5 , 3D bhδ 5 , 2D , g 2 =                  π 2 ad 2 δ 9 8 , 3D 8ad 2 h 2 δ 7 9 , 2D . ( 33
)
By taking the expressions of b and d from ( 15) and ( 16), the absolute value of the ratio g 2 g 1 is given by:

| g 2 g 1 | =                              | 27(4ν-1) 256(1-2ν) |, 3D | 8(3ν-1) 27π(1-ν) |, plane stress | 8(4ν-1) 27π(1-2ν) |, plane strain . ( 34
)
It is seen that the value of | g 2 g 1 | depends on Poisson's ratio ν, as shown in Figure 4. One can see that the value of | g 2 g 1 | remains much smaller than unit for the range of ν from 0.15 to 0.35. Based on (9), that means that the energy caused by the distortional deformation (controlled by g 1 ) is the dominating part against the energy related to the volumetric deformation (controlled by g 2 ). Based on this fact, it is possible to take

(1 + 1 k 1 ) (1 + 2 k 1 )
. Accordingly, the value of s 0 can now be calculated by:

s 0 =                          G c ( 2πbδ 6 5 + π 2 ad 2 δ 9 8 )(1 + 2 k 1 ) , 3D G c (bhδ 5 + 8ad 2 h 2 δ 7 9 )(1 + 2 k 1 )
, 2D

. By replacing again the expressions of b and d, the elastic limit of bond stretch s 0 is expressed as a function of macroscopic elastic properties:

s 0 =                        G c (3µ + ( 3 4 ) 4 (k -5µ 3 ) (1 + 2 k 1 )δ , 3D G c 6 π µ + 16 9π 2 (k -2µ) (1 + 2 k 1 )δ , 2D . (36)
Hence, with a given value of the horizon radius δ, the value of s 0 can be calculated from the macroscopic fracture energy G c once the value of k 1 is known. The parameter k 1 mainly controls the mechanical response in the post-peak regime and can be fitted from experimental results. It is worth noticing that the value of s 0 is here determined by setting k 2 = 0. Starting from this primary value, an improved value of s 0 can be iteratively identified by putting another value of k 2 .

Adaptive coupling method

The non-local numerical method based on the state-based PD theory is suitable to dealing with progressive damage and cracking process in solid materials. For large scale problems, the computer time can be considerably high. On the other hand, for many engineering problems, cracking generally occurs inside some small zones. For example, during excavation of an underground cavity, damage and cracking are generated only in the close zone to the excavated cavity wall. There is no need to use the damage model with the PD theory in far field. In this case, it is generally more efficient to use the classical finite element method for solving elastic (even plastic without softening and localization) problems in far field. Therefore, it is very convenient to combine the non-local PD theory for dealing with cracking process and the finite element method (FEM) for dealing with classical elastic and plastic problems. Further, it is also needed to consider that the size of cracking zone progresses during loading history. The boundary between the PD and the FEM zones is then not fixe but evolves. For this purpose, an adaptive coupling algorithm is here developed. In order to achieve such an adaptive coupling problem with moving PD-FEM boundary, two issues should be addressed: when the coupling is needed and how the coupling is realized.

In the present study, we shall develop an adaptive switching strategy. It is based on the relative elongation of the bond linking two adjacent FEM nodes (exactly the same as the stretch calculation defined in ( 17)). More precisely, the proposed coupling method is illustrated in Figure 5. One considers two initial FEM nodes (the gray square nodes) x k and x j . When the relative elongation between these nodes reaches a switching threshold value [START_REF] Zaccariotto | Coupling of fem meshes with peridynamic grids[END_REF], they are converted to PD points together with all other points inside their horizon zones (the blue circle nodes). With the new bond damage model presented above, the switching threshold value is set as the bond stretch elastic limit s 0 . In this way, the switching from FEM nodes to PD points is consistent with the crack propagation process. After the switching is turned on, the computational domain is divided into three regions (as shown in Figure 5b For the effective realization of FEM-PD coupling, the static equilibrium equation of a PD point is first written as follows by substituting ( 22) for ( 5), it is:

0 = N j=1 ζ k j         4bδs k j + 2ad 2 δ 2 |x j -x k | ( N m=1 ζ km s km V m + N n=1 ζ jn s jn V n )         V j x j -x k |x j -x k | + b(x k ) (37)
Referring to the definition of stretch s of two pairwise points given in (17), and by means of the multiplier V k representing the volume of material point x k , the equilibrium equation at x k can be rewritten as: 38), it is seen that the equilibrium equation of a PD material point has a similar expression to that of FEM nodes:

N j=1 ζ k j         4bδV j V k |x j -x k | (u k -u j ) + N m=1 ζ km 2ad 2 δ 2 V m V j V k |x j -x k ||x m -x k | (u k -u m ) + N n=1 ζ jn 2ad 2 δ 2 V n V j V k |x j -x k ||x n -x j | (u j -u n )         = b(x k )V k (38) Observing (
K P k U = f k (39)
in which K P k is treated as the stiffness components of the PD material point x k , it is defined as:

K P k = N j=1 • • • ζ k j p 1 + N m=1 ζ km p 2 • • • ζ k j (-ζ km p 2 ) • • • ζ k j -p 1 + N n=1 ζ jn p 3 • • • ζ k j (-ζ jn p 3 ) • • • (40) where p 1 = 4bδV j V k |x j -x k | , p 2 = 2ad 2 δ 2 V m V j V k |x j -x k ||x m -x k | , p 3 = 2ad 2 δ 2 V n V j V k |x j -x k ||x n -x j |
. Correspondingly, U represents the node displacement vector, which is given as:

U T = • • • u k • • • u m • • • u j • • • u n • • • (41)
Finally f k = b(x k )V k denotes the external load acting on x k . Then owing to the expression in (39), the coupling can be realized by assembling the stiffness components respectively from FEM nodes and PD points into the global stiffness matrix [START_REF] Galvanetto | An effective way to couple fem meshes and peridynamics grids for the solution of static equilibrium problems[END_REF][START_REF] Zaccariotto | Coupling of fem meshes with peridynamic grids[END_REF]:

K = K F + K P (42)
where K denotes the global stiffness matrix. K P denotes the stiffness matrix from the PD points as described in (40). K F represents the stiffness matrix from FEM nodes, which is obtained by the classical assemblage of the elementary stiffness matrices K F e as follows [START_REF] Zienkiewicz | The finite element method[END_REF]:

K F e = B T DBdV ( 43
)
where B is the elementary displacement dradient matrix, and D denotes the elastic stiffness matrix for the FEM zones. It is worth noticing that the internal forces evaluated with the FEM only act on the FEM nodes, while the internal forces calculated using the PD theory are applied on the PD nodes. In other words, the internal force acting on a node is of the same nature as the node [START_REF] Zaccariotto | Coupling of fem meshes with peridynamic grids[END_REF]. Subsequently, the total system of equilibrium equations can be expressed as follows [START_REF] Galvanetto | An effective way to couple fem meshes and peridynamics grids for the solution of static equilibrium problems[END_REF][START_REF] Zaccariotto | Coupling of fem meshes with peridynamic grids[END_REF]:

K U = F (44)
F represents the total nodal force vector, which is assembled from the elementary nodal force vectors F e given in the classical FEM framework [START_REF] Zienkiewicz | The finite element method[END_REF]:

F e = N T p v dV + N T p s dS ( 45
)
where N is the matrix of shape functions, p v denotes the body force vector, p s the surface traction on the external boundary. With this system in hand, the nodal displacements can be determined by:

U = K -1 F (46)
Finally, the numerical method for modeling the progressive failure process in cohesive materials by using the new bond damage model and the adaptive FEM-PD coupling is established and illustrated in Algorithm 1.

Algorithm 1: Flowchart of the adaptive PD-FEM coupling algorithm

Input: E, ν, G c , δ, ∆load, k 1 , k 2 , s 0 Output: U n , f n , ϕ n Initialize U n = U n-1 (U 1 = 0);

Numerical assessment and experimental validation

The proposed numerical method is implemented with the MATLAB software. In this section, a few linear elastic examples with a fixed PD region are first considered to assess the accuracy of the proposed coupling method. Then several typical experimental tests on concrete structures are investigated to verify the efficiency of the new bond damage model implemented with the adaptive PD-FEM coupling method for modeling the progressive failure process in cohesive materials. In those cases, the PD region is no more fixed a priori but depends on the evolution of cracking process. All the numerical calculations are performed in 2D plane stress conditions using uniform meshes composed of quadrilateral elements. The unit of all length variables is millimeter (mm).

Elastic response verification

We consider first two linear elastic examples to verify the accuracy of the proposed FEM-PD coupling method. The first example is a uni-dimensional bar subjected to tensile force. As shown in Figure 6, the bar with a length of 50mm is constrained at the left end (which is set to the origin of the coordinates frame), and is stretched by a unit force F(F = 1N) at the right end. The cross section area is A = 1mm 2 . For the mesh, the bar is divided by a uniform grid space ∆. The selected value of Young's nodulus is E = 1GPa. For the FEM-PD coupling, the region of 20 ≤ x ≤ 30 (the area filled with pink color) is chosen for the PD modeling.

And the horizon radius δ is specified as δ = 3∆, as widely used in previous studies [START_REF] Zaccariotto | An enhanced coupling of pd grids to fe meshes[END_REF][START_REF] Madenci | Ordinary state-based peridynamics for thermoviscoelastic deformation[END_REF][START_REF] Gao | Ordinary state-based peridynamic modelling for fully coupled thermoelastic problems[END_REF]. In order to explore the robustness sensitivity of the adopted coupling method with respect to the grid size, three cases respectively with ∆ = 2mm, ∆ = 1mm and ∆ = 0.5mm are considered. For each case, the solutions obtained by the coupled FEM-PD method and from the pure FEM modeling are compared with the analytical ones. From Figure 7, it can be observed that the displacements obtained by the FEM calculation and by the coupled method perfectly coincide with the analytical solutions. It seems to confirm the effectiveness of the developed state-based PD method as well as the FEM-PD coupling strategy. Moreover, the numerical results are not dependent on the grid size. Poisson's ratio ν = 0.2. The two plates are uniformly meshed by the quadrilateral elements with a size of ∆x = ∆y = ∆. For the plate without notch, the filled region with a range of 20 ≤ x ≤ 30, 0 ≤ y ≤ 20 is set as the PD region. And for the notched plate, the PD nodes are in the region of 20 ≤ x ≤ 30, 0 ≤ y ≤ 15. The horizon radius δ is also taken as three times of the grid length ∆. Again, the different grid size are chosen as ∆ = 1mm, ∆ = 0.5mm and ∆ = 0.25mm. The displacement solutions given by pure FEM calculation and the coupling FEM-PD method are compared along the geometric centerlines (line 1 and line 2 for the plate without notch, and along the lines line 1 and line 2 for the notched plate). The obtained results are present in Figure 9 and Figure 10 respectively. It can be seen that the results based on the coupling method are well consistent with those by the FEM calculation for all the cases. However, it is undeniable that the grid size does have a certain impact on the results. In general, the relative error between the coupling method and the FEM calculation decreases with the decrease of the grid size. For the two cases studied here, the relative error between the two calculations is less than 2%. On the other hand, the maximum relative error between 305 three grid sizes for the results obtained by the FEM-PD coupling method is about 3%.

Wedge splitting test 307

We shall now assess the ability of the proposed FEM-PD coupling method for modeling the progressive 308 failure process.

309

The first example considered is the wedge splitting test of a concrete specimen CP250 reported in 310 [START_REF] Trunk | Einfluss der bauteilgrösse auf die bruchenergie von beton[END_REF], and widely used to investigate mode-I fracture process. The geometric and boundary con-311 ditions of the test are presented in Figure 11. The notched specimen is subjected to the symmetrically 312 prescribed displacement (u) at the middle point of the upper notch. The material properties are taken from 313 the previous study [START_REF] Trunk | Einfluss der bauteilgrösse auf die bruchenergie von beton[END_REF]): Young's modulus E = 28.3GPa, Poisson's ratio ν = 0.18 and fracture 314 energy G c = 0.3N/mm. The grid size of the adopted mesh is ∆ = 5mm. And the horizon radius δ is equal to

Sensitivity analysis of bond damage parameters

The sensitivity of numerical results to the key parameters involved in the new bond damage model is first evaluated. As described in Section 3, there are three key parameters introduced in the proposed bond damage model, which are k 1 , k 2 and s 0 . s 0 can be determined from k 1 and k 2 according to (30) and (31). In consequence, the parameters that play the major role are k 1 and k 2 . Numerical calculations are then performed with different values of those parameters in order to verify their impacts on the obtained numerical solutions.

The influence of k 1 is first studied by taking k 1 = 0.02, k 1 = 0.01 and k 1 = 0.005, while keeping k 2 to 0.

From Figure 12, one can see that k 1 not only affects the peak value of structure strength but also controls the reduction rate of force after the peak value. The peak strength decreases when the value of k 1 is smaller.

This is due to the fact that the bond stretch elastic limit s 0 decreases when k 1 is smaller as predicted by (36). As a consequence, the cracking process starts earlier. On the other, the parameter k 1 does not affect the residual strength. Accordingly, when its value is higher, as the peak strength is higher, one gets a more rapid drops of force after the peak strength.

Three different values of k 2 (k 2 =0, k 2 = 0.02, k 2 = 0.05) are also considered by taking k 1 = 0.005. The obtained force-displacement curves are depicted in Figure 13. One observes a quasi uniform increase of both the peak and residual strengths when the value of k 2 increases. It does not affect the reduction rate of force in the post-peak regime. According to ( 26) and ( 27), the bond damage rate is lowered by a higher value of k 2 . This leads to the increase of the peak and residual strengths.

As mentioned above, the value of s 0 is theoretically calculated by suing (36). However, several assumptions are used to reach the simplified form of (36). Therefore, it appears interesting to evaluate the sensitivity of numerical results to any perturbation of s 0 . For this purpose, the reference value of s 0 is first calculated by (36). Then, two other disturbed values are selected as s 01 = 1.10s 0 and s 02 = 1.05s 0 . The values of k 1 and k 2 are hold constant in three cases (k 1 = 0.005, k 2 =0). The obtained results are shown in Figure 14. One can see that the influence of s 0 is a kind of mixture of those of k 1 and k 2 . The values of s 0 341 affects significantly the peak strength and slightly the residual strength.

Sensitivity analysis of loading increment size

Crack propagation is a strongly nonlinear problem. The whole loading history is divided into a number of loading steps. Increments of forces and displacements are prescribed at each loading step. Due to the strong nonlinearity, numerical solutions can be significantly influenced by the size of load increment. For this purpose, the influence of the prescribed displacement increment size (∆u) is here evaluated by taking four different values. As presented in Figure 15, the size of (∆u) has a significant impact on the overall force-displacement responses. Large sizes of ∆u, such as ∆u = 5 • 10 -3 mm, may delay and under-estimate cracking process and generates numerical oscillations. However, the numerical results are convergent when the size of ∆u is small enough. For instance, the numerical results seem to become stable when the size of ∆u is about 5 • 10 -4 mm or less.

Experimental verification

The numerical results are now compared with the experimental data reported in [START_REF] Trunk | Einfluss der bauteilgrösse auf die bruchenergie von beton[END_REF] in terms of overall force-displacement curves. For the parameters involved in the bond damage model, the value of k 1 is adjusted to 0.005. Two values of k 2 are chosen as 0, 0.02. The adopted value of s 0 is 10% higher than the calculated one by ( 36). This is consistent with the fact that the value of s 0 is slightly under-estimated by The crack propagation process in the wedge splitting test is also investigated. The overall damage distributions are calculated and presented in Figures 17 and 18 at four subsequent loading steps (u = 0.1mm, u = 0.2mm, u = 0.4mm, u = 0.8mm), and respectively for two values of k 2 (k 2 = 0, k 2 = 0.02). One can see that in both cases the cracking emerges and propagates vertically along the symmetric line during the loading history. This is consistent with the experimental observations reported in [START_REF] Trunk | Einfluss der bauteilgrösse auf die bruchenergie von beton[END_REF]. Comparing the two cases, it is concluded that the residual strength parameter k 2 hardly affects the crack propagation pattern, but slightly reduces the maximum damage value.

L-shape test

The second case is a mixed-mode test of an L-shaped concrete structure illustrated in Figure 19. This test was carried out in [START_REF] Winkler | Experimental verification of a constitutive model for concrete cracking[END_REF]. As set in the experiment, a vertical displacement is applied to the point at a distance of 30mm from the right edge while the fixed constraints are prescribed on the bottom edge. Based on previous studies [START_REF] Winkler | Experimental verification of a constitutive model for concrete cracking[END_REF][START_REF] Le | Localised failure mechanism as the basis for constitutive modelling of geomaterials[END_REF], the mechanical parameters chosen are Young's modulus E = 25.85GPa, Poisson's ratio ν = 0.18 and critical fracture energy G c = 0.065N/mm.

For the mesh, the grid size is set as ∆ = 5mm. The horizon radius δ is again equal to 3∆.

In the numerical calculation, the size of applied displacement increment is chosen as 10 -3 mm. The damage parameter k 1 is taken as 0.01. The value of s 0 is adjusted to 1.5 times the value given by (36). The The crack propagation patterns are illustrated in Figure 21 in terms of the overall damage contour at different loading steps for the case k 2 = 0.02. As it can be found, the crack initiates from the corner and then curvedly propagates toward the left side. This coincides well with the experimentally observed results in [START_REF] Winkler | Experimental verification of a constitutive model for concrete cracking[END_REF], particularly the curved cracking path. Moreover, in Figure 22, one presents the cracking paths predicted by the bond and state based peridynamics theories and that observed in the experiment. On the whole, both numerical predictions are consistent with the experimental observation.

But the smoothly curvilinear cracking path is better captured by the state-based peridynamics theory. It is noticed that by using the adaptive coupling method proposed here, the boundary between the PD and FEM zones is not fixed and evolves during the cracking process. For instance, for this example, at the initial state, all the domain belongs to the FEM nodes. With the nucleation and propagation of induced cracks, the PD zone is introduced. In Figure 23, one can see the current PD points distributions at u=1mm, which are surrounded by the FEM nodes.

Tension-shear test

Finally, the combined tension-shear test performed in [START_REF] Nooru-Mohamed | Experimental and numerical study on the behavior of concrete subjected to biaxial tension and shear[END_REF] is investigated.

The geometric parameters and boundary conditions of this test are shown in Figure 24. The double-edge notched specimen is under combined tension and shear. The specimen is firstly subjected to the lateral force P s = 5kN along the left edge above the notch, and then subjected to the axial tensile force at the top edge while P s is maintained unchanged. The tensile force is applied in the form of displacement in order to capture the post-peak response. Based on previous studies [START_REF] Wu | Extended embedded finite elements with continuous displacement jumps for the modeling of localized failure in solids[END_REF][START_REF] Zhang | Strong discontinuity embedded approach with standard sos formulation: Element formulation, energy-based crack-tracking strategy, and validations[END_REF], the mechanical parameters are taken as: Young's modulus E = 32GPa, Poisson's ratio ν = 0.2 and critical fracture energy G c = 0.11N/mm. The grid size of adopted mesh is ∆ = 1.25mm and the value of PD horizon radius δ is again chosen as 3∆.

After preliminary calculations, the prescribed displacement increment is determined as 2 • 10 -4 mm.

Regarding the parameters in the new bond damage model, k 1 has a value of 0.002, two values (0 and 0.01) are assigned to k 2 , and s 0 is 20% higher than the value calculated with (36).

The predicted force-displacement curves are presented in In addition, it seems that the state-based PD theory is more suitable to described the curvilinear cracking path than the bond-based PD theory due to the simultaneous consideration of shear and volumetric strains.

As future prospectives, the proposed method will be applied to numerical analysis of cracking process in large scale structures in order to highlight the efficiency of adaptive FEM-PD coupling strategy.

Figure 1 :

 1 Figure 1: Illustration of the ordinary stated-based peridynamics theory

Figure 2 :

 2 Figure 2: Evolutions of bond force and status function respectively in the PMB model (a) and new cohesive bond damage model (b)

Figure 3 :

 3 Figure 3: Sketch map of the interaction between pairwise points crossing a crack surface

Figure 4 :

 4 Figure 4: Variation of | g 2 g 1 | with ν

  ): the original FEM region, referring to the region composed of the quadrilateral elements with four gray square FEM nodes connected by the gray straight solid lines; the pure PD region, referring to the region including the blue circle PD points interacting with other PD points through the bonds indicated by the magenta curved solid lines; and the coupling region, referring to the region consisting of the nominal quadrilateral elements with FEM nodes and PD points connected by the gray straight dotted lines and solid lines, and the nominal bonds interacting PD points and FEM nodes indicated by the magenta curved dotted lines.

Figure 5 :

 5 Figure 5: Diagram of adaptive coupling strategy: gray square nodes are FEM nodes and blue circle nodes are PD points; the gray straight solid lines represent the edges of the quadrilateral element in FEM, while the gray straight dotted lines represent the edges of the nominal elements in coupling region; the magenta curved solid lines indicate the PD bonds (the short magenta lines represent the bonds omitted for brevity), and the magenta curved dotted lines indicate the nominal bonds in coupling region.

Figure 7 :

 7 Figure 6: Geometry and boundary conditions of the 1D bar problem

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: Geometry and boundary conditions of the 2D plate problems

Figure 11 :

 11 Figure 11: Geometrical and boundary conditions of the wedge splitting test

Figure 12 :Figure 13 :Figure 14 :

 121314 Figure 12: Sensitivity analysis of k 1 with k 2 = 0

Figure 15 :

 15 Figure 15: Convergence study of overall response for wedge splitting test with respect to displacement increment size ∆u (k 1 = 0.005, k 2 = 0)

( 36 )Figure 16 :

 3616 Figure16: Comparisons of force-displacement curves of wedge splitting test between different bond damage models and experimental data ([START_REF] Trunk | Einfluss der bauteilgrösse auf die bruchenergie von beton[END_REF]

Figure 17 :Figure 18 :

 1718 Figure 17: Global damage value contours in wedge splitting test (k 2 = 0)

  obtained force-displacement curves for k 2 = 0 and k 2 = 0.02 are shown in Figure 20a. In an overall manner, the numerical results agree well with the experimental data, especially for the response in the post-peak region. Again, the results based on the PMB model are presented in Figure 20b. Once more, large scatters with experimental data are observed. The PMB model predicts a too high peak strength and a too stiff (or brittle) behavior in the post-peak region. The advantage of the new bond damage model in modelling the 389 progressive failure of cohesive materials is again clearly demonstrated.

Figure 20 :Figure 21 :

 2021 Figure 19: Sketch map of the L-shape test

Figure 22 :Figure 23 :

 2223 Figure 22: Comparisons of cracking path between bond-based PD, state-based PD and experiment observation ((Winkler et al., 2001))

  Figure 25. Both the peak strength value and the post-peak response are in good concordance with the experimental data. The use of the residual force parameter (k 2 = 0.01) improves the numerical prediction. The results obtained by the classical PMB model are once more significantly in disagreement with the experimental data.In Figure26, one illustrates the crack propagation modes by the overall damage contour at different loading steps for the case k 2 = 0.01. The experimental observations reported in[START_REF] Nooru-Mohamed | Experimental and numerical study on the behavior of concrete subjected to biaxial tension and shear[END_REF] are well reproduced. The crack initiates at the corner of the notch and propagates along a curvilinear path owing to the combined effect of lateral compression induced shear and axial tension.

Figure 24 :Figure 26 :

 2426 Figure 24: Geometrical parameters and boundary conditions of the tension-shear test

  Solve the system of equilibrium equations and obtain the node displacements U n with (46); Update the PD bond stretch value s k j with (17); Update the PD bond status value ζ k j with (27); Calculate the global damage value ϕ n of PD points with (25).

	Calculate the relative elongation s of any two adjacent FEM nodes with (17);
	if s < s 0 then
	Assemble the global stiffness matrix with (42) but K P = 0 ;
	else
	Update the PD nodes set based on the adaptive switching strategy;
	Calculate the bond stretch value s k j with (17);
	Calculate the bond status value ζ k j with (27);
	Assemble the global stiffness matrix with (40), (42) and (43);
	end

Acknowledgments

This work has been partially supported by the I-site program of University of Lille through the project Multifrac.