L Farge

Sandrine Hoppe

V Daujat

F Tournilhac

Stéphane André

Solid rheological properties of PBT-based vitrimers

published or not. The documents may come

Introduction

Vitrimers were introduced in the beginning of the 2010s and are chemically crosslinked polymers characterized by the presence of exchangeable links. [START_REF] Montarnal | Silica-like malleable materials from permanent organic networks[END_REF][START_REF] Capelot | Metal-catalyzed transesterification for healing and assembling of thermosets[END_REF][START_REF] Capelot | Catalytic control of the vitrimer glass transition[END_REF] This feature allows topology changes while maintaining the integrity of the network (see Denissen et al. and Van Zee et al. for review articles [START_REF] Denissen | Vitrimers: permanent organic networks with glass-like fluidity[END_REF][START_REF] Van Zee | Vitrimers: Permanently crosslinked polymers with dynamic network topology[END_REF] ). The bond exchange is thermally activated and takes places according to an "associative" scheme corresponding to a substitution reaction between a pendent reactive group within the network and an existing cross-link. [START_REF] Kloxin | Covalent adaptable networks: smart, reconfigurable and responsive network systems[END_REF][START_REF] Scheutz | Adaptable crosslinks in polymeric materials: resolving the intersection of thermoplastics and thermosets[END_REF] This exchange mechanism implies that the number of covalent bonds (or cross-links) remains constant after the synthesis step. The general concept of vitrimers implies exchange reactions are at work but this is in no way bound to any particular chemistry or even to a particular type of network. Prototypes were based on bond exchange through metal catalyzed transesterifications and synthesized from multifunctional monomers, with relevant functions to ensure gelation and establish exchangeable links. Since then a number of catalyzed and catalyst-free exchange reactions have been tested. [START_REF] Denissen | Vinylogous urethane vitrimers[END_REF][START_REF] Arnebold | Network dynamics in cationically polymerized, crosslinked epoxy resins and its influence on crystallinity and toughness[END_REF][START_REF] Azcune | Aromatic disulfide crosslinks in polymer systems: Self-healing, reprocessability, recyclability and more[END_REF][START_REF] Denissen | Chemical control of the viscoelastic properties of vinylogous urethane vitrimers[END_REF][START_REF] Zhang | Effect of sterics and degree of cross-linking on the mechanical properties of dynamic poly (alkylurea-urethane) networks[END_REF][START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF][START_REF] Hendriks | Poly (thioether) vitrimers via transalkylation of trialkylsulfonium salts[END_REF][START_REF] Tang | Strong and efficient self-healing adhesives based on dynamic quaternization cross-links[END_REF][START_REF] Obadia | Tuning the Viscosity Profile of Ionic Vitrimers Incorporating 1, 2, 3-Triazolium Cross-Links[END_REF][START_REF] Stukenbroeker | Polydimethylsiloxane quenchable vitrimers[END_REF][START_REF] Zhang | Dynamic polysulfide shape memory networks derived from elemental sulfur and their dual thermo-/photo-induced solid-state plasticity[END_REF][START_REF] Snyder | Reprocessable aciddegradable polycarbonate vitrimers[END_REF][START_REF] Lessard | Catalyst-free vitrimers from vinyl polymers[END_REF][START_REF] Self | Brønsted-Acid-Catalyzed Exchange in Polyester Dynamic Covalent Networks[END_REF][START_REF] Spiesschaert | Filler reinforced polydimethylsiloxane-based vitrimers[END_REF] Furthermore, new vitrimers have been synthesized from polymer precursors by introducing exchangeable chemical cross-links between their chains. [START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF][START_REF] Lu | Making insoluble polymer networks malleable via olefin metathesis[END_REF][START_REF] Brutman | Polylactide vitrimers[END_REF][START_REF] Imbernon | Stress relaxation and self-adhesion of rubbers with exchangeable links[END_REF][START_REF] Röttger | Highperformance vitrimers from commodity thermoplastics through dioxaborolane metathesis[END_REF][START_REF] Caffy | Transformation of polyethylene into a vitrimer by nitroxide radical coupling of a bis-dioxaborolane[END_REF] This second route is particularly appealing because it enables the transformation of commercial thermoplastics into new materials with minimal engineering effort and additionally so because this kind of modification can preserve crystallinity. [START_REF] Leibler | Processable semi-crystalline polymer networks European patent EP[END_REF] The presence of dynamic (as opposed to static) covalent networks whose characteristics depend on the synthesis process parameters gives vitrimers specific and potentially tunable properties [START_REF] Zhao | Catechol-Mediated Glycidylation toward Epoxy Vitrimers/Polymers with Tunable Properties[END_REF][START_REF] Zhou | Tuning PBT vitrimer properties by controlling the dynamics of the adaptable network[END_REF][START_REF] Erice | Effect of Regioisomerism on Processability and Mechanical Properties of Amine/Urea Exchange Based Poly (ureaurethane) Vitrimers[END_REF] . This is the main reason why these new materials are considered to be viable solutions to overcome certain weaknesses of existing polymers.

