P Jonnard

J P Zhang

C Pendenque

K Le Guen

R Delaunay

I Vickridge

D Schmaus

Q G Fu

Kossel interferences of proton-induced X-ray emission lines to study thin film waveguides

Keywords: X-ray waveguide, thin film, PIXE, Kossel diffraction, X-ray reflectometry

  L'archive ouverte pluridisciplinaire

Introduction

Kossel interferences, evidenced by Walther Kossel in 1935 [START_REF] Kossel | Die Richtungsverteilung der in einem Kristall entstandenen charakteristischen Röntgenstrahlung[END_REF], are the typical diffraction patterns of the characteristic X-rays emitted within a crystalline structure and sensitive to the structure itself. Since 1970s, it has been applied as a tool to investigate the crystallographic orientations and lattice parameters [START_REF] Tixier | Kossel patterns[END_REF][START_REF] Martínez-González | Directed self-assembly of liquid crystalline blue-phases into ideal single-crystals[END_REF][START_REF] Langer | Lattice constant determination from Kossel patterns observed by CCD camera[END_REF]. The observation of Kossel interferences requires two steps. Firstly, to produce the characteristic X-ray emissions, an ionization source is necessary to generate the electron transitions. Then the produced X-ray emissions are diffracted by the periodic structure itself.

The ionization source can be energetic electrons [START_REF] Jonnard | Soft-x-ray Kossel structures from W/C multilayers under various electron ionization conditions[END_REF][START_REF] Aiken | Elemental analysis of organic species with electron ionization high-resolution mass spectrometry[END_REF][START_REF] Goldstein | Scanning electron microscopy and X-ray microanalysis[END_REF], X-ray photons [START_REF] Rousse | Production of a keV X-ray beam from synchrotron radiation in relativistic laser-plasma interaction[END_REF][START_REF] Jonnard | Spontaneous soft x-ray fluorescence from a superlattice under Kossel diffraction conditions[END_REF], charged particles including protons or ions [START_REF] Lucarelli | The upgraded external-beam PIXE/PIGE set-up at LABEC for very fast measurements on aerosol samples[END_REF][START_REF] Subercaze | High energy PIXE: A tool to characterize multi-layer thick samples[END_REF]. Kossel diffraction associated to X-rays is related to the grazing exit X-ray fluorescence (GEXRF) technique [START_REF] Tixier | Kossel patterns[END_REF][START_REF] Guen | Kossel Effect in Periodic Multilayers[END_REF].

GEXRF measurement representing the angular resolved X-ray fluorescence technique, is sensitive to the elemental depth distributions in the nanometer range [START_REF] Baumann | Laboratory Setup for Scanning-Free Grazing Emission X-ray Fluorescence[END_REF][START_REF] Szwedowski-Rammert | Laboratory based GIXRF and GEXRF spectrometers for multilayer structure investigations[END_REF]. It can be used as a complementary characterization method of X-ray reflectometry (XRR).

In XRR measurements, the spectra are acquired by varying the incident angle in the grazing incidence regime while detecting the specular reflected X-ray beam. This technique uses the average electron density of a layer and thus cannot provide information on the elemental composition. By the combination of XRR and XRF analysis, we could overcome the limitations of XRR when dealing with elemental distributions and lateral inhomogeneities of multilayers [START_REF] Gil | Limitations of X-ray reflectometry in the presence of surface contamination[END_REF][START_REF] Spear | Advances in X-ray Reflectivity (XRR) and X-ray Fluorescence (XRF) Measurements Provide Unique Advantages for Semiconductor Applications[END_REF][START_REF] Ingerle | Combination of grazing incidence x-ray fluorescence with x-ray reflectivity in one table-top spectrometer for improved characterization of thin layer and implants on/in silicon wafers[END_REF][START_REF] Yang | Structural characterization of multilayer using the analysis combining GIXRF with GIXRR method[END_REF][START_REF] Hönicke | Characterization of ultra-shallow aluminum implants in silicon by grazing incidence and grazing emission X-ray fluorescence spectroscopy[END_REF]. For example, P.

Hönicke et al. characterized the ultra-shallow Al-implanted silicon wafers by the combination of XRR and GIXRF [START_REF] Hönicke | Characterization of ultra-shallow aluminum implants in silicon by grazing incidence and grazing emission X-ray fluorescence spectroscopy[END_REF].

In our previous works [START_REF] Guen | Kossel Effect in Periodic Multilayers[END_REF][START_REF] Wu | Kossel interferences of proton-induced X-ray emission lines in periodic multilayers[END_REF][START_REF] Wu | Kossel diffraction observed with X-ray color camera during PIXE of nano-scale periodic multilayer[END_REF], we develop a methodology, that is particle-induced X-ray emission (PIXE) combined with Kossel diffraction, to study the structure of nano-scale Cr/Sc and Pd/Y periodic multilayers. PIXE analysis possesses a good sensitivity to measure the radiations emitted by electron state changes and distinguishes the elements of the sample in a wide range of atomic numbers, without destruction of the sample [START_REF] Subercaze | High energy PIXE: A tool to characterize multi-layer thick samples[END_REF][START_REF] Midy | Application of a new algorithm to depth profiling by PIXE[END_REF]. In Ref. [START_REF] Guen | Kossel Effect in Periodic Multilayers[END_REF][START_REF] Wu | Kossel interferences of proton-induced X-ray emission lines in periodic multilayers[END_REF][START_REF] Wu | Kossel diffraction observed with X-ray color camera during PIXE of nano-scale periodic multilayer[END_REF], the generated X-ray emission is recorded in Kossel geometry to obtain the so-called Kossel curve, which is the intensity of characteristic X-ray emission as a function of the detection angle. Kossel curves are convenient to study the interfacial environment in multilayers made of a periodic alternation of nanometer-thick films.

Compared with other ionization excitation sources, electron excitation for example, the main advantage of proton excitation is the enhanced sensitivity due to substantial reduction of the primary beam bremsstrahlung [START_REF] Johansson | Particle-induced X-ray emission spectrometry (PIXE)[END_REF]. This makes it possible to obtain reliable Kossel curves from small amounts of material: either very thin layers or from elements present in low concentrations. A second advantage is that, within the thickness of typical structures investigated by PIXE-Kossel, proton trajectories within the structure are rectilinear/straight. The target atoms may be assumed to be uniformly ionized throughout the thickness of the films, which could possess well-known ionization cross sections and x-ray emission probabilities. The sensitivity of the Kossel curve to sample structure (depth resolution) depends only on the X-ray transport properties of the structure: there is no intrinsic influence of the excitation source here. However, according to the comparison of measured/simulated Kossel curves, the noisiness of the observed curve will influence the confidence with which different possible sample structures can be differentiated/discriminated. In this case, compared with electron excitation, the lower primary beam bremsstrahlung and the "cleaner" X-ray spectra obtained upon proton irradiation are also clear advantages.

The above merits significantly contribute to simplify the analysis of Kossel curves.

Let us note that X-rays can also be used as a source with an excitation depth that can be up to 1 mm. Such an excitation also produces a low bremsstrahlung radiation. Moreover, as in this case the technique is photon-in photon-out, it can be applied on insulating sample, for example on multilayers deposited on silica substrates. Instead of multilayers, generally made of a few tens of layers, the aim of this study is to check the feasibility of this method when used in non-periodic stacks consisting in only a few layers. To observe well-contrasted Kossel curves, the stack is an x-ray planar waveguide formed of three layers. Waveguides are designed to locate and concentrate the electric field of the emitted radiation in a specific layer of the stack, thus being sensitive to the depth distribution of the emitting atoms. Here, Pt/Fe/Pt and Ta/Cr/Pt samples are designed, deposited on Si and then investigated by the PIXE-Kossel method. By comparing the experimental and simulated Kossel curves in order to obtain the description of the actual stack, the possible application of the PIXE-Kossel method is discussed.

Experimental

Design and fabrication of the thin film waveguides

Planar X-ray waveguides can be used as a focusing tool to produce nanometer-sized X-ray beams [START_REF] Jarre | Two-dimensional hard x-ray beam compression by combined focusing and waveguide optics[END_REF][START_REF] Salditt | Compound focusing mirror and X-ray waveguide optics for coherent imaging and nano-diffraction[END_REF][START_REF] Hoffmann-Urlaub | Advances in fabrication of X-ray waveguides[END_REF]. The structure generally consists of a low-Z layer sandwiched between two high-Z layers [START_REF] Hayashi | Monochromatization of characteristic X-rays using stepped X-ray waveguide[END_REF]. In the waveguide, the X-rays propagate with total internal reflection at the interfaces between the high-and low-Z layers. So, it is promising to study the feasibility of PIXE-Kossel for the structural characterization of planar X-ray waveguides.

In this study, in the process for designing planar X-ray waveguide, we first performed electric field depth distribution simulations to optimize the structural parameters (including elements and thickness) using IMD software [START_REF] Windt | IMD-Software for modeling the optical properties of multilayer films[END_REF]. Interfacial and substrate roughness were assumed to be zero for the simulations. Two thin films waveguides were designed: Pt(4 nm)/Fe(6 nm)/Pt(10 nm)/Si and Ta(4 nm)/Cr(10 nm)/Pt(12 nm)/Si, noted PFP and TCP, respectively. As an example, Fig. 1 The designed thin films were deposited at room temperature on a Si (100) substrate by magnetron sputtering. The power applied on the sputtering targets was 10 W, the base pressure was 10 -8 mbar, and the sputtering gas was argon at the working pressure of 5×10 -2 mbar. The sputtering rates of Pt, Fe, Cr and Ta were 0.066, 0.024, 0.024 and 0.035 nm/s, respectively, as measured by a quartz microbalance.

Proton-induced X-ray emission

The schematic diagram of the PIXE-Kossel experiment is shown in Fig. 2. A 1.6 MeV proton beam generated by the Van de Graaff accelerator of the Système d'Analyse par Faisceaux d'Ions Rapides (SAFIR) platform, incident normal to the sample surface, was used to analyse the sample. An ANDOR iKon-M energy dispersive CCD camera equipped with a 1024×1024 sensor array and 13×13 µm pixels, was used to record the generated X-ray emissions. This CCD camera produces an X-ray spectrum for each of its pixels. Each column of pixels corresponds to a detection angle. To facilitate the data acquisition and subsequent data processing, a 4×4 binning is selected to obtain 256×256 pixel images (as shown in Fig. 2). The angular acceptance of the camera is 2.71°, as determined in Ref. [START_REF] Wu | Kossel diffraction observed with X-ray color camera during PIXE of nano-scale periodic multilayer[END_REF]. The measurement time is 2 h with 24 collected frames. Each frame acquisition time was 5 min.

