Yun Kai

Somayajulu Dhulipala

Rachel Sun

Jet Lem

Washington Delima

Thomas Pezeril

Carlos M Portela

Dynamic Diagnosis of Metamaterials via Laser-Induced Vibrational Signatures

 or lack of existing testing protocols. Here we demonstrate a high-throughput non-contact framework that employs MHz-wave propagation signatures within a metamaterial to nondestructively extract dynamic linear properties, omnidirectional elastic information, damping properties, and defect quantification. Using rod-like tessellations of microscopic metamaterials, we report up to 94% directionand rate-dependent dynamic stiffening at strain rates approaching 10 2 s -1 , in addition to damping properties 3 times higher than their constituent materials. We also show that frequency shifts in the vibrational response allow for characterization of invisible defects within the metamaterials, and that selective probing allows for construction of experimental elastic surfaces, previously only possible computationally. Our work provides a route for accelerated data-driven discovery of materials and microdevices for dynamic applications such as protective structures, medical ultrasound, or vibration isolation.

 , or high-rate deformation 11 have been explored theoretically and computationally, but experimental demonstrations have remained scarce. While

interest in dynamic properties of metamaterials has outpaced current characterization routes, recent studies have shown promising advancements in some dynamic regimes, from piezo-excitation [START_REF] Frenzel | Ultrasound experiments on acoustical activity in chiral mechanical metamaterials[END_REF][START_REF] Iglesias Martínez | Experimental observation of roton-like dispersion relations in metamaterials[END_REF] to underwater ultrasound propagation [START_REF] Krödel | Microlattice Metamaterials for Tailoring Ultrasonic Transmission with Elastoacoustic Hybridization[END_REF] , elucidating routes to induce polarization and form band gaps for mechanical waves. However, these pioneering methods have been limited to extracting wave-propagation properties, while extraction of dynamic elastic properties remains challenging. Their contact-based and fluid-coupled measurements pose limitations towards performing highthroughput characterization-a requirement for accelerated progress in understanding metamaterials under dynamic conditions. Microscale metamaterial properties have historically been extracted from stress-strain curves captured by nanoindentation tools, 1-3, 18, 19, 28 which face a trade-off between deformation amplitude and increases in strain rate. As a result, dynamic mechanical properties of microscale metamaterials [START_REF] Hawreliak | Dynamic Behavior of Engineered Lattice Materials[END_REF][START_REF] Lind | In situ dynamic compression wave behavior in additively manufactured lattice materials[END_REF] and full elasticity-tensor approximations [START_REF] Deshpande | Effective properties of the octet-truss lattice material[END_REF][START_REF] Bastek | Inverting the structure-property map of truss metamaterials by deep learning[END_REF] have only been possible using dynamic numerical frameworks or computational homogenization approaches, respectively. These limitations have left a void in understanding anisotropy, damping properties, and wave propagation in advanced metamaterials. While dynamic characterization techniques at large scales like the split Hopkinson bar method [START_REF] Weeks | High strain-rate compression behavior of polymeric rod and plate Kelvin lattice structures[END_REF][START_REF] Weeks | Shock compression behavior of stainless steel 316L octet-truss lattice structures[END_REF][START_REF] Tancogne-Dejean | Additively-manufactured metallic microlattice materials for high specific energy absorption under static and dynamic loading[END_REF] have served as a proxy to understand microscopic metamaterial behaviour, the contact-based, destructive nature of these methods results in low-throughput, onetime measurements with considerable variability. With advances in data-driven design tools for metamaterials, [START_REF] Bastek | Inverting the structure-property map of truss metamaterials by deep learning[END_REF][START_REF] Gongora | Designing lattices for impact protection using transfer learning[END_REF][START_REF] Mao | Designing complex architectured materials with generative adversarial networks[END_REF] experimental throughput has failed to keep up with vast design parameters spaces, particularly in dynamic regimes. A high-throughput method and underlying principles to uncover dynamic responses of metamaterials at small scales could enable characterization and data collection necessary for accelerating design and discovery of new types of metamaterials for various multifunctional applications.

Here, we introduce a framework for characterizing the dynamic response of 3D metamaterials at the microscale, enabled by a non-destructive and non-contact optical framework that converts vibrational signatures to dynamic elastic properties and wave propagation measurements. We demonstrate our high-throughput framework on several morphologies of 3D microscale metamaterials by leveraging pulsed-laser-based photoacoustic excitation [START_REF] Abi Ghanem | Longitudinal eigenvibration of multilayer colloidal crystals and the effect of nanoscale contact bridges[END_REF] of three types of elastic waves within the materials. Given low-amplitude elastic deformation and short timescales of excitation, we employed 10 4 measurements per sample to extract precise damping-property measurements, omnidirectional stiffness, and defect quantification within complex microscale samples. This research paves the way for employing nondestructive vibrational signatures as an alternative to slower microscale contact-based approaches, enabling understanding and discovery of metamaterials in regimes closer to real-world applications.

Laser-induced elastic waves

To enable high-throughput non-destructive characterization of metamaterials at the microscale, we designed a pulsed-laser-based technique to leverage photoacoustic excitation of elastic waves on the surface of the material sample, analogous to photoacoustic thin-film and polycrystalline metal characterization studies [START_REF] Akimov | Coherent hypersonic closed-pipe organ like modes in supported polymer films[END_REF][START_REF] Dryburgh | Measurement of the single crystal elasticity matrix of polycrystalline materials[END_REF] . Our method-laser-induced resonant acoustic spectroscopy (LIRAS)-is composed of two separate pump and probe modules (Fig. 1a). The pump module consists of a picosecond pulsed laser applied to a three-dimensional sample of interest with a thin chromium coating which partially absorbs the laser pulse and induces broadband acoustic waves (Methods). Resulting sample surface displacements are then measured by the probe module, consisting of a phase-mask interferometer. To facilitate extraction of metamaterial mechanical properties, we fabricated polymeric microlattices out of IP-Dip photoresist in slender aspect ratios using a two-photon lithography system (Nanoscribe GmbH), which approached the limit of rod-like 1D-wave propagation upon excitation. Specifically, we fabricated 5×5×n tessellations of 9.0±0.1 µm octet and tetrakaidecahedron unit cells, with the height of the tessellation n ranging from 4 to 22. Octet and tetrakaidecahedron morphologies were selected as representatives of kinematically rigid and non-rigid architectures, respectively. At relative densities of ρ = 17% for the octet and 14% for the tetrakaidecahedron, these translated to stretching-and bending-dominated mechanical responses [START_REF] Meza | Reexamining the mechanical property space of threedimensional lattice architectures[END_REF] . Exploiting cubic symmetry in these classical geometries enabled metamaterial anisotropy detection (insets in Fig. 1a) by fabricating metamaterial samples along lattice orientations such as [100], [110], and [111], and monolithic IP-Dip polymer samples for comparison.

To investigate different types of elastic metamaterial responses, our LIRAS technique enabled excitation and measurement of three types of elastic waves: (i) longitudinal, (ii) torsional, and (iii) flexural waves. Primary excitation of each of these waves was achieved by tuning the pump-continuous-wave (CW) probe scheme of the lasers (Fig. 1b). Specifically, coincident pump-CW probe lasers at the center of the sample induced a primarily longitudinal wave in the metamaterial, while lateral excitation at the top of a modified sample with separate pump-probe locations induced a torsional wave. Similarly, off-center excitation from the top of the sample with different pump-probe sites induced a primarily flexural wave in the metamaterial. As shown in Fig. 1c, characteristic response waveforms attenuated due to inherent vibrational decay but lasted up to ∼200 µs depending on sample type and pump scheme (Methods), and experimental calibration enabled translation from photovoltage to surface displacement (Supplementary Information section I).

Metamaterial fabrication onto a glass substrate resulted in clamped and free boundary conditions at the base and top, respectively, inducing standing waves of wavelength λ = 4H/m upon laser excitation, where H denotes sample height and m denotes harmonic order (i.e., m = 1 indicates the fundamental first harmonic). Consequently, varying sample heights allowed for wavevector k tunability, with the magnitude defined as wavenumber |k| = 1/λ. Fabricating and testing samples with systematically increasing heights enabled extracting their resonant frequencies for different wavevectors (Fig. 2a).

Figure 2 LIRAS dynamic measurements on metamaterials. a, SEM micrographs of selected octet [100] samples of different n tessellations, corresponding to different wavenumbers, scale bar 50 µm. b, Fast-Fourier transform (FFT) spectra of detected signals from the surface of the metamaterials using the center pump-CW probe scheme, depicting longitudinal mode (L) and its higher harmonics (e.g., 3L). Spectra corresponding to different n-tessellations are plotted with an offset for illustration purposes. c, FFT spectra from off-center scheme, depicting appearance of flexural modes (F) and higher harmonics; and d, spectra from the transverse scheme depicting torsional modes (T) and several flexural harmonics. e, Dispersion relations f (k) extracted from FFT spectra along specific edges of the octet's irreducible Brillouin zone (top left, along with Bravais-and reciprocal-lattice vectors, a i and b i , respectively).