When vitrimers are obtained with thermoset precursors at a temperature high enough to allow the activation of the covalent bond migrations, this makes topological modifications of the polymer network possible thus rendering the material malleable. This in turn means it can be shaped and paves the way for recycling processes. [START_REF] Montarnal | Silica-like malleable materials from permanent organic networks[END_REF][START_REF] Wang | Facile in situ preparation of highperformance epoxy vitrimer from renewable resources and its application in nondestructive recyclable carbon fiber composite[END_REF] In the case of thermoplastic-based vitrimers, if the exchange reaction rate is low at service temperature this gives an enhanced dimensional stability due to crosslinking. Conversely, at processing temperatures a higher rate of exchange reactions means the material remains easy to shape. [START_REF] Lu | Making insoluble polymer networks malleable via olefin metathesis[END_REF][START_REF] Zhou | Tuning PBT vitrimer properties by controlling the dynamics of the adaptable network[END_REF] With the exception of shape-memory vitrimer formulations which include a volume fraction of semi-crystalline phase, [START_REF] Arnebold | Network dynamics in cationically polymerized, crosslinked epoxy resins and its influence on crystallinity and toughness[END_REF][START_REF] Zhao | Shape memory polymer network with thermally distinct elasticity and plasticity[END_REF][START_REF] Zheng | Thermoset Shape-Memory Polyurethane with Intrinsic Plasticity Enabled by Transcarbamoylation[END_REF][START_REF] Pei | Regional shape control of strategically assembled multishape memory vitrimers[END_REF][START_REF] Arnebold | Fast switchable, epoxy based shape-memory polymers with high strength and toughness[END_REF] the semi-crystalline vitrimers that have been reported so far have been based on two types of semicrystalline polymer (SCP) precursors: 1) PBT (Poly(butylene terephthalate)) [START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF][START_REF] Zhou | Tuning PBT vitrimer properties by controlling the dynamics of the adaptable network[END_REF][START_REF] Zhou | Poly (butylene terephthalate)/glycerol-based vitrimers via solid-state polymerization[END_REF] , the polymer we used in the present study, and 2) PE (Polyethylene) [START_REF] Röttger | Highperformance vitrimers from commodity thermoplastics through dioxaborolane metathesis[END_REF][START_REF] Caffy | Transformation of polyethylene into a vitrimer by nitroxide radical coupling of a bis-dioxaborolane[END_REF][START_REF] Tellers | Reprocessable vinylogous urethane cross-linked polyethylene via reactive extrusion[END_REF][START_REF] Zych | Polyethylene vitrimers via silyl ether exchange reaction[END_REF] . The influence of the dynamic permanent covalent network on the rheological properties of SCP-based vitrimers has mainly been studied in the "molten state", i.e. with measurements obtained above the crystal melting temperature ( m T ). With DMA (Dynamic Mechanical Analysis) measurements carried out with a temperature sweep, the storage modulus has been found to exhibit a rubbery plateau above m T . [START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF][START_REF] Röttger | Highperformance vitrimers from commodity thermoplastics through dioxaborolane metathesis[END_REF][START_REF] Caffy | Transformation of polyethylene into a vitrimer by nitroxide radical coupling of a bis-dioxaborolane[END_REF][START_REF] Zhou | Tuning PBT vitrimer properties by controlling the dynamics of the adaptable network[END_REF][START_REF] Zhou | Poly (butylene terephthalate)/glycerol-based vitrimers via solid-state polymerization[END_REF][START_REF] Tellers | Reprocessable vinylogous urethane cross-linked polyethylene via reactive extrusion[END_REF] The value of the storage modulus on the plateau is mainly governed by the content of cross-linkers added during the vitrimer synthesis. [START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF][START_REF] Caffy | Transformation of polyethylene into a vitrimer by nitroxide radical coupling of a bis-dioxaborolane[END_REF][START_REF] Zhou | Tuning PBT vitrimer properties by controlling the dynamics of the adaptable network[END_REF][START_REF] Zhou | Poly (butylene terephthalate)/glycerol-based vitrimers via solid-state polymerization[END_REF] By contrast, in the case of a neat semi-crystalline polymer, the storage modulus measured by DMA drastically collapses above m T due to the "liquid-like" behaviour of the material in the molten state. The dynamic character of the permanent covalent network can be demonstrated using relaxation experiments which again need to be carried out above m T . At such temperatures, the covalent bond exchange process is likely to be activated, which makes full stress relaxation possible even though the material is crosslinked. [START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF][START_REF] Caffy | Transformation of polyethylene into a vitrimer by nitroxide radical coupling of a bis-dioxaborolane[END_REF][START_REF] Zhou | Poly (butylene terephthalate)/glycerol-based vitrimers via solid-state polymerization[END_REF][START_REF] Tellers | Reprocessable vinylogous urethane cross-linked polyethylene via reactive extrusion[END_REF] Furthermore, it is possible to extract the relaxation times [START_REF] Zhou | Tuning PBT vitrimer properties by controlling the dynamics of the adaptable network[END_REF][START_REF] Zhou | Poly (butylene terephthalate)/glycerol-based vitrimers via solid-state polymerization[END_REF] or zero shear viscosities [START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF] from the stress relaxation curves and to check that these quantities follow Arrhenius laws which proves that the flow is mainly controlled by the covalent bond exchange kinetics. [START_REF] Denissen | Vitrimers: permanent organic networks with glass-like fluidity[END_REF] The rheological behaviour of SCP-based vitrimers in the "molten state" can also be analyzed from the variations of the storage ( ' G ) and loss ( '' G ) moduli measured by rheometry with temperature or frequency sweeps. [START_REF] Kloxin | Covalent adaptable networks: smart, reconfigurable and responsive network systems[END_REF][START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF][START_REF] Caffy | Transformation of polyethylene into a vitrimer by nitroxide radical coupling of a bis-dioxaborolane[END_REF][START_REF] Zhou | Tuning PBT vitrimer properties by controlling the dynamics of the adaptable network[END_REF][START_REF] Zhou | Poly (butylene terephthalate)/glycerol-based vitrimers via solid-state polymerization[END_REF][START_REF] Chabert | Multiple welding of long fiber epoxy vitrimer composites[END_REF][START_REF] Snijkers | Curing and viscoelasticity of vitrimers[END_REF] In particular, it is possible to roughly estimate whether a vitrimer has predominantly "rubber" ( ' '' GG  ) or "viscous" behaviour ( '' ' GG  ) for a given temperature and time scale. During vitrimer synthesis, the measurement of the temporal increase of ' G and '' G also makes it possible to monitor how the crosslinking reaction progresses and find the gelation point, which occurs approximately when ' '' GG  . [START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF][START_REF] Chabert | Multiple welding of long fiber epoxy vitrimer composites[END_REF][START_REF] Snijkers | Curing and viscoelasticity of vitrimers[END_REF] After the gelation point, the crosslinking reaction ends and the storage and loss moduli reach constant values with ' '' GG  . The molten polymer remains then in the "rubbery" state except when destruction of cross-links takes place by thermal degradation. [START_REF] Caffy | Transformation of polyethylene into a vitrimer by nitroxide radical coupling of a bis-dioxaborolane[END_REF] At service temperatures, i.e. in the solid state or below m T , the published results concern mainly the viscoelastic characteristics of vitrimers, or more specifically the creep properties of these materials. Röttger et al. and Caffy et al. showed that PE-based vitrimers have a better creep resistance than their HDPE precursors, [START_REF] Röttger | Highperformance vitrimers from commodity thermoplastics through dioxaborolane metathesis[END_REF][START_REF] Caffy | Transformation of polyethylene into a vitrimer by nitroxide radical coupling of a bis-dioxaborolane[END_REF] which can naturally be attributed to crosslinking of these materials. This also implies that the rate of covalent bond exchange remains low at the measurement temperature because, if this were not the case, the material topology modifications of the dynamic covalent network would result in material flowing. In the case of a PBT/glycerolbased vitrimer with Zn(II) catalyst, Zhou et al. have shown that creep properties can be tuned by controlling the glycerol content and the Zn(II)/glycerol ratio. [START_REF] Zhou | Tuning PBT vitrimer properties by controlling the dynamics of the adaptable network[END_REF] The cross-link density is directly related to the former parameter and the dynamics of the exchangeable network to the latter. Regarding the mechanical behaviour of highly deformed solid vitrimers, it was shown using tensile test measurements performed at room temperature that elongation at break of PE vitrimers is smaller than that of their pristine HDPE precursor. [START_REF] Caffy | Transformation of polyethylene into a vitrimer by nitroxide radical coupling of a bis-dioxaborolane[END_REF] However, there are no published studies which describe the details of the plasticity development in semi-crystalline based vitrimers. This requires at least measuring the true strain-true stress curve which is intrinsic to the material, and performing strain field measurements at the specimen scale in order to study the characteristics of the possible formation of a neck.

To date, despite their promising properties for industrial applications, vitrimers have only been produced in relatively small quantities in chemistry labs. It is therefore a challenge to find processes that would enable the mass production of vitrimers before envisaging the dissemination of these materials. Reactive extrusion is a widely-used process in the polymer industry and has already been used to produce vitrimers from a PBT precursor [START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF] and PE-based vitrimers. [START_REF] Caffy | Transformation of polyethylene into a vitrimer by nitroxide radical coupling of a bis-dioxaborolane[END_REF][START_REF] Zych | Polyethylene vitrimers via silyl ether exchange reaction[END_REF] However the laboratory-scale extruders that were used were equipped with a recirculation channel to control the residence time and production cannot be carried out continuously as is required in industry. Furthermore, it has been shown that progressive material crosslinking causes the viscosity of the melt to increase drastically during the reaction which can lead to blockage of the extruder. [START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF] Taplan et al. demonstrated that highly crosslinked vinylogous urethanes vitrimers can be extruded in their fully cured form. [START_REF] Taplan | Fast processing of highly crosslinked, lowviscosity vitrimers[END_REF] These low viscosity vitrimers were obtained by modifying the distribution of reactive moieties and by selecting an amine functional cross-linker with the aim of increasing the chemical exchange rate. It was possible to achieve stress relaxation times below 1s at 150°C and cross-over between ' G and '' G at time scales and temperatures that are compatible with continuous processes (approximately 130°C and 1 rad/s).

Demongeot et al. [START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF] have already demonstrated that the mechanical properties of PBT-based vitrimers in the "molten state" were superior to those of their pristine PBT precursor. In this study, our main objective was to assess whether the improvement of vitrimer mechanical behaviour remains observable in the "solid state", between g T and m T and in particular in the range approaching the Heat Deflection Temperature (HDT ≈ 160°C for unreinforced PBT [START_REF] Radusch | Poly(Butylene Terephthalate) in Handbook of Thermoplastic Polymers: Homopolymers, copolymers, Blends, and Composites[END_REF] ). The paper ends with a brief perspective paragraph where we suggest a simple two-stage process allowing significant production of PBT-based vitrimers.

Experimental section 1.Vitrimer synthesis by non continuous reactive extrusion

We synthesized vitrimers with a PBT precursor using an experimental procedure developed by Demongeot et al. [START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF] The synthesis is based on a catalyzed addition and transesterification chemistry. Crosslinking occurs during the reactive mixture between PBT, a diepoxide and the Zn(II) catalyst. An antioxidant (Irganox 1010) was also added to the mixture. The PBT precursor was a commercial Crastin® grade from DuPont in pellet form (Material reference: CRASTIN FGS600F40, Melting temperature 225°C, Glass transition temperature 55°C, dried 24h at 160°C before injection). As diepoxide, we used the commercial bisphenol-A epichlorhydrin resin DER332, (Olin epoxy) with an [OH]/epoxy ratio lower than 0.02, further on referred as diglycidyl ether of bisphenol A (DGEBA).

We produced three different polymers: Pristine PBT and two vitrimers denoted 1

V and 2 V . The subscript indicates the weight percentage of crosslinking agent (DGEBA). The detail of the prepared blends is presented in We introduced the mixture in a 12 cm 3 DSM Xplore twin-screw micro-compounder under a nitrogen flux. The micro-compounder was equipped with a recirculation channel to control the residence time. We set the temperature and the rotor speed to 270°C and 100 rpm respectively. We also recorded the axial force exercised on the barrel of the compounder in order to monitor the melt viscosity changes. For 1 V (1% DGEBA) and 2 V (2% DGEBA), the force signal increased up to 5000 and 7000 N respectively while remaining lower than 2000 N for the pristine PBT. This proves that crosslinking actually did occur when the diepoxide and the catalyst were incorporated in the mixture. When the force reached its maximum value after about 5 minutes, we transferred the melt to a shooting pot and injected it into a tensile specimen shape mould at 200° C (figure 1c). The polymer was next cured for 4 minutes before we removed the specimen from the mould. The volume of the microcompounder means that only one tensile specimen can be obtained at a time (see figure 1b for the specimen dimensions) and the production of this unique specimen takes approximately half an hour.

Unlike the pristine PBT which is white, the two vitrimers each have a pale brown colour which is more pronounced for 2 V .

Mechanical testing

We studied the solid rheological properties of the polymers produced by performing tensile tests using a Bose® 3000 machine equipped with a thermal chamber and a 3kN load cell (see figure 1a). The tensile force is denoted F . The relative speed displacement of the tensile machine grips was 0.02 ms -1 . The maximum displacement of the grips is 18 mm. Dog-Bone shaped specimens were designed with a 6mm long central cubic part to control the localization and development of plastic instability exactly in the middle of specimens (figure 1b). In the absence of strain localization (necking), the strain in the specimen central part remains approximately constant.

We measured the strain on the specimen surface using an Aramis 3 D Digital Image Correlation (3D DIC) device provided by the Gom® company (see the cameras in figure 1a). To achieve image correlation, a random black speckle pattern was applied on the specimen surface with a paint airbrush. The measurement was performed through a glass window placed on the front door of the temperature chamber. The longitudinal strain is expressed in terms of logarithmic or true strain (see section 1 of the Supporting Information for more details):

0 ln l l      . 0
l is a small length element taken on the undeformed specimen along the x tensile axis and l is its current size during deformation. The strain in the specimen central cross-section is denoted M  . As previously mentioned,

M

 is the maximum strain along the specimen length. The true (or Cauchy) stress is defined as

F S   .
S is the current area of the specimen central cross-section. Working with the assumption of an isochoric deformation process, the true stress in the specimen central section can be shown to be given by (see section 1 of the Supporting Information): TC  , see section 2 of the SI). For a pristine PBT, the former condition is well known to insure that a plastic flow will occur and that the material will not exhibit brittle behaviour.

0 M FF e SS   

Rheometry

Dynamic shear measurements were carried out on a strain-controlled ARES G2 rheometer from TA Instruments at 270°C under a nitrogen flux. Parallel plate geometry was used. The pulsation was set to

1/ rad s   ( 0.16 f Hz 
) and the imposed strain to 0.5% .

Thermal analysis

The details about the experiments of Differential Scanning Calorimetry (DSC) and ThermoGravimetric Analysis (TGA) are given in section 2 and 3 of the SI.

X ray-diffraction

The description of the X ray-diffraction measurement is given in section 4 of the SI.

2 Solid rheology of PBT-based vitrimers

The semi-crystalline structure of solid PBT-based vitrimers

Crosslinking PBT using epoxy vitrimer chemistry in batch extrusion has already been described by Demongeot et al. [START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF] The reagents are depicted in Table 1. In short, the reactive functions (either -OH or COOH) present at the chain ends of PBT react with the epoxy functions of DGEBA, leading to the creation of resp. hydroxy ether and hydroxy ester links. The result of these addition reactions is a chain extension of PBT and the appearance of lateral hydroxide groups (in red in scheme 1), in addition to those already present (in minute amounts) in DGEBA (in red in Table 1).

Scheme 1

Then, the Zn(II) transesterification catalyst comes into action and promotes crosslinking through exchange reactions between main chain ester links and the newly formed hydroxide groups (scheme 2). At temperatures above 230°C, chain extension and crosslinking occur within a few minutes.

Scheme 2

Following this, the exchanges by transesterification continue, provided that they are thermally activated. The important feature, related to the characteristic associative behavior of crosslinking in vitrimers, is that new bonds form before existing ones are broken, thus allowing malleability at high temperature while maintaining connectivity at any time.

After extrusion, Demongeot et al. checked that the material was insoluble but remained heat processable and thereby showed the characteristic features of a vitrimer. [START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF] Crystallization was preserved (and even accelerated) in the vitrimer and the overall degree of crystallinity slightly decreased according to the amount of crosslinker.

According to the manufacturer, [START_REF][END_REF] here, the starting PBT has a density of  p = 1.310 g / cm [START_REF] Capelot | Catalytic control of the vitrimer glass transition[END_REF] , intermediate between that of the amorphous ( a = 1.265 to 1.268 g / cm [START_REF] Capelot | Catalytic control of the vitrimer glass transition[END_REF] ) and that of the  crystalline form ( c = 1.403 g / cm 3 ) [START_REF] Radusch | Poly(Butylene Terephthalate) in Handbook of Thermoplastic Polymers: Homopolymers, copolymers, Blends, and Composites[END_REF] . Based on density values, the volume fraction of crystal is equal to:

% 34 1 1 1 1                     c a p a c     
On the diffraction spectrum of unstretched injected tensile bars of PBT and 2 V (figures S2 and S3), the spectrum of the  form superimposed on an amorphous signal can be recognized. The integration of the Bragg peaks between 0.4 and 2.4 Å -1 compared to the area of the amorphous signal in I (q) representation gives a ratio of 0.43 for 2 V against 0.52 for the PBT, showing that despite a decrease in crystallinity of about 15% the structure remains the same after vitrimerization. The measurements by DSC on injected specimens confirm that the crystallinity is present in all cases and that the degree of crystallinity goes from approximately 34% for the PBT to 32% for 1 V and 30% for 2 V (see table S1). Concerning the TGA results (see figure S1), they are basically the same for PBT and the vitrimers, showing a single degradation step between 350 and 400 °C. Vitrimers only differ by a slightly higher char residue, consistent with their crosslinked structure. [START_REF] Van Krevelen | Some basic aspects of flame resistance of polymeric materials[END_REF] As PBT is a hardly entangled polymer, exhibiting very little melt viscosity, its rheological properties in the solid phase between g T and m T are mainly governed by crystallinity. In this temperature range, the cohesion of the uncross-linked PBT is therefore ensured by the presence of tie chains uniting the crystalline lamellae to each other and of entangled loops which add physical nodes (figure 2a). In this context, vitrimerization brings about two things -firstly chemical cross-linking of the tie chains (figure 2b) which increases the cohesion and secondly the possibility of relaxing stresses by rearrangement of the chemical bonds (figure 2c). As these two effects go in opposite directions, the net result on solid phase properties such as tensile strength or creep resistance is not easy to predict.

Plastic strain localization and strain hardening

The strain profiles   x  along the specimen length ( x -axis) are represented in figure 3 for the pristine PBT as well as for the 1 V and 2 V vitrimers. These three curves were measured when the maximum strain in the central cross-section was the same:

0.4 M  
. The experiments were performed at 160 TC  , which practically corresponds to the maximum temperature at which unreinforced PBT retains good mechanical properties. We made two distinct tests for each material to check repeatability. V , the strain is nearly constant along the specimen length.

Plasticity is associated there with a uniform development of the strain over the entire specimen. The slight strain decrease that can be observed from either side of the curves can be attributed to the initial geometry of our dogbone specimen (see figure 1b). The

2

V strain profile appears noisy. This is due to the lack of contrast between the random black pattern and the brown surface of the , the total length of the specimen (and therefore the relative displacement of the tensile machine grips) is much longer than that of the pristine PBT specimen. This means that the maximum strain corresponding to a given displacement of the grips differs greatly between the pristine PBT and 2

V .

The strain localization phenomenon can also be observed for the 1 V vitrimer (Figure 3) but it is much less pronounced than for the pristine PBT. This shows that this material has intermediate behaviour between that of the pristine PBT (strong strain localization) and In order to analyse the difference in behaviour between the three materials further, we also plotted the curves, intrinsic to each material, which relate the true stress  to the true stain M  both measured in the central cross- section of the specimen (figure 5).   ) shows a maximum and a decreasing part (see figure 4).

From the   M  curves, it is easy to obtain the evolution of the strain-hardening coefficient defined as 6). This quantity makes it possible to establish a link between the possible occurring of necking observable at the specimen scale, and the material's intrinsic properties. Assuming there is an isochoric deformation process, a simple calculation (section 1 of the SI) shows that the force curve is only likely to exhibit a maximum (yield point) if the intrinsic quantity  becomes smaller than 1. [START_REF] Vincent | The necking and cold-drawing of rigid plastics[END_REF][START_REF] Séguéla | On the Natural Draw atio of Semi-Crystalline Polymers: eview of the Mechanical, Physical and Molecular Aspects[END_REF] Figure 6 confirms that  remains larger than 1 throughout the deformation process for the In figure 5, the first part of the curves roughly corresponds to the elastic behaviour and subsequently the mechanical properties of the vitrimers and pristine PBT are strongly different. The

    ln M M d d     (figure
). This confirms that the vitrimer's solid mechanical properties can be tuned using the concentration of the crosslinking agent DGEBA.

At the microstructure level, the strain hardening of SCPs has been proven to mainly depend on chain entanglement density. [START_REF] Seguela | Critical review of the molecular topology of semicrystalline polymers: The origin and assessment of intercrystalline tie molecules and chain entanglements[END_REF][START_REF] Schrauwen | Intrinsic deformation behavior of semicrystalline polymers[END_REF][START_REF] Pawlak | Plastic deformation of crystalline polymers: the role of cavitation and crystal plasticity[END_REF][START_REF] Bartczak | Effect of chain entanglements on plastic deformation behavior of linear polyethylene[END_REF] Logically, the higher the entanglement density, the less extensible the macromolecular network and at macro-level this results in strain hardening enhancement. In addition to the "physical nodes" associated with the entanglements, a SCP vitrimer includes "chemical nodes" created during material crosslinking (see figure 2). In the case of the vitrimer, the semi-crystalline morphology is preserved; the additional presence of "chemical nodes" means the distance between the nodes is smaller and the macromolecular network of the SCPbased vitrimer is less extensible than that of its pristine SCP precursor. Consequently, strain hardening increases and necking may be suppressed.

This interpretation to explain necking suppression does not take into account the dynamic character of the vitrimer covalent networks. Temperature of 160°C is both well above give an exchange half time of about 10 min and complete equilibration after about 1h whereas stress relaxation data of an amorphous systems in the rubbery state give relaxation times of about 10 4 s. [START_REF] Montarnal | Silica-like malleable materials from permanent organic networks[END_REF] Here, temperature (160 °C) and the Zn molar concentration (about 12 mol%) are both a little higher, therefore chemical exchanges cannot be excluded, within the time scale of a tensile experiment (about 900s). However, we can see that on this time scale, the situation does not look at all like in the absence of cross-links. It therefore appears that although exchanges can take place, the permanence of the network associated with a high degree of crystallization in the semi-crystalline vitrimer gives the material sufficient connectivity to overcome the low degree of entanglement of the PBT.

Few examples of necking suppression in SCPs related to network properties have been reported in the literature.

Unwin et al. found that a branched polyethylene exhibits no necking if the tensile force is applied along the material extrusion direction. [START_REF] Unwin | Suppression of necking in polyethylene[END_REF] This was attributed to a residual anisotropy of the macromolecular network remaining after the extrusion process. Although this is not discussed in the article, the results obtained by Khonakdar et al. (see figure 7 of the cited article) show that necking can be suppressed for HDPE by chemically crosslinking the polymer with a peroxide agent. [START_REF] Khonakdar | An investigation of chemical crosslinking effect on properties of high-density polyethylene[END_REF]

Temperature dependency of vitrimers' mechanical properties

The evolution of the tensile force F at 80°C, 120°C and 160°C is represented in figure 7 for the 1 V vitrimer. As expected, the higher the temperature, the lower the force as is the case for all polymer materials. At 80°C, a maximum is again observable on the curve, which suggests that the strain localization phenomenon is accentuated when the temperature decreases. This can also be seen in the evolution of the strain-hardening coefficient which becomes smaller than 1 during the deformation process (see figure S4). Analysis of the strain (figure S8). For the 2

V vitrimer, at 120°C the strain profile looks like a broad hump which is slightly narrower than that obtained at 160°C (figure S9). ).

Figure 9 shows the true strain-true stress curves obtained for the three studied materials at 80°C (figure 9a) and 120°C (figure 9b). At 120°C as at 160°C, we can see that the V vitrimer shows the best resistance to deformation and also that the 1 V vitrimer's resistance to deformation is better than that of the pristine PBT.

For the vitrimer with the highest DGEBA concentration ( (or engineering strain 84%) at 120°C. The Reduction of elongation at break is a well-known characteristic of crosslinked materials. [START_REF] Nielsen | Cross-linking-effect on physical properties of polymers[END_REF] Recent studies have confirmed the increase in brittleness with the crosslinking degree for solid vitrimers (PE-based vitrimers [START_REF] Röttger | Highperformance vitrimers from commodity thermoplastics through dioxaborolane metathesis[END_REF][START_REF] Caffy | Transformation of polyethylene into a vitrimer by nitroxide radical coupling of a bis-dioxaborolane[END_REF][START_REF] Zych | Polyethylene vitrimers via silyl ether exchange reaction[END_REF] and polydimethylsiloxane-based vitrimers [START_REF] Stukenbroeker | Polydimethylsiloxane quenchable vitrimers[END_REF] ). Furthermore, the PBT elongation at break has already been observed to sometimes strongly decrease with temperature in a 50°C temperature interval above G T . 58 This could provide an insight as to why the elongation at break of 2 V is much smaller at 80°C than at 120°C. Rupture was not observed for the 1 V vitrimer or for the pristine PBT. The mechanical properties of the pristine PBT as well as those of , figure 10b). These experiments were performed with the same experimental setup as that shown in figure 1. A constant nominal stress was applied for 3600s and the strain in the specimen centre M  was measured through 3D DIC. Necking was never observed during the creep experiments and the strain can be considered as constant in the specimen 6mm long central part (see figure 1b). Three creep characteristics are indicated in table 3: the instantaneous strain

  0 M   ,
  0 M   ) (see table 3).
At the end of the tests shown in figure 10a, the increase in length of the specimen central part where the stresses and strains are approximately constant can be calculated from the values indicated in table 3. It is 5.8%

( exp(0.056) ) for the neat PBT and 4.5% ( exp(0.044) ) for the 2 V vitrimer. This proves that specimen dimensional stability is better for the vitrimer material. The improvement in dimensional stability becomes very significant when tests were carried out with the highest nominal stress level (

, figure 10b). In this case, for a one-hour duration test the increase in length of the specimen central part was 29% for the pristine PBT and only 10% for the 2 V vitrimer. After 4 hours (see Supporting Information, figure S10), the increase in length was 35% for the pristine PBT and 11% for the 2 V vitrimer. Overall, these results prove that the dimensional stability of a PBT-based vitrimer part can be greatly improved compared to that of a identical part made of the PBT precursor. This is largely because the vitrimer materials exhibit smaller instantaneous strains compared to the pristine PBT. The creep rates of the vitrimer materials were also found to be smaller than that of the pristine PBTabout four times smaller in the case of the 2 V vitrimer for the test carried out at the highest nominal stress level.

However, in the case of the curves of figure 10b, the three materials were subjected to significantly different true stresses during the test. It is difficult to estimate the extent to which the reduction of the vitrimer creep rates are really related to the intrinsic behaviour of these materials or whether this is simply because the true stresses are smaller for the vitrimers than for the pristine PBT. Answering this question would require creep tests to be carried out at constant true stress levels which is very challenging on the experimental level. However, this point is of little interest if we are to consider the real case of a polymer part in service subjected to a constant force which possibly results from an accidental overload. Nonetheless, it should be born in mind by researchers working on modelling material behaviour laws. The experiments reported in section 2 show that the plastic deformation mechanisms of vitrimers differ drastically from those of their pristine PBT precursor. Furthermore, our creep experiments clearly suggest the potential benefits of vitrimers in terms of dimensional stability. These promising results highlight the need for systematic and thorough studies of the mechanical properties of these new materials to be carried out. This kind of objective requires experiments to be run using several variable parameters namely DGEBA concentration, catalyst concentration, temperature and strain rate. In addition to monotonic tensile tests, a wide range of mechanical tests is required such as tensile tests including loading and unloading stages, fatigue tests, multi-axial tests, relaxation tests or creep tests at various constant stresses. Overall, this would also require a great amount of tensile specimens which cannot be produced using the non-continuous reactive extrusion process described in section 1.1. It is therefore desirable to develop an alternative procedure for the production of significant quantities of vitrimer materials. In the following perspective section, we briefly describe a process based on CRE (Continuous Reactive Extrusion) which could respond to the aforementioned requirement and could be transferred to the industrial scale.