As shown in Fig. 2, when the proton beam irradiates the sample, we first obtain a picture on the CCD camera, where the intensity of a pixel depends on the total number of X-rays collected in this pixel. After the energy calibration of the spectrum, a region of interest (ROI) is selected around a characteristic line, to obtain a 256×256 filtered spectral image. Then, this image is integrated in the vertical direction so as to obtain a 256 channel scan, representing the intensity of the given line as a function of the detection angle (grazing exit angle θ), i.e. the Kossel curve. Fig. 2. Scheme of the experimental setup for the PIXE-Kossel experiment. The angular acceptance of the camera is 2.71°.

X-ray reflectivity

To know the actual thickness and roughness of the deposited stacks, X-ray reflectivity (XRR) measurement is conducted. XRR is widely employed to study thin films and multilayers. In this study, the XRR measurements were performed at the 0.154 nm wavelength on a Rigaku five-circle diffractometer. Each measurement time was about 40 min. Then with the help of the IMD software [START_REF] Windt | IMD-Software for modeling the optical properties of multilayer films[END_REF] and the optical constants of Ref. [START_REF] Henke | CXRO database for x-ray attenuation length[END_REF], the measured XRR curves were fitted to determine the thickness, roughness and density of the various layers.

Results and Discussions

Analysis of Pt (4 nm)/Fe (6 nm)/Pt (10 nm) film waveguide

The XRR curve of the PFP sample is shown in Fig. 3. It displays the scattered X-ray intensity (in logarithmic scale) as a function of the glancing angle. For the fitting procedure, no interfacial compound is firstly considered. The values of thickness and roughness of each layer obtained from the fitting results are listed in Table 1. The densities, determined with an uncertainty of 10%, are close to those of the bulks. The above structural parameters (density, thickness and roughness of the layers) are then applied to the calculation of Kossel curves. The simulation is performed with the FLUORT code [START_REF] Chauvineau | Analyse angulaire de la fluorescence du fer dans une multicouche périodique Fe/C[END_REF]. Details of the Kossel curve simulation have been described in Ref. [START_REF] Wu | Kossel diffraction observed with X-ray color camera during PIXE of nano-scale periodic multilayer[END_REF]. The simulation method is originally designed for X-ray fluorescence calculation under Kossel effect. It calculates the electric field generated within the stack by the emitted radiation according to the reciprocity theorem. In the calculation, all parameters (thicknesses of the layers and roughness of each interface) determined from XRR fitting are applied. The simulation results will then be compared with the experimental results. We then carry out the PIXE experiment and obtain the overall X-ray emission spectrum (intensity in logarithmic scale), as shown in Fig. 4. The spectrum is the sum of the spectra obtained over all the pixels of the CCD camera. The energy calibration of the spectrum is performed from the position of the Fe Kα, Fe Kβ, Pt Lα and Pt Lβ emissions. Only Fe Kα and Pt Lα Kossel curves are discussed here. To extract the Fe Kα and Pt Lα Kossel curves, we select ROI on the overall spectrum (shown in Fig. 5).

The ROI range is chosen to minimize the influence of neighboring lines. The data processing is conducted using the PyMCA software [START_REF] Solé | A multiplatform code for the analysis of energy-dispersive X-ray fluorescence spectra[END_REF]. As an angular calibration is missing, we shift one experimental Kossel curve so that its main peak is at the same angle as the simulated one. The same shift is applied for all the Kossel curves. The simulated Kossel curves of the as-designed structure Pt(4 nm)/Fe(6 nm)/Pt(10 nm), and of that having the structure given in Table 1 are compared in Fig. 6. For Fe Kα, the curves (Fig. 6 appear at about 0.6° and 0.4° respectively, which are corresponding to the waveguide modes (in agreement with the results of Fig. 1). Small shifts of the main peaks are observed between the two simulations. We find that the main experimental features are almost reproduced by the simulated XRR curves. However, for the Fe Kα Kossel curve, the intensity toward higher angles has obvious difference, indicating that the description of the stack should be modified. The above results demonstrate that it is not enough to characterize the stack only according to the XRR fitting result.

Fig. 6. Comparison of the experimental (blue line) and simulated Fe Kα (a) and Pt Lα (b) Kossel curves of the PFP stack. Simulated curves are obtained with the parameters both from the as-designed structure and from the structure fitting the XRR measurements (Table 1).

Composition variation of the stack is then considered to improve the fit. In

Ref. [START_REF] Sato | Magnetooptical Spectra in Pt/Co and Pt/Fe Multilayers[END_REF], magneto-optical spectra of Pt/Fe multilayers were investigated. The 2).

Analysis of Ta(4 nm)/Cr(10 nm)/Pt (12 nm) waveguide

Fig. 9 shows the measured and fitted XRR curves of the TCP stack. It is difficult to obtain a reliable fitting result if no oxidation of the designed structure (Ta/Cr/Pt) is considered. Owing to the protective ability of Ta 2 O 5 compound against oxidization [START_REF] Díaz | Chromium and tantalum oxide nanocoatings prepared by filtered cathodic arc deposition for corrosion protection of carbon steel[END_REF]34], the Ta 2 O 5 /Cr/Pt structure is firstly considered. Again, fitting remains difficult if no oxidation of Cr is considered. It is reported that Cr 2 O 3 (the oxidation product of Cr) also possesses good corrosion resistance, which is usually used for the protection of components operating in harsh chemical environments [START_REF] Asteman | A comparison of the oxidation behaviours of Al 2 O 3 formers and Cr 2 O 3 formers at 700 °C-Oxide solid solutions acting as a template for nucleation[END_REF].

So the Ta 2 O 5 /Cr 2 O 3 /Cr/Pt structure is then taken into consideration. It is observed that although there exist some intensity discrepancies, the fit is well aligned with the measured XRR data. The obtained structural parameters of the stack are listed in Table 3. The fitted thickness of Ta 2 O 5 layer is almost two times higher than that of the original Ta layer. This could be explained by its lower density (bulk densities of Ta and Ta 2 O 5 are 16.65 [START_REF] Henke | CXRO database for x-ray attenuation length[END_REF] and 8.2 g/cm 3 [36] respectively). Compared with Cr (bulk density of 7.19 g/cm 3 ), the bulk density Cr 2 O 3 (5.21 g/cm 3 [START_REF] Henke | CXRO database for x-ray attenuation length[END_REF]) does not change too much. As a result, the sum of the fitted thicknesses of Cr 2 O 3 and Cr layers is almost equal to the designed Cr layer thickness. From the fitted results, we note that the oxidation of the stack is major (total for Ta and partial for Cr), so that the studied stack differs substantially from the designed structure. Based on the structure obtained from XRR (Table 3), Kossel curve simulations are performed (shown in Fig. 11), which are also compared with the as-designed structure and the experiment for the Cr Kα, Pt Lα and Ta Lα emissions. It is clearly seen that the as-designed structure cannot reproduce the experiment since large angular shifts and intensity discrepancies exist. Whereas the experimental Pt Lα curve (Fig. 11b) shows a reasonable agreement with the simulated one deduced from the XRR fit including oxidization effect, for the Cr Kα and Ta Lα curves (Fig. 11a andc) the relative intensities of the main features still exhibit obvious differences.

As a partial conclusion, since only introducing the oxidation of the stack does not allow fitting the experimental Kossel curves, we now study the possibility of compound formation. Due to the negative standard Gibbs energy of the reaction calculated with the HSC software [START_REF] Roine | HSC Chemistry® for Windows, Chemical Reaction and Equilibrium Software with Extensive Thermochemical Database[END_REF] (Ta (s) + 2Cr (s) = TaCr 2 (s), ΔG = -26.83 kJ/mol at 25°C [START_REF] Chase | NIST-JANAF thermochemical tables[END_REF][START_REF] Europe | Thermodynamic properties of inorganic materials, Landolt-Boernstein New Series, Group IV[END_REF]), we then assume that TaCr 2 is possibly formed between the Ta and Cr layers during the deposition process. Oxidation of the stack is then assumed when the sample is transferred in air. It has been reported that CrTaO 4 possesses good protective ability, which could decrease the oxidation rate of the underlying alloys [START_REF] Li | Oxidation behavior of a single-crystal Ni-base superalloy in air. I: at 800 and 900 °C[END_REF][START_REF] Ren | The influence of CrTaO 4 layer on the oxidation behavior of a directionally-solidified nickel-based superalloy at 850-900° C[END_REF]. So the new stack description, CrTaO 4 /TaCr 2 /Cr/Pt, is applied in the XRR fitting process. The corresponding reflectivity curve and its fit are displayed in Fig. 12 while the structural parameters of the stack extracted from the fit are collected in Table 4. The obtained structure differs greatly from the as-designed structure and also from the previously considered oxidized stack (Table 3). Compared with the as-designed structure, the stack thickness change (Table 4) can be explained by the density variation of the layers. Based on the conservation of mass, when the designed Ta layer is consumed totally for the TaCr 2 formation, the thickness of the TaCr 2 layer will be 13.1 nm and the residual thickness of the Cr layer will be 4.7 nm.

So the total thickness of TaCr 2 and Cr will be 17.8 nm (if no oxidation is considered).

When oxidation occurs, as the fitted density difference of CrTaO 4 and TaCr 2 is small, the sum of the fitted thickness of CrTaO 4 , TaCr 2 and Cr only changes a little, which is 18.73 nm. 13 shows the revised PIXE-Kossel simulations, which improve significantly.

Compared with the simulation curves in Fig. 11, especially for the Cr Kα, the relative intensities are much closer to the experimental ones. Hence, according to the XRR (Fig. 12) and the Kossel (Fig. 13) results, it is anticipated that the new designed structure is close to the real structure of the TCP stack. However, a refinement of this structure is still in need to reach an even better fit of the Kossel curves. Particularly the large discrepancy for the Ta curve (Fig. 13b) probably originates from a wrong number and depth distribution of the Ta atoms. The non-uniform Ta depth distribution makes it more difficult to improve the result. The uncertainty could be reduced if the optical constants used in simulation were better known. In addition, precise definition of the angular calibration would be also helpful to decrease the difference.

Conclusions

We have studied the feasibility of observing and exploiting proton-induced X-ray emission in Kossel geometry in stacks consisting of only a few layers, many less than the number in periodic multilayers designed for X-ray optics. Pt/Fe/Pt and Ta/Cr/Pt nano-scale tri-layers were designed to be used as X-ray planar waveguides. In the waveguide configuration, owing to the high sensitivity of the electric field versus the detection direction of the emitted radiation, PIXE-Kossel curves are sensitive to the depth distribution of the emitting atoms. In addition, the prepared stacks were also investigated by grazing incident X-ray reflectometry.