Following excitation by 10 4 pump-laser pulses per sample (over 1 to <5 minutes, Supplementary Information section II), frequency content of time-domain dynamic responses for the metamaterials was extracted by taking their fast Fourier transform (FFT). Using the center pump-CW probe scheme, we determined the presence of a fundamental longitudinal-mode peak for all 5×5×n samples (L in monolithic and octet data in Fig. 2b. Besides this dominant peak, higher-order harmonics were also measurable for most samples (e.g., 3L and 5L, Fig. 2b). We additionally measured a height-independent Poisson-effect transverse mode in monolithic samples, which was not measurable in metamaterials due to prominent microstructural scattering in this shorter-wavelength regime. Using the off-center pump-CW probe scheme, we primarily excited a flexural mode (and higher-order harmonics) in addition to a weaker longitudinal mode (Fig. 2c). To measure torsional modes, we implemented a transverse pump-CW probe scheme, where metamaterial samples were fabricated with a flat reflective surface on the sides and were pumped and probed at two different transverse locations on that surface (see Fig. 2d). This scheme excited primarily flexural modes, but enabled measurement of a torsional mode. To correlate experimental spectra to corresponding vibrational modes, we confirmed a long-wavelength-limit linear relation between wavenumber and frequency for longitudinal and torsional waves, while a quadratic relation was identified for flexural waves. We further confirmed modal categorization via numerical eigenfrequency analyses on monolithic and metamaterial geometries, which provided qualitative and quantitative agreement for longitudinal, flexural, and torsional modes.

After identifying frequency peaks and their corresponding modes for all samples, each frequency was correlated to a wavenumber based on sample height, enabling construction of experimental dispersion relations for monolithic and metamaterial samples as shown in Fig. 2e. Cubic symmetry enabled probing along edges of the irreducible Brillouin zone, corresponding to Γ-X for [100], Γ-M for [110], and Γ-R for [111]. For longitudinal modes in all samples, dispersion was observed for larger wavenumbers, primarily since the samples in this regime violated assumptions for 1D-wave propagation due to their non-slender aspect ratio. Nondispersive longitudinal waves were observed for metamaterial tessellations of n ≥ 10 when considering only the fundamental mode. Additionally, torsional and flexural waves in monolithic and [100] metamaterial samples were measured, but were excluded in other metamaterial-orientation analyses due to redundancy.

High-throughput dynamic mechanical characterization

While dispersion relations obtained via the LIRAS technique provided insight into different types of elastic waves supported by these metamaterials, we aimed to leverage this information to determine their dynamic mechanical properties. Using measurements in the long-wavelength limit (|k| ≤ 0.0025 µm -1 , tessellations of n ≥ 10), we calculated group velocity of each wave type from experimental dispersion relations as c g = ∂ f ∂k . Following theory for 1D-longitudinal-wave propagation in a rod-like medium, we used longitudinal group velocities c g,L to calculate the effective stiffness for each sample as E * = ρc 2 g,L , where ρ is the sample density (Methods). From displacement-calibrated LIRAS experiments, which enabled effective engineering strain measure-ments ranging from 0.01% to 0.1%, we determined characteristic root mean square (RMS) strain rates of ε = 38 s -1 for monolithic polymer and ε = 70 s -1 for octet metamaterials-corresponding to the MHz ultrasound dynamic regime (Methods, Supplementary Information Sections III & IV). As expected for octets and tetrakaidecahedra, measured effective stiffness (Fig. 3a) exhibited direction-dependent values, with [111] and [100] directions yielding the highest stiffness for octet and tetrakaidecahedron samples, respectively. [START_REF] Garrett | Understanding Acoustics: An Experimentalist's View of Sound and Vibration[END_REF] , where β( f ) is the temporal damping coefficient and c p,L is the longitudinal phase velocity, depicting higher and anisotropic damping over our measurable frequency range as a power law α( f ) = α 0 f η with constant a 0 and tunable exponent η ranging from 1.35 to 1.51 (slope of the log-log plot) for octet and 1.39 to 2.02 for tetrakaidecahedron, compared to 1.58 for monolithic polymer.

( f ) = β( f )/c p,L
To quantify dynamic effects on the responses of our samples, we performed quasi-static ( ε = 10 -3 s -1 ) uniaxial compression experiments using a displacement-controlled nanoindenter (Alemnis AG) on identical samples as those tested via LIRAS. Comparing the effective stiffness, E * (= ∂σ/∂ε), obtained through these two techniques (Fig. 3a) showed dynamic stiffening effects 8 ACCEPTED MANUSCRIPT / CLEAN COPY in all samples. From experiments on the monolithic polymer samples, we measured a viscoelastic response in the base material with a dynamic stiffness of 4740 ± 30 MPa, corresponding to a 74% stiffening compared to the measured quasi-static effective Young's modulus, consistent with other dynamic compressions on this material [START_REF] Rohbeck | Effect of high strain rates and temperature on the micromechanical properties of 3D-printed polymer structures made by two-photon lithography[END_REF] . This stiffening response was also observed in the metamaterial samples, with an on-average stiffening of 46% for the octets and 82% for the tetrakaidecahedra. Despite this direction-dependent dynamic stiffening, the anisotropic response of the metamaterials remained qualitatively unchanged.

Additional uniaxial compression experiments to strain rates of up to 10 2 s -1 on both monolithic polymer micropillars and 5 × 5 × 5 tessellations of the [100] octet and tetrakaidecahedron metamaterials (see Methods), confirmed the rate-dependent response observed in LIRAS measurements. Despite the different nature of both characterization techniques, i.e., the contact-based method targets a single strain rate while a range of strain rates are accessed in the free-vibration case, agreement to within 10% was observed across all sample types (Fig. 3b, Supplementary Information Section IV). These measurements confirmed that achieving high-rate nanometer-scale displacements via photoacoustic excitation effectively allows for non-destructive dynamic characterization, presenting an alternative to single-measurement destructive experiments and their concomitant signal-to-noise ratio challenges.

In addition to enabling high-throughput linear mechanical property measurements, LIRAS frequency-spectrum analysis provided metamaterial damping properties. [START_REF] Hussein | An emergent phenomenon in dissipative metamaterials[END_REF][START_REF] Salari-Sharif | Damping of selectively bonded 3D woven lattice materials[END_REF] Analogous to a classical damped harmonic oscillator, whose oscillation amplitude is described by u(t) = u 0 e -βt cos(2π fδ), the exponential temporal decay of the LIRAS experiments from an initial amplitude u 0 at a frequency f and phase shift δ is dictated by the decay parameter β (Supplementary Information Section V). Using center-pump FFT spectra, we calculated the more informative spatial attenuation coefficient 42 α = β/c p,L (units m -1 ), where β scales with the half-width-half-max of the longitudinal peak, and c p,L = f /|k| is phase velocity of the longitudinal waves. Plotting the experimental spatial attenuation coefficient against the dominant longitudinal frequency peak, we measured direction-dependent damping in our metamaterials (Fig. 3c). Specifically, higher damping was observed for the more compliant direction [100] in octets, while the stiffer [110] and [111] directions exhibited an average of 30% lower damping at a frequency of 3 MHz. In tetrakaidecahedra, the compliant [111] direction also exhibited higher damping at a given frequency than its stiffer [100] direction. Remarkably, the octet metamaterials exhibited up to a 3.5× increase in damping ability compared to their constituent monolithic polymer. Fitting spatial attenuation values to the classical frequency-dependent attenuation power law 43 α( f ) = α 0 f η yielded a scaling exponent η = 1.51 (slope in Fig. 3c) for monolithic samples, in alignment with other polymer measurements [START_REF] Szabo | A model for longitudinal and shear wave propagation in viscoelastic media[END_REF] . In contrast, the metamaterials exhibited direction-dependent scaling, with slopes of 1.51, 1.44, and 1.35 for the octet [100], [110] and [111] directions, respectively, while the tetrakaidecahedron had slopes of 1.39, 2.02, and 1.56 for the same orientations. The measurements indicate that architecture and anisotropy can tune damping properties at a wide range of ultrasonic frequencies. Both metamaterials indicate higher damping capacity in compliant directions and bending-dominated architectures, over stiff and stretching-dominated responses.

Omnidirectional elastic properties

Established characterization methods for metamaterials have primarily been limited to reporting effective Young's moduli only along selected orientations [START_REF] Zheng | ultrastiff mechanical metamaterials[END_REF][START_REF] Crook | Plate-nanolattices at the theoretical limit of stiffness and strength[END_REF][START_REF] Meza | Reexamining the mechanical property space of threedimensional lattice architectures[END_REF] , providing only partial information of material anisotropy and lacking measurement of properties such as effective shear moduli. Numerical methods like computational homogenization have been the preferred technique to obtain full elastic-tensor information of metamaterial geometries, but provide metrics that have not been experimentally validated [START_REF] Portela | Extreme mechanical resilience of self-assembled nanolabyrinthine materials[END_REF][START_REF] Bastek | Inverting the structure-property map of truss metamaterials by deep learning[END_REF][START_REF] Patil | Effective property evaluation and analysis of three-dimensional periodic lattices and composites through Bloch-wave homogenization[END_REF] . To fill this experimental gap, we harnessed the capabilities of the LIRAS technique to provide fully experimental elastic-tensor extraction for metamaterials under dynamic conditions.

To obtain preliminary validation of properties obtained from LIRAS, pertaining both to the measured frequencies and apparent anisotropy, we employed three computational approaches that together predict expected anisotropy and wave propagation characteristics within the metamaterials. For consistency with the viscoelastic response in the constituent polymer, the material properties used in dynamics simulations corresponded to the RMS strain rate in our experiments, while upper-and lower-bound simulations employed properties corresponding to the maximum strain rate and quasi-static conditions, respectively (Methods, Supplementary Information section IV). The first comparison approach consisted of finite eigenfrequency calculations on a finite-element representation of discrete 5 × 5 × n tessellations, with boundary conditions identical to those in our LIRAS experiments. Identification of longitudinal, torsional, and flexural resonant modes and associated eigenfrequencies agreed with experimental measurements especially in the longwavelength limit, as shown in Fig. 4a. As a second validation method, we performed Bloch-wave analysis on metamaterial unit cells, which provided dispersion relations that enabled longitudinal and shear mode identification. Using these two computational methods, we calculated the predicted [100] stiffness for the metamaterials and confirmed an agreement to within 5% when comparing the stiffness from eigenfrequency simulations to the effective stiffness in the LIRAS experiments. The third computational approach consisted of linear static homogenization, which extracted the metamaterials' elasticity tensor C, and enabled calculation of their effective stiffness in all directions [START_REF] Portela | Extreme mechanical resilience of self-assembled nanolabyrinthine materials[END_REF][START_REF] Meza | Reexamining the mechanical property space of threedimensional lattice architectures[END_REF] (Methods) to be represented as an elastic surface (Fig. 4b, top). After confirming wave propagation response and resulting dynamic stiffness along specific directions, we combined all properties extracted from LIRAS measurements to produce an experimental dynamic elastic surface for our metamaterials (Fig. 4b, Supplementary Information Section VI). Given their cubic symmetry, we used the experimental effective stiffness along the [100], [110], and [111] directions, and the torsional response of the [100] samples, to construct dynamic elasticity and compliance tensors Ĉ and Ŝ = Ĉ-1 , respectively. Through the relation Êd = 1/( Ŝi jkl d i d j d k d l ) for a given probing direction d, we extracted non-shear components of the compliance tensor Ŝ. To obtain the shear components, we leveraged torsional waves to relate the torsional wave velocity c T to the effective dynamic shear stiffness Ĝ = Ŝ -1 44 for these cubic metamaterials via the expression c T = K Ĝ/ρ, where K is a constant defined by the cross-section of the body undergoing torsion [START_REF] Graff | Wave Motion in Elastic Solids[END_REF] and ρ again denotes density of the metamaterial (Supplemen-tary Information Section VII). The experimental dynamic elastic surfaces preserved the qualitative anisotropy for the metamaterials, with [111] and [100] remaining as the stiffest directions in the octet and tetrakaidecahedron, respectively, and exhibited dynamic stiffening (Fig. 4b, center) which was computed as the difference between static and dynamic elastic surfaces. Dynamic stiffening for both metamaterials was higher for the stiffest direction, at ∼49% and ∼95% for the octets and tetrakaidecahedron, respectively. Despite the octet exhibiting omnidirectional stiffening greater than 40%, the tetrakaidecahedron exhibited higher overall dynamic stiffening of 94%, hinting to higher inertial contribution in this bending-dominated architecture (Supplementary Information section VIII).

10

ACCEPTED MANUSCRIPT / CLEAN COPY

Non-destructive defect identification

Despite multiple routes to fabricate mechanical metamaterials across scales, various forms of defects emerge from fabrication or processing issues and ultimately degrade their mechanical properties [START_REF] Meza | Reexamining the mechanical property space of threedimensional lattice architectures[END_REF][START_REF] Gross | Correlation between topology and elastic properties of imperfect truss-lattice materials[END_REF][START_REF] Liu | Elastic and failure response of imperfect three-dimensional metallic lattices: the role of geometric defects induced by Selective Laser Melting[END_REF][START_REF] Glaesener | Predicting the influence of geometric imperfections on the mechanical response of 2D and 3D periodic trusses[END_REF] . At nano-to-microscales, metamaterials still rely on post-fabrication characterization via SEM or 3D X-ray computer tomography (XCT) [START_REF] Crook | Plate-nanolattices at the theoretical limit of stiffness and strength[END_REF] , with inevitable limitations of resolution and time-consuming scans. Other non-destructive techniques employing contact-based ultrasound have been widely employed in structural health monitoring, but they remain inapplicable to small length scales pertinent to metamaterials or microscale electro-mechanical devices.

As an alternative route, we leveraged LIRAS to provide high-throughput defect identification within microscopic samples by measuring changes in their wave propagation response. Specifically, following the center-pump scheme, we identified changes in resonant frequencies as a function of defect density. To validate these capabilities, we studied solid inclusions and missing features, which are common defects in microfabrication techniques such as two-photon lithography [START_REF] Meza | Reexamining the mechanical property space of threedimensional lattice architectures[END_REF][START_REF] Gross | Correlation between topology and elastic properties of imperfect truss-lattice materials[END_REF] or larger-scale light-based printing techniques. To mimic solid-inclusion defects commonly caused by unpolymerized resin trapped within the sample, we designed and fabricated 5 × 5 × 10 octet metamaterials with varying degrees of cubic solid inclusions per vertical tessellation layer-undetectable via standard SEM characterization (Fig. 5a). Matching inclusion sizes to those of ∼10 µm octet unit cells, we fabricated samples spanning from zero to four defects per layer, measuring their primary resonant mode upon excitation. These experiments identified a ∼2% to 8% increase in resonant frequency when the number of defects per layer ranged from 1 to 3, with respect to an inclusion-free sample (Fig. 5b). A 23% increase in frequency was observed at 4 inclusions per layer, likely due to cumulative solid inclusions resulting in solid domains across several layers which contributed to a significant stiffness increase.

To mimic the second defect type, we designed 5 × 5 × 10 octahedron metamaterials (unit cell size of 10 µm) with intentionally removed beams throughout the tessellation. We selected the octahedron morphology since it retained the octet's stretching-dominated response [START_REF] Meza | Reexamining the mechanical property space of threedimensional lattice architectures[END_REF] , but its lower beam connectivity enabled easier defect visualization (see Supplementary Information section IX for other morphologies). Starting with a perfect tessellation, we fabricated samples with various percentages of missing beams, with a maximum of 40% missing features to ensure structural integrity (Fig. 5c). LIRAS measurements indicated a gradual decrease in resonant frequency for longitudinal and flexural modes (Fig. 5d, Supplementary Information Section IX). Assuming a linear decrease in longitudinal frequency as a function of missing beams, we identified a 27 kHz decrease per percent of missing beams, compared to a 7 kHz decrease per missing-beam percent for flexural modes. Both defect explorations indicated that a resonant frequency shift provided sufficient resolution to determine defect densities that were otherwise impossible to identify or quantify efficiently.

Summary and Outlook

In summary, we have designed and demonstrated a laser-based framework that enables non-contact high-throughput characterization of microscopic metamaterials. Leveraging elastic wave propagation within the materials, we extracted the full dynamic elasticity tensor and characterized mechanical-wave propagation, damping properties, and quantified the presence of defects. Fur-

13

ACCEPTED MANUSCRIPT / CLEAN COPY thermore, we uncovered the interplay between viscoelastic material properties and 3D architecture, demonstrating direction-and architecture-dependent dynamic stiffening and damping properties. With growing scientific interest in the dynamics of metamaterials, experimental validation to predicted responses is needed, opening the door to high-throughput and non-destructive characterization under extreme conditions. [START_REF] Dattelbaum | Shockwave dissipation by interface-dominated porous structures[END_REF] High-throughput characterization techniques (refer to tabular comparison between different techniques in Supplementary Information Section X) in these regimes will also enable rapid creation of experimental datasets for data-driven design, accelerating metamaterial applications in thriving fields relevant to medical ultrasound [START_REF] Wang | Bioadhesive ultrasound for long-term continuous imaging of diverse organs[END_REF] and nonlinear wave propagation. [START_REF] Jin | Guided transition waves in multistable mechanical metamaterials[END_REF] This framework paves the way for high-throughput discovery of novel nanoand micro-architected materials, on-the-fly defect identification in advanced microscale devices, and microscopic acoustic metamaterial design for medical ultrasound devices.

Methods

Laser-Induced Resonant Acoustic Spectroscopy (LIRAS) method

The pump laser was generated from the uncompressed output of a Ti:Sapphire regenerative amplifier (Coherent Astrella) with a central wavelength of 800 nm, 1 kHz repetition rate, and a pulse duration of 300 ps. A pulse energy of 20 to 50 µJ was applied to induce sudden thermal expansion on the metal-coated (50 nm thick Cr or Au/Pd) sample surface. The probe laser was comprised of a phase mask interferometer driven by a continuous-wave (CW) laser (Coherent Verdi) with 532 nm wavelength and 10 mW maximum output power, while only ∼2.5 mW was shone on the sample to avoid damage, monitored by a CCD camera (Hamamatshu Orca C13440). This phase-mask interferometer measured the surface displacement of a given sample as a function of time. Specifically, two CW-laser beams-formed as a result of beam-splitting by the phase mask-were focused on (i) a vibrating sample, and (ii) a static reference sample, respectively. The reflected laser beams from the two samples were recombined by the same phase mask, and the resulting interferometric beam was registered by an avalanche photodetector (c5658, Hamamatsu) with 1 GHz bandwidth and 50 kHz low-frequency cut-off. While vibrating, the sample surface displacement u z (t) induced an optical path difference ξ = 2u z (t) between both arms of the interferometer, which was translated into a change of intensity of the signal beam and a measurable change in photovoltage by the photodetector.

Since the pump laser repetition rate was 1 kHz, it excited a resonant vibration 1,000 times per second in a given sample. Consequently, our measurement had a high signal-to-noise ratio attained by averaging of detected waveforms over 10,000 repeated measurements. Due to the acquisition speed in the waveform averaging mode of our oscilloscope (Tektronix TDS 7404), the actual measurement time ranged from 1 to ∼3 minutes depending on the signal length and sampling rate.

The interferometric signal V(t) registered by the photodetector was translated to the sample surface displacement u z (t) by u z (t) = λ/(4π) sin -1 (2V/V 0 ) with λ = 532 nm as the wavelength of the CW laser and V 0 as the operational range of the photovoltage. Note that the measurement range was limited to a maximum surface displacement of 66.5 nm when the interferometer was operating at the highest sensitivity. However, that range could be extended by tuning the initial phase offset of the interferometer, achieved by a horizontal translation of the phase mask. Because the pump laser induced not only acoustic waves, but also thermal expansion, the surface displacement of the sample showed a combined effect at two different time scales. To eliminate the slow-varying thermal signal (typically µs-periodicity), we fit the displacement waveform with a polynomial function of up to fifth order and subsequently subtracted it from the signal.

Fabrication of polymeric samples

Polymeric samples were fabricated out of IP-Dip resin using a two-photon polymerization system (PPGT2, Nanoscribe GmbH). Octet and tetrakaidecahedron unit cells of nominal dimensions of 10×10×10 µm 3 and strut diameters of 1 µm were the standard metamaterial unit cells. Tessellations of 5 × 5 × n unit cells, with n ranging from 4 to 22, comprised the rod-like metamaterial samples in this work. Samples were fabricated with a writing speed of ∼10 mm s -1 and laser power of 30 mW for lattices and 25 mW for monolithic samples. The samples were printed on a glass chip coated with an indium tin oxide (ITO) thin film as the substrate; the coating facilitated better interface finding between the resin and the glass. All samples were developed for 20 min in propylene glycol monomethyl ether acetate (PGMEA), followed by a 5 min. rinse in isopropanol and subsequent drying in a critical point dryer (Autosamdri 931, Tousimis). After fabrication and drying, the chips were coated with 50 nm of Cr by a sputter deposition tool (Orion, AJA International Inc.), yielding a reflective surface for LIRAS excitation. Relative density of each lattice was determined by measuring geometric feature dimensions of each unit cell via scanning electron microscopy (SEM), including unit cell size and strut dimensions, then modeling the unit cell in SolidWorks. The volume of the unit cell was then tabulated and divided by the unit cell bounding-box volume to obtain the relative density. The IP-Dip polymer density was 1170 kg/m 3 (Nanoscribe GmbH), and the relative density for the as-fabricated octet metamaterials was 17 %, yielding 199 kg/m 3 as the effective octet lattice density. The relative density of the as-fabricated tetrakaidecahedron metamaterials was 14%, yielding 162 kg/m 3 as the effective tetrakaidecahedron lattice density.

Fabrication of defect-containing samples

The defect-containing 5 × 5 × 10 samples with bulk inclusions or missing beams were fabricated using the same parameters and process as the octet and tetrakaidecahedron samples. The inclusions were distributed inside the samples, within the inner 3 × 3 × 9 core, such that they were not identifiable via SEM characterization. The missing-beam locations in the octahedron (or braced cubic, Supplementary Information IX) geometries were randomly assigned, resulting in a random distribution of defects within a given sample. A custom code was employed to create missing-beam geometries with various percentages of missing struts, which produced STL files for fabrication via two-photon lithography.

Experimental dispersion relation calculation

The LIRAS voltage-time signals (converted to displacement-time signals via interferometer calibration) were analyzed in MATLAB to remove the thermal background by subtracting the raw data with a fifth-order polynomial fit function. The resulting waveform signal was conditioned with a moving-average filter spanning 5 data points. We extracted the resonant frequencies in the signal by taking the fast Fourier transform (FFT) of the time derivative of the waveform, followed by identification of frequency peaks. The wavenumber for a given peak was calculated as k = 1/(4H) with H denoting the sample height, which was measured via laser confocal microscopy and SEM characterization. Using frequencies extracted from the FFT and calculated wavenumbers from each sample, we constructed partial experimental dispersion relations by plotting frequency f as a function of wavenumber k. To correctly assign each vibration mode to its corresponding dispersion relation, we selectively excited resonant modes and distinguished dispersion trends accordingly. As demonstrated in Fig. 2b-d, a center-pump scheme mainly induced a longitudinal mode, which appeared as the dominant peak in the frequency spectra. Correspondingly, the offcenter-pump scheme led to a dominant spectrum peak related to a bending or flexural mode. The transverse-pump scheme excited both flexural and torsional modes. Dispersion analysis [START_REF] Graff | Wave Motion in Elastic Solids[END_REF] was used to distinguish flexural and torsional modes. Since the fundamental mode was the most dominant peak in the spectra, it was the first to be determined. Subsequently, the expected 3rd and 5th-order harmonics of the modes were calculated by multiplying the measured fundamental frequencies by factors of 3 and 5. All the other measured spectrum peaks were plotted on the same dispersion figure to compare with the higher harmonics of the modes. The modes that matched the expected values were confirmed as higher harmonics. Torsional frequencies were found to be distinctly different to those of other modes.

LIRAS strain rate estimation

Once we obtain the thermal-signal-free time-displacement curves from a LIRAS measurement, dividing by the sample height approximates the effective strain as a function of time ε(t) (see Supplementary Information Section III). Subsequently, the corresponding strain rates ε(t) can be obtained by taking the time derivative of the strain. The characteristic strain rate for each structure is then calculated by taking the average of the RMS strain rate of the tallest 5 samples (longwavelength limit). Since free-vibration and subsequent attenuation of the samples encompassed a range of strain rates, the bounds for these ranges were calculated by considering the RMS strain rates for the first and last wavelengths of the signal as discussed in Supplementary Information Section IV.

Experimental dynamic elastic surface calculation

The directional stiffness was obtained from the longitudinal and torsional wave speeds for 5×5 × n octet and tetrakaidecahedron lattices in different crystallographic orientations. The components of the dynamic compliance tensor, Ŝi jkl were obtained by considering their relation to the effective directional stiffness in prescribed directions (see Supplementary Information Section V). The elastic surface can be computed by calculating the directional stiffness, which for a given direction d is given by E d = Ŝ i jkl d i d j d k d l .

Damping quantification

The damping coefficients are extracted using only the longitudinal mode due to its strong signal. Using MATLAB to analyze the FFT spectra, we extracted the central frequency of the fundamental longitudinal spectral peak in addition to its full width at half maximum (FWHM), as described 16 ACCEPTED MANUSCRIPT / CLEAN COPY in Supplementary Information Section V. Applying a damped harmonic oscillator model to our metamaterial, the FWHM in the frequency domain is related to the decay parameter β( f ) through β( f ) =FWHM/2. The spatial attenuation coefficient α( f ) can be thus derived from the measured acoustic wave velocities c p,L through α( f ) = β( f )/c p,L . The decay parameter of a classical damped harmonic oscillator can be found by fitting the envelope of the waveform with an exponential function Ae -βt . However, due to the multi-modal signal with multiple frequency components (Supplementary Information Section V), we elect to perform the damping analysis only in the frequency domain using solely the longitudinal-mode peak.

Quasi-static and dynamic contact-based uniaxial compression experiments

Uniaxial compressions were performed using an Alemnis ASA nanoindenter with a 400 µmdiameter diamond flat-punch tip at strain rates ranging from 10 -3 s -1 (quasi-static) to 10 2 s -1 . Quasi-static experiments were performed on IP-Dip pillars (162.8 µm tall and 50 µm-wide square cross-section) and 5×5×10 octet and tetrakaidecahedron [100], [110], and [111] lattices. The IP-Dip pillars and octet lattices were printed on a silicon substrate using two-photon polymerization at a laser power of 22.5 mW and 30 mW, respectively. A standard load cell with a load range of 2.5 N and a piezostack with a displacement range of 100 µm were used for compression. The tip was operated in a displacement-controlled mode. The high-strain-rate experiments were performed on IP-Dip pillars (diameter of 21.9 µm, height of 43.3 µm) and 5×5×5 octet and tetrakaidecahedron [100] lattices (for consistency, these samples were also compressed at a quasi-static strain rate of 10 -3 s -1 ). For high-rate compressions, we employed the Alemnis high strain rate (HSR) module, comprised of a mini load cell (MLC) connected to a signal amplifier. Displacement rates from 50 nm/s to 5 mm/s were applied while maintaining close to 2 × 10 3 points for sampling the data at all strain rates. The displacement rates were set by assuming a constant sensitivity (the ratio of the change in displacement to the change in voltage of the piezo stack) determined through an air indent and then setting up a voltage ramp with smoothing. For strain rates of 10 -3 and 10 -2 s -1 , a maximum displacement of 10 µm was applied. For a strain rate of 10 0 s -1 , a displacement of 30 µm was applied and the samples were compressed in the last 15 µm of the compression. This procedure ensured that the indenter tip accelerated and reached a constant displacement rate before coming in contact with the sample. Similarly for strain rates above 10 0 s -1 , a total displacement of 80 µm was applied with the sample coming in contact in the last 40% of the displacement range. The effective stiffness of the samples was measured by calculating the slope of a least-squares linear fit to the linear loading region of the stress-strain data. The strain range for the linear fit was restricted to 3% strain-well within the yield point of the samples. The strain rate for the samples was determined by measuring the total change in strain over the range of the linear fit for the stiffness and dividing that by the time taken for the change in strain. These are all found to be within 5% of the strain rate set for each study.

Computational Analysis Dynamic explicit simulations

We conducted dynamic explicit simulations in ABAQUS/Explicit on an octet and tetrakaidecahedron unit cells with unit cell size of 7.5 mm and relative densities of 14% and 17% respectively. The unit cells were meshed with C3D10M elements. The compression was enforced with top and bottom planes tied and coupled kinematically to the top and bottom of the unit cell, respectively. The top plate was set to compress such that a maximum 5% engineering strain was enforced at a 60 s -1 strain rate with general contact properties, and an encastre boundary condition for all degrees of freedom on the bottom plate. Material properties were set such that Young's modulus of the constituent material was 2.7 GPa, density was 1,200 kg/m 3 , and Poisson's ratio was 0.49.

Eigenfrequency analysis

We conducted the eigenfrequency simulations using the solid mechanics module in COMSOL. As in the experiments, the bottom of each metamaterial sample was fixed while all other parts were free to vibrate. The modeled geometry reflected measured values of fabricated samples, with octet unit cell sizes of 9 µm and strut radius of 0.58 µm. The models output the eigenfrequency f corresponding to each vibration mode for different sample heights, i.e., different wavenumbers k. Thereby, we constructed a simulated acoustic dispersion relation f (k) for finite lattice tessellations of a cross-section of 5 × 5 unit cells in x and y directions with varying numbers of unit cells in the z direction. Monolithic samples were also simulated with a similar aspect ratios, with a fixed cross-section of 50 µm by 50 µm with different heights such that a rod-wave approximation was achieved. Since different strain-rate compressions of monolithic samples confirmed a ratedependent stiffness of the IP-Dip polymer, the scaling was determined to be E * = 4.216 ε0.0445 GPa. The strain rate ε was estimated via the method described in the Supplemental. We simulated lower, middle, and upper eigenfrequencies by simulating with Young's moduli corresponding to the quasi-static effective Young's Modulus of 2.7 GPa, mean RMS dynamic effective stiffness of 4.522 GPa (strain rate of 4.84 s -1 ), and maximum dynamic effective stiffness of 5.168 GPa (strain rate of 97 s -1 ), respectively. Other material properties were held constant with Poisson's ratio of 0.49 and density of 1,170 kg/m 3 . Mesh size was determined by keeping a constant maximum and minimum element size, leading to ∼ 85, 000 tetrahedral elements per tessellation layer, as this was the finest possible mesh size due to computational constraints.

Bloch-wave analysis

In addition to the eigenfrequency simulation, we applied Bloch-wave analysis within a finite element method framework to obtain a dispersion relation for octet lattices. We simulated Bloch-wave analysis with experimental geometric values of octet unit cell size 9 µm, elliptical strut width of 0.9 µm, elliptical strut height of 1.96 µm, and material properties of Young's modulus of 4.522 GPa (corresponding to strain rate of 4.84 1/s), Poisson's ratio of 0.49, and density of 1170 kg/m 3 . We used Floquet-Bloch periodic boundary conditions in the COMSOL Solid Mechanics module, using symmetric fine physics-controlled mesh. We extracted the first three modes corresponding to the lowest longitudinal and two transverse waves of the lattice structure and calculated a best fit using least-squares regression to extract the longitudinal and shear wave velocities, which were used to populate the stiffness tensor; we then obtained Young's modulus and shear modulus from the compliance tensor [START_REF] Patil | Effective property evaluation and analysis of three-dimensional periodic lattices and composites through Bloch-wave homogenization[END_REF] . We used volume-averaged modal displacement along principal directions to identify the mode shapes.

Static homogenization

To obtain the directional stiffness, computational homogenization for the octet and tetrakaidecahedron morphologies was performed on ABAQUS/Standard (Simulia). Octet and tetrakaidecahedron unit cells were built using the ABAQUS CAD module and meshed with quadratic tetrahedral (C3D10) elements with full integration. A Young's Modulus of 2.7 GPa and a Poisson's ratio of 0.49 were used for the simulations. A mesh-convergence study was performed to give a final mesh size of 38,660 elements and 152,682 elements for the octet and tetrakaidecahedron unit cell meshes, respectively. Periodic Boundary Conditions (PBC) [START_REF] Portela | Extreme mechanical resilience of self-assembled nanolabyrinthine materials[END_REF] were applied such that the relative displacement between opposite faces was constrained and an average strain was imposed on the unit cell. Linear perturbation was performed on the unit cells for 6 linearly independent strain vectors. The corresponding stresses for the applied strains were obtained and the compliance tensor Ŝ was constructed from these.

Figure 1

 1 Figure 1 Non-contact laser-induced excitation and probing of elastic waves in metamaterials. a, Schematic of the custom laser-induced resonant acoustic spectroscopy (LIRAS) method, employing a picosecond pump laser and a continuous-wave (CW) laser for interferometry displacement signals. SEM micrographs (insets) depict metamaterial architectures (octet, tetrakaidecahedron, and monolithic) used in our experiments. b, Schematic of pump-CW probe schemes to excite longitudinal, torsional, and flexural modes in metamaterials, accompanied by a depiction of resulting resonant modes. c, Representative waveforms obtained from longitudinal excitation of octet, tetrakaidecahedron, and monolithic samples. Scale bars in SEM insets, 20 µm.

Figure 3

 3 Figure 3 Dynamic elastic and damping properties of metamaterials. a, LIRAS (color) and quasistatically obtained (gray, uniaxial compression) effective stiffness E * depicting dynamic stiffening in samples. b, Rate-dependent stiffness of monolithic samples (circular IP-Dip pillars), tetrakaidecahedron, and octet [100] metamaterials depicting agreement of destructive uniaxial compression method with highthroughput LIRAS experiments. c, Spatial attenuation coefficient α as a function of frequency f , that follows α( f ) = β( f )/c p,L[START_REF] Garrett | Understanding Acoustics: An Experimentalist's View of Sound and Vibration[END_REF] , where β( f ) is the temporal damping coefficient and c p,L is the longitudinal phase velocity, depicting higher and anisotropic damping over our measurable frequency range as a power law α( f ) = α 0 f η with constant a 0 and tunable exponent η ranging from 1.35 to 1.51 (slope of the log-log plot) for octet and 1.39 to 2.02 for tetrakaidecahedron, compared to 1.58 for monolithic polymer.

Figure 4

 4 Figure 4 Wave propagation validation and omnidirectional elastic properties. a, Eigenfrequency analysis (hollow data points) on octet [100] 5×5×n tessellations presenting modal and frequency agreement with LIRAS experiments (solid data points). Shaded region depicts upper and lower bounds for material properties, hollow data points correspond to RMS-strain-rate material properties. b, Computational static elastic surfaces (top) compared to LIRAS-determined dynamic elastic surfaces (bottom) for the metamaterials, depicting anisotropic and architecture-dependent stiffening (center) under dynamic conditions, calculated as the difference between the static and dynamic surfaces.

Figure 5

 5 Figure 5 Vibrational defect identification. a, Octet 5 × 5 × 10 metamaterials with internal solid cubic inclusions (CAD and SEM), increasing from 0 to 4 inclusions per n tessellation. b, Measured shift in resonant frequencies as a function of inclusion-defect density using the center-pump LIRAS scheme. c, Octahedron 5 × 5 × 10 metamaterials with missing beams (CAD and SEM), increasing from 0 to 40% missing features. d, Measured shift in longitudinal and flexural resonant frequencies, denoting a linear decrease in frequencies as a function of defect density. Black scale bars, 20 µm; white scale bars, 10 µm.

ACCEPTED MANUSCRIPT / CLEAN COPY

ACCEPTED MANUSCRIPT / CLEAN COPY

Acknowledgements We acknowledge financial support from KCNSC (PDRD# 705288). C.M.P. acknowledges support from the National Science Foundation (NSF) CAREER Award [CMMI-2142460], partial support from DEVCOM ARL Army Research Office through the MIT Institute for Soldier Nanotechnologies (ISN) under Cooperative Agreement number W911NF-23-2-0121, and the MIT MechE MathWorks Seed Fund. This work was carried out in part through the use of MIT.nano's facilities. The authors thank A. Maznev, D.W. Yee, and K.A. Nelson for helpful discussions, T. Butruille for experimental support, R. Glaesener for imperfect-lattice design, and B. Aymon for explicit dynamics simulation support.

Data Availability

The data sets generated during and/or analysed during the current study are available within the main text or Supplementary Information.

Correspondence Correspondence should be addressed to C.M.P. or T.P. (emails: cportela@mit.edu, pezeril@mit.edu).

Supplementary information is available online at Nature journal.

Competing interests

The authors declare that they have no competing financial interests.


U Krupp
email: u.krupp@fh-osnabrueck.de

R Orosz

H.-J Christ
email: cchrist@ifwt.mb.uni-siegen.de

D Monceau

On the Mutual Interaction between Mechanical Stresses and Internal Corrosion during Isothermal and Cyclic Oxidation of Nickel-Base Superalloys

Keywords: Nickel-base alloys, creep, cyclic oxidation, internal oxidation, internal nitridation, embrittlement

Thermal cycling has been observed to cause a transition from superficial alumina formation to internal oxidation and nitridation, an effect that was shown to depend on the specimen thickness and geometry, which can be described by a spalling-probability model. Once protection by a dense and adherent alumina scale got lost, the internal-corrosion rate is determined by the diffusivity and solubility of nitrogen and oxygen in the alloy. These parameters seem to depend not only on the temperature and the alloy composition but also on the applied mechanical stress. Internal nitridation under a superimposed creep loading was found to follow a higher rate constant than under just isothermal exposure. This effect can probably be attributed to dislocation-pipe diffusion, a mechanism which has been claimed also to be relevant for outward solvent diffusion during internal corrosion, a phenomenon, which was observed as a stress-relief mechanism during various internal-reaction processes

Introduction

Nickel-based superalloys are designed to withstand loading conditions at high temperatures where static or cyclic mechanical stresses are superimposed by corrosive attack from the atmosphere. Structural-integrity concepts for systems that have to survive several ten thousands hours in service, e.g. land-based gas turbines or heat exchangers, require the extrapolation of creep, fatigue and hightemperature corrosion data, which are obtained usually from short-term tests under more severe conditions. Here, the focus is mostly being placed either on the mechanical response or on the surface degradation due to aggressive atmospheres. This may cause a non-conservative design for conditions, where mechanical stresses and the high-temperature-corrosion mechanism interact.

Even though, the general implications of mechanical stresses during the growth of oxide scales are known and reviewed, e.g., by Evans [1], the complex relationships between the phase stability, internal and external corrosion mechanisms and applied mechanical loads are hardly to predict in a quantitative manner. Prediction is made even more difficult, when taking into account that variations in the chronology of the occurrence of mechanical loading and corrosion attack might completely alter the overall damage process. At very high temperatures, solid-state diffusion of both the corrosive species and the metallic alloying elements is the key issue of high-temperature degradation. It was shown by O 18 tracer measurements on creep-loaded Ni and Zr (Zircaloy 4) [2,3] that tensile stresses lead to an increase in the diffusivities . The consequence is a change in the oxide-scale-growth mechanism as it was shown, e.g., by Calvarin-Almiri et al. [4] for chromiaforming Ni-20Cr foils. On the other hand, high-temperature-corrosion attack strongly deteriorates the mechanical properties of Ni-based superalloys. This was, e.g., shown by Dryepondt et al. [5] for a single-crystalline Ni-base superalloy. They attribute this degradation not only to the reduction in the load-carrying cross section by superficial oxide-scale formation and internal corrosion but also to an enhanced dissolution of the γ´ phase and the injection of vacancies into the surface.

Beside the effect of applied creep loading, the mechanical integrity of the oxide scale itself is of substantial significance for the long-term performance of high-temperature components. Often, thermal-cycling loading conditions lead to repeated cracking and/or spalling of superficial oxide scales. This kind of damage may result in a depletion of the oxide-forming element, being particularly pronounced in the case of thin-walled components or cracks [6,7]. Once the concentration of the oxide-forming element has fallen below a critical value, a strong acceleration of high-temperature corrosion can be observed due to the formation of an oxide of higher k p and the onset of massive internal oxidation and/or nitridation. This process is often referred to as "breakaway oxidation", cf. [8]. As it is also shown in the present paper, the relationship between the high-temperature-corrosion behavior and thermal-cycling conditions depends strongly on the specimen geometry and thickness, the cycle duration, and the kind of oxide being formed. This was the reason why a draft of a standard for cyclic oxidation testing was developed within the framework of the European project COTEST (http://cotest.dechema.de). However, reliable prediction of cyclic oxidation is hardly possible; probably the most promising approaches were suggested by Poquillon and Monceau [9] (probabilistic spalling model) and by Smialek [10] (deterministic interfacial oxidation spalling program, DICOSM) taking partial spalling of the superficial oxide into account. The situation becomes even more complex, if thermal cycling in combination with creep loading occurs. This was recently shown to strongly increase the steadystate creep rate in a superalloy, where a significant dissolution of the γ´ precipitates occurred during the high-temperature dwell times [11].

Once internal corrosion has set in, additional factors have to be taken into consideration: The higher specific volume of the internal precipitates as compared to the substrate causes the generation of residual compressive stresses in the surface layer. These stresses may (i) induce tensile stresses in the superficial oxide [12] or (ii) act as driving force for outward diffusion of the non-reacting solvent atoms by Nabarro-Herring creep and/or dislocation pipe diffusion [13,14].

For the example of two commercial Ni-base superalloys the present paper illustrates the complexity of mutual interactions between oxide-scale formation, internal corrosion and the occurrence of mechanical stresses and provides a summary of the aspects to be considered for physically-based concepts of service-life assessment.

Experimental

Studies on oxide-scale integrity and internal corrosion were carried out on two commercial Ni-base superalloys, the chromia-forming polycrystalline Alloy 80A and the single-crystalline aluminaforming alloy CMSX-4. The chemical compositions of them are given in Table 1. Both isothermal/cyclic thermogravimetric (TGA) measurements as well as creep tests were carried out in temperature ranges between T=850-900°C (creep) and T=1000°C-1100°C using selfdesigned systems, which are schematically represented in Fig. 1. To establish a low-p(O 2 ) nitriding atmosphere an electronically-controlled gas mixture (N 2 -45 Vol.-% He-5 Vol.-% H 2 is fed through a package of porous titanium sponge, which acts as an oxygen getter.

The superimposition of thermal cycles during continuous TGA is established by means of a specimen-lift device, by which the specimens (small polished discs or wedge-type specimens) are removed for defined time intervals from the hot furnace area into a cooling zone, and vice versa (cf. Fig. 1a).

Some of the creep tests were carried out on polished specimens with stepwise modified cross sections (cf. Fig. 1b) to allow a correlation between different values of creep strain and the corresponding internal-nitridation attack on one single specimen.

Nitridation and Mechanical Stresses during Isothermal Exposure with and without Applied Creep Loading

Exposure of Ti-and/or Al-containing Ni-base superalloys in nitrogen-based atmospheres (air or at very low p(O 2 )) results in the formation of internal nitrides. The introduction of additional phases with new interfaces gives rise to changes in the diffusion kinetics. Due to the higher specific volume of the nitrides as compared to the Ni-base matrix, compressive stresses are generated within the surface layer acting as driving force for an outward diffusive flux of the non-reacting matrix atoms. This kind of diffusion manifests itself in the formation of surface protrusions with the chemical composition of the non-reacting solvent on top of the surface, as it is shown in Fig. and as it was observed in several studies on internal oxidation and nitridation, e.g. [13,14]. Nabarro-Herring creep, dislocation-pipe diffusion and grain-boundary sliding were suggested to be most relevant for the diffusive flux but it turned out to be extremely difficult to distinguish between the contributions of bulk, interfacial or pipe diffusion and to estimate the flux-blocking effect of the precipitates themselves.

The stresses generated by the internal precipitation process can be superimposed by applied creep loading. For the case of internal nitridation it was shown that the higher the creep strain or stress the higher is the internal nitridation depth (cf. Fig. 3) [15,16]. A pronounced coverage of the specimens by surface protrusions was never observed, hinting to the fact that the residual compressive stresses are compensated by the tensile creep stress. With respect to creep damage, so far no beneficial effect of the nitrides was observed. Rather, the preferential precipitation of nitrides along the alloy grain boundaries promotes the occurrence of brittle intergranular cracking.

It is worth mentioning that a specimen of Alloy 80A crept under vacuum conditions at T=900°C at 30MPa reached a strain value of 14.6%, while the specimen crept under the same conditions but in nitriding atmosphere failed in a more brittle manner (N 2 -45 Vol.-% He-5 Vol.-% H 2 + Ti sponge) at a strain value of 8.6%. In agreement with Pieraggi [17] it can be concluded that the deteriorating effect of high-temperature corrosion attack (in this case internal nitridation) is not restricted to the corrosion-affected surface layer.

Formation of surface protrusions as a consequence of internal precipitation of TiN in Ni-

10Cr-2Ti during exposure to nitriding atmosphere for 150h at 1100°C (ref. [13]).

Fig. 3. TiN penetration depth vs. creep strain after creep testing of alloy 80A (Nicrofer 7520 Ti) in oxygen-free nitrogen atmosphere (ref. [15]).

Geometrically-Induced Mechanical Stresses during Cyclic Oxidation

High temperature exposure under cyclic-loading conditions causes the formation of thermal stresses in the oxide scale/substrate surface layer system due to strongly different coefficients of thermal expansion (CTE). The response of the system depends on (i) the scale adherence on the substrate surface, (ii) the scale coherence, and (iii) the ability of the scale/substrate to accommodate the stresses by creep. Furthermore, the component geometry and thickness as well as the cycle duration and the cooling rate are of high significance for the long-term behaviour of the oxide scale.

Figure 4a shows that in the case of Cr 2 O 3 -scale formation on Alloy 80A thermal cycling leads to an increased mass gain, which corresponds to an oxidation process that is accelerated by the presence of a crack network [START_REF] Krupp | Proc. EFC-Workshop Life Time Modelling of High-Temperature Corrosion Processes[END_REF]. The situation is different in the case of Al 2 O 3 formation on CMSX-4 (Fig. 4b): Here, strong cyclic spalling was observed. One would expect that the extent of spalled oxide depends strongly on the cycle duration; the longer the high-temperature dwell time the thicker is the oxide scale and the higher are the stresses within the scale/substrate interface. Surprisingly, the spalled gross mass is more or less a function of the number of cycles, the mass change for 20 cycles with 1h dwell at 1100°C is almost the same than for 20 cycles with 5h dwell at 1100°C (Fig. 4b). Finally, Fig. 4c shows the effect of the specimen geometry. TGA under thermal-cycling conditions on wedge-type specimens (cf. Fig. 1a and [19,20]) results in an enhanced spalling behaviour as compared to a conventional rectangular specimen. The effect of the geometry on the extent of cyclic spalling behaviour becomes obvious from the cross section micrographs in Fig. 5. Isothermal exposure results in the formation of an adherent alumina scale without showing any indication of spalling, cracking, Al depletion, and internal oxidation, even at the very sharp tip of the wedge-type specimen (Fig. 5a). However, exposure under thermal-cycling conditions causes a transition from external to internal oxidation and nitridation. Obviously, cyclic spalling leads to a depletion in the Al concentration in the specimen below the critical value required for superficial scale formation. Eventually, the formation of internal oxides and nitrides causes stresses in the thin-walled area which are accommodated by cracking.

The results show clearly that components or materials of low wall thickness and/or sharp edges exhibit a particular susceptibility to "breakaway oxidation" in combination with massive internal attack. further set of thermogravimetric experiments were carried out on specimens with different thicknesses ranging from t=100µm to t=500µm (cf. [21]). The results are summarized in Fig. 6. There is a clear trend for thick specimens to exhibit stronger oxide spallation (Al 2 O 3 ). Any kind of "breakaway oxidation" or massive internal corrosion could not be identified. Also, the depth of the dissolution of the γ´ phase (Ni 3 Al) due to oxidation-induced Al consumption was more or less the same for all specimens. Hence, the thickness effect on the spalling behaviour must be attributed to in the stress-state within the scale-substrate system. In agreement to [22], thin specimens can accommodate stress peaks in the oxide substrate interface by creep with the consequence that the spalling rate is strongly reduced. The thickness effect can be extrapolated by means of the probabilistic spalling model (p-k p model) suggested by Poquillon and Monceau [9]. The results are shown by the solid lines in Fig. 6. The model accounts for a change in the amount of spalled oxide volume, which was also found experimentally by evaluating the specimen surfaces after certain thermal-cycling intervals, assuming that the nature of the oxide scale remains unaffected. The results summarized in the preceding sections illustrate strong interactions between mechanical stresses and the external and internal high-temperature corrosion behaviour. To implement these effects into service-life assessment methods, it is required to use (i) stress terms in the diffusion differential equations to account for superimposed creep loading and for residual stresses arising from internal precipitation, (ii) quantitative failure criteria for the different types of superficial oxide (cracking/spalling), (iii) the mechanical-stress and the alloying-element distributions in thin walled or sharp-edged specimens to account for depletion or stress-relief effects, and (iv) a quantitative description of phase stability within the complete system (substrate and corrosion products) by means of the CALPHAD approach in combination with a numerical treatment of multi-element diffusion.

Fig. 1 .

 1 Fig. 1. Schematic representation of (a) the thermobalance used for TGA under isothermal and thermal cycling conditions on cylindrical and wedge-type specimens and (b) the creep machine for tests in defined gas atmospheres.

Fig. 4 .

 4 Fig. 4. Thermogravimetrically measured mass change vs. exposure time for (a) oxidation of Alloy 80A in air at 1000°C (b) oxidation of CMSX-4 in air at 1100°C and (c) oxidation of rectangular and wedge-type specimens of CMSX-4 in air at 1100°C (thermal cycling: intervals of 1h/5h dwell at high temperature followed by 15min. cooling at 50°C).

Fig. 5 .

 5 Fig. 5. Sharp tip of wedge-type specimens of CMSX-4: (a) after 100h isothermal exposure at 1100°C to air and (b) after 100h exposure under thermal-cycling conditions at 1100°C to air.

Fig 6 .

 6 Fig 6. Measured (data points) and simulated (solid lines, p-k p model [9]) net mass change vs. number of thermal cycles for exposure of CMSX-4 specimens different thickness in air at 1100°C (1h dwell times, 15min cooling at 50°C).

Table . 1

 . Chemical composition of the Ni-base superalloys Alloy 80A (Nicrofer 7520 Ti) and CMSX-4 (in wt.%)

	Ni	Cr	Al	Ti	Co	Ta Mo	W	Re Hf	S	C
	Alloy 80A (75.2) 20.2	1.6	2.7 0.01	--	--	--	--		20ppm	0.045
	CMSX-4 (67.1) 6.0	5.6	1.0 10.0 6.0	0.6	6.0	3.0	0.1 <12ppm

Acknowledgments

The financial support by Deutsche Forschungsgemeinschaft (DFG) and the supply with the test materials Alloy 80A (Nicrofer 7520 Ti) by Krupp VDM, Altena, Germany, and CMSX-4 by Alstom, Baden, Switzerland, is gratefully acknowledged.


N Di Miceli

S Elgue

L Prat

P Cognet

M Cabassud

Optimisation of operating conditions in batch for more sustainable continuous process transposition

Keywords: Optimisation, exothermic and fast reactions, process intensification

Implementation of exothermic and fast reactions requires a perfect temperature control to avoid thermal runaway and in most cases to limit by-products production. In order to fit with the heat removal capacity of common devices, expensive strategies are currently used to slow down this kind of reactions in order to avoid a strong temperature increase such as reactants dilution. Within the concept of process intensification, industries could move towards more sustainable process by reducing technology constraints to the benefit of chemistry. For that purpose, a two-step methodology is implemented. The first step consists in the optimisation of the operating conditions only based on stoichio-chemical scheme and kinetic laws. This is carried out by adjusting temperature profile and feeding rate strategy in a batch operation. Then a design for a continuous process is proposed, trying to approach the optimal batch operating conditions. This methodology is applied to the linear alkylbenzene sulfonation.

Introduction

Fast and exothermic reactions are difficult to carry out at industrial scale because their processing requires a perfect temperature control. In conventional devices, the heat transfer capacity is often limited as shown in Table 1: while the jacketed batch reactor offers a very poor thermal efficiency, it is still used in many manufactures to carry out fast and exothermic reactions. In such reactors, the temperature control is insured by diluting the reactants, or by implementing a feeding strategy which generally does not favour the selectivity of the reactions and lead to an increase in the residence time. This procedure tends to adapt the chemistry to the technology constraints. In agreement with the concept of process intensification [START_REF] Stankiewicz | Process intensification: transforming chemical engineering[END_REF], the current trend is to adapt the technology to the chemistry. The emergence of highly efficient technologies in terms of heat removal (heat exchanger reactors, microreactors) allows a better fitting of the technology with the chemical constraints.

Table 1. Heat exchange capability for different reactors [START_REF] Ferrouillat | Open loop thermal control of exothermal chemical reactions in multifunctional heat exchangers[END_REF]. This can be illustrated by the sulfonation, which is a key reaction in the detergent production. Sulfonates have a white colour and are mainly used as surfactants in laundry, while by-products are black: consequently, sulfonation manufacture aims at the purest product to satisfy consumers' requirements. Processes using gaseous sulfur trioxide (SO 3 ) as sulfonating agent are predominant since there are very versatile (a wide variety of feedstocks can be sulfonated or sulfated) and they produce high quality products with a purity of more than 95 % [START_REF] Foster | Sulfonation and sulfation processes[END_REF]. Because by-products formation is mainly favoured by high temperature, temperature has to be perfectly controlled in the reactor. This makes the sulfonation processing complex since the main reaction and the secondary ones are very fast and exothermic. Moreover, the viscosity of the sulfonate solution is between 15 and 300 times higher than the feedstock one, leading to a drastic decrease in the local heat transfer coefficient of the liquid phase as the reaction proceeds. This reaction is generally carried out in a falling film reactor coupled with a cooling water jacket where the cooling area is around 3 000 m 2 /m 3 and the global heat transfer coefficient around 800 W/(m 2 K) (Chemithon annular falling film reactor, Ballestra multi-tube film sulfonator [START_REF] Foster | Sulfonation and sulfation processes[END_REF]. Therefore this technology is characterized by a high thermal efficiency (2 400 kW/(m 3 K)). Nevertheless, this thermal efficiency added to the small thickness of the liquid are not sufficient to efficiently dissipate the heat generated by the reactions, the heat of reaction being of the order of -170 kJ per mole of sulfonate produced. Consequently, the increase in temperature is currently limited by diluting the gaseous reactant SO 3 with air (volumetric fraction in SO 3 between 3 % and 7 % (de Groot, 1991)). However, air needs to be dried before being in contact with SO 3 to avoid the formation of sulfuric acid which can corrode the process devices. Therefore, since an important flow rate of air is used, the air/SO 3 process is highly energy consuming. In order to avoid air consumption, other reactor technologies have to be considered to enhance heat transfer between the reaction mixture and the cooling system.

In this article a methodology to adapt process designs to chemical constraints is presented, illustrated by the linear alkylbenzene (LAB) sulfonation with SO 3 . In order to suppress the impact of the transfer of SO 3 from the gas phase to the liquid mixture, this study focuses on LAB sulfonation with liquid SO 3 . This way, the present work concentrates only on the chemistry of the system. This methodology is developed in two steps: first, the optimal operating conditions are determined only considering the thermo-kinetic parameters of the reaction in batch conditions. Then a reactor design is proposed trying to approach the optimal operating conditions in a continuous reactor.

A fast and exothermic reaction: the LAB sulfonation

Reaction scheme

Since sulfonation is very fast, very few data on the reaction kinetic are available. Indeed, previous works used to model this process assuming the kinetic rate instantaneous compared to the rate of SO 3 mass transfer from the gas phase to the organic medium [START_REF] Johnson | Modeling of a thin-film sulfur trioxide sulfonation reactor[END_REF][START_REF] Davis | An analysis of the falling film gasliquid reactor[END_REF]. This approach has notably been validated in falling film reactors by [START_REF] Talens-Alesson | The modelling of a falling film chemical reactors[END_REF]. However, this work aims to propose a method of process design optimisation to the benefit of the chemistry. For that purpose, chemical kinetic is the only contribution that we consider in this work to model the sulfonation process. As illustrated in Figure 1, a reaction scheme is proposed to simulate and optimise the LAB sulfonation, inspired from Roberts's works (2003):

(1) The first reaction R1 allows the production of a small amount of pyrosulfonic acid (Pyro) required to prime the main mechanism of sulfonation.

(2) The main mechanism of the LAB sulfonation is composed of two reactions (R2 and R3) and leads to the formation of the desired product: the linear alkylbenzene sulfonate (LAS). (3) Numerous competitive side reactions were suggested by [START_REF] Roberts | Optimization of the linear alkyl benzene sulfonation process for surfactant manufacture[END_REF]; to simplify the simulation of the sulfonation, only two side reactions have been considered: the production of the linear alkyl benzene sulfonate anhydride LASA (R4) and of the pyroanhydride LASP (R5). Both reactions also lead to the formation of sulphuric acid H 2 SO 4 . This scheme suggests that both reactants should be introduced in the reactor in equimolar ratio. Indeed, in one hand an excess in SO 3 would favour the consumption of the desired product, LAS, to form pyrosulfonic acid through R3. On the other hand, an excess of LAB does not seem to influence the selectivity of the sulfonation. Therefore, it would lead to a poor conversion of this reactant while the reaction scheme does not justify it.

Thermo-kinetic parameters

The thermo-kinetic parameters used to represent the LAB sulfonation are based on the following considerations (Table 2): (1) [START_REF] Roberts | Optimization of the linear alkyl benzene sulfonation process for surfactant manufacture[END_REF] suggested to represent the self-priming reaction R1 as first order in LAB and second order in SO 3. For the other reactions, we assumed that the order of the reactions could be deduced from the stoichiometric coefficients of each reactant as described in Table 2.

(2) The rate constant k can be written in term of Arrhenius law. Roberts estimated the pre-exponential factor k 0 regarding to the number of reactant molecules and product in agreement with the kinetic theory of gases [START_REF] Frost | Kinetics and Mechanism[END_REF]. For instance, R2 being a reaction in which two reacting molecules form two product molecules, he proposed a pre-exponential factor of 10 11 L/(mol s).

The same assumption allows the estimation of k 0 for reactions R4 and R5. For reactions where two reacting molecules form one product molecule such as R3, Roberts suggested a pre-exponential value of 10 8 L/(mol s). Finally, for the self-priming reaction R1 which is a third-order reaction, he estimated k 0 around 10 4 L 2 /(mol 2 s).

(3) Roberts gave 54 300 J/mol as the estimated activation energy of R2 (production of LAS from LAB and the pyrosulfonic acid). Then, he suggested that the activation energy of R3 should be very low since the chemical mechanism to form the purosulfonic acid from LAS and SO 3 might be rather "facile". Therefore, we arbitrary set activation energy lower than the main sulfonation reaction one. Finally, Roberts estimated the activation energy of the self-priming reaction R1 by considering the enthalpy of formation of the probable intermediate products that lead to the formation of the pyrosulfonic acid from LAB and SO 3 . He assumed it ranged from 8 500 J/mol to 12 500 J/mol. (4) By-products production is favoured by high temperature; therefore the activation energy of the side reactions R4 and R5 must be higher than the others. We considered that E a of R4 and R5 were 20 % higher that the one estimated for R2. (5) Roberts estimated the enthalpies of reaction for the main reactions from the enthalpies of formation of the different products by means of calculations of molecular orbital energies. (6) The reaction enthalpies of the side reactions should be very small assuming that these data can be deduced from the dissociation enthalpies of the broken and formed bonds during the reaction. In both R4 and R5, two bonds Sulphur -Oxygen are broken and then formed. The evolution over time of the variables that describe the reactor state such as temperature and composition is calculated by integrating mass and energy balances in respect with constraints. To solve the problem, the residence time t is divided into a finite number of intervals n int equals to 5 in this work. The optimisation algorith allows the determination of the (n int + 1) temperatures at each bond of each interval. The temperature profile is then assumed to be linear in each interval:

        - - ⋅ - + = - - - - 1 j j 1 j j 1 j 1 j ) ( t t T T t t T T with 1 < j < (n int + 1) (1)
For the feeding profile, the program determines the mass that should be introduced in each interval. The flow rate is then obtained by assuming it remains constant over each interval. [START_REF] Garcia | Constrained optimization for the fine chemicals production in batch reactors[END_REF] described in more details the numerical formulation of the problem and validated the method with general exemples of reaction schemes (consecutive-competitive and parallel reactions).

In the present work, the criterion to be optimised is the final amount ratio of LAS x LAS using pure reactants (LAB and liquid SO 3 ). In our study, the molecule of LAB is composed of a radical R of molecular formula C 12 H 25 (dodecylbenzene). The residence time is set to one second. The operating conditions to be defined are the temperature profile in the reactor and the SO 3 feeding profile (c SO3 = 10.6 mol/L). We decided to feed the reactor with SO 3 rather than LAB since an excess of SO 3 would dramatically affect the selectivity of the sulfonation regarding the reaction scheme considered in this study (Figure 1). Therefore, at initial time, the reactor contains only LAB (c LAB = 3.46 mol/L). It is assumed that SO 3 is totally and instantaneously mixed in the organic matter. The physico-chemical properties of the mixture in terms of density and heat capacity are supposed constant during the reaction, calculated at the initial time. Finally, the temperature in the reactor is constrained: it must not be lower than 25 °C and higher than 150 °C to maintain the solution in liquid phase (molecules that contain long carbon chain may freeze at 20-25 °C; the decomposition temperature of the products is arbitrary set to 150 °C).

Optimisation results

First, we determined the optimal operating conditions for a discontinuous operation (batch or semi-batch) without constraining the heat removal capability of the reactor. This optimisation supposes that the temperature in the reactor can be perfectly controlled whatever the heat that needs to be removed or provided to the reactive system. The result of the optimisation computation is a semi-batch reactor at minimum temperature to achieve 95% of the final amount of LAS, which corresponds to a residence time of 0.4 s (Figure 2Erreur ! Source du renvoi introuvable.). Then a higher conversion is reached with an increase of the temperature (x LAS = 92.06 %). As expected, the total amount of SO 3 introduced in the reactor is in equimolar ratio with the initial amount of LAB. However, this simulation shows that a lot of heat needs to be removed from the reactor in a hundred seconds to reach such a conversion as illustrated on Figure 3 (2•10 6 kW/m 3 of reactor). It corresponds to the heat flow rate produced by the reaction since the temperature of the mixture is constant during this period. We can compare the maximal heat flow rate necessary to optimise the sulfonation with the thermal powers of different devices:

(1) Considering a new concept of continuous heat exchanger reactor, the "open plate reactor" designed by Alfa Laval Vicarb (OPR), [START_REF] Prat | Performance evaluation of a novel concept "open plate reactor" applied to highly exothermic reactions[END_REF] and [START_REF] Benaissa | Développement d'une méthodologie pour la conduite en sécurité d'un réacteur continu intensifié[END_REF] proved its high thermal efficiency (heat transfer coefficient up to 2500 W/(m 2 K) for a volumic area of 240 m 2 /m 3 ). With this reactor of 0.4 L, they decreased the temperature of hot water of 45 °C by exchange with cold water (15 °C). Therefore, this device can be characterized by a thermal power of 9 000 kW/m 3 which is very far from the 2•10 6 kW/m 3 of reactor required to obtain maximal conversion. (2) To approach so high thermal power, micro-heat exchangers seem to be appropriate apparatus to carry out the reaction since they offer a high exchange area for a very small reactive volume. Indeed, McGovern et al. (2006) calculated that microreactors can offer an area per volume unit ratio of more than 20 000 m 2 /m 3 compared to roughly 500 m 2 /m 3 for usual device. [START_REF] Schubert | Microstructure devices for applications in thermal and chemical process engineering[END_REF] obtained thermal powers up to 20 kW with microstructured devices having 1 cm 3 active volume, i.e. up to 2.10 7 kW/m 3 .

The previous results shows that technology constraints to carry out the sulfonation reaction can be reduced to the benefit of the chemistry (use of pure reactants) by using highly efficient devices in terms of thermal power such as microreactors.

Optimal operating conditions were determined in semi-batch reactor (temperature and feeding profiles): the transposition to continuous operation is now required to fit with the use of miniaturised structures.

4. Process design optimisation for a continuous operation 4.1. Geometry of the reactor

The previous part demonstrates that sulfonation manufacture is strongly constrained by the capability of the reactor to dissipate the heat generated by the reactions. Since micro-heat exchangers are considered as highly efficient apparatus [START_REF] Qu | Heat transfer for water flow in trapezoidal silicon microchannels[END_REF][START_REF] Gao | Scale effects on hydrodynamics and heat transfer in two-dimensional mini and microchannels[END_REF][START_REF] Gamrat | Conduction and entrance effects on laminar liquid flow and heat transfer in rectangular microchannels[END_REF], sulfonation is simulated in a microchannel coupled with a cooling system as illustrated in Figure 4: the cooling of the reactive media is insured by thermal exchange with an utility fluid that flows around the process channel. The simulated structure is in steel and the size of the channels and the walls between the process and the utility channels are given in Table 3. This geometry is inspired by a microstructured reactor designed by the Forschungszentrum Karlsruhe (FZK, Germany), where the process passage consists of 20 foils with 34 microchannels of such geometry. The length of the reactor is estimated from the desired residence time which is set to one second to compare the results of the simulations with those obtained with Optibatch. The channel has roughly 25 000 m 2 /m 3 of volumic area. The inlet flow rates are estimated so that the fluid velocity in the channel is around 3 cm/s. Pure LAB is introduced at the inlet of the reactor (c LAB = 3.46 mol/L ) while SO 3 is injected at several points in the main channel to approach the semi-batch reactor behaviour (c SO3 = 10.6 mol/L). In order to approach the results obtained with Optibatch, the total amount of SO 3 injected is in equi-molar ratio with the LAB introduced at the inlet of the channel. The cooling fluid is water, with a volume flow rate that is approximately twice the process fluid one. In such conditions, the temperature of the utility fluid remains roughly constant. It assumes that no heat is exchanged with the ambient environment, at room temperature. Table 4 summarizes the properties of the inlet streams. Different multi-injection strategies are compared in terms of final amount of LAS.

Simulator description

The simulator is based on the integration of a dynamic model. It has been described and validated by [START_REF] Elgue | Optimal operations of a novel intensified reactor based on dynamic modelling[END_REF] through the simulation of sodium thiosulfate oxidation in the OPR. The different parts of the reactor (process and utility channels) are represented by series of perfectly stirred tank reactors called cells. Since mixing is supposed to be very fast compared to the reaction time, we considered that the flow in the channel can be represented by a perfect plug flow, which is equivalent to a serie of a high number of stirred tanks, set to 100 in this work. The modeling of a cell is based on the resolution of mass and energy balances and constraint equations (geometrical characteristic of the reactor, physical properties of the fluids...). This resolution allows the calculation of the parameters characterizing the flows: temperature, flow rates, composition, pressure.... The balance equations and the correlations used to model the system require the estimation of the physical properties of each fluid. For that purpose, the model implemented has been linked to a commercial thermodynamic and equilibrium database: Bibphy32 (Prosim ®). The local heat transfer coefficients on the utility and process sides are calculated considering that Nusselt number is constant in microchannels, equal to 2.181 [START_REF] Shen | Flow and heat transfer in microchannels with rough wall surface[END_REF].

The mathematical model that represents the operation is an extremely complex model due to the number and also the type of equations that it involves: differential and algebraic equations (global and component mass balances, energy balances within the process and utility fluids, pressure balance, …). The resolution of the model calls for a specific strategy and solving method. The system is solved using a global approach (i.e. simultaneous treatment of all equations with no discrimination according to the type of equation) by means of the specifically designed solving tool based on the Gear method: DISCo [START_REF] Sargousse | DISCo : un nouvel environnement orienté objet[END_REF].

Results and discussions

Results of the simulations done with different injection configurations are summed up in Table 5. Figure 5 shows that the more numerous the injection points are, the best the yield of LAS is. However, the number of injection channels of the optimal device will be limited by manufacture constraint. Figure 6 illustrates that the injections should occur during 0.4 seconds in terms of residence time: these results fit with those previously obtained despite the discretisation of feeding.

Conclusion

The methodology presented aims to optimise a process design in two steps according to the process intensification principle. First, the optimal operating conditions are determined in a batch or semi-batch reactor only considering kinetic and thermal constraints. The second step consists in defining an optimal geometry, suggested by the results previously obtained, to carry out the reaction in a continuous reactor. All the parameters that could influence the optimisation have not been thoroughly studied (utility flow rate, channel dimensions, wall properties) but this work aims to demonstrate the capability of such a methodology of optimisation, and the performances of the associated tools we used.

In the case of the LAB sulfonation the application of this methodology reveals that the reaction should be carried out in a micro-exchanger. Indeed this technology allows a good heat removal and so a good temperature control while implementing a multi-injections strategy of pure SO 3 in the process channels. Mixing limitations have not been brought up in this study but [START_REF] Hessel | Micromixers -a review on passive and active mixing principles[END_REF] reviewed numerous ways to enhance mixing in micro-device. [START_REF] Kockmann | Convective mixing and chemical reactions in microchannels with high flow rates[END_REF] tested different mixers based on successions of T-mixers and obtained a characteristic mixing time between 0.1 and 10 ms depending on the Reynolds number in the channel. For reaction processing, additionally to mixing elements, the size of the main and injection channels could be designed in order to create turbulences where the fluids get in contact. The reactants flow rates could also have a major effect on mixing and should be properly chosen to obtain a mixing time very small compared to the characteristic reaction time.

Figure 1 .

 1 Figure 1. Reaction scheme proposed to model the LAB sulfonation.

Figure 2 .

 2 Figure 2. Result of the optimisation: (a) amount ratio profile of LAB, SO 3 and LAS; (b) optimal temperature and feeding profiles obtained for the production of one mole of LAS.

Figure 3 .

 3 Figure 3. Optimal heat flow rate per volume of reactor.

Figure 4 .

 4 Figure 4. Scheme of the microchannel used for the simulation of the LAB sulfonation in a continuous reactor.

Figure 5 .

 5 Figure 5. Comparison between different process designs varying the number of injections in 0.4 s: (a) yield of LAS profile; (b) temperature profile.

Figure 6 .

 6 Figure 6. Comparison between different process designs varying the space interval between the injections: (a) yield of LAS profile; (b) temperature profile.

Table 2 .

 2 Thermo-kinetic parameters used to model the LAB sulfonation. LicenseINPT, 1998) is a software developed for the optimisation of operating conditions in a batch or semi-batch reactor only based on the chemistry of a process (no mass transfer considerations). The objective of this tool is to optimise a criterion defined in terms of yield, selectivity, time of reaction, …, while taking into account different constraints specified by the operator: technology and/or safety constraints, productivity and/or selectivity constraints.

	Reaction	Reaction rate r (mol/(L s))	Pre-exponential factor k 0	Activation energy E a (J/mol)	Reaction enthalpy ∆H (kJ/mol)
	R1 R2 R3 R4 R5	2 k c LAB c SO3 k c LAB c Pyro k c LAS c SO3 k c LAS c Pyro k c Pyro 2	10 4 L 2 /(mol 2 s) 10 11 L/(mol s) 10 8 L/(mol s) 10 11 L/(mol s) 10 11 L/(mol s)	10 500 54 300 45 000 65 000 65 000	-70 -140 -70 0 0
	3. Operating conditions optimisation
	3.1. Description of the optimiser Optibatch

Table 3 .

 3 Geometric characteristics of the simulated microchannel.

		Process channel	Utility channel	Wall
	Length L (mm)	33	33	33
	Depth d (µm)	110	55	90
	Width w (µm)	320	320	320

Table 4 .

 4 Properties of the inlet streams in the simulated microchannel.

		Process channel	Injections	Utility channels
	Entity	LAB	SO 3	water
	Flow rate q V (µL/s)	1.11	0.16 (in total)	1.39 (per channel)
	Temperature T (°C)	25	25	25
	Pressure p (bar)	5	5	5

Table 5 .

 5 Result of the simulations with different multi-injection strategies (residence time=1s).

	Injections	Hot spot	Yield of LAS
	1 in 0.4 seconds	123.64°C	86.91 %
	4 in 0.4 seconds	53.15°C	92.37 %
	10 in 0.4 seconds	43.85°C	92.62 %
	10 in 0.65 seconds	40.99°C	92.08 %
	10 in 0.9 seconds	40.04°C	91.11 %
	10 in 0.2 seconds	51.08°C	92.50 %