N MPa   Pristine PBT 1 V 2 V Pristine PBT 1 V 2 V Instantaneous strain   0 M   0.

Perspective: production of PBT pre-vitrimers by Continuous Reactive Extrusion (CRE)

The production of fully cross-linked vitrimers by CRE is a very challenging task because of the high viscosity of these materials in the molten state. To our knowledge, this has only been achieved for vitrimers with specific properties which were deliberately designed to have a very high rate of chemical exchange. [START_REF] Taplan | Fast processing of highly crosslinked, lowviscosity vitrimers[END_REF] In section 3.1, we describe how CRE can be used to produce materials for which the reaction time in the extruder is reduced so that the crosslinking reaction is not finished at the extruder outlet. In the following, these materials are called "pre-vitrimers". In section 3.2, we show that the crosslinking reaction can resume by a simple curing process above m T leading to the final vitrimer. Our overall approach is rather straightforward but to our knowledge this is the first time the validity of such a simple CRE-based process to produce vitrimers has been demonstrated.

Pre-vitrimer production by Continuous Reactive Extrusion

The samples of PBT-based pre-vitrimers were prepared using a single-step reactive extrusion. Well-dried PBT The extruder used to carry out the crosslinking reaction was a modular intermeshing co-rotating twin-screw extruder TSE 24 from Thermo Scientific™ with a screw diameter, D, of 24 mm and an L/D ratio of 40 (where L is the length of the screw). The screw profile was designed to consecutively allow the melting of the PBT and the further mixing of the melted PBT with the reactants (DGEBA+catalyst+antioxidant) thanks to one or two zones of kneading blocks. The temperature of the barrel was constant at 270°C all along the portion of the screws where the reaction effectively took place. During the preliminary experiments, the screw rotation speed was varied from 50 rpm to 300 rpm and the locations of the feeding of PBT pellets and solution were modified. The objective of this optimization work was to reduce the reaction time to avoid the screw blockage caused when the melt viscosity is too high (figure S11). Figure S12 shows the final screw and temperature profiles along with the feeding locations of PBT and the solution of reactant (pre-mixed with the catalyst and antioxidant powders) determined after optimization of the operating process parameters. The reaction time in the extruder is about 200 s.

The flow rate of the PBT was 1 kg/h. Three "pre-vitrimers" were produced by CRE. The syringe pump flow rates respectively corresponded to 1%, 1.3% and 1.6% for the pre-vitrimers DGEBA mass fraction. Table 4 Operating conditions for the preparation of PBT pre-vitrimer samples.

Gelation and curing of the pre-vitrimers obtained by CRE

After solidification at room temperature, the three pre-vitrimers ( PV ) were ground to a characteristic size very similar to PBT pellets (a few millimetres) and introduced in the rheometer at 270°C. The test was launched when melting was observed. The storage ( ' G ) and loss ( '' G ) shear modulus evolutions are shown in figure 12. The three curves for ' G as well as the three curves for '' G show the same overall trends but are shifted -the higher the DGEBA concentration, the higher the storage and loss moduli. At the beginning of the test, ' G is smaller than '' G for the three pre-vitrimers, which indicates that the materials still have a predominantly "liquid-behaviour" and are therefore again likely to flow. The materials are approximately in the same state as at the extruder outlet. Within 10 40 s  '' G simultaneously become nearly constant, which strongly suggests that the crosslinking reaction leading to the final vitrimer (that could be called i V as in section 2) is approximately over. However, at the end of the test, the '' G modulus decreases slightly. The decrease of the storage and loss moduli can be attributed to polymer degradation. [START_REF] Caffy | Transformation of polyethylene into a vitrimer by nitroxide radical coupling of a bis-dioxaborolane[END_REF] The results of the DSC experiments obtained from the pre-vitrimers are summarized in the section 2 of the SI. As for the vitrimers produced with the micro-compounder (see section 2.1), the crystallinity degree is close to that of the pristine PBT.

The overall conclusion of the experiments reported in section 3 is that the CRE process can be considered as a preliminary mixing stage. The residence time in the extruder and the reactive agents input flow rate have to be short enough to avoid screw blocking but there is a possible range of variability to set a desired "pre-vitrimer".

The reported observations from the rheometric curves show that raising the temperature of the pre-vitrimer to 270°C restarts the crosslinking reaction. This is an important point for the development of an industrial process.

The pre-vitrimers could be initially produced by CRE in the form of storable and portable pellets. Next, the material could be shaped using fast injection molding processes. Vitrimer parts with optimal advancement in the crosslinking reaction could finally be obtained by curing at 270°C.

Conclusion

In the main part of the paper, we presented an experimental investigation of the plastic deformation mechanisms in vitrimers obtained from a PBT precursor. With monotonic tensile tests carried out at temperatures ranging from the glass transition to the melting point, the development of plasticity in the studied vitrimers was found to differ drastically from that of both their PBT precursor and of the other semi-crystalline polymers. The vitrimer specimens were found to be able to deform plastically without necking as the strain remains uniform over the entire specimen throughout tests up to high strain levels. The suppression of necking is a consequence of the strain-hardening increase caused by vitrimer crosslinking. With creep tests, it was also possible to demonstrate the significant improvement of the vitrimers' dimensional stabilities compared to their PBT precursor.

These results are of significant practical interest. Previous studies have shown that a part in the PBT vitrimer retains its original shape in the event of temporary overheating above the melting point (for example 30 min at 250 ° C) while the same part in unmodified PBT drips. [START_REF] Demongeot | Cross-linking of poly (butylene terephthalate) by reactive extrusion using Zn (II) epoxy-vitrimer chemistry[END_REF] Poly(butylene terephthalate) is present in moulded parts such as connector sheaths, kitchen utensils, cogwheels, household appliance housings, etc and such parts must withstand heat even in normal use conditions. Our study has shown that even below the melting point the PBT vitrimer offers enhanced dimensional stability in comparison to pristine PBT.

Our results are very promising and highlight the need for comprehensive systematic studies of vitrimers' solid mechanical properties. This objective requires a way to be found to achieve rapid mass production of this type of material. With this perspective in mind, we have proposed a synthesis process in which materials we named "previtrimer" are first produced by continuous reactive extrusion. In our feasibility study, the precursor was a PBT but the procedure principle is not specific to this polymer. At the extruder outlet, the crosslinking reaction is not over which makes it possible to avoid extruder blockage. During the second stage of the process, the crosslinking reaction is simply restarted by curing above the melting temperature until the final vitrimer is obtained. However, transferring the process to industry will require careful control of reaction kinetics. The gelation time must be sufficiently long for the "pre-vitrimers" to retain the ability to flow at the extruder outlet but, on the other hand, a gelation time which is too long may excessively lengthen the second stage of the process during which the crosslinking reaction ends. Furthermore, above the melting temperature, the bond migration mechanism in highly crosslinked vitrimers must be sufficiently fast to allow the easy shaping and recycling of the material. However, in normal use, high bond migration speeds could be detrimental to the material's dimensional stability. Overall, the reaction kinetics can be controlled by adjusting the chemical composition. Controlling the catalyst activity is of paramount importance in view of developing processable semi-crystalline-based vitrimers with fast crosslinking in the melt, compatible with the extrusion time and sufficiently low exchange rate at lower temperatures to get net improvement of mechanical properties in the range of PBT usage.

.

  is the initial area of the central cross-section. The nominal stress (or engineering stress) is given by tensile tests at temperatures between the glass transition ( 50 g TC  ) and the melting point ( 225 m

Figure 1

 1 Figure 1 Experimental set-up a) Photography showing the specimen in the tensile machine grips and the 3D DIC

Figure 2

 2 Figure 2 Schematic representation of the crystalline and amorphous parts of a neat SCP polymer a) and of a SCP-

Figure 3 1 V and 2 V 2 V

 3122 Figure 3 Strain profiles measured along the specimen central axis ( x -axis) at 160 TC  for the pristine PBT and

2 V vitrimer. The difference in behaviour between 2 V 2 V

 222 and the pristine PBT is effectively illustrated by the specimen photography and by the strain fields obtained when the strain in the specimen centre was insert in the figure3right part). In the case of the vitrimer, because the strain is uniform and nearly equal to the maximum value0.4 M  

2 V 2 V 1 Vof the 1 V

 2211 (quasi-uniform deformation). Compared to , the noise is reduced for the 3 D DIC measurements that were performed on the specimen. As previously mentioned, the colour specimen is lighter and the contrast with the black speckle pattern is better.

Figure 4 for 1 V , 2 V 2 Vto the 1 V 1 V and 2 V

 4122112 Figure 4 shows the force variations as a function of M  for

Figure 4 1 V and 2 V

 412 Figure 4 Evolution of the force signal ( F ) at 160 TC  for the pristine PBT and the 1 V and

Figure 5 1 V and 2 V

 512 Figure 5 True stress ( ) against true strain ( M  ) curves for the pristine PBT and the

1 V and 2 V

 12 vitrimers for tests performed at 160 TC  . For the pristine PBT,  becomes smaller than 1 at the yield point. Next, for 0.9 M   ,  again becomes larger than 1 which corresponds to neck stabilization.[START_REF] Séguéla | On the Natural Draw atio of Semi-Crystalline Polymers: eview of the Mechanical, Physical and Molecular Aspects[END_REF][START_REF] Farge | Use of Digital Image Correlation to study the effect of temperature on the development of plastic instabilities in a semi-crystalline polymer[END_REF][START_REF] Ye | Kinematic study of necking in a semi-crystalline polymer through 3D Digital Image Correlation[END_REF]

Figure 6 TC

 6 Figure 6 Evolution of the strain-hardening coefficient

2 V 2 V 1 V

 221 vitrimer clearly shows better resistance to deformation than the pristine PBT: at a given  level, the strain is much smaller for than for PBT (2V PBT   in figure5). Again the behaviour of is intermediate when compared to the pristine PBT and 2

50 gTC

 50  and well below 225 m TC  . At 150 TC  , studies of Zn(II) catalyzed transesterification on model molecules 54 with a [Zn]/[OH] ratio of 5%

(figure 8 ) 2 V 1 V and 2 V 1 V

 82121 confirmed this assumption -the smaller the temperature the narrower the strain peaks. The evolution of the tensile force F for and the pristine PBT are shown at 80°C, 120°C and 160°C in the supporting information (figureS6and figureS7). Both for , no yield point can be seen on the force signals at120 C  .The pristine PBT is similar to in that the higher the temperature the narrower the strain peaks measured at

Figure 7 1 V 18 Figure 8

 71188 Figure 7 Evolution of the force signal during the test for the

2

 2