For all systems, the as-designed parameters (thickness, roughness and density of the layers) fail to generate both XRR and PIXE-Kossel simulated curves reproducing the experimental ones. It is always possible to find a structure able to reproduce the XRR curve, by considering some oxidation or the formation of some interfacial compound for example, but then the PIXE-Kossel curve simulations still differ markedly from the experimental ones. To get a satisfying agreement between experiment and simulation for both XRR and PIXE-Kossel curves, we had to consider that interlayers form at the interfaces of the Pt/Fe/Pt waveguide, while in the Ta/Cr/Pt system, some intermixing takes place between the Ta and Cr layers to form a compound, which is subsequently partly oxidized. This work demonstrates the ability of the PIXE-Kossel methodology to characterize the stacks of nano-scale thin films.

  (a) and (b) show the simulated depth and angular distribution of the electric field in the PFP stack at the photon energies of the Fe Kα and Pt Lα emissions, respectively. It is clearly seen that waveguide modes are formed in the Fe layer (guiding layer) at glancing angles close to 0.6° and 0.4° for the Fe Kα and Pt Lα emissions, respectively. This means that when the angle of detection is equal to these angles, an intensity maximum is expected for the corresponding emission. Any shift of the positions of the Fe or Pt atoms from the designed layers should give different position and amplitude of their features in the PIXE-Kossel curve.

Fig. 1 .

 1 Fig. 1. Simulated depth and angular distribution of the electric field in the planar PFP X-ray waveguide at the photon energy corresponding to (a) the Fe Kα emission (6403 eV) and (b) the Pt Lα emission (9442 eV). The vertical color bar gives the intensity scale of the electric field.

Fig. 3 .

 3 Fig. 3. Measured (black line) and fitted (red line) XRR curves of the PFP stack.

Fig. 4 .

 4 Fig. 4. X-ray emission spectrum of the PFP stack summed over all the pixels of the CCD camera.

Fig. 5 .

 5 Fig. 5. Fe Kα (a) and Pt Lα (b) ROI selected on the X-ray emission spectrum of the PFP stack.

  (a)) present one main peak while the simulated curves of Pt Lα (Fig. 6(b)) show two main peaks. Due to the current fluctuation of PIXE apparatus (no current regulation), the generated Kossel curves are in part noisy. It is noticed that the main peaks of Fe Kα and Pt Lα Kossel curves

  comparison of the simulation and the experiment confirmed the formation of Pt 69 Fe 31 alloy layer at the Pt-Fe interface. As a result, to improve the fitting result for the PFP stack, the presence of the interfacial compound Pt 69 Fe 31 is taken into account at both Pt-on-Fe and Fe-on-Pt interfaces. The reflectivity curve and corresponding fit are presented in Fig.7. After the consideration of Pt 69 Fe 31 interlayers, the χ² value decreases, indicating the fitting result improves. The geometrical information of the stack obtained from the modified fit is listed in Table2. It is also noticed that the thickness and roughness of Pt 69 Fe 31 at the Pt-on-Fe interface are different to those at the Fe-on-Pt interface, indicating that the interfaces are asymmetrical. Based on the obtained structure information (Table2), the Fe Kα and Pt Lα Kossel lines are simulated and then compared with the experimental ones. Fig.8shows the revised PIXE-Kossel simulations. The simulated XRR curves with introduction of Pt 69 Fe 31 interlayers show a better agreement with the experimental ones than the previous simulations (Fig.6). The difference between Fe Kα experimental and simulated Kossel curves at higher angles derives from the density and composition of Pt x Fe y intermixing layer. In this study, Pt 69 Fe 31 interlayers are considered at both interfaces, which could efficiently reduce the discrepancy between experiment and simulation curve at higher angle. It is difficult to define the exact composition and density of the formed Pt x Fe y intermixing layer. But by the combination of XRR and PIXE-Kossel, it is convenient to know whether intermixing occurs.

Fig. 7 .

 7 Fig. 7. Measured (black line) and fitted XRR curves of the PFP stack. Fitting result-1 (blue line) and Fitting result-2 (red line) represent the fitting results without and with the consideration of the presence of Pt 69 Fe 31 interlayers at both interfaces.

Fig. 8 .

 8 Fig. 8. Comparison of the experimental (blue line) and simulated Fe Kα (a) and Pt Lα (b) Kossel curves of the PFP stack. Simulated curves are obtained with the parameters both from the as-designed structure and from the structure fitting the XRR measurements (Table2).

Fig. 9 .

 9 Fig. 9. Measured (black line) and fitted (red line) XRR curves of the TCP stack.

Fig. 10 shows

 10 Fig.10shows the overall X-ray emission spectrum. The energy is calibrated from the positions of the Cr Kα, Ta Lα and Pt Lα emissions. Other peaks can be ascribed to the fluorescence of impurities inside the sample and/or of the experimental setup.Then Cr Kα, Ta Lα and Pt Lα Kossel curves are extracted similarly to those characteristic from the previous stack.

Fig. 10 .

 10 Fig. 10. X-ray emission spectrum of the TCP stack summed over all the pixels of the CCD camera.

Fig. 11 .

 11 Fig. 11. Comparison of the experimental and simulated (see Table 3) Cr Kα (a), Pt Lα (b) and Ta Lα (c) Kossel curves.

Fig. 12 .

 12 Fig. 12. Measured (black line) and fitted (red line) XRR curves of the TCP stack accounting for the compound formation and subsequent oxidation.

Fig. 13 .

 13 Fig. 13. Comparison of the experimental and revised simulated (see Table 4) Cr Kα (a), Pt Lα (b) and Ta Lα (c) Kossel curves.

Table 1 .

 1 Fitted values from XRR measurements of thickness, roughness and density of the different layers in the PFP stack (listed from the top layer to the substrate).

		Thickness/Roughness (nm)	Density
		(±0.05 nm)
	Pt	4.62/0.60	bulk
	Fe	6.65/0.80	bulk
	Pt	11.06/0.77	bulk
	Si substrate	bulk/0.42	bulk

Table 2 .

 2 Fitted values from XRR measurements of thickness, roughness and density of the different layers in the PFP stack (listed from the top layer to the substrate) accounting for the presence of Pt 69 Fe 31 interlayers at both interfaces.

		Thickness/Roughness (nm)	Density
		(±0.05 nm)	(g/cm 3 )
	Pt	3.83/0.57	bulk
	Pt 69 Fe 31	1.62/0.90	16 ± 2
	Fe	5.19/1.20	bulk
	Pt 69 Fe 31	1.41/0.79	16 ± 2
	Pt	10.25/0.39	bulk
	Si substrate	bulk/0.50	bulk

Table 3 .

 3

		Thickness/Roughness (nm)	Density
		(±0.05 nm)
	Ta 2 O 5	7.87/1.09	bulk
	Cr 2 O 3	6.55/1.30	bulk
	Cr	4.54/0.80	bulk
	Pt	13.79/0.54	bulk
	Si substrate	bulk/0.67	bulk

Fitted values from XRR measurements of thickness, roughness and density of the different layers in the TCP stack (listed from the top layer to the substrate).

Table 4 .

 4 Fitted values from XRR measurements of thickness, roughness and density of the different layers in the TCP stack (listed from the top layer to the substrate) accounting for the compound formation and subsequent oxidation.

		Thickness/Roughness (nm)	Density
		(±0.05 nm)	(g/cm 3 )
	CrTaO 4	10.75/1.10	6 ± 1
	TaCr 2	3.96/1.32	8 ± 1
	Cr	4.02/0.53	bulk
	Pt	13.79/0.54	bulk
	Si substrate	bulk/0.69	bulk
	Fig.

Acknowledgement

This work has been supported by the mobility scholarship within the framework of the Erasmus+ International Credit Mobility Program and the China Scholarship Council (CSC) scholarship within a CSC-Sorbonne University (SU) doctoral program. N. Casaretto and Y. Zheng from INSP are thanked for performing the XRR experiments.


Noël Lahellec
email: lahellec@lma.cnrs-mrs.fr

Pierre Suquet
email: suquet@lma.cnrs-mrs.fr

On the effective behavior of nonlinear inelastic composites: II. A second-order procedure

Keywords: Composites, homogenization, nonlinear viscoelasticity, elasto-viscoplasticity, variational method, second-order method, periodic micro-structure, random micro-structure, Hashin-Shtrikman estimates

A new method for determining the overall behavior of composite materials comprised of nonlinear viscoelastic and elasto-viscoplastic constituents is presented. Part I of this work showed that upon use of an implicit time-discretization scheme, the evolution equations describing the constitutive behavior of the phases can be reduced to the minimization of an incremental energy function. This minimization problem is rigorously equivalent to a nonlinear thermoelastic problem with a transformation strain which is a nonuniform field (not even uniform within the phases). In part I of this paper the nonlinearity was handled using a variational (or secant) technique. In this second part of the study, a proper modification of the secondorder procedure of Ponte Castañeda is proposed and leads to replacing, at each time-step, the actual nonlinear viscoelastic composite by a linear viscoelastic one. The linearized problem is even further simplified by using an "effective internal variable" in each individual phase. The resulting predictions are in good agreement with exact results and improve on the predictions of the secant model proposed in part I of this paper.

Introduction