2 V

 2 ), rupture was observed before the end of the tensile test for strain levels corresponding to

figure 9

 9 figure 9 True stress ( ) against the maximum strain ( M  ) curves for the pristine PBT and the 1 V and 2 V

1 V and 2 V

 12 be significantly smaller for the vitrimers than for the pristine PBT. This confirms that the materials' rigidity is increased during the vitrimerization process. During test corresponding to 10 three materials remained nearly equal (difference smaller than 1.2%).The creep rate measurement is illustrated in the insert in the lower part of figure10a. The pristine PBT's CR was found to be respectively 33% and 40% larger than those of the vitrimers. The total creep strain ( two vitrimers was 0.014, which is close to that of the pristine PBT (0.016). The differences observed between the three curves of figure10aincreased weakly during the creep phase although the creep rates are indeed significantly different. The curve shifts are mainly due to the difference in instantaneous strain values (

Figure 10 2 V

 102 Figure 10 Creep experiments for 160 TC  , 10 N

  pellets were introduced trough the feed hopper of the extruder via a volumetric feeder (K-Tron). A syringe pump (PS KDS-100) was used to add a controllable volume flow in the form of a steady flow of drops of the DGEBA/catalyst/antioxidant solution (figure 11) farther in the extruder. The mass proportions of the DGEBA/catalyst/antioxidant solutions were 10/1/1.

Figure 11

 11 Figure 11 Scheme of the process with the consecutive feedings of the PBT pellets and solution of reactants

Figure 12 6 PV

 126 Figure 12 Time dependence of the storage ' G (open symbols) and loss shear moduli '' G (filled symbols) for

6 PV 1 PV

 61 . In this time interval and for this viscosity range, polymer injection remains very easy. Next, within ). We can therefore estimate that gelation takes place in this time range. Thereafter (

table 1 .

 1

	Compounds
	Pristine PBT	V 1	V 2

Table 1 Prepared blends

 1

Table 2 .

 2 The Young modulus values (Table2, last line) are in accordance with our previous comments regarding the materials' rigidities: at a given temperature,[START_REF] Self | Brønsted-Acid-Catalyzed Exchange in Polyester Dynamic Covalent Networks[END_REF]

	V and	V at 80°C, 120°C and 160° are
	1	2

Table 2

 2 Mechanical properties of the pristine PBT as well as of theTo estimate the potential benefits of vitrimers in terms of dimensional stability, creep experiments were carried

	V and	V vitrimers. * Test ends at 18 mm
	1	2

Table 3 .

 3 Pristine PBT,

	V and	V creep characteristics.
	1	2

  Table 4 also gives the DGEBA solution mass flow rate as well as the torque and pressure measured at the die during the preparation of the three samples. In the following, these three Pre-Vitrimers are denoted

	1 PV 1.3 PV PV	Flow rate of DGEBA solution (g/min) 0.17 0.22 0.28	DGEBA mass fraction in pre-vitrimers 1% 1.3% 1.6%	Torque (Nm) 12 19 22	Pressure at the die (bars) 11 14 19
	1.6

i PV ( i = DGEBA mass fraction).

Acknowledgements

We thank Mickael Pomes-Hadda for technical assistance with TGA. F. T. acknowledges financial support from the ANR through the MATVIT project (ANR-18-CE06-0026-01). This research used resources of the Advanced Photon Source, a U.S. Department of Energy (DOE) Office of Science User Facility operated for the DOE Office of Science by Argonne National Laboratory under Contract DE-AC02-06CH11357. We thank Steven J. Weigand of the Advanced Photon Source for performing our measurements at beamline DND-CAT 5ID-D.

Associated Content

Supporting Information

Further information on true stresses and strains as well as on the strain-hardening coefficient. Additional results obtained during tensile tests. DSC, TGA and X-ray diffraction data.


J A Sanz

M Bonnet
email: bonnet@lms.polytechnique.fr

J Dominguez

Fast multipole method applied to 3D frequency domain elastodynamics 1

Keywords: Fast Multipole Method, Boundary Element Method, 3D Elastodynamics

This article is concerned with the formulation and implementation of a fast multipoleaccelerated BEM for 3-D elastodynamics in the frequency domain, based on the so-called diagonal form for the expansion of the elastodynamic fundamental solution, a multi-level strategy. As usual with the FM-BEM, the linear system of BEM equations is solved by GMRES, and the matrix is never explicitly formed. The truncation parameter in the multipole expansion is adjusted to the level, a feature known from recent published studies for the Maxwell equations. A preconditioning strategy based on the concept of sparse approximate inverse (SPAI) is presented and implemented. The proposed formulation is assessed on numerical examples involving O(10 5 ) BEM unknowns, which show in particular that, as expected, the proposed FM-BEM is much faster than the traditional BEM, and that the GMRES iteration count is significantly reduced when the SPAI preconditioner is used.

INTRODUCTION

The boundary element method (BEM) is a mesh reduction method. However, in traditional BEM implementations, the resulting dimensional advantage with respect to domain discretization methods is offset by the fully-populated nature of the governing discretized system of equations, with set-up and solution times rapidly increasing with the number N of unknown BEM degrees of freedon (DOFs). It is therefore important to formulate alternative strategies which allow to take advantage of the proven efficiency of BEMs for e.g. analyses involving unbounded media, fracture mechanics, or inverse problems when large numbers of degrees of freedom (DOFs) are considered. Approaches such as symmetric Galerkin BEMs and BEM-FEM coupling, when developed within traditional set-up and solution algorithms, do not fully address these concerns because they do not lead to lowered complexities (i.e. CPU and memory requirements grow with N at the same rate).

Fast BEMs, i.e. BEMs having a complexity lower than that of traditional BEMs, appeared around 1985. Rokhlin [START_REF] Rokhlin | Rapid solution of integral equations of classical potential theory[END_REF] presented a O(N) per iteration integral-equation approach for solving 2-D Laplace problems. The Fast Multipole Method (FMM) terminology was introduced first by Greengard [START_REF] Greengard | A fast algorithm for particle simulations[END_REF], in the context of many-particle simulations. Thereafter, the FMM led to fast multipole-accelerated boundary element methods (FM-BEMs), whose scope and capabilities have since rapidly progressed, especially in connection with application in electromagnetics [START_REF] Lu | Fast algorithm for solving hybrid integral equations[END_REF][START_REF] Lu | A multilevel algorithm for solving a boundary integral equation of wave scattering[END_REF], but also in other fields including computational mechanics. Many of these investigations are summarized in a review article by Nishimura [START_REF] Nishimura | Fast multipole accelerated boundary integral equation methods[END_REF].

This article is concerned with the formulation and implementation of a FM-BEM for 3-D elastodynamics in the frequency domain. Very few references address this particular area of application. FM-BEM formulation for elastodynamics in the time domain and the frequency domain have been proposed by Takahashi et al. [START_REF] Takahashi | A fast BIEM for three-dimensional elastodynamics in time domain[END_REF] and Fujiwara [START_REF] Fujiwara | The fast multipole method for solving integral equations of threedimensional topography and basin problems[END_REF], respectively. This work is based on the same general approach as that of [START_REF] Fujiwara | The fast multipole method for solving integral equations of threedimensional topography and basin problems[END_REF], which uses the so-called diagonal form for the expansion of the elastodynamic fundamental solution. This approach is suitable for moderate to medium frequencies. The upper limit stems from the fact that the size N becomes intractable at high frequencies, while the diagonal form breaks down at very low frequencies and must be replaced with other types of expansions, based on e.g. the Wigner-3j symbols [START_REF] Messiah | Quantum mechanics[END_REF] or the so-called stable plane wave expansion [START_REF] Darve | A fast multipole method for Maxwell equations stable at all frequencies[END_REF]. The present work differs from [START_REF] Fujiwara | The fast multipole method for solving integral equations of threedimensional topography and basin problems[END_REF] in details of the formulation that substantially impact the overall efficiency, and especially the fact that the truncation parameter in the multipole expansion should be adjusted to the subdivision level. This feature, known from previous studies for the Maxwell equations (e.g. [START_REF] Darve | The fast multipole method: numerical implementation[END_REF]), is incorporated in the present implementation but not in that of [START_REF] Fujiwara | The fast multipole method for solving integral equations of threedimensional topography and basin problems[END_REF]. This article is organized as follows. The necessary background concepts and notations are reviewed in section 2. Then, some details of the elastodynamic FM-BEM are discussed in section 3. A preconditioning strategy based on the concept of sparse approximate inverse (SPAI) is presented in section 4. Finally, the present formulation is assessed on numerical examples in section 5, and the paper closes with concluding remarks.

PRELIMINARIES

The main idea behind the FMM is to reformulate the fundamental solution (or kernel) in terms of products of functions of the integration and collocation points. Together with the use of iterative techniques for solving linear systems, this provides algorithms of complexity O(N log α N) per iteration (see [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF] for a mathematical justification), where N is the number of BEM DOFs. Such kernel reformulations hold for sets of integration and collocation points each clustered around a origin point (called pole); moreover, the separation distance between the two poles must be such that the two clusters are well-separated. Hence, many such poles must be introduced.

There are basically two kinds of FMMs: the one-level FMM and the multilevel FMM. Both are based on dividing the spatial region of interest (i.e. that containing the boundary over which the integral equation is formulated) into cubic regions called cells. Pairs of well-separated (i.e. completely disjoint) cells are then used for defining the above-mentioned clusters of integration and collocation points. This article is concerned with the multi-level approach, known to have a lower complexity than the one-level approach, where the spatial region of interest is first divided into 4 × 4 × 4 cubic cells. Each nonempty cell is then divided into 2 × 2 × 2 smaller cells. This subdivision process is repeated until some criterion (such as a preset minimum number of degrees of freedom per cell, or a minimum cell size) is reached.

Classical BEM formulation

Let Ω denote the (bounded or unbounded) region of space occupied by an elastic medium with isotropic constitutive properties defined by µ (shear modulus), ν (Poisson's ratio) and ρ (mass density). Assuming a time-harmonic motion with circular frequency ω and the absence of body forces, the classical elastodynamic BEM formulation is based on the well-known integral identity (see e.g. [START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF])

C ip (x)u p (x) = Γ U ⋆ ip (x,y)t n p (y)dΓ y - Γ T ⋆n ip (x,y)u p (y)dΓ y (x ∈ Γ) (1)
where t n p and u p , respectively, denote the Cartesian components of the displacements and tractions on the boundary Γ of Ω, C ip are the coefficients of the free term, which depend on the local geometry of Γ at x and is equal to (1/2)δ ij in the usual case where Γ is smooth at x (with δ ij denoting the Kronecker symbol). In [START_REF] Rokhlin | Rapid solution of integral equations of classical potential theory[END_REF] and thereafter, the usual Einstein convention of implicit summation over repeated indices is used. The second integral term in Eq. ( 1) is a Cauchy principal value singular integral, which can be dealt with by means of either a regularization procedure [START_REF] Mantic | A new formula for the C-matrix in the Somigliana identity[END_REF] or a specialized direct quadrature approach [START_REF] Giuggiani | A general algorithm for multidimensional Cauchy principal value integrals in the boundary element method[END_REF]. The fundamental solution of 3D frequency domain isotropic elastodynamics is given, in terms of Cartesian components of displacements U ⋆ ip and tractions T ⋆n ip , by

U ⋆ ip (x,y) = 1 4πµk
T ⋆n ip (x,y) = E jkpl ∂ ∂y l U ⋆ ij (x,y)n k (y) (2)
where k S , k P are the transversal and longitudinal (S and P) wave numbers, respectively, i = √ -1 is the imaginary unit, e abc are the components of the permutation tensor, n is the unit normal to Γ exterior to Ω, and E ijkl are the components of the isotropic elasticity tensor, i.e.:

E ijkl = µ 2ν 1 -2ν δ ij δ kl + δ ik δ jl + δ il δ jk (4)

Multipole expansion of the elastodynamic fundamental solution

The function exp(ikr)/(4πr) appearing in [START_REF] Greengard | A fast algorithm for particle simulations[END_REF], which is the fundamental solution for the Helmholtz equation with wavenumber k, admits the plane-wave expansion [START_REF] Abramowitz | Handbookof mathematical functions with formulas, graphs, and mathematical tables[END_REF]:

exp(ik |x-y|) 4π |x-y| = ik 16π 2 lim L→∞ S exp ik s, O x -x T L (s; k, O y -O x ) exp ik s, y-O y ds (5)
where S = {s ∈ R 3 , |s = 1|} denotes the unit sphere, the brackets •, • indicate the inner (dot) product in R 3 , and the transfer function T L is given by

T L (s; k, z) = 0≤l≤L (2l + 1) i l h (1) l (k |z|) P l (cos (s, z)) (6)
in terms of the spherical Hankel functions of the first kind h

(1) l and the Legendre polynomials P l . The poles O x and O y must be close to the collocation point x and the integration point y, respectively. Then, on substituting (5) and (6) into expressions (2) and (3) of the elastodynamic fundamental solution, the multipole expansion of U ⋆ ip and T ⋆n ip is readily obtained as

U ⋆ ip (x,y) = 1 16π 2 µ Q q=1 ω q i k S D ri (O x -x, s q )T L (s; k S , O y -O x )D rp (y-O y , s q ) + k 3 P k 2 S D i (O x -x, s q )T L (s; k P , O y -O x )D p (y-O y , s q ) (7)
T ⋆n ip (x,y) = -

1 16π 2 µ Q q=1 ω q E jkpl n k (s q ) l k 2 S D ri (O x -x, s q )T L (s; k S , O y -O x )D rj (y-O y , s q ) + k 4 P k 2 S D i (O x -x, s q )T L (s; k P , O y -O x )D j (y-O y , s q ) ( 8
)
where the integration on the unit sphere of Eq. ( 5) has been replaced with a Qpoint numerical quadrature, with s q and ω q denoting the quadrature points on S and the associated weights. Note that Eqs. ( 7) and ( 8) feature, through the transfer function T L , the two wavenumbers k S and k P , as expected. Moreover, the functions D ri (z, s) and D i (z, s) are defined by the formulae

D ri (z, s) def = e rti (s q ) t exp ik S s q , z D i (z, s) def = (s q ) i exp ik P s q , z 2.
3 Numerical quadrature over the unit sphere

On parameterizing unit vectors s ∈ S using spherical angular coordinates (θ, φ), a frequently chosen [START_REF] Darve | The fast multipole method: numerical implementation[END_REF][START_REF] Sylvand | La méthode multipôle rapide en électromagnétisme: performances, parallélisation, applications[END_REF] numerical quadrature scheme is based on Q = (L + 1)(2L + 1) points s q of the form

s q = sin(θ a ) cos(φ b ) sin(θ a ) sin(φ b ) cos(θ a ) ,
where L is the truncation level used in [START_REF] Takahashi | A fast BIEM for three-dimensional elastodynamics in time domain[END_REF], cos θ a (1 ≤ a ≤ L + 1) are the abscissae for the (L+1)-point Gauss-Legendre 1-D quadrature rule over [-1, 1] (with θ a such that 0 < θ a < π) and

φ b = 2πb/(2L + 1) (0 ≤ b ≤ 2L)
are uniformly-spaced abscissae on [0, 2π]. The associated weights are w q = 2πw θ a /(2L + 1), where w θ a are the Gauss-Legendre weights for the (L + 1)point one-dimensional rule. The chosen dependence of Q in L stems from the fact that the transfer function [START_REF] Takahashi | A fast BIEM for three-dimensional elastodynamics in time domain[END_REF] involves Legendre polynomials of degree up to L.

Moreover, it is known [START_REF] Darve | The fast multipole method: numerical implementation[END_REF][START_REF] Sylvand | La méthode multipôle rapide en électromagnétisme: performances, parallélisation, applications[END_REF][START_REF] Carpentieri | Combining fast multipole techniques and an approximate inverse preconditioner for large electromagnetism calculations[END_REF] that to ensure a given accuracy in [START_REF] Takahashi | A fast BIEM for three-dimensional elastodynamics in time domain[END_REF], the truncation parameter L must be adjusted to the cell size d. A often-used empirical formula for that purpose is [START_REF] Darve | The fast multipole method: numerical implementation[END_REF][START_REF] Sylvand | La méthode multipôle rapide en électromagnétisme: performances, parallélisation, applications[END_REF]

L(k) = √ 3kd + C log 10 √ 3kd + π (9)
where the constant C is determined on the basis of numerical experiments (C = 4 -8 being typical values). A similar approach, based on fine-tuning the coefficient C of another empirical truncation formula with a O((kd) 1/3 ) term replacing the logarithmic term of ( 9), is reported in [START_REF] Ohnuki | Numerical accuracy of multipole expansion for 2D MLFMA[END_REF]. Using either approach, L is roughly proportional to the cell size. In practice, the truncation parameter in ( 6) is set according to

L α = Max(5, L(k α )) (α = P, S) (10)
depending on whether k = k P or k = k S is used in the expansion. On noting that L S ≥ L P , the number of quadrature points Q is set for simplicity to Q = (L S + 1)(2L S + 1) in all cases.

ELASTODYNAMIC FMM

To exploit optimally the separation of variables x and y achieved by expansions ( 7), ( 8) and the clustering of influence terms associated with the introduction of poles O x , O y , a hierarchical oct-tree structure of elements is introduced. For that purpose, a cube of linear size d 0 containing the boundary Γ, called 'level-0 cell', is divided into eight cubes (level-1 cells), each of which is divided in the same fashion, and so on. The linear size of a level-ℓ cell is thus d(ℓ) = 2 -ℓ d 0 . A level-ℓ cell is divided into level-(ℓ+1) cells unless a subdivision-stopping criterion is satisfied. Such criterion is usually based on either a preset minimum cell size D (defined relative to the shear wavelength) or a preset minimum number E of boundary elements in the cell. Such "terminal" cells are termed leaves As they do not feature disjoint cells, and hence do not permit the use of multipole expansions, the level-0 and level-1 subdivisions are not actually used in the computations. The level-2 subdivision, containing 4×4×4 cells, is thus the highest level used in practice. A given cell may be used as a cluster of either collocation points or integration points, its centroid serving as pole O x or O y .

The truncation parameters L P , L S depend upon the cell size through equation [START_REF] Darve | A fast multipole method for Maxwell equations stable at all frequencies[END_REF]. Hence, in the multi-level framework, L P , L S are level-dependent: when applied at subdivision level ℓ, they are evaluated according to [START_REF] Darve | A fast multipole method for Maxwell equations stable at all frequencies[END_REF] with

d = d(ℓ) = 2 -ℓ d 0 .

Computation of the outgoing rays

Let C y denote a generic leaf cell (i.e. a cell which is not subdivided into smaller cells). The radiation functions associated with C y are defined by

F r (s; C y ) def = F U r -F T r F (s; C y ) def = F U -F T (11)
with

F U r (s; C y ) def = k S i Γ∩Cy D rp (y-O y , s q )t n p (y)dy F U(s; C y ) def = k P i Γ∩Cy D p (y-O y , s q )t n p (y)dy F T r (s; C y ) def = -k 2 S (s q ) l E jkpl Γ∩Cy n k (y)D rj (y-O y , s q )u p (y)dy F T (s; C y ) def = -k 2 P (s q ) l E jkpl Γ∩Cy n k (y)D j (y-O y , s q )u p (y)dy
With reference to expansions ( 7), ( 8), the radiation functions [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF] are to be evaluated at quadrature points s = s q on the unit sphere.

Then, an upward pass is performed where radiation functions [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF] are first computed for the lowest-level cells and then recursively aggregated by moving upward in the tree until level 2 (for which there are 4 × 4 × 4 cells) is reached, a phase known as the outer-to-outer translation. The translation from a level-(ℓ+1) cell to its parent level-ℓ cell consists in evaluating the radiation functions at the pole O ℓ y through (see Fig. 1)

F r (s; C ℓ y ) = C ℓ+1 y ∈C ℓ y F r (s; C ℓ+1 y ) exp ik S s, O ℓ y -O ℓ+1 y (12)
F (s; C ℓ y ) = C ℓ+1 y ∈C ℓ y F (s; C ℓ+1 y ) exp ik P s, O ℓ y -O ℓ+1 y (13)
However, the previously mentioned level-dependency through Eq. ( 9) of the truncation parameters L S and L P implies that the number and location of quadrature points on the unit sphere are also level-dependent. Application of the multipole expansions ( 7) and ( 8) at level ℓ+1 requires that radiation functions F r (s; C ℓ+1 y ), F (s; C ℓ+1 y ) be known at the level-ℓ+1 quadrature points s ℓ+1 q . A subsequent outer-to-outer translation via (12) or (13) therefore yields the values F r (s ℓ+1 q ; C ℓ y ), F (s ℓ+1 q ; C ℓ y ), whereas F r (s ℓ q ; C ℓ y ), F (s ℓ q ; C ℓ y ) are needed. Evaluating the latter from the former relies on a rapid interpolation technique described in [START_REF] Sylvand | La méthode multipôle rapide en électromagnétisme: performances, parallélisation, applications[END_REF], which takes advantage of the uniform distribution of quadrature points along φ (φ j = 2πj/K, 0 ≤ j ≤ K). This rapid interpolation takes advantage of the fact that radiation functions have a finite bandwidth in the Fourier space.

Let the grids of quadrature points for levels ℓ and ℓ+1 be of the form

(θ ℓ i ′ , φ ℓ j ′ ) 0 ≤ i ′ ≤ L ℓ S , 0 ≤ j ′ ≤ 2L ℓ S (θ ℓ+1 i , φ ℓ+1 j ) 0 ≤ i ≤ L ℓ+1 S , 0 ≤ j ≤ 2L ℓ+1 S
with the angular coordinates chosen as explained in Section 2.3. Let F (θ, φ) denote a generic function defined on the unit sphere and having a bandwidth of L ℓ+1 S , i.e. of the form

F (θ, φ) = L ℓ+1 S l=0 l m=-l F lm Y lm (θ, φ), (14)
having put L = L ℓ+1 S and where the functions Y lm (θ, φ), given by

Y lm (θ, φ) = Q m l (cos(θ)) exp(imφ), Q m l (x) = 2l + 1 4π (l-m)! (l+m)! P m l (x) (15)
in terms of the Legendre functions P m l , define the L 2 (S)-orthonormal basis of spherical harmonics. The coefficients F lm in ( 14) are therefore obtained by taking the scalar product (in the L 2 (S) sense) of F (θ, φ) with the Y lm (θ, φ). On using the numerical quadrature on S associated with level ℓ+1, one therefore has

F lm ≈ Q(ℓ+1) q=1 ω ℓ+1 q Ȳlm (θ ℓ+1 i , φ ℓ+1 j )F (θ ℓ+1 i , φ ℓ+1 j )
with indices q and (i, j) related through s q = s(θ ℓ+1 i , φ ℓ+1 j ) and the overbar indicating complex conjugation. On substituting the above expression into [START_REF] Giuggiani | A general algorithm for multidimensional Cauchy principal value integrals in the boundary element method[END_REF], the following linear relationship between the values of F at the quadrature points for levels ℓ and ℓ+1 is arrived at:

F (θ ℓ i ′ , φ ℓ j ′ ) = Q(ℓ+1) q=1 ω ℓ+1 q L l=0 l m=-l Y lm (θ ℓ i ′ , φ ℓ j ′ ) Ȳlm (θ ℓ+1 i , φ ℓ+1 j )F (θ ℓ+1 i , φ ℓ+1 j ) ( 16
)
This relation provides the basis for computing radiation functions at the levelℓ quadrature points from values at the level-(ℓ + 1) quadrature points. On inserting (15) into ( 16), the following three-step procedure is obtained for computing the F (θ ℓ i ′ , φ ℓ j ′ ) from the F (θ ℓ+1 i , φ ℓ+1 j ), where in practice F is one of the radiation functions [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF].

• Step 1: direct discrete Fourier transform. Compute the transforms F ℓ+1 im de-fined by:

F ℓ+1 im = 2L j=0 exp(-imφ ℓ+1 j )ω ℓ+1 ij F (θ ℓ+1 i , φ ℓ+1 j )
• Step 2: matrix-vector products. Compute the F ℓ i ′ m defined by

F ℓ i ′ m = L i=1 B m i ′ i F ℓ+1 im
where the entries B m i ′ i of matrix B m are given by

B m i ′ i = |m|≤l≤L Q m l (x i ′ )Q m l (x i ) = (L + 1) 2 -m 2 4(L + 1) 2 -1 × Q m L+1 (x i ′ )Q m L (x i ) x i ′ -x i - Q m L+1 (x i )Q m L (x i ′ ) x i ′ -x i having put x i = cos(θ ℓ+1 i ) and x i ′ = cos(θ ℓ i ′ )
, and with the second equality stemming from the Christoffel-Darboux formula.

• Step 3: inverse Fourier transform. The values F (θ ℓ i ′ , φ ℓ j ′ ) at the quadrature points for level ℓ are finally recovered by means of a discrete inverse Fourier transform:

F (θ ℓ i ′ , φ ℓ j ′ ) = L m=-L exp(imφ ℓ j ′ ) F ℓ i ′ m
Note that the above summation is in fact a (2L ℓ S +1)-point discrete transform, with the missing transformed coefficients

F ℓ i ′ m for L = L ℓ+1 S < |m| ≤ L ℓ S
implicitly set to zero.

The above-described scheme for transitions from a generic level ℓ+1 to the next upper level ℓ, combining evaluation of radiation functions at level-ℓ quadrature points and the translation formulae ( 12) and [START_REF] Mantic | A new formula for the C-matrix in the Somigliana identity[END_REF], is applied until level ℓ = 2 is reached.

Computation of the incoming rays

The radiation functions [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF] synthesize, for a given cell C y , the integral contributions of all boundary elements contained in that cell. Once computed for all non-empty cells C y at all levels, they have to be converted into quantities defined at collocation points. The structure of the multipole expansion [START_REF] Nishimura | Fast multipole accelerated boundary integral equation methods[END_REF] indicates that this conversion entails a transfer from a cell C y which supports element integrations to a (disjoint) cell C x which supports collocation points (sometimes referred to as a outer to inner translation).

Let I(C) denotes the interaction list of cell C, i.e. the list of cells of the same level as C that (a) are well-separated from C, i.e. have no point in common (not even a single vertex) with C, and (b) are such that their parent cell is adjacent to the parent cell of C. For a cell C of level ℓ = 2, I(C) lists all cells which are well-separated from C.

The incoming rays at level ℓ = 2 are defined, for a generic level-2 cell C (2) x , by applying the transfer functions to the level-2 radiation functions associated with all cells in I(C (2) x ) and aggregating the result, i.e.:

G r (s; C (2) x ) = C (2) y ∈I(C (2)
x )

F r (s; C (2) y )T L (s; k S , O (2) y -O (2) x ) (17)
G(s; C (2) x ) = C (2) y ∈I(C (2)
x )

F (s; C (2) y )T L (s; k P , O (2) y -O (2) x ) (18)
Then, a downward pass is performed whereby, for ℓ ≥ 2, the incoming rays G r (s; C ℓ+1 x ), G r (s; C ℓ+1 x ) associated with a level-(ℓ+1) cell C ℓ+1

x aggregate the contributions from (a) the level-(ℓ+1) cells C ℓ+1 y belonging to the interaction list of C ℓ+1

x , and (b) the incoming rays associated with the parent cell C ℓ x with its pole translated from O ℓ x to O ℓ+1

x (an operation sometimes referred to as inner to inner translation), i.e.:

G r (s; C ℓ+1 x ) = C ℓ+1 y ∈I(C ℓ+1 x ) F r (s; C ℓ+1 y )T L (s; k S , O ℓ+1 y -O ℓ+1 x ) + G r (s; C ℓ x ) exp ik S s q , O ℓ+1 x -O ℓ x (19)
G(s; C ℓ+1 x ) = C ℓ+1 y ∈I(C ℓ+1 x ) F (s; C ℓ+1 y )T L (s; k P , O ℓ+1 y -O ℓ+1 x )) + G(s; C ℓ x ) exp ik S s q , O ℓ+1 x -O ℓ x ( 20
)
Like with the recursive (ascending) evaluation of the radiation functions, application of equations ( 19) and ( 20) must be preceded with a procedure (known as anterpolation) whereby G r (s ℓ q ; C ℓ x ), G(s ℓ q ; C ℓ x ) are evaluated from G r (s ℓ+1 q ; C ℓ x ), G(s ℓ+1 q ; C ℓ x ). This step simply consists of a transposition of the interpolation scheme of section 3.1.

Final expression of BIE

Let C x denote a generic leaf cell with center O x . On performing the complete upward and downward passes described in Sections 3.1 and 3.2, the contribution to the residual of integral equation ( 1) at all collocation points in C x of boundary elements located in leaf cells not adjacent to C x are accounted for, as one has

y∈Γ\A(Cx) U ⋆ ip (x,y)t n p (y) -T ⋆n ip (x,y)u p (y) dy = 1 16π 2 µ Q q=1 ω q D ri (O x -x, s q )G r (s q ; C x ) + k 2 P k 2 S D i (O x -x, s q )G(s q ; C x ) (x ∈ C x ) (21)
where A(C) denotes the list set of cells of the same level as C that are adjacent to C (including C itself). Hence, integral equation ( 1) at all collocation points in C x is recast in the form

C ip (x)u p (x) - y∈Γ∩A(Cx) U ⋆ ip (x,y)t n p (y) -T ⋆n ip (x,y)u p (y) dy - 1 16π 2 µ Q q=1 ω q D ri (O x -x, s q )G r (s q ; C x ) + k 2 P k 2 S D i (O x -x, s q )G(s q ; C x ) = 0 (x ∈ C x ) (22)
where the near contributions, coming from integrals over elements in Γ∩A(C x ), are evaluated using the traditional numerical integration techniques developed for the BEM and include all singular integrals. The latter are here computed using the direct approach of [START_REF] Giuggiani | A general algorithm for multidimensional Cauchy principal value integrals in the boundary element method[END_REF].

Choice of FMM parameters

The accuracy and efficiency of the FMM is predicated on careful choice of two parameters. One is the truncation parameter L in the transfer function [START_REF] Takahashi | A fast BIEM for three-dimensional elastodynamics in time domain[END_REF], which in turn influences the choice of quadrature points on the unit sphere; this has been discussed in Sec. 2.3. The other is the linear size d of a leaf cell, and is influenced by the fact that the spherical Hankel functions h

l (v) featured in Eq. ( 6) behave in the large-index limit as

h (1) l (v) ∼ C v 2l + 1 v l (l ≫ v)
(where C is a constant) and therefore grow more than exponentially fast with l for l > |v|. This behavior is illustrated in Fig. 2 for several values of l. As a consequence, as discussed in [START_REF] Darve | A fast multipole method for Maxwell equations stable at all frequencies[END_REF][START_REF] Darve | The fast multipole method: numerical implementation[END_REF], the FMM algorithm should avoid any evaluation of h

(1) l (v) with l significantly larger than |v|.

In expansion [START_REF] Takahashi | A fast BIEM for three-dimensional elastodynamics in time domain[END_REF], the value v = k |z| taken by the argument of h

(1) l is such that v ≥ min k|O y -O x |
where O y and O x are centers of disjoint same-level cells. Hence, at each level, one has |O y -O x | ≥ 2d ℓ , where d ℓ is the size of a level-ℓ cell. So, in the worst case, one has

2l + 1 v ≤ 2L(k) + 1 2kd ℓ = √ 3 + 1 + C log 10 √ 3kd ℓ + π 2kd ℓ
and that upper bound should not be too large compared to unity. In particular, the size d of a leaf cell should not be a too small fraction of the wavelength, a typical range of recommended values [START_REF] Sylvand | La méthode multipôle rapide en électromagnétisme: performances, parallélisation, applications[END_REF][START_REF] Carpentieri | Combining fast multipole techniques and an approximate inverse preconditioner for large electromagnetism calculations[END_REF] being kd = O(1 -2).

SYSTEM SOLUTION

The BIE formulation ( 22) leads to a linear system of equations of the form

A near x + A far x = b ( 23
)
where the vector x ∈ C N ×1 collects all displacement and traction DOFs that are left unknown after accounting for the boundary conditions, the vector b ∈ C N ×1 results from the contributions of all known values of displacements and tractions in [START_REF] Leung | Iterative solution of large three-dimensional BEM elastostatics analyses using the GMRES technique[END_REF], and A near , A far ∈ C N ×N are the square matrices associated with the near and far contributions involving unknown DOFs. For BEM models of large size, A near is a sparse matrix while A far is almost fully populated.

The FMM is inherently based on using an iterative solution technique for the set of linear equations arising from the BEM discretization. The GMRES (Generalized Minimum RESiduals [START_REF] Saad | Gmres: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]) algorithm is the natural, and usual, choice for this purpose. It relies on the ability to evaluate, for a given solution candidate x (k) , the residual r (k) = b -A near x (k) -A far x (k) (which may be seen as the BEM-discretized version of ( 22)), a task requiring in turn the computation of matrix-vector products A near x and A far x. As usual with iterative solvers, matrix storage is not required. In practice, storing A far entails close to O(N 2 ) and is thus out of the question, whereas storing A near is a reasonable (but not mandatory) option which avoids costly reevaluations of singular and near-singular element integrals.

An important issue for algorithms such as GMRES is preconditioning. Numerical tests [START_REF] Mansur | Solution of BEM systems of equations via iterative techniques[END_REF][START_REF] Kane | Iterative solution techniques in boundary element analysis[END_REF][START_REF] Leung | Iterative solution of large three-dimensional BEM elastostatics analyses using the GMRES technique[END_REF] establish that the rate of convergence of the iterative solver is strongly dependent on the spectral properties of the matrix system.

The preconditioning aims at transforming the original system of equations into another, formally equivalent, one with improved spectral properties. Here, an explicit left preconditioning approach is followed, where the original sys-tem ( 23) is left-multiplied by a preconditioning matrix P, to obtain

PAx = Pb (24)
with A = A near +A far . The matrix P should ideally be such that the governing matrix of the preconditioned system [START_REF] Alléon | Sparse approximate inverse preconditioning for dense linear system arising in computational electromagnetics[END_REF] is the identity matrix. Of course, this would amount to solve the original system by means of a direct solver, and therefore be of little practical value. In practice, a good preconditioning strategy is such that the matrix of the preconditioned system (24) is as close as possible to the identity matrix, with computing time and storage requirements for the preconditioning significantly lower, and growing slower with N, than those expected for the original system of equations. Here, following [START_REF] Margonari | Fast multipole method applied to elastostatic BEM-FEM coupling[END_REF] and earlier works from the computational electromagnetism community [START_REF] Alléon | Sparse approximate inverse preconditioning for dense linear system arising in computational electromagnetics[END_REF], the preconditioning matrix P is chosen as the SParse Approximate Inverse (SPAI) of A. A complete description of this approach is given in [START_REF] Grote | Parallel preconditioning with sparse approximate inverses[END_REF][START_REF] Benzi | A sparse approximate inverse preconditioner for nonsymmetric linear sytem[END_REF], where the nonsingular character and other theoretical properties of the SPAI preconditioner are discussed.

In its simplest version, a SPAI P of A is defined by choosing a priori a sparsity pattern for P, and then finding P as the solution to the minimization problem

min I -PA 2 (25)
(where • denotes the Frobenius matrix norm) which satisfies the preset sparsity pattern. The minimization problem ( 25) is then shown to take the form min I -

P l A 2 = n k=1 min e k -P l k A 2 (26)
where e k and P l k are the k-th row of the identity matrix and matrix P l respectively. Hence, each row P l k of P can be found by solving an independent least-squares problem whose size is equal to the number of nonzero entries in P l predicted by the preset sparsity pattern. The sparsity pattern of P l should ideally be chosen such that the entries of A -1 with largest moduli are retained. However, the locations of these entries are not known a priori. In [START_REF] Grote | Parallel preconditioning with sparse approximate inverses[END_REF] the sparsity pattern is found for every row (or column) with an adaptative strategy. The positions of the non-zero terms in P l are iteratively modified to minimize the residual provided by Eq. [START_REF] Benzi | A sparse approximate inverse preconditioner for nonsymmetric linear sytem[END_REF]. In [START_REF] Alléon | Sparse approximate inverse preconditioning for dense linear system arising in computational electromagnetics[END_REF] a different approach is adopted, where four heuristics are proposed for choosing a fixed sparsity pattern for P.

Here, the sparsity pattern of P is defined as corresponding to the location of the m largest entries in each row of A (with m ≪ n). On defining the sparse matrix  from A according to this same sparsity pattern, the minimization problem [START_REF] Benzi | A sparse approximate inverse preconditioner for nonsymmetric linear sytem[END_REF] can be simplified by dropping out the zero rows, to obtain:

n k=1 min êk -P l k  2 (27)
where (•) denote entities reduced according to the sparsity pattern. Each least squares problem in ( 27) is therefore of small size m ≪ n (and hence can be solved inexpensively). Well-posedness of the m × m least squares problems is assured by including all diagonal entries of A in the sparsity pattern.

This type of preconditioners seems to be quite effective also for non-symmetric matrices even if a sparsity of around 97 -99% is imposed for P l .

The stopping criterion adopted for the GMRES algorithm is based on the norm of the residual of the original system of equations, according to

b -Ax b ≤ Tol
where the tolerance Tol is chosen a priori.

NUMERICAL EXAMPLES

Numerical results obtained on the basis of the present FM-BEM are now presented. The convergence rate of the iterative system solution GMRES was accelerated by using the left-preconditioning strategy of Section 4. The efficiency (in CPU time sense) and the accuracy of results is compared to the classical BEM, for which GMRES with left-preconditioning was also used. The truncation parameters were set according to [START_REF] Darve | The fast multipole method: numerical implementation[END_REF], with C = 7.5 in (9) found to provide the best compromise between efficiency and accuracy. The recursive cell subdivision is stopped at the level for which the cell size d is such that

k S d min ≤ 2π 3
The GMRES algorithm is used with a stopping criterion defined by

b -Ax b ≤ 10 -3
and a restart parameter of 600 iterations. The problems were run on a work station with 4 Processors 64bits Opteron 2.2 Ghz and 16 GB of core memory.

Cubic box under traction

A cube with dimensions H × H × H (Fig. 3) is considered. Normal displacements and in-plane tractions are set to zero on faces x, z = ±H/2 and y = -H, while a time-harmonic normal traction p is imposed on the face y = 0. The relevant parameters are set as follows: H = 1 m, µ = 10 11 N/m 2 , ν = 0.25, p = 100 N/m 2 , ρ = 8 10 3 Kg/m 3 .

BE meshes made of four-noded bilinear elements are set up, so that each face features a uniform grid of M × M four-noded bilinear elements of size a × a. The linear element size a is chosen according to the frequency so that the shear wavelength is λ T = 20a. Four cases with increasing frequencies and problem sizes have been considered, with meshes such that each face of the cube features M × M square elements, with M = 20, 34, 48, 86 respectively. The resulting numbers of BE DOFs are indicated in Table 1, and the mesh for problem 1 shown in Fig. 4. The oct-trees featured in problems 1, 2, 3 and 4 feature 2, 3, 3 and 4 levels, respectively.

CPU times and GMRES iteration counts recorded for the various phases of the computation are shown in Table 2, where the heading "other" indicates the cumulative CPU time entailed by all tasks other than the computation of incoming and outgoing rays. Figure 5 graphically shows that the CPU times recorded for the present FM-BEM approach are significantly lower than those needed for the traditional BEM on the same problems. Problem 4, involving O(10 5 ) DOFs, could not be solved using traditional BEM because the size of the resulting governing fully-populated BEM matrix (about 250 GB) largely exceeds the RAM capacity of the computer on which the computations were run. Moreover, the CPU times per GMRES iteration are consistent with the theoretical estimate O(N ln N) [START_REF] Darve | The fast multipole method: numerical implementation[END_REF], as shown in Fig. 6.

In Fig. 7, the y-component of displacement results along the s-coordinate defined in Fig. 3), obtained with the present FM-BEM and the traditional BEM are seen to agree well for the four above-defined problems. The number of GMRES iterations is seen to increase with the frequency (to which the problem size is adjusted), as expected since the spectral properties of the linear system being solved is frequency-dependent. The preconditioning strategy presented in section 4 leads to a substantial reduction of the GMRES iteration count, as can be seen in Table 2, where the iteration counts recorded with and without preconditioning are given. For completeness and comparison purposes, the GMRES iteration counts and overall CPU times recorded using the traditional BEM with preconditioning are given in Table 3.

Ellipsoidal cavity in a rectangular box

A traction-free ellipsoidal cavity (with principal semiaxes (a 1 , a 2 , a 3 ) = (3, 2, 2)a, a being an arbitrary reference length, and axis a 1 oriented along direction (1, 1, 1)) embedded in a rectangular box (with dimensions (15a, 15a, 20a) in directions (x, y, z), respectively, see Fig. 8) is considered. The material param-eters and the boundary conditions on the external surface are defined as in the first example. The frequency is set so that λ T = 5a. The cavity center and the box center coincide. Four-noded elements are used for the cavity boundary (Fig. 8) and the box surface, the mesh densities being adjusted so as to feature 20 and 30 elements per shear wavelength on the external and cavity surface, respectively. Overall, the problem involves N = 81054 BE DOFs. Its solution required 441 iterations and 29714 seconds. The y-component of the displacement on the box surface is shown in Fig. 9. This example corresponds to a generic forward simulation for a defect identification problem.

CONCLUSIONS

A fast multipole accelerated BEM formulation for solving 3-D frequency domain elastodynamics has been presented in this paper. The expansion of the elastodynamic kernels is based on the diagonal-form, plane-vave, expansion of the scalar acoustic full-space fundamental solution. The divergent behaviour in the large-index limit of the spherical Hankel function of the first kind featured in that expansion led to the imposition of a minimum size (relative to the wavelength) for the leaf cells. Another important parameter, namely the truncation parameter of the transfer function, is determined as a function of the cell size (and hence the level) by means of a rule previously introduced for time-harmonic acoustic or electromagnetic computations.