This article is the second part of a study devoted to the overall response of nonlinear composites comprised of phases which have a partly reversible and partly irreversible mechanical behavior [START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I incremental variational principles[END_REF]. The individual constituents of the composites under consideration are generalized standard materials. As such, their mechanical state is characterized by their deformation ε at a given material point and a collection of internal variables α (Germain et al., 1983) and their constitutive relations derive from two thermodynamic potentials, a free-energy density w and a dissipation potential ϕ. The driving forces available in the system to modify its state are obtained by derivation of the free-energy density with respect to the state variables :

σ = ∂w ∂ε (ε, α), A = - ∂w ∂α (ε, α). ( 1
)
The evolution of the internal variables is obtained by derivation of the dual of the dissipation potential ϕ * : α = ∂ϕ * ∂A (A).

(2)

A representative volume element V of a nonlinear composite comprised of N different homogeneous phases, with free-energy w (r) and dissipation potential ϕ (r) is considered and subjected to a loading history. In part I of this study [START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I incremental variational principles[END_REF], an incremental variational principle has been derived following the lines of [START_REF] Mialon | Eléments d'analyse et de résolution numérique des relations de l'élasto-plasticité[END_REF], [START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF] and [START_REF] Miehe | Strain-driven homogenization of inelastic micro-structures and composites based on an incremental variational formulation[END_REF] (see part I of this study for more details). After time-discretization of the equations governing the quasi-static evolution of this representative volume element, it can be shown that the local fields of deformation and internal variables in the composite have a variational property which is used to reduce the problem of composites with two potentials to a problem for a nonlinear composite with a single potential. More specifically the overall stress Σ at time t = t n+1 can be deduced from the strain E at time t = t n+1 by derivation of a single condensed effective potential :

Σ = ∂ w ∆ ∂E (E) , (3)
where w∆ has the variational characterization :

w∆ (E) = Inf ε =E w ∆ (ε) , w ∆ (ε) = Inf α J(ε, α), (4)
and

J(ε, α, x) = N r=1 w (r) (ε, α) + ∆t ϕ (r) α -α n (x) ∆t χ (r) (x). ( 5
)
w ∆ is the local condensed potential, whereas J is the incremental potential.

For reasons exposed in [START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I incremental variational principles[END_REF], it is more convenient to work with J rather than with w ∆ . A first difference with the case of nonlinear composites with a single potential lies in the fact that the incremental potential J not only depends on the strain ε but also depends on the internal variables α. However, given the variational properties expressed by the relations (4), this additional dependence can be overcome by considering more variables in the minimization. Second, as noted in the first part of this study [START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I incremental variational principles[END_REF], two salient features of the incremental potential J complicate the analysis and require suitable approximations :

a) The incremental potential J, just like the dissipation potential ϕ (r) , is a nonquadratic function of α, even when the elasticity of the constituents is linear (w (r) is quadratic). Therefore J has to be approximated by a potential which is quadratic with respect to ε and α. This "linearization" procedure, which is usual in nonlinear composites with a single potential, amounts in the present setting to replacing the actual nonlinear viscoelastic (or elastoviscoplastic) composite by a linear viscoelastic one. b) A second difficulty stemming from the dependence of J on the field α n (x) requires a further approximation to approach J by a piecewise uniform potential (as a function of ε and α). This is done, as in part I of this study, by introducing an effective internal variable (EIV) α (r) n which replaces, as accurately as possible, the field α n (x) in phase r.

In the first part of this study, these two steps were taken care of simultaneously by a procedure derived from the variational procedure. But the two approximations are distinct in general and, in this part of the study, they will be handled in separate steps. The decoupling between the two approximations is illustrated by the case of linear viscoelasticity where no linearization is required but where the "effective internal variable" approximation is nevertheless needed (Lahellec and Suquet, 2006a).

Regarding the linearization strategy (point a), the procedure developed in part I of this study was based on the so-called variational method (in a form which has no real bounding property and which was therefore called a secant method). However, having in mind the improvements brought by second-order methods [START_REF] Ponte Castañeda | Exact second-order estimates for the effective mechanical properties of nonlinear composite materials[END_REF], 2002) over the variational method when the constitutive relations of the individual constituents are governed by a single potential, the use of one of these second-order estimates (SOE) is a natural alternative to the approach followed in the first part of the study. A first attempt in this direction is made in section 2, where a procedure based on the original second-order procedure is explored. However the resulting predictions are found to exhibit unrealistic features, such as a strong dependence of the predicted response on the time-step. After elucidating the origin of the problem, which lies in the gap between the overall stress field in the linear thermoelastic comparison composite and the derivative of the energy, use is made in section 3 of a modification due to [START_REF] Lahellec | Nonlinear composites: a linearization procedure, exact to second-order in contrast and for which the strain-energy and the affine formulations coincide[END_REF] of the original SOE. The resulting predictions are compared in section 5 with fullfield solutions and show no abnormal dependence on the time-step. They also improve significantly on the secant estimates derived in the first part of this study.

2 A procedure based on the original second-order method 2.1 The original second-order procedure of Ponte Castañeda

The idea of the original second-order method (SOE for second-order estimates) of Ponte [START_REF] Ponte Castañeda | Exact second-order estimates for the effective mechanical properties of nonlinear composite materials[END_REF] for nonlinear composites with a single potential is to replace the potential by its second-order Taylor expansion in the neighborhood of the average strain in each individual phase. In doing so, the initial nonlinear composite is replaced by a linear thermoelastic composite with anisotropic elastic moduli. The resulting prediction has no bounding property and yields only an estimate of the effective potential. On the other hand, this estimate is shown [START_REF] Ponte Castañeda | Exact second-order estimates for the effective mechanical properties of nonlinear composite materials[END_REF] to be exact to second-order in the contrast for nonlinear composites comprised of phases having a small contrast between their material properties (exact results for this case being available from the work of [START_REF] Suquet | Small-contrast perturbation expansions for the effective properties of nonlinear composites[END_REF].

The SOE can be transposed in the present context of composites with two potentials (one of them, w (r) being quadratic), by replacing the nonquadratic potential ϕ (r) by its second-order Taylor expansion :

ϕ (r) α -α n ∆t ϕ (r) soe α -α n ∆t = ϕ (r) α(r) + ∂ϕ (r) ∂ α α(r) : ∆ α + 1 2 ∆ α : ∂ 2 ϕ (r) ∂ α2 α(r) : ∆ α (6)
where α(r) is an effective rate for the internal variable in phase r and

∆ α = α -α n ∆t -α(r) .
As initially proposed in Ponte Castañeda (1996), α(r) is taken to be the phaseaverage of the rate of the internal variables :

α(r) = α -α n ∆t r = α r . (7)
Making the approximation (6) in the variational problem (4), amounts to replacing the incremental potential J by :

J soe (ε, α) = N r=1 w (r) (ε, α) + ∆tϕ (r) soe α -α n ∆t χ (r) (x) = N r=1 1 2 (ε -α) : L (r) : (ε -α)+ ∆t ϕ (r) ( α r ) + ∂ϕ (r) ∂ α ( α r ) : ∆ α + 1 2 ∆ α : L (r) v : ∆ α χ (r) (x), (8)
where χ (r) denotes the characteristic function of the domain V r occupied by the r-th phase and where :

L (r) v = ∂ 2 ϕ (r) ∂ α2 ( α r ) . (9)
Under this approximation, the variational problem (4) becomes :

w ∆ (E) ≈ w ∆soe (E) = Inf ε =E Inf α J soe (ε, α) . ( 10
)
The associated Euler-Lagrange equations read :

       σ = L (r) : (ε -α) = ∂ϕ (r) ∂ α ( α r ) + L (r) v : ∆ α in V r , div σ = 0 in V, ε = E, + boundary conditions on ∂V. (11)
The first line in (11) can be developed into

σ = L (r) : (ε -α) = L (r) v : α -α n ∆t + τ (r) , τ (r) = ∂ϕ (r) ∂ α ( α r ) -L (r) v : α(r) .
This last relation shows that the approximation (10) amounts to replacing the initial nonlinear composite by a linear thermo-viscoelastic composite with an anisotropic viscosity tensor L (r) v , an initial (at time t n ) internal variable field α n (x) and a piecewise uniform eigenstress τ (r) . Note that α n depends on x so that the thermo-viscoelastic problem (11) has to be solved numerically in general.

Once the equations ( 11) are solved for the local fields σ, ε and α the macroscopic stress-strain relations is obtained by derivation of the overall potential

Σ soe = d w ∆soe dE (E). ( 12
)
In practice, the potential w ∆soe is not known explicitly as a function of E even when analytical schemes (such as the Hashin-Shtrikman bounds) are used to solve the thermoelastic problem (11). Therefore the derivative involved in (12) has to be computed numerically by taking difference quotients in E. However, unexpected and inconsistent results, such as a significant dependence of the predictions on the time-step ∆t were obtained, as shown by the following example.

An enlighting example

The above comments are illustrated in figures 1 and 2. We consider a periodic two-dimensional fiber-reinforced composite whose microstructure is characterized by the unit-cell shown in figure 6a. It is composed of an isotropic powerlaw nonlinear viscoelastic matrix reinforced by elastic fibers. The free-energy density in both phases is given by :

w (r) (ε, α) = 1 2 (ε -α) : L (r) : (ε -α), L (r) = 3k (r) J + 2µ (r) K, k (r) = E (r) 3(1 -2ν (r) ) , µ (r) = E (r) 2(1 + ν (r) ) , (13)
where J and K are the two orthogonal projectors on spherical and deviatoric tensors respectively.

The dissipation potential is identically 0 in the fibers and reads in the matrix as :

ϕ( α) = σ 0 ε0 m + 1 αeq ε0 m+1 , (14)
The material characteristics of the phases are :

• Fibers : E = 400 GPa, ν = 0.2 • Matrix : E = 70 GPa, ν = 0.3, σ 0 = 480 MPa, ε0 = 10 -2 s -1 and m = 0.2.

The fiber volume-fraction is c (1) = 0.3. The loading is an extension in the first direction :

E(t) = E 11 (t)e 1 ⊗ e 1
, where E 11 is increased monotonically with constant strain-rate E 11 = 5 10 -3 t.

The thermo(visco)elastic problem ( 11) is solved exactly by means of the Finite Element Method (FEM). Different time-steps have been used. The "exact" result is shown in figure 1 as a solid line, whereas the different predictions of the second-order procedure for different time-steps are shown by dashed lines. As can be seen, the second-order procedure shows a significant sensitivity of the predictions to the time-step. In particular, the asymptotic overall stress (as t tends to +∞) as predicted by the second-order procedure for different time-steps varies between to limits denoted as Σ ∞ aff and Σ ∞ soe (more details in section 2.4). It is close to Σ ∞ soe for large time-steps (∆t = 2 s) and close to Σ ∞ aff for small time-steps (∆t = 2 10 -2 s). Another way of evaluating the dependence of the predictions on the time-step is to consider the asymptotic predictions (at large times) of the stress obtained for different time-steps. This is done in figure 2 where the strong dependence of the asymptotic stress on the time-step can be seen.

Second-order and affine procedures

The original second-order method of Ponte Castañeda (1996), developed for composites governed by one nonquadratic potential, proceeds by approximating the initial potential by its Taylor expansion up to order 2. This expansion is a polynomial of degree 2 and the Euler-Lagrange equations associated with its minimization are those for a linear, anisotropic, thermoelastic composite.

The effective energy of this thermoelastic composite is taken as an approximation for the actual effective potential and the overall stress is obtained by taking the derivative of this potential with respect to the overall strain as indicated by (12).

The affine method proposed by [START_REF] Masson | Self-consistent estimates for the rate-dependent elastoplastic behaviour of polycrystalline materials[END_REF] (see also [START_REF] Masson | An affine formulation for the prediction of the effective properties of nonlinear composites and polycrystals[END_REF][START_REF] Brenner | Improved affine estimates for nonlinear viscoelastic composites[END_REF][START_REF] Brenner | Improved affine estimates for nonlinear viscoelastic composites[END_REF], for subsequent developments) considers the Taylor expansion of the stress-strain relations in the phases up to order one (thus the method is called the affine method) and leads again to a linear thermoelastic problem, identical to that encountered in the secondorder procedure (expanding the potential to order 2 is equivalent to expanding the constitutive relations to order one). However there is a difference in the way in which the overall stress is deduced from the resolution of the thermoelastic problem. In the affine method, the overall stress is approximated by the average of the stress field in the thermoelastic composite and not by the derivative of the effective energy. The two results for the overall stress differ, although the affine and the second-order methods consider the same thermoelastic composite.

More specifically the affine prediction for the overall stress reads as :

Σ aff = σ , (15)
which is different from (12) (Ponte Castañeda and [START_REF] Ponte Castañeda | Nonlinear composites[END_REF]. When applied to (11) the affine method showed none of the inconsistencies of the second-order method.

To understand the origin of the different sensitivity of the SOE and the affine methods, the difference between the two predictions can be explicitly evaluated. The overall stress predicted by the affine procedure reads as :

Σ aff = N r=1 c (r) L (r) : ( ε r -α r ) = N r=1 c (r) ∂ϕ (r) ∂ α ( α r ) , (16)
whereas the stress Σ soe is obtained by taking the derivative with respect to E of the homogenized incremental potential

w ∆soe (E) = N r=1 c (r) 1 2 (ε -α) : L (r) : (ε -α) +∆t ϕ (r) ( α r ) + ∂ϕ (r) ∂ α ( α r ) : ∆ α + 1 2 ∆ α : L (r) v : ∆ α r . (17)
Its derivative reads as :

d w ∆soe (E) dE = N r=1 c (r) (ε -α) : L (r) : dε dE + ∆t 1 2 ∆ α : dL (r) v d α : d α r dE : ∆ α r , (18)
and finally

Σ soe = Σ aff + ∆t N r=1 c (r) 1 2 ∆ α : dL (r) v d α : d α r dE : ∆ α r . ( 19
)
2.4 Why there should be no gap between the energy formulation and the affine procedure

Consider the case where the loading history consists of an imposed average strain history with a constant (in time) strain-rate Ė. As already noticed in part I of this study [START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I incremental variational principles[END_REF], the overall response of the composite is governed initially (for small times t 1) by its effective linear elastic properties (in other words by the effective potential w alone) whereas its response at large times (t → +∞) is governed by its viscous properties (in other words by the effective potential φ alone). In the second limit the local fields σ and ε tend to stationary values σ ∞ and ε∞ as t tends to +∞, the stress-rate σ vanishing. It follows from the elastic relations that the difference εα tends to 0 and therefore that α can be asymptotically replaced by ε∞ in the constitutive relations. The local problem solved by the local fields therefore reads asymptotically

σ ∞ = ∂ϕ (r) ∂ ε ( ε∞ ) in V, div σ ∞ = 0 in V, ε∞ = Ė + boundary conditions.
(20) The standard SOE can now be applied to this problem with a single potential. The potential ϕ r) in each phase is consequently replaced by its second-order Taylor expansion. This leads to approximating the fields σ ∞ and ε∞ by the solutions of a thermoelastic problem "tangent" to the nonlinear problem (20) (for simplicity, the corresponding local fields will still be denoted by σ ∞ and ε∞ ). Let Σ ∞ aff denote the overall stress corresponding to the affine estimate for the problem (20) :

Σ ∞ aff = N r=1 c (r) ∂ϕ (r) ∂ ε ( ε∞ r ) . (21)
The affine estimate for the transient problem (finite t) obtained by the incremental procedure is given by ( 16) and is to be compared to (21). In the limit as t tends to +∞, one has ε = α and the estimates ( 16) and ( 21) coincide.

The second-order estimate of the same problem (20), denoted by Σ ∞ soe , is given by (Ponte Castañeda and [START_REF] Ponte Castañeda | Nonlinear composites[END_REF] :

Σ ∞ soe = Σ ∞ aff + N r=1 c (r) 1 2 ( ε∞ -ε∞ r ) : dL (r) v d ε : d ε∞ r d Ė : ( ε∞ -ε∞ r ) r , (22)
and is to be compared with ( 19). The two estimates Σ aff and Σ ∞ aff coincide asymptotically and it is sufficient to compare the last terms in ( 19) and ( 22). The following asymptotic results are worthy of notice :

lim t→+∞ ∆ α = ε∞ -ε∞ r and lim t→+∞ dL (r) v d α ( α r ) = dL (r) v d ε ( ε∞ r ).
It remains to compare the terms ∆t(d α r /dE) in Σ soe and d ε∞ r /d Ė in Σ ∞ soe . To this end, the relation (11a) can be averaged in phase r :

L (r) : ( ε r -α r ) = ∂ϕ (r) ∂ α ( α r ) = ∂ϕ (r) ∂ α α -α n ∆t r ,
and this relation can be differentiated with respect to E to give :

d α r dE = L (r) + 1 ∆t L (r) v -1 : L (r) : d ε r dE ,
which implies, after time derivation :

∆t d α r dE = L (r) + 1 ∆t L (r) v -1 : L (r) : d ε r d Ė . ( 23
)
As a consequence :

lim t→+∞ ∆t d α r dE = L (r) + 1 ∆t L (r) v -1 : L (r) : d ε∞ r d Ė . ( 24
)
This last relation sheds light on the problems encountered with the original second-order method:

-The second-order estimate (19) does not converge, as t → +∞, towards ( 22), as it should. -The asymptotic limit of (19) depends on the elastic moduli of the phases (the L (r) 's), on their dissipative properties (the L (r) v 's) and on the time-step ∆t. This last point is illustrated in figure 2.

-It is only in the extreme case where L (r) 1 ∆t L (r) v that the correct limit is retrieved.

-In the other extreme case L (r) 1 ∆t L (r) v , the second-order estimate (19) converges to the affine estimate (21).

Remark 1 : Since the original second-order procedure of Ponte Castañeda (1996) which is followed here, other propositions have been made which lead to more accurate results in the context of composites governed by a single potential (see for instance Idiart et al., 2006a). In all these improved methods, the overall stress is computed as the derivative of the effective energy with respect to the overall strain and not as the average of the local stress field in the linear comparison composite (doing so would be in line with the above affine method). The gap between the two possible ways of evaluating the overall stress (derivative of the energy and average of the stress field in the linear comparison composite) is significantly reduced by these improved methods. However it is not clear that the dependence on the time-step would be removed by adopting one of these methods. Our guess is that the problem will persist as long as a gap between the two approaches, derivative of the energy and average of the stress field in the comparison composite, remains, but this statement deserves a deeper investigation which is left for future work.

Orientation for the next section

To avoid the above inconsistency in the prediction of the second-order method ( lim t→+∞ Σ soe = Σ ∞ soe ) it would be quite natural to make use of the affine method rather than of the second-order original approach. On the other hand the affine approach is known to be too stiff and less accurate than the secondorder procedure (for instance it is not exact to second-order in the contrast of the phases [START_REF] Masson | An affine formulation for the prediction of the effective properties of nonlinear composites and polycrystals[END_REF]. These limitations led us to develop a method having some of the features of the second-order method (based on energy considerations and exact to second-order in the contrast), combined with some of the features of the affine method (no dependence on the timestep). The direction that we have pursued is to enforce the stationarity of the effective potential w∆ with respect to all (adjustable) variables defining the approach. This stationarity prescription, which was not present in the original second-order procedure of Ponte [START_REF] Ponte Castañeda | Exact second-order estimates for the effective mechanical properties of nonlinear composite materials[END_REF], ensures that the affine and energy-based approaches coincide

σ = d w ∆ dE (E),
where σ is the stress field in the linear thermoelastic comparison composite. When this equality is satisfied, we can expect that there is no dependence on the time-step in such an approach (and the numerical results show that this guess is correct).

3 A procedure based on a new second-order method

Reduction to a linearly viscoelastic problem

Ponte Castañeda and [START_REF] Ponte Castañeda | Nonlinear composites[END_REF] have shown that, in the case of a single potential, the affine and the energy-based predictions differ by a term involving the third derivative of the potential. The relation (19) shows that in the present context they differ by a term involving the first derivative of L (r) v with respect to α or, equivalently, the third derivative of the dissipation potential ϕ (r) . Therefore, in order to close the gap between the affine and the second-order (or energy-based) method, the initial potential ϕ (r) is expanded to third-order and additional adjustable parameters are introduced in the third-order term in such a way that the resulting approximation remains a polynomial of degree 2. More specifically the initial potential ϕ (r) is replaced by the expression

ϕ (r) α -α n ∆t ϕ (r) toe α -α n ∆t ,
where

ϕ (r) toe α -α n ∆t = ϕ (r) α(r) + ∂ϕ (r) ∂ α α(r) : α -α n ∆t -α(r) + 1 2 ∂ 2 ϕ (r) ∂ α2 α(r) :: α -α n ∆t -α(r) ⊗ α -α n ∆t -α(r) + 1 6 ∂ 3 ϕ (r) ∂ α3 α(r) ::: α -α n ∆t -α(r) ⊗ α (r) -α(r) ⊗ α (r) -α(r) , (25)
where α(r) (to be specified) is the uniform reference value for the rate of the internal variables in phase r and where the last term, which looks like the third-order term in the Taylor expansion of the potential, except that it has been transformed into a first-order term by adding an additional adjustable variable α (r) which remains to be specified.

Substitution of the expansion (25) into the variational problem (4) leads to

w ∆ (E) = Inf ε =E Inf α J toe (ε, α) , ( 26
)
where J toe is defined as:

J toe (ε, α) = N r=1 w (r) (ε, α) + ∆t ϕ (r) toe α -α n ∆t χ (r) (x).
The variables α(r) and α (r) are chosen to achieve stationarity in (26). In the proposed procedure, the final form of the variational problem (4) is :

w ∆ (E) = Stat α(r) Stat α (r) Inf ε =E Inf α J toe (ε, α) . ( 27
)
The stationarity condition with respect to α (r) reads as:

∂ 3 ϕ (r) ∂ α3 α(r) : α -α n ∆t -α(r) ⊗ α (r) -α(r) r = 0. ( 28
)
Several choices for α(r) are possible in order to satisfy (28). The simplest one is to set α(r) equal to the average of the time derivative of the internal variables in each individual phase :

α(r) = α -α n ∆t r = α r . ( 29
)
When this choice is made, the stationarity condition with respect to α(r) reduces to :

1 3 N (r) v ( α r ) :: α (r) -α r ⊗ α (r) -α r = N (r) v ( α r ) :: C (r) α , (30)
with :

N (r) v ( α r ) = ∂ 3 ϕ (r) ∂ α3 ( α r ) , C (r) α = ( α -α r ) ⊗ ( α -α r ) r . (31)
Again several choices for α (r) are possible in order to satisfy (30) depending on the specific form of the potentials ϕ (r) . In the next section devoted to powerlaw materials, a specific choice is detailed for specific loading conditions.

Once the two stationarity conditions ( 28) and ( 30) are satisfied, the final expression for J toe (ε, α) is:

J toe (ε, α) = N r=1 1 2 (ε -α) : L (r) : (ε -α) + J (r) ref +τ (r) : (α -α n ) + 1 2 (α -α n ) : 1 ∆t L (r) v : (α -α n ) χ (r) (x), (32)
where

L (r) v = ∂ 2 ϕ (r) ∂ α2 α(r) , J (r) ref = ∆t ϕ (r) α(r) -τ (r) : α(r) -1 2 α(r) : L (r) v : α(r) , τ (r) = ∂ϕ (r) ∂ α α(r) -L (r) v : α(r) + 1 6 N (r) v ( α r ) :: α (r) -α r ⊗ α (r) -α r .                       
The variational problem ( 27) is then reduced to :

w ∆ (E) Inf ε =E Inf α J toe (ε, α) , (33)
where J toe is given by (32). The associated local problem reads as :

σ = L (r) : (ε -α) = τ (r) + L (r) v : α -α n ∆t in V r , div σ = 0 in V, ε = E, periodicity conditions on ∂V.        (34)
One can recognize in (34) a problem for a generalized linear thermo-viscoelastic composite where the r-th individual phase is characterized by its elastic moduli L (r) , its viscous moduli L (r) v , an eigenstress τ (r) and a field of internal variables α n at time t n .

In conclusion, the linearization rule (25), in association with the stationarity conditions over α(r) and α (r) , has allowed us to reduce the minimization of the incremental potential to the minimization of a potential for an anisotropic linear viscoelastic problem. However, as in the case of one potential, the problem remains nonlinear because of the closing conditions ( 29) and ( 30). The coefficients of the linear viscoelastic composite depend on the solutions of these nonlinear equations which themselves depend on the local strain field in the linear thermo(visco)elastic composite.

Remark 2 : It is essential to note that, in this approach, the derivative of the effective energy w ∆ with respect to the overall strain E coincides with the average of the stress in the linear thermo-viscoelastic composite. Indeed one has

d w ∆ dE (E) = ∂J toe ∂ε : ∂ε ∂E + ∂J toe ∂α : ∂α ∂E + ∂J toe ∂ α(r) : ∂ α(r) ∂E + ∂J toe ∂ α(r) : ∂ α(r) ∂E .
However, thanks to the stationarity conditions on α, α(r) and α (r) in ( 27), the three last terms in these relations vanish. The first term can be transformed by Hill's lemma and the resulting relation reads as :

d w ∆ dE (E) = ∂J toe ∂ε : ∂ε ∂E = ∂J toe ∂ε : I = σ , (35)
where σ = ∂J toe ∂ε (ε, α) is the stress field in the thermo-viscoelastic composite.

Resolution of the linear viscoelasticity problem

The variational problem ( 26) with (32) involves the non piecewise uniform field of internal variables α n (x) arising from the previous time-step. This difficulty can be overcome by introducing an effective internal variable (EIV) α (r) n at time t n , following the same procedure as in part I of this study [START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I incremental variational principles[END_REF].

The potential J toe (ε, α) is compared to a reference potential J 0 (ε, α) given by :

J 0 (ε, α) = N r=1 J (r) 0 (ε, α)χ (r) (x),
where 36) Note the two differences between J 0 and J toe given by (32). First, the field α n (x) of internal variables at time t n has been replaced in each individual constituent by an internal variable α (r) n which is now uniform in phase r. Second, the viscous moduli L (r) v have been multiplied by a scalar factor θ (r) which is closely linked to the replacement of the field α n (x) by α (r) n as shown in Lahellec and Suquet (2006a) or [START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I incremental variational principles[END_REF] to which the reader is referred for more details. It should be noted that if θ (r) is taken equal to 1 (no modification of the tangent moduli), then, by virtue of (40), this amounts to considering only fields α n (x) which have no fluctuations in the phases (piecewise uniform fields). As shown in the appendix A of the first part of this study [START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I incremental variational principles[END_REF], the resulting approach then reduces to the Transformation Field Analysis which is known to be too stiff.

J (r) 0 (ε, α) = w (r) (ε, α)+J (r) ref +τ (r) : (α -α n )+ 1 2 (α-α (r) n ) : θ (r) ∆t L (r) v : (α-α (r) n ). (
With these two modifications in mind, J toe can be written as J 0 + (J toe -J 0 ), which leads to replacing w ∆ in (33) by :

w ∆ (E) = Inf ε =E Inf α J 0 (ε, α) + Stat α [J toe (ε, α) -J 0 (ε, α)] . (37)
The last stationarity problem, which is local in character, reads at each point x in phase r :

Stat α 1 2 (α -α n ) : 1 ∆t L (r) v : (α -α n ) - 1 2 (α -α (r) n ) : θ (r) ∆t L (r) v : (α -α (r) n ),
whose solution is :

α(x) = N r=1 α n (x) -θ (r) α (r) n 1 -θ (r) χ (r) (x). ( 38
)
With this relation, the estimate (37) becomes :

Inf ε =E Inf α J 0 (ε, α) + N r=1 c (r) 2 (α n -α (r) n ) : θ (r) ∆t(θ (r) -1) L (r) v : (α n -α (r) n ) r .
(39) This estimate can again be optimized with respect to θ (r) and α (r) n which have not yet been specified. Stationarity of (39) with respect to these variables yields

∂ ∂θ (r) J (r) 0 (ε, α) + 1 2 (α n -α (r) n ) : θ (r) ∆t(θ (r) -1) L (r) v : (α n -α (r) n ) r = 0 ⇒ θ (r) = 1 ± (α n -α (r) n ) : L (r) v : (α n -α (r) n ) r (α -α (r) n ) : L (r) v : (α -α (r) n ) r , (40)
and

∂ ∂α (r) n J (r) 0 (ε, α) + 1 2 (α n -α (r) n ) : θ (r) ∆t(θ (r) -1) L (r) v : (α n -α (r) n ) r = 0 ⇒ α (r) n = α n r + (θ (r) -1) α r θ (r) , (41)
The expression (39) can be further simplified by eliminating α. The field α solution of the infimum problem in ( 39) is given in phase r by :

α(x) = L (r) + θ (r) ∆t L (r) v -1 : L (r) : ε(x) + L (r) + θ (r) ∆t L (r) v -1 : θ (r) ∆t L (r) v : α (r) n -τ (r) . ( 42
)
Upon substitution of this expression in ( 39), the condensed effective potential can be estimated as :

w 0 (E) + N r=1 c (r) 2 (α n -α (r) n ) : θ (r) ∆t(θ (r) -1) L (r) v : (α n -α (r) n ) r . (43)
In this expression w0 is the effective energy of an auxiliary thermoelastic problem :

w0 (E) = Inf ε =E w 0 (ε) , (44)
with

w 0 (ε) = N r=1 w (r) 0 (ε) χ (r) (x), w (r) 0 (ε) = 1 2 ε : L (r) 0 : ε + ρ (r) 0 : ε + f (r) 0 ,
where the tensors L (r) 0 and ρ (r) 0 and the energy f (r) 0 are piecewise uniform and defined as :

L (r) 0 = L (r) : L (r) + θ (r) ∆t L (r) v -1 : θ (r) ∆t L (r) v ρ (r) 0 = L (r) : L (r) + θ (r) ∆t L (r) v -1 : τ (r) - θ (r) ∆t L (r) v : α (r) n f (r) 0 = J (r) ref -τ (r) : α n r + 1 2 α (r) n : θ (r) ∆t L (r) v : α (r) n - 1 2 τ (r) + θ (r) ∆t L (r) v : α (r) n : L (r) + θ (r) ∆t L (r) v -1 : τ (r) + θ (r) ∆t L (r) v : α (r) n                               
(45) The first and second-order moments of the strain field (in the thermoelastic body) are given by (see for instance appendix C of part I of this study) :

ε r = A (r) 0 : E + a (r) 0 and ε ⊗ ε r = 2 c (r) ∂ w0 ∂L (r) 0 (46)
where

A (r) 0 , a (r)
0 are the classical localization tensors in a thermoelastic composite with the same microstructure as the actual nonlinear composite. They can be obtained either numerically or analytically for specific microstructures.

Remark 3 : Using the same arguments (stationarity conditions) as in remark 2, one can prove that

d w ∆ dE (E) = ∂J 0 ∂ε (ε, α) = σ , ( 47
)
where σ is now the stress field in the thermoelastic composite defined by the relations (45).

Remark 4 : As in the first part of this study, two different expressions ( 38) and ( 42) for the field of internal variables α have been derived. Let α * and α * * denote respectively these two different expressions. The first one solves the stationarity condition in the evaluation of the gap between J and J 0 whereas the second realizes the infimum condition which serves to define the thermoelastic potential w 0 from the potential J 0 .

A straightforward calculation starting from (37) shows that the stationarity conditions ( 40) and ( 41) express equality between the first moment of α * and α * * in each individual phase and equality of a scalar quantity related to the second moment of these fields :

α r = α * r and L (r) v :: α ⊗ α r = L (r) v :: α * ⊗ α * r (48)

Effective response of the composite

Let us first note that four unknowns are a priori to be determined, three of them being tensorial α(r) , α (r) n , α (r) and one unknown being scalar, θ (r) . However, the relations ( 29) and (41) show that these unknowns are related in such a way that one of them can be eliminated. It was chosen to work with α r , α (r) and θ (r) as the main unknowns.

The procedure for determining the effective response of the composite then goes as follows:

1. At time t n , the first moment α n r and the second moment α n : α n r of the internal variables are known in each individual phase r.

2. The nonlinear equations consisting of the average of (42) over each phase r, (30) and ( 40), are solved for α r , α (r) and θ (r) by a modification of the hybrid method of [START_REF] Powell | A hybrid method for nonlinear algebraic equation[END_REF]. In the present study use has been made of the MINPACK library (HYBRD1 routine).

At each step of this iterative procedure the residues of these equations are obtained by the following procedure :

2.1. The first and second moment of the strain are deduced from the relations (46). These moments depend on the thermoelastic constants (45).

2.2. The first and second moment per phase of the unknown field of internal variable α at time t n+1 are obtained by the relation (42) from the first and second moments of the strain field ε computed at step 1.

2.3. The residues are computed using the first moment of the field of internal variables α computed at step 2.2 and the relations (30) and (40).

3. Once convergence is reached, the macroscopic stress Σ n+1 is obtained by taking the average of the microscopic stress field :

Σ n+1 = σ n+1 = N r=1 c (r) L (r) 0 : ε r + ρ (r) 0 . ( 49
)
Details about the convergence-rate of this algorithm are given in section 5.2.2.

The implementation of the TOE method is, in its principle, very similar to that of the secant method proposed in the first part of this study. The main differences lie first in the elastic moduli of the linear thermoelastic composite which are now anisotropic and second in the nonlinear equation for α(r) . Even the writing of this nonlinear equation requires the explicit knowledge of the third derivative of the dissipation potential which is not required in the secant method. Therefore, despite its similarities with the secant method, the TOE approach is more complex.

3.4

The TOE procedure applied to composites with a single potential Although the TOE procedure was not initially meant for nonlinear composites governed by a single potential, it has also been applied in this case (Lahellec andSuquet, 2004, Rekik et al., 2005). For completeness, the main results of these authors are summarized here.

(1) The TOE procedure is exact to second-order in the contrast, as is the original SOE. Incidentally it provides an affine method which is exact to second-order in the contrast. (2) In most cases which have been explored, the effective properties predicted by the TOE are close to the predictions of the SOE and exhibit similar trends as functions of the nonlinearity exponent or volume-fraction of the phases. Therefore it is likely that the limitations of the SOE, such as bounds violation in extreme cases, also apply to the TOE. (3) In particular the field fluctuations in the phases do not influence the elastic moduli of the linear comparison composite, but affect the eigenstress. (4) The SOE is known to exhibit a "duality gap" between the primal and dual formulations, in the sense that the predictions obtained by taking a Taylor expansion of the strain potential or of the stress potential are different. This duality gap is also present in the TOE. (5) By construction, the TOE has no gap between the energy-based formulation and the affine (or stress-averaging) formulation. Therefore, in the TOE, the stress and strain fields in the linear comparison composite (thermoelastic) can be consistently taken as approximations of the actual stress and strain fields in the nonlinear composite. This is not the case in the SOE approach.

Power-law materials

In order to implement the TOE procedure, and in particular in order to solve the equation (30) for α(r) , one has to specify further the form of the potentials. Attention will be restricted here to isotropic nonlinear viscoelastic constituents with a power-law dissipation potential. The above theory applies also in principle to elasto-viscoplastic constituents (with a yield surface delimiting an elasticity domain), but we have not checked its practical implementation (we recall that it has been checked, in an unpublished study, that the secant method of part I works for elasto-viscoplastic materials as well).

Detailed form of the operators

The potentials are taken in the form :

w (r) (ε, α) = 1 2 (ε -α) : L (r) : (ε -α), L (r) = 3k (r) J + 2µ (r) K, ϕ (r) ( α) = σ (r) 0 ε0 m + 1 αeq ε0 m+1 .        (50)
The first two derivatives of ϕ (r) read as :

∂ϕ (r) ∂ α α(r) = 2µ (r) v α(r) , L (r) v = ∂ 2 ϕ (r) ∂ α2 α(r) = +∞J+2λ (r) v E (r) +2µ (r) v F (r) = +∞J+2µ (r) v mE (r) + F (r) , (51)
where +∞ refers to the fact that ϕ (r) is finite only for incompressible α, and :

µ (r) v = σ (r) 0 3 ε0 α(r) eq ε0 m-1 , λ (r) v = mµ (r) v , E (r) = 2 3 α(r) ⊗ α(r) ( α(r) eq ) 2 , F (r) = K-E (r) .
(52) Finally, the third derivative of ϕ (r) reads as :

N (r) v = 4 9 (m-1) σ (r) 0 ε2 0 α(r) eq ε0 m-2    3 K ⊗ α(r) 3s α(r) eq + 2 3 (m -3) α(r) ⊗ α(r) ⊗ α(r) ( α(r) eq ) 3    , (53)
with the convention that T 3s ijkhmn = 1 3 (T ijkhmn + T khmnij + T mnijkh ).

4.2 Determination of α (r) .

The equations ( 30) are solved for each phase separately and the unknowns α (r) can be decomposed into two terms respectively parallel and orthogonal to α(r) :

α (r) = β α(r) + α⊥(r) , with α⊥(r) : α(r) = 0. ( 54
)
Upon use of this decomposition and of the expression (53) of N , the equation ( 30) reduces to :

α⊥(r) (β -1) α(r) eq + 1 3 α(r) α(r) eq 3 2 (β -1) 2 ( α(r) eq ) 2 m + α⊥(r) : α⊥(r) = 2C (r) α : α(r) α(r) eq + α(r)   K :: C (r) α α(r) eq + 2 3 (m -3) C (r) α :: α(r) ⊗ α(r) ( α(r) eq ) 3   . ( 55
)
The left-hand side of this equation has a component parallel to α⊥(r) and another component parallel to α(r) . The right-hand side has a term parallel to α(r) and another term parallel to C (r) α : α(r) .

In all cases that we have investigated, due to symmetries in the microstructure and in the loading, the term C (r) α : α(r) has been observed numerically to be always parallel to α(r) . This was in particular the case for the two-dimensional periodic microstructures with two orthogonal planes of symmetry, subjected to a loading having the same axes of symmetry, such as the example studied in section 5.2.

Therefore, we will not discuss in general the resolution of equation ( 55) but, instead, work under the assumption that

C (r) α : α(r) is parallel to α(r) . ( 56
)
Under this assumption, the solution α (r) of equation ( 55 (57)

The above equation does not allow for a complete determination of α⊥(r) which is, in fact, not needed. The only expression involving α(r) is the eigenstress τ (r) given by ( 33). In the case of power-law materials, and thanks to the result β = 1, this expression reduces to :

τ (r) = 2(m -1)µ (r) v -1 + α⊥(r) : α⊥(r) 9( α(r) eq ) 2 α(r) , (58)
where only ( 57) is required.

The EIV approximation for power-law materials

The EIV approximation introduces two unknowns α (r) n and θ (r) given by equations ( 41) and ( 40) respectively. The first relation ( 41) does not depend on the material under consideration (power-law or else) whereas the second relation (40) becomes for power-law materials, using the expression (51) for L (r) v :

θ (r) = 1 - mE (r) + F (r) :: (α n -α (r) n ) ⊗ (α n -α (r) n ) r mE (r) + F (r) :: (α -α (r) n ) ⊗ (α -α (r) n ) r (59)
The sign in front of the square root, which was left undetermined in ( 40), has been taken to be a minus, for it was observed by Lahellec and Suquet (2006a) that the plus sign leads to predictions which are excessively stiff.

Note that the relations ( 57) and ( 59) require the evaluation of the contracted product (over four indices) of (α -

α (r) n ) ⊗ (α -α (r) n ) r
with the fourth-order tensor mE (r) + F (r) which is itself proportional to L (r) v according to (51). It is clearly the case for (59). Regarding (57), we note that

C (r) α = 1 ∆t 2 (α -α n ) ⊗ (α -α n ) r -α r ⊗ α r , = θ (r) ∆t 2 2 α -α (r) n ⊗ α -α (r) n r -α r ⊗ α r .
Therefore the product of C (r) α with L (r) v , which is needed in (57), can be directly deduced from the product of (αα

(r) n ) ⊗ (α -α (r) n ) r
with the fourth-order tensors L (r) v and from the product of α r ⊗ α r with the same tensor. This last product is straightforward, given the relation (51).

In conclusion, the equations to be solved involve only the first moment α r in each phase of the field of internal variables α and the projection of the second-order moment of this field on L (r) v .

The thermoelastic problem

The projection of interest for α is obtained from the same projection of the strain field ε in the linear thermoelastic comparison composite with effective energy given by ( 44) with, for power-law materials :

L (r) 0 = 3k (r) J + 2λ (r) 0 E (r) + 2µ (r) 0 F (r) , (60)
with : r) and µ (r) refer to the linear elastic bulk and shear moduli of the phases. The eigenstress reads as :

λ (r) 0 = µ (r) θ (r) λ (r) v /∆t µ (r) + θ (r) λ (r) v /∆t and µ (r) 0 = µ (r) θ (r) µ (r) v /∆t µ (r) + θ (r) µ (r) v /∆t . k (
ρ (r) 0 = µ (r) µ (r) + θ (r) λ (r) v /∆t   2µ (r) v (1 -m) α(r)   1 - α⊥(r) : α⊥(r) 9( α(r) eq ) 2   - θ (r) λ (r) v α (r) n ∆t . (61)
The detailed expression of f (r) 0

is not needed.

The first moment of ε is given by the relation (46 a) and its projections over r) and F (r) are obtained by derivation of the effective energy with respect to λ (r) 0 and µ (r)

E (
0 :

E (r) :: ε ⊗ ε r = 1 c (r) ∂ w 0 ∂λ (r) 0 and F (r) :: ε ⊗ ε r = 1 c (r) ∂ w 0 ∂µ (r) 0 . ( 62
)
5 Examples

Orientation

In this section, the accuracy of the TOE model derived in section 3 is assessed by comparing its predictions with the exact response (obtained numerically) of composite with specific microstructures. Comparisons with the secant procedure derived in the first part of this study also show the improvements brought by the present second-order procedure (TOE).

As in the first part of this study, it should be noted that there are two levels of approximation involved in the practical implementation of the present procedure. The first approximation (which is in fact two-fold) consists in the linearization procedure (TOE) and in the introduction of an effective internal variable (EIV). The second level of approximation pertains to the resolution of the thermoelastic problem. Our objective is to check the discrepancy introduced by the first set of approximations only. In order to minimize the error due to the resolution of the thermoelastic problem we have used either an exact (numerical) solution of the thermoelastic problem (section 5.2) or a class of microstructures for which an (almost) exact analytical solution of the thermoelastic problem is available (section 5.3).

The same five results are shown for each microstructure or material data.

(1) The first set of results, identified as "Exact (FEM)" or "Exact (FFT)", is obtained by solving numerically the nonlinear nonlinear viscoelastic problem, either by the Finite Element Method (FEM) for the simple periodic microstructure of section 5.2, or by a method based on Fast Fourier Transforms for the cylinder assemblage in section 5.3. (2) The second set of results, identified as "TOE (FEM)", is obtained by approximating the nonlinear incremental potential (4) by the quadratic potential J

(r)

toe given by the relation (32). The corresponding thermoelasticity problem (34) involves nonuniform eigenstrain fields and is solved by the FEM.

(3) The third set of results, identified as "TOE (EIV+HS)", is obtained by approximating, in the linear incremental potential, the nonuniform fields of internal variables by piecewise uniform "effective internal variables", leading to the thermoelastic problem (44). In addition, the effective energy of this thermoelastic problem is evaluated either by one of the classical Hashin-Shtrikman bounds (in the case of the composite cylinder assemblage) or by a generalized Hashin-Shtrikman bound due to [START_REF] Suquet | Une méthode simplifiée pour le calcul des propriétés élastiques de matériaux hétérogènes à structure périodique[END_REF] (see also Lahellec and Suquet, 2006a) in the case of the periodic microstructure of section 5.2. In all cases, the Hashin-Shtrikman bounds are obtained with the matrix as the reference medium. (4) The fourth set of results, identified as "SEC (FEM)" refers to the intermediate secant model (secant linearization without the further approximation of effective internal variables). More details about this model were given in the first part of this study. (5) The fifth set of results, identified as "SEC (EIV+HS)", corresponds to the complete secant model of part I of this paper. The Hashin-Shtrikman bounds are obtained following the same lines as in the third set of results.

The focus will be on two-phase fiber composites. The fibers, aligned with the direction of e 3 , will be identified as phase 1, whereas the continuous phase, called the matrix, will be identified as phase 2. Two microstructures will be considered, a simple periodic microstructure consisting of a circular fiber in a square unit-cell repeated periodically (figure 6a) and the assemblage of composite cylinders already used in part I [START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I incremental variational principles[END_REF]) of this study. The composites are subjected to an imposed macroscopic shear strain :

E(t) = E 11 (t) (e 1 ⊗ e 1 -e 2 ⊗ e 2 ) . (63)
The strain E 11 (t) is increased monotonically with a constant strain-rate.

Square array of circular fibers

Material data

We consider here an elastically-reinforced matrix with a fiber volume-fraction ranging from c (1) = 0.05 to c (1) = 0.5. The fibers are isotropic and linearly elastic whereas the matrix is nonlinear viscoelastic with potentials in the form (50). The material data are :

-Fibers : E = 400 GPa, ν = 0.2.

-Matrix : E = 70 GPa, ν = 0.3, σ 0 = 480 MPa, ε0 = 10 -2 s -1 .

In addition, the matrix rate-sensitivity exponent m varies from 0.1 to 1. The loading is given by ( 63) with a constant strain-rate Ė11 = 5 10 -3 s -1 .

Implementation of the TOE approach

The TOE approach has been implemented along the lines of section 3.3. In the present case of purely elastic fibers (phase 1), the unknowns of the problem are α 2 , α (2) and θ (2) . In all examples investigated in the present study, the iterative algorithm of section 3.3 has always converged to a solution. We believe that this solution is unique (at least in rate-dependent problems) but we have no proof of this assertion. Convergence is reached when the relative errors between two iterates for θ (2 ), α (2) and < α > 2 is less than a given threshold (typically 10 -8 ). With this criterion and for a typical example of the present section corresponding to m = 0.2 and c (1) = 0.3, convergence was attained in 20 to 40 iterations (typically), the slowest convergence being observed in the transient regime where elastic and viscous effects are of the same order of magnitude. For the same problem, the secant method of part I took between 15 and 30 iterations to converge.

Results and discussion

As already noted in the first part of this study, the main trends of the stressstrain curve of these composites under monotone loading at constant strainrate are typically an initial linear-elastic response for small strains, then a transient part which finally leads to a plateau corresponding to the purely viscous response of the composite. This plateau is characterized by a stress Σ ∞ 11 . These trends are shown in figure 3 where the five sets of results described in section 5.1 are plotted for different fiber volume-fraction or rate-sensitivity exponents. Figures 3 (a),(b) and (c) correspond to the same exponent m = 0.2 and to different volume fractions c (1) = 0.2, 0.3 and 0.4. Figures 3 (d),(e) and (f) correspond to the same fiber volume-fraction c (1) = 0.4 with different rate-sensitivity exponents m = 0.1, 0.5 et 1.

The following comments can be made :

(1) The error introduced by the EIV approximation and by the HS bound is measured for both models, SEC and TOE, by the distance between SEC(FEM) and SEC(EIV+HS) (same for TOE). As can be seen this error is small and depends on the size of the transient part of the response. (2) Except in the linear case where no linearization is needed, the secant model overestimates the stiffness of the composite. Although it has not been shown to be a rigorous upper bound, it tends to lie systematically above the exact results. The discrepancy increases with the fiber volumefraction and with the nonlinearity (the smaller m, the larger the nonlinearity).

(3) By contrast with the SEC model, the TOE model underestimates the composite stiffness. The discrepancy, which is in all cases less than in the SEC model, increases with the fiber volume-fraction. However, there is a specific fiber volume-fraction around c (1) = 0.4 where it seems that the nonlinearity does not affect the error. (4) In all cases, the largest error is observed in the asymptotic regime. Therefore the error introduced by the model can be measured by comparing the asymptotic stresses Σ ∞ 11 . This is done in Figure 4 as a function of the fiber volume-fraction for m = 1 and m = 0.2.

(5) Comparison between the strain fields in the unit-cell can be performed from the full-field computations, either in the actual nonlinear composite or in the thermoelastic comparison composites introduced in the SEC and TOE procedures. The corresponding snapshots of the ε 11 component of the strain are shown in figure 6, corresponding to the exact result (figure 6b), the prediction of SEC(FEM) (figure 6c) and the prediction of TOE(FEM) (figure 6d). White or light grey correspond to high deformations, whereas black corresponds to low deformation. The field estimated by the TOE model is in good agreement with the actual strain field. The better agreement of the strain fields compensates a higher complexity of the TOE model as compared to the SEC model.

Composite cylinder assemblage

The microstructures considered in the present section have been described in the first part of this study and in [START_REF] Moulinec | Intraphase strain heterogeneity in nonlinear composites: a computational approach[END_REF] and Idiart et al. (2006b). They approach, as closely as possible within computational limitations, an ideal composite cylinder assemblage for which one of the Hashin-Shtrikman bounds is a sharp estimate for isotropic phases. Two fiber volumefractions have been considered : c (1) = 0.21 and c (1) = 0.41. Both phases are nonlinear viscoelastic with potentials in the form (50) with the same ratesensitivity exponent m. The material data were :

-Fibers : E = 100 GPa, ν = 0.45, ε0 = 1 s -1 .

-Matrix : E = 100 GPa, ν = 0.45, σ 0 = 1 GPa, ε0 = 1 s -1 . Two different cases for the flow-stress in the fibers were considered, σ

(1) 0 = 5 and σ

(1) 0 = 0.2 GPa, corresponding to fiber-reinforced or fiber-weakened composites. The rate-sensitivity exponent m, identical for both phases, was varied from 0.1 to 1. The applied loading was a pure shear (63) with a constant strain-rate Ė11 = 5 10 -1 s -1 . Three predictions for the overall stress-strain response of the composites are presented in figure 7 in the fiber-reinforced case and in figure 8 for the fiberweakened case. The exact result corresponds to the full-field computation performed by FFT. The TOE(EIV+HS) prediction corresponds to the TOE linearization together with the EIV approximation. In addition the thermoelas- tic is solved using the relevant Hashin-Shtrikman bound for two-dimensional composites (lower bound in the fiber-reinforced case, upper bound in the fiberweakened case).

When m = 1 (linear viscoelasticity), the secant and second-order predictions coincide, since there is no linearization needed. In the fiber-reinforced case, the predictions are softer than the exact result because of the use of the lower HS bound (conversely in the fiber-weakened case, the model is stiffer than the exact results). These bounds are accurate rigorous bounds for the microstructure under consideration, but not exact results.

When m < 1, the fiber-reinforced case confirms the trends of the preceeding section: for high enough nonlinearity, the SEC predictions are stiffer than the exact results whereas the TOE predictions tend to lie on the compliant side. In the fiber-weakened case, the SEC prediction remains too stiff whereas the TOE prediction lies now on the stiff side, but with a better agreement than the SEC model. In all cases the error is maximal in the asymptotic regime and the different models can be compared by inspecting the asymptotic stresses Σ ∞ 11 . This is done in figure 9 where the three predictions for the asymptotic stress are compared when the rate-sensitivity exponent m is varied. The overall stress has been normalized by the similar asymptotic stress in the matrix alone under the same loading conditions. The TOE model improves on SEC model in all cases. The improvement is especially significant for the fiber-weakened composite where the trends of the SEC model are opposite to the trends of the exact results (as m varies) whereas the TOE model captures well the influence of m.

Closing remarks

This study is devoted to nonlinear composite materials comprised of nonlinear viscoelastic constituents whose constitutive relations can be derived from two thermodynamic potentials. The present study is based on an incremental variational principle derived in the first part of this paper [START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I incremental variational principles[END_REF]). In the present paper the application of a new second-order procedure, where a modification of the initial second-order of Ponte Castañeda (1996) has been proposed to close the gap between the so-called affine and energy approach, has been investigated.

Our main findings can be summarized as follows :

(1) The prediction obtained by a procedure following strictly the original second-order method showed inconsistent features due to the difference between the stress obtained by derivation of the energy and the average stress in the linear thermoelastic composite (gap between the affine and strain-energy based approaches).

(2) By adding to the second order expansion of the dissipation potential a term depending on the third derivative of the potential to the eigenstress, it has been possible to close this gap and to recover consistent results.

The resulting scheme remains exact to second-order in the contrast of the phases. . Comparison between the exact result "Exact (FEM)", the linearization scheme solved with the EIV approximation associated with the Hashin-Shtrikman method "TOE (EIV+HS)" and "SEC (EIV+HS)".(a) : σ

32

(3) In addition to this more consistent linearization, an "effective internal variable approach" has been proposed in order to reduce the problem to that for a N-phase thermoelastic composite. Examples shown in section section 5 show that this approximation is indeed very accurate. (4) Comparison with the secant model (SEC) proposed in the first part of this study [START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I incremental variational principles[END_REF] shows on one hand that the predictions of this second-order procedure (TOE) are more accurate than those obtained with the SEC model. On the other hand, they are obtained at the expense of a more complicated implementation, making use of an anisotropic thermoelastic composite and also relying on assumptions which hold true in the two-dimensional examples that have been investigated here but which remain to be checked in more general (3-dimensional) analyses.

There remains room for significant improvement of the proposed schemes.

(1) First the modification of the SOE proposed by [START_REF] Lahellec | Nonlinear composites: a linearization procedure, exact to second-order in contrast and for which the strain-energy and the affine formulations coincide[END_REF] is not optimal. Although it achieves its initial goal, which was to close the gap between the overall stress evaluated either from the energy formulation or from the average stress in the linear thermoelastic composite, its predictions remain close to that of the original SOE. As such, its limitations are similar: there remains a significant gap between the primal and the dual formulations and its predictions can be inaccurate in extreme cases. A first improvement could be brought by using a more refined second-order procedure such as those recently developed by Idiart et al. (2006a).

(2) Second the procedure used to derive the "effective internal variable" is directly inspired from the incompressible isotropic case of the first part of this study. It should be improved by considering directly the EIV in the thermoelastic potential, instead of treating this step separately.

Fig. 3 .

 3 Fig. 3. Square array of circular fibers. Elastic fibers in a nonlinear viscoelastic matrix. Response of the composite under in-plane shear for different fiber volume-fractions and for different rate-sensitivity exponents.(a): c (1) = 0.2 and m = 0.2. (b): c (1) = 0.3 and m = 0.2.(c): c (1) = 0.4 and m = 0.2. (d): c (1) = 0.4 and m = 0.1. (e): c (1) = 0.4 and m = 0.5. (f ): c (1) = 0.4 and m = 1.

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig. 4. Square array of circular fibers. Elastic fibers in a nonlinear viscoelastic matrix. In-plane macroscopic shear. Asymptotic macroscopic stress Σ ∞ 11 as a function of the fiber volume-fraction c (1) for two different rate-sensitivity exponents. a) m = 0.2, b) m = 1. The asymptotic stress is normalized by the stress in the unreinforced matrix under the same loading conditions, Σ (2),∞ 11.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Composite cylinder assemblage. Fiber-reinforced case σ (1) 0 /σ (2) 0 = 5. (a): c (1) = 0.41 and m = 1.(b): c (1) = 0.41 and m = 0.2.

Fig. 9 .

 9 Fig. 9. Composite cylinder assemblage. Asymptotic macroscopic stress Σ ∞ 11 as a function of the nonlinearity exponent m. The asymptotic stress is normalized by the stress in the unreinforced matrix under the same loading conditions, Σ (2),∞ 11

Acknowledgements

This study is part of program supported by the Centre National de la Recherche Scientifique under grant CNRS/NSF n • 14555.