Numerical results show that the numerical solution obtained by the present elastodynamic FM-BEM is, in practical terms, identical to that obtained by the traditional BEM (the two solutions differ by at most 3%) while requiring about 10 times less CPU time for the examples run (which were of size N = O(10 4 -10 5 ). The convergence rate of the GMRES iterative solver was found to be significantly accelerated by left-preconditioning the linear system with the sparse approximate inverse (SPAI) technique.

Fig. 1 .Fig. 3 .

 13 Fig. 1. Pole translations.

Fig. 4 .

 4 Fig. 4. Cube under traction: mesh for problem 1 (2400 elements).

Fig. 5 .Fig. 6 .

 56 Fig. 5. Cube under traction: overall CPU times for the FM-BEM and classical BEM computations.

Fig. 7 .

 7 Fig. 7. Cube under traction: results.

Fig. 8 .

 8 Fig. 8. Ellipsoidal cavity in a rectangular box: mesh of the cavity boundary. Overall problem size: N = 81054.

Fig. 9 .

 9 Fig. 9. Ellipsoidal cavity in a rectangular box: y-component of displacement [m] on the box surface.

Table 2

 2 Cube under traction: CPU times (seconds) and GMRES iteration counts both with and without SPAI preconditioning (Pre and No-Pre, respectively). Problem Out-rays In-rays Other Iters (Pre) Iters (No-Pre) Total time (Pre)

	1	2.3	1.2	0.03	26	32	100.5
	2	4.4	6.5	0.05	53	79	608.9
	3	13.2	12.8	0.2	78	130	2098.1
	4	46.3	74.8	1.9	279	685	34708.1

Table 3

 3 Cube under traction: CPU times (seconds) and GMRES iteration counts using the traditional BEM approach.

	Problem Iters (Pre) Total time (Pre)
	1	24	359.1
	2	44	4270.2
	3	69	22135.8

ACKNOWLEDGEMENTS

This work was supported by the Ministerio de Educación y Ciencia of Spain. (Project: DPI2004-08147-C02-02). The financial support is gratefully acknowledged. We also acknowledge Dr. Massimiliano Margonari for having kindly provided his elastostatic FM-BEM source code.

List of Figures

List of Tables