H Choukri

A Fischer

S Forget

S Chenais

M-C Castex

B Geffroy

D Adès

A Siove

Doped and non-doped organic light-emitting diodes based on a yellow carbazole emitter into a blue-emitting matrix

Keywords: carbazole, OLED, doped and non-doped devices, white emission ___________________________________________________________________________

A new carbazole derivative with a 3,3'-bicarbazyl core 6,6'-substituted by dicyanovinylene groups (6,6'-bis(1-(2,2'-dicyano)vinyl)-N,N'dioctyl-3,3'-bicarbazyl; named (OcCz2CN) 2 , was synthesized by carbonyl-methylene Knovenagel condensation, characterized and used as a component of multilayer organic light-emitting diodes (OLEDs). Due to its π-donor-acceptor type structure, (OcCz2CN) 2 was found to emit a yellow light at λ max =590 nm (with the CIE coordinates x=0.51; y = 0.47) and was used either as a dopant or as an ultra-thin layer in a blue-emitting matrix of 4,4'-bis(2,2'-diphenylvinyl)-1,1'-biphenyl (DPVBi). DPVBi (OcCz2CN) 2 -doped structure exhibited, at doping ratio of 1.5 weight %, a yellowish-green light with the CIE coordinates (x = 0.31; y = 0.51), an electroluminescence efficiency η EL =1.3 cd/A, an external quantum efficiency η ext = 0.4 % and a luminance L= 127 cd/m 2 (at 10 mA/cm 2 ) whereas for non-doped devices utilizing the carbazolic fluorophore as a thin neat layer, a warm white with CIE coordinates (x = 0.40; y= 0.43), η EL = 2.0 cd/A, η ext = 0.7 %, L = 197 cd/m 2 (at 10 mA/cm 2 ) and a color rendering index (CRI) of 74, were obtained. Electroluminescence performances of both the doped and non-doped devices were compared with those obtained with 5,6,11,12-tetraphenylnaphtacene (rubrene) taken as a reference of highly efficient yellow emitter.

.

Introduction

Organic light emitting diodes (OLEDs) have attracted great attention for full-colour flat-panel displays since the demonstration of efficient electroluminescent devices [1]. Moreover, the tremendous progress in device luminous efficiency [2] and the emerging viability as a commercial display open new markets for white OLEDs (WOLEDs), like backlight for liquid crystal displays or lighting sources. WOLEDs can be considered as potential large area next generation of solid state lighting sources to replace traditional incandescent white light sources, thanks to their potential of energy saving, their high efficiencies and their possibility to fabricate thin and flexible devices [3]. To achieve white emission, various methods have been used such as e.g. multilayer diodes using the combination of the three primary colours Red, Green and Blue (RGB), multiple-quantum-well architecture [4], or exciplex formation [5]. Among these various devices, numerous doped-type WOLEDs using two emission colours to produce white have been fabricated [6,7] but it remains technologically difficult to attain an accurate control of the dopant concentration at low doping ratio through the co-evaporation process. We recently utilized the delta-doping technique [8], a way to colour-control, including balanced white emission, in a multilayer non-doped OLED [9][10][11] based on blue matrices, in which we inserted an ultrathin layer of rubrene as a yellow emitter [12]. Hence, we intended to test another yellow material with a different structure and coming from the carbazolic molecular derivatives.

__________________ * Corresponding author

Indeed, these materials have been widely used in optoelectronics [13] and particularly as emitters in OLEDs by several of us [14][15][16][17]. The main goal of our work is turned towards synthetic chemistry of luminescent organic materials belonging to the carbazole family. Unsusbstituted carbazolic members are well-known for their good hole transporting properties and for their blue-light emitting capabilities (see e.g. review ref. [14]).

With such a goal, we report, herein, on the synthesis of a novel carbazole derivative with an electron-donating 3,3'-bicarbazyl core linked, at the 6,6'-positions, to electron-accepting dicyanovinylene groups i.e. 6,6'-bis(1-(2,2'-dicyano)vinyl)-N,N'-dioctyl-3,3'-bicarbazyl (named (OcCz2CN) 2 ). This compound, was synthesized by carbonyl-methylene Knovenagel condensation between bicarbazyl dialdehyde and dicyanomethane and was characterized by 1 H and 13 C NMR, electronic absorption and emission spectroscopies and cyclic voltammetry. Molecular structure of (OcCz2CN) 2 , is shown in Figure 1. Due to its π-donor-acceptor type structure, (OcCz2CN) 2 exhibits an energy band-gap of Eg = 2.4 eV and was found to emit a yellow light with Commission Internationale de l'Eclairage (CIE) coordinates of (x=0.51; y = 0.47). To the best of our knowledge, it is the first carbazolic derivative exhibiting such properties.

It was used in multilayer OLEDs either, as a dopant of a blue-emitting matrix of 4,4'-bis(2,2'diphenylvinyl)-1,1'-biphenyl (DPVBi) or as an ultrathin layer inserted in the latter. For a doping ratio of 1.5 weight %, devices exhibited modest electroluminescence (EL) performances e.g. an electroluminescence efficiency η EL =1.3 cd/A, an external quantum efficiency η ext = 0.4 % and a luminance L= 127 cd/m 2 (at 10 mA/cm 2 ), a yellowishgreen light being emitted. For the non-doped diodes utilizing a thin layer of (OcCz2CN) 2 within the exciton diffusion zone in DPVBi, diodes exhibited better EL performances and an emission colour varying from yellow to white depending on the position of the carbazole film within the DPVBi layer.

For such delta-doped devices utilizing a carbazolic thin neat layer, white emission with CIE coordinates (x = 0.40; y= 0.43), η EL = 2.0 cd/A, η ext = 0.7 %, L = 197 cd/m 2 , were obtained. To the best of our knowledge, it is the first report of a white OLED using a carbazolic material as a thin emitting layer. EL performances of both devices were compared with those obtained with rubrene which is by far the best organic yellow fluorophore and thereby which is taken as a yellow quantum efficiency reference in the literature and is known to enhance multilayer device efficiency and stability [18,19].

Experimental

Materials

Aluminium tris(8-hydroxyquinolinate) (AlQ 3 ), copper pthalocyanine (CuPc), 4,4'-bis[N-(1-naphtyl)-N-phenylamino]biphenyl, (NPB), 4,4'-bis(2,2'diphenyl vinyl)-1,1'-biphenyl (DPVBi), 5,6,11,12tetraphenylnaphtacene (rubrene) were purchased from H.W. Sands Company. Molecular structures are shown in Figure 1.

The synthetic route of (OcCz2CN) 2 is shown in Figure 2. All the solvents were purified by the usual methods prior to utilization. (OcCz2CN) 2 was synthesized by carbonyl-methylene Knovenagel condensation between N,N'-dioctyl-3,3'-bicarbazyl-6,6'-dicarbaldehyde and dicyanomethane as the active methylene compound. N,N'-dioctyl-3,3'-bicarbazyl-6,6'-dicarbaldehyde was prepared according to the method we previously described [20]. We briefly recall, that the dialdehyde was obtained by diformylation of N,N'-dioctyl-3,3'-bicarbazyl (itself prepared by oxidative dimerization of N-octylcarbazole by iron trichloride in chloroform) by N,Ndimethylformamide (DMF) in o-dichlorobenzene (ODCB) (first reaction step of Figure 2). After column chromatography with chloroform/hexane (3:1) and recrystallization from diethylether, yellow crystals of the dialdehyde were obtained in 55% yield.

Knovenagel condensation is carried out as follows. Equimolar quantities of the dialdehyde and dicyanomethane (2 mmol), were dissolved in 20 ml of dry tetrahydrofuran (THF) and introduced into a threenecked flask fitted with a condenser and a nitrogen inlet. A 2 mmol solution of tetrabutylammonium hydroxide (Fluka) in an isopropanol-methanol mixture (10:1), was then added to the reaction medium under nitrogen. After refluxing for 1 hour, reaction was stopped by cooling the reaction mixture at 0°C, and poured into ice-water (150 ml). The resulting brownish solid obtained was filtered off and purified by column chromatography with chloroform/hexane (3:1). Orange powder was obtained in 72 % yield.

Physico-chemical characterization

1 H and 13 C NMR spectra were recorded in solution in CDCl 3 with a Brucker AC200 spectrometer at 200 and 50 MHz respectively; all chemical shifts (δ in ppm) were referenced to tetramethylsilane (TMS). UV-visible spectroscopy was carried out using a Perkin Elmer spectrometer (model Lambda 12) with 0.5 cm path quartz cells. For solid-state measurements, (OcCz2CN) 2 chloroform solutions were cast onto quartz plates. Optical band gap E g was estimated from the absorption threshold of the spectra. Photoluminescence spectra were obtained with a Safas model flx-Xenius spectrofluorometer with a Xenon 150 W lamp. Cyclic voltammetry (cv) experiments of solution-cast films of (OcCz2CN) 2 onto conducting support were carried out in acetonitrile (HPLC grade, Fisher Scientific) containing 0.3 M of LiClO 4 (Aldrich) as a supporting electrolyte in a onecompartment cell equipped with three electrodes. The working anode was glassy carbon, the counterelectrode a platinum wire, the reference electrode being a saturated calomel electrode (SCE) at 4.75 eV.

OLEDs fabrication and characterization

Figure 3 details the layer structures of the two types of devices with the detail of the organic compounds and the thickness of each layer of the two types of devices a)-type (DPVBi layer doped with yellow emitter) and b)-type (non-doped type, with a thin layer of yellow emitter in DPVBi). Indium tin oxyde (ITO)-coated glass with sheet resistance of about 15 Ω/ was obtained from Asahi. Prior using, the glass was cleaned by sonication in a detergent solution and then rinsed in de-ionized water. The organic compounds are deposited onto the ITO anode by sublimation under high vacuum (10 -7 Torr) with a rate of 0.2 -0.3 nm/s. An in-situ quartz crystal was used to monitor the thickness of the vacuum depositions. The active area of each OLED was 0.3 cm 2 . In the doped configuration, we sequentially deposited onto the ITO anode: a thin (10 nm) layer of copper phtalocyanine (CuPc) as a hole-injecting layer, then a 50 nm thick NPB layer as a hole-transporting layer (HTL). A 60 nm thick DPVBi layer doped with the yellow emitter (OcCz2CN) 2 is then deposited and acts as a blue matrix and an electron-transporting layer (ETL). The doping ratio is controlled by coevaporation of host and dopant and was optimized to be 1.5% (in weight percent) of (OcCz2CN) 2 , the accuracy of the doping rate being around 10%. Finally, a thin (10 nm) AlQ 3 layer was used as an electron-injecting layer onto the top cathode consisting in 1.2 nm of LiF capped with 100nm of Al. In the nondoped configuration, the same sequential process was followed, keeping constant the thicknesses of the organic layers, the only difference from the doped configuration being that an ultrathin layer of (OcCz2CN) 2 was inserted within the DPVBI layer, at a position d from the NPB/DPVBi interface in the range 0 ≤ d ≤10 nm. Such a range corresponds to the exciton diffusion domain within the DPVBi side [12]. Electroluminescence spectra and chromaticity coordinates were recorded with a PR 650 SpectraScan spectrophotometer at a constant current density of 10mA/cm 2 . The Current-Voltage and Current-Luminance (I-V and I-L) characteristics of the diodes were measured with a regulated power supply (ACT 10 Fontaine) combined with a multimeter and a 1 cm 2 area silicon calibrated photodiode (Hamamatsu). All the measurements were performed at room temperature and under ambient atmosphere, without any encapsulation. An important quality of white light emitters for illumination purpose is the Color Rendering Index (CRI). The CRI defines how well colors are rendered under different illumination conditions in comparaison to a standard (i.e. a thermal radiator or daylight). As defined by the Commission Internationale de l'Eclairage [START_REF]Method of measuring and specifying color rendering properties of light sources[END_REF], a set of eight calibrated samples with different absorption spectra and a reference illuminant are used to define the CRI. The reference source used for the CRI calculation is not given in the CIE definition. To choose the appropriate illuminant, we first calculate the Correlated Color Temperature (CCT) (according to reference [START_REF] Hernandez-Andres | Calulating correlated color temperatures across the entire gamut of daylight and skylight chromaticities[END_REF]) of our White OLED. The reference source used for the CRI calculation is then the black body at this temperature if the CCT is less than 5500K, and the D65 Daylight illuminant otherwise.

Results and discussion

Synthesis, characterization and physicochemical properties

As illustrated in Figure 2 the synthesis of (OcCz2CN) 2 proceeds in two steps. At first, bicarbazyldialdehyde was prepared by formylation of bicarbazyl according to the method that we previously described [20]. In a second step, (OcCz2CN) 2 was obtained by Knovenagel condensation of the bicarbazyldialdehyde as the carbonyl reactive compound, with dicyanomethane as the methylene active derivative. Both the 1 H and 13 C NMR spectra and elemental analysis are in agreement with the chemical structure of (OcCz2CN) 2 shown in Figure 1.

Absorption and photoluminescence spectra of 50nm thin film of the tetracyanobicarbazyl derivative deposited onto a quartz plate are shown in Figure 4. Absorption spectrum is obtained by exciting the thin film in the range of 300-700 nm wavelengths. The photoluminescence spectrum is so determined by exciting the film around the absorption maxima.

As already reported for several donor-acceptor carbazole-based molecules, substitution in 3-position of the electron-donor heterocycle by cyanovinylene electron-withdrawing groups leads to an internal charge transfer which shows up in the visible domain of the optical spectrum [23]. In the case of (OcCz2CN) 2 , i.e 3,3'-dimer of OcCz2CN, the charge-transfer band is peaking at λ max = 450 nm. From the onset of the absorption threshold of the spectrum the energy gap value Eg was found to be close to 2.4 eV. Furthermore, it is noteworthy that π-conjugation within the bicarbazyl core provokes, as expected, a red-shift of the visible band, since the corresponding signal appears at λ max = 428 nm for films of OcCz2CN monomer (Eg=2.7 eV) [23]. Hence, a yellow-light emission is observed for (OcCz2CN) 2 with a maximum at λ max = 590 nm (Figure 4) whereas a green one is observed for the monomer derivative (λ max = 528 nm). Such an effect of the extension of π-conjugation on the shift of the emission color appeared in red-emitting Narylcarbazoles substituted in 3,(6)-position(s) by three conjugated double bonds ended by a dicyano group, recently reported by Fu et al. [24].

Cyclic voltammetry of (OcCz2CN) 2 (at a scan rate of 10 mV/s), allowed an estimation of its ionization potential IP (HOMO level). An oxidation process at E ox = 0.95 V/SCE (HOMO level at 5.7eV) is observed corresponding to the formation of carbazolic radicalcations [25]. This value is 0.10 V higher than that observed for 3,3' unsubstituted bicarbazyl, as expected from the effect of the electron-acceptor groups. Indeed, it has been shown that such substituents make the carbazole core less oxidizable and lead to an anodic shift of the oxidation potential. On the opposite, introduction of acceptor groups is expected to increase the electronic affinity (EA). Nevertheless, no clear reduction signal was observed in the electrochemical reduction domain of cv and LUMO level value was deduced from the optical energy bandgap Eg and HOMO values according to the expression E HOMO -E LUMO = Eg. Such estimation leads to a value of 3.3 eV for the LUMO.

Electroluminescence

Doped and non-doped OLED structures are depicted Figure 3.

Doped-OLEDs

EL spectrum of doped devices in configuration (a) i.e. with the following structure: ITO/CuPc/NPB/DPVBi:(OcCz2CN) 2 /AlQ 3 /LiF/Al (DPVBi doped with 1.5% and 3.3% of the yellow emitter is shown on Figure 5).

Although NPB is also a blue fluorescent emitter, the position of the EL blue luminescence peak in our devices match perfectly (after small correction due to calculated microcavity effects) the DPVBi PL peak, which was measured on a separate substrate and differs from about 20 nm from the NPB one. We consequently consider that the blue emission is due to DPVBi in our devices.

A good overlapping between (OcCz2CN) 2 absorption and DPVBi emission (inset Figure 5) leads to an efficient Förster transfer process. With increased (OcCz2CN) 2 concentration from 1.4 to 3.3 wt %, the wavelength of the peak red shifts from 520 to 540nm. This red-shift of the EL peak could be attributed to polarization phenomena [26]. Meanwhile, external quantum efficiency η ext decreases from 0.4 to 0.3 %.

Such an effect of the dopant concentration on η ext may arise from a quenching process. EL parameters of doped devices are presented in Table 1.

Doping DPVBi with 1.5% wt of rubrene, taken as a reference of highly efficient yellow emitter, led to an important improvement of all the EL parameters as e.g. a η ext = 2.5 % (Table 1), a greenish-yellow light being observed. These results show that rubrene is a much more efficient dopant than (OcCzCN 2 ) 2 although having close LUMO and HOMO levels (3.2 eV against 3.3 eV and 5.4 eV against 5.7 eV respectively). This is not surprising as it has been reported that rubrene dispersed in a host layer acts as a very efficient trap charge carrier and recombination center [27]. In order to attain a color-control we tested another type of OLED configuration i.e. a non-doped device with a thin film of (OcCz2CN) 2 inserted in the DPVBi layer (Figure 3b).

Non-doped-OLEDs

A 1 nm-thick layer of the carbazolic dimer was inserted within the DPVBi blue-emitting matrice at a position d from the NPB/DPVBi interface in the range 0 ≤ d ≤10 nm. By using an ultrathin rubrene sensing layer we recently found for the same structure that such a d range corresponds to the exciton diffusion domain within the DPVBi layer (diffusion length around 9nm) [12]. Typical emission spectra for different values of d (5, 7 and 10 nm) are shown on Figure 6 with corresponding CIE coordinates in Inset. In these configurations, the emission spectrum consists in two blue and yellow peaks. The relative intensity of these peaks is sensitive to the position of the carbazolic layer from the interface. With increased d from 5 to 10nm, the amount of blue is increased by the exciton diffusion in the DPVBi matrix. Influence of the position d of the (OcCz2CN) 2 layer on performances and CIE coordinates for several devices are shown on Table 2. In absence of yellow emitter, a deep blue luminescence with CIE coordinates (0.17; 0.15), a quantum efficiency η ext = 2.8 % typical for DPVBi were observed (see EL performances of undoped DPVBi Table 1). Pure yellow emission with CIE coordinates (0.51; 0.47) corresponding to the emission light of (OcCz2CN) 2 molecules is achieved when the carbazolic compound lies exactly at the interface (d = 0). Here, the energy transfer between DPVBi excitons and (OcCz2CN) 2 is complete, as a consequence of the efficient Förster process. Nevertheless, the quantum efficiency is η ext = 0.3% around one order of magnitude lower than that of the blue OLED, showing that important quenching limits the efficiency. A warm white standard emission closed to -illuminant A-with CIE coordinates (0.40; 0.43), a η ext = 0.7 %, a η EL = 0.7 cd/A, a luminance L = 932 cd/m 2 at 60mA/cm 2 (197 cd/m 2 at 10mA/cm 2 ) and a quite good Color Rendering Index of 74, were obtained for a position d = 10 nm. To the best of our knowledge, it is the first report of a white OLED using a carbazolic molecular compound as a thin emitting layer. Intermediary yellowish-green light (CIE: 0.47; 0.48) and yellow (CIE: 0.43; 0.46) lights were obtained for d = 5 nm and d = 7 nm respectively.

For comparison, devices utilizing 1nm of rubrene instead of (OcCz2CN) 2 have found to produce white light with similar CIE coordinates [12] and a quantum efficiency η ext = 2.6 %, at a closier d position from the emitting interface i.e. d = 5 nm from the emitting interface. We can also observe that when increasing d from 0 to 10nm, CIE coordinates of rubrene-doped DPVBi decrease from yellow (0.51; 0.48) to bluewhite (0.25; 0.23). In the case of the (OcCz2CN) 2 , the decrease is weaker, from (0.51; 0.47) to (0.40; 0.43). At d = 10nm, blue emission is more important in the rubrene configuration. Comparison between EL performances of devices utilizing our carbazole molecule with those based on the highly efficient rubrene reference show that the former are only 4 times lower. Increasing the quantum efficiency of our yellow light-emitting carbazoles would imply theoretical investigations coupled with subsequent experimental chemical engineering, which is beyond the scope of this paper.

Conclusion

3 Figure 1 .

 31 Figure 1. Molecular structures of the organics used for doped and non-doped OLEDs with (OcCz2CN) 2 as a yellow fluorophore in blue-emitting DPVBi layer.

Figure 3 .

 3 Figure 3. Configuration of the multilayer OLEDs; (a): (OcCz2CN) 2 -doped DPVBi, (b): non-doped, with a 1 nm thick (OcCz2CN) 2 film inserted in DPVBi layer. HOMO level of (OcCz2CN) 2 was obtained from cyclic voltammetry. LUMO level was deduced from the optical energy band-gap Eg and HOMO values.

Figure 4 .

 4 Figure 4. Absorption (line + square) and photoluminescence (solid line) spectra of the 50nm thick film of (OcCz2CN) 2 .

Figure 5 .

 5 Figure 5. Electroluminescence spectra at a current density of 10mA/cm² of doped OLEDs with 1.4% and 3.3% wt of (OcCz2CN) 2 in DPVBi. Inset: proposed (OcCz2CN) 2 absorption and DPVBi emission spectra.

Figure 6 .

 6 Figure 6. EL spectra for different positions d of a 1nm thick layer of (OcCz2CN) 2 from the NPB/DPVBi interface. The spectra are normalized to the yellow peak at 568nm. Inset: corresponding CIE diagram (the solid line corresponds to the edge of the spectrum locus).

A

  new carbazole derivative 6,6'-bis(1-(2,2'dicyano)vinyl)-N,N'-dioctyl-3,3'-bicarbazyl ((OcCz2CN) 2 ), was synthetised by Knovenagel condensation and used as an active molecule in multilayer organic light-emitting diodes (OLEDs). Due to the π-donor-acceptor conjugation developed between the electron-donor carbazolic moieties and the cyanovinylene electron-withdrawing groups, (OcCz2CN) 2 was found to emit a yellow light at λ max =590 nm with CIE coordinates (x = 0.51; y = 0.47) and was used either as a dopant or as an ultrathin layer in a blue-emitting matrice of 4,4'-bis(2,2'diphenylvinyl)-1,1'-biphenyl (DPVBi). DPVBi (OcCz2CN) 2 -doped structure exhibited, at doping ratio of 1.5 weight %, a yellowish-green light with CIE coordinates (x = 0.31; y = 0.51), an electroluminescence efficiency η EL =1.3 cd/A, an external quantum efficiency η ext = 0.4 % and a luminance L= 127 cd/m 2 (at 10 mA/cm 2 ) whereas for non-doped devices utilizing a thin carbazole yellow-emitting layer inserted in a blue-emitting DPVBi matrix, white OLEDs carbazole-based have been achieved for the first time. A warm white light with CIE coordinates (x = 0.40; y= 0.43), η EL = 2.0 cd/A, η ext = 0.7 %, L = 197 cd/m 2 (at 10 mA/cm 2 ) and a quite good Color Rendering index (CRI) of 74, were obtained.

  Anal. Calc. for C 48 H 48 N 6 : C 81.31, H 6.82, N 11.87; Found: C 81.79, H 6.65, N 11.53. 1 H NMR (CDCl 3 , 200MHz, δ ppm) : 8.4 (H 5,5' ), 8.1 (H 4,4' ), 7.8 (H 2,2' ), 7.6 (H 7,7' ), 7.5 (2 H vinylene), 7.3 (H 1,1' ), 7.2 (H 8,8' ), 4.1-1.1 (28H CH 2 ), 0.7 (6H CH 3 ). 13 C NMR (CDCl 3 , 50MHz, δ ppm) : 159.2 (CH-CN vinylene), 143.8 (C 8a,8'a ), 140.1 (C 9a,9'a ), 134.3 (C 3,3' ), 128.3 (C 5,5' ), 126.6 (C 7,7' ), 124.4 (C 2,2' ), 123.4 (C 6,6' ), 122.9 (C 4b,4b' ), 122.1 (C 4a,4a' ), 118.9 (C 4,4' ), 114.5 and 113.6 (CN(s)

	cis, trans), 109.6 (C 8,8' ), 109.3 (C 1,1' ), 76.3 (C-CN
	vinylene), 43.5-22.2 (CH 2 (s) octyl), 13.8 (CH 3 ).
																C 8 H 17
				C 8 H 17											N
				N		POCl 3 / DMF OHC										CHO
							ODCB		N
		N								C 8 H 17
		C 8 H 17
			C 8 H 17															C 8 H 17
	OHC	N	N	CHO +	2	H 2 C	CN CN	N(Bu) 4 OH iPrOH/MeOH	NC	NC	6 7	5 8	4b 8a	N 9	4a 9a	4 1	3 2	N	CN	CN
		C 8 H 17												C 8 H 17
																( O c C z 2 C N ) 2
				Figure 2. Synthetic route of bicarbazyldialdehyde and (OcCz2CN) 2


David Dureisseix

Béatrice Pascale Royer

Béatrice Faverjon

Numerization of a memory effect for an homogenized composite material with a large contrast in the phase thermal conductivities *

Keywords: computational material, asymptotic expansion, 2-field model, periodic media, thermal conduction problems PACS: 44.05.+e, 44.10.+, 44.35.+c MSC: 35B27, 35K05

teaching and research institutions in France or abroad, or from public or private research centers.

Numerization of a memory effect for an homogenized composite material with a large contrast in the phase thermal conductivities

Introduction

Models with memory effects (or heredity models) have been used quite extensively in structural mechanics for viscoelasticity behaviors. Especially when the structures are subjected to forced vibration loading, their response may be easily found, analytically in simple cases, or numerically. For transient responses, concurrent models, mainly with internal variables have proved their ability to be used numerically [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF].

Several justifications can be stated for the physical meaning of memory effects, one being the presence of a microstructure in the material, leading to a delayed response at the macroscopic scale, therefore related to upscaling [START_REF] Allaire | Homogenization and two-scale convergence[END_REF] and a kind of time non-locality. This is for instance the case for scattering of elastic waves for hyperbolic problems [START_REF] Fish | Uniformly valid multiple spatial-temporal scale modeling for wave propagation in heterogeneous media[END_REF], for composite materials where the constituents themselves exhibit a memory effect [START_REF] Suquet | Effective properties of time dependent heterogeneous materials[END_REF], or for particular non-linear behaviors or time-dependent material properties [START_REF] Abdessamad | Memory effect in homogenization of a viscoelastic Kelvin-Voigt model with time-dependent coefficients[END_REF].

In this article we focus on the case of evolution (elliptic) problems, even with linear constituents, but where the contrast of constitutive properties is large. The delayed effect is therefore directly provided by the homogenization procedure [START_REF] Jiang | Two-scale homogenization and memory effects of a first order differential equation[END_REF][START_REF] Auriault | Effective macroscopic description for heat conduction in periodic composites[END_REF][START_REF] Tartar | Memory effects and homogenization[END_REF][START_REF] Royer | Time analysis of the three characteristic behaviours of dual-porosity media. I: Fluid flow and solute transport[END_REF]. We proposed to study a thermal problem for illustration purpose, though the same framework also applies to different physics [START_REF] Auriault | Upscaling heterogeneous media by asymptotic expansions[END_REF].

The aim is to be able to compute numerically the memory effect function as a macroscopic material characteristic from simulations on the microstructure (i.e. a periodic cell), as it is classically done for the homogenized other parameters, such as the capacity and the conductivity.

In this article, section 2 recalls the results of the periodic homogenization method for the kind of problem that is under concern. Section 3 proposes a numerical procedure for the memory effect function determination. Finally, section 4 provides two numerical examples: (i) a 1D problem with an analytical solution, allowing to test the numerical procedure, and whose study provides approximations for determining the memory function at small and large times; (ii) a 3D problem for testing the discretization choices and for validation of the previous approximations.

Asymptotic expansion and periodic homogenization

This homogenization technique is based on the spatial description of the various fields at two scales, i.e. with two coordinate systems: one is the slow coordinates x at the scale of the whole structure (typically x ∈ [0, L] where L is a characteristic size of the whole structure), the other one is the fast coordinate y, related to the scaled microstructure cell (typically y ∈ [0, l/ε] where l is the characteristic size of the cell). The scaling parameter is ε = l/L; it is expected to be small, traducing the assumption that the scales are separated. The spatial derivative then reads ∂ ∂x + 1 ε ∂ ∂y . The temperature field is expanded with respect to the (small) parameter ε:

θ = ∞ α=0 ε α θ α (x, y, t)
where all the temperature fields θ α (x, y, t) are Ω-periodic for the variable y, where Ω is the scaled basic cell that allows the reconstruction of the whole microstructure by periodic duplication in each space dimension. For sake of simplicity, we assume that this microstructure is composed of two phases, denoted by F and S occupying the domains Ω F and Ω S ; Γ = ∂Ω S ∩ ∂Ω F is the interface. Equating all the problem equations for the different powers of ε allows to successively provide the models for the different temperature fields [START_REF] Sanchez-Palencia | Non homogeneous media and vibration theory[END_REF][START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF].

The large contrast in material characteristics arises when the ratio between the capacities or the conductivities is of the order of ε n with n = 0. A special case arises when the ratios between conductivity is of order ε 2 while the capacities are of the same order, i.e. ρ F c F ≈ ρ S c S and

k F = ε 2 k F with k F ≈ k S .
In such a case, the emergent macroscopic thermal problem is of a different kind [START_REF] Auriault | Effective macroscopic description for heat conduction in periodic composites[END_REF]. This case is usefully interpreted with characteristic times [START_REF] Royer | Time analysis of the three characteristic behaviours of dual-porosity media. I: Fluid flow and solute transport[END_REF][START_REF] Néron | A computational strategy for thermo-poroelastic structures with a timespace interface coupling[END_REF] since the characteristic macroscopic time is expected to be τ M = (C M /k M )(L/π) 2 (if C M and k M are the homogenized total capacity and conductivity), while their counterparts for the microscopic S and F phases are τ S = (ρ S c S /k S )(l S /π) 2 and

τ F = ρ F c F l 2 F k F π 2 ≈ ρ F c F L 2 F k F π 2 ≈ τ S /ε 2 ≈ τ M
since the characteristic lengths associated to each phase are l S ≈ l F ≈ l and

L S = l S / ≈ L F = l F / ≈ L.
Therefore, the transient phenomena at the micro scale for the F -phase may well arise at macro scale. Details on the successive derivations of the micro and macro problems can be found in [START_REF] Auriault | Effective macroscopic description for heat conduction in periodic composites[END_REF], and the main results are recalled in the following.

Macroscopic property of θ S0

The first result for the S-phase is that the temperature at order 0, i.e. multiplied by ε 0 , does not depend on the fast variable:

θ S0 = θ S0 (x, t)
it is therefore a macroscopic temperature field, and its gradient with respect to the fast variable y is null. Its gradient with respect to the slow variable x is a macroscopic temperature gradient

Z S0 (x, t) = grad x θ S0

Localization property of θ S1

The next order equations give the problem that θ S1 should satisfy:

div y q S0 = 0 in Ω S q S0 • n S = 0 on Γ q S0 = -k S (Z S0 + grad y θ S1 ) in Ω S θ S1 is Ω-periodic
It is a steady-state thermal problem whose variational formulation is

Ω S grad y θ • k S grad y θ S1 dΩ = - Ω S grad y θ • k S Z S0 dΩ (1)
for all Ω-periodic test functions θ (y). It is usually discretized with finite elements, and its solution is linear with respect to the macroscopic loading Z S0 . Therefore this solution, and its gradient Z S1 = grad y θ S1 can be expressed using a linear operator which is a characteristic of the S-microstructure. For instance,

Z S1 (x, y, t) = -L S (y)Z S0 (x, t)
The numerical determination of L S requires: (i) as many resolutions of the previous linear micro problem [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF] on Ω S as there are independant components in Z S0 , i.e. 2 for 2D problems, and 3 for 3D problems, and (ii) the storage of the same number of temperature gradient fields on Ω S (the columns of the matrix L S ).

This operator is a localization operator for the temperature gradient. Indeed, it allows to recover the first two terms of the full temperature gradient, obtained from the asymptotic expansion of θ S , once the macroscopic gradient is known:

Z S ≈ Z S0 (x, t) + Z S1 (x, y, t) = (1 -L S )Z S0
and of the thermal flux q S = -k S Z S ≈ q S0 = -k S (1 -L S )Z S0

Problem to solve for θ F 0

The zero-order terms for the F -phase provide the following problem:

ρ F c F θF 0 = -div y q F 1 in Ω F θ F 0 = θ S0 on Γ q F 1 = -k F grad y θ F 0 in Ω F θ F 0 is Ω-periodic
This problem couples the micro and the macro scales due to the fact that θ F 0 still depends on the y coordinate, and due to the presence of the θ S0 (x, t) term as a boundary condition on Γ. On the F -cell, this term is uniform due to the scale separation. A change in variable may be useful: if we denote

θ F 0 (x, y, t) = θ S0 (x, t) + w(x, y, t) (2)
the new problem consists in finding w which is Ω-periodic, null on Γ, such that

ρ F c F θF 0 + ρ F c F ẇ = div y (k F grad y w) in Ω F
Using a test function which is Ω-periodic and null on Γ, the variational formulation reads:

Ω F w ρ F c F ẇdΩ + Ω F grad y w • k F grad y wdΩ = = - Ω F w ρ F c F θF 0 dΩ (3)
2.4 A macro model with a 2-field temperature, or with a memory function

The macroscopic model is obtained by ensuring that the next higher order problem for θ S2 is well-posed. It reads:

-div y q S1 -div x q S0 = ρ S c S θS0 in Ω S q S1 • n S + q F 1 • n F = 0 on Γ q S1 = -k S (grad x θ S1 + grad y θ S2 ) in Ω S θ S2 is Ω-periodic
The partial averaging operations are defined as

• S = 1 vol(Ω) Ω S • dΩ and • F = 1 vol(Ω) Ω F • dΩ, so that 1 vol(Ω) Ω • dΩ = • S + • F .
The compatibility condition is obtained by averaging the first equation (and using periodicity), leading to [START_REF] Auriault | Effective macroscopic description for heat conduction in periodic composites[END_REF]:

- 1 vol(Ω) Γ q S1 • n S dΓ + div x ( k S (1 -L S ) S Z S0 ) = ρ S c S S θS0
with the macroscopic parameters

C M S = ρ S c S S and k M S = k S (1 -L S ) S (4)
The average of the previous problem giving θ F 0 reads:

- 1 vol(Ω) Γ q F 1 • n F dΓ = ρ F c F θF 0 F = 1 vol(Ω) Γ q S1 • n S dΓ
Finally, the macroscopic problem reads:

C M S θS0 = div x (k M S Z S0 ) -ρ F c F θF 0 F (5)
It also involves the second thermal field rate θF 0 defined on the microstructure. This problem can be seen as a macroscopic thermal problem with additional internal variables at each macro material point which are the underlying values of θF 0 on a full micro cell. Once discretized by finite elements, this involves at each macro integration point the storage of a full field θF 0 at each nodes of the underlying discretized cell, as used in the so-called FE 2 approach [START_REF] Feyel | FE 2 multiscale approach for modelling the elastoviscoplastic behaviour of long fiber SiC/Ti composite materials[END_REF]. In our context, this approach is too storage consuming. Note that in the past, memory effect models were designed for visco-elasticity problems, but have been superseded by models with internal variables [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF]. The converse way is used herein, since we promote an available memory effect model [START_REF] Auriault | Effective macroscopic description for heat conduction in periodic composites[END_REF]:

C M S θS0 = div x (k M S Z S0 ) + τ =t τ =-∞ θS0 (τ )β(t -τ )dτ (6)
where β(t) is a memory function, that should depend only on the macro material (independent of boundary conditions) and that should vanish for large t (evanescent memory). Concerning practical initial value for this problem, we may consider that a stationary state is given for t ≤ 0, i.e. θS0 (t ≤ 0) = 0 and θF 0 (t ≤ 0) = 0, so that the back-in-time can start for τ = 0, and we may define the searched temperature differences as the differences with respect to these initial temperatures, so that one also has:

θ S0 (t ≤ 0) = 0 and θ F 0 (t ≤ 0) = 0.
The results given in [START_REF] Auriault | Effective macroscopic description for heat conduction in periodic composites[END_REF] and reused previously are established using a Fourier transform method. This approach also enabled the author to derive the memory function as an infinite sum of a series, for a simple 1D problem. In this article, we aim to establish the possibility of computing the memory function in a general way, using only transient problems on a unique cell, in a similar fashion as the two other homogenized parameters C M S and k M S that are derived from stationary simulations on the same unique cell.

Numerical determination of macro material parameters

As mentioned before, the parameters C M S and k M S can be computed in a standard fashion, solving simple steady state problems on a single cell with periodic boundary conditions. The last parameter is now the memory function β to be determined.

Discretization of the memory effect

The first step is the discretized-in-time version of the memory effect. We propose to use a simple implicit Euler integration scheme. Note that the same procedure can be used for any other scheme. We denote t k the time instant number k, and h the time step (assumed constant for sake of simplicity), and we use the superscript or the subscript k to denote the value of a quantity at time t k . Moreover the considered instant for unknown quantities is t i+1 , and the past instants are t k for k = 0, ..., i; by convention, t 0 = 0.

The time integration scheme relates the values for the derivative of a quantity w to the values of the same quantity. Here it reads: w i+1 = w i + h ẇi+1 or equivalently ẇi+1 = 1 h (w i+1 -w i ). Concerning the heredity integral with the second derivative of the micro temperature field, dealing with first order derivatives is easier for a time integration scheme dedicated to parabolic equations. We therefore discretize the integral using only values for the first order derivative of the temperature. To do so, we integrate by parts to formally use the derivative of the memory function γ(t) = β(t) , and the time integration scheme:

t 0 θS0 (τ )β(t -τ )dτ = θS0 (t)β(0) -θS0 (0)β(t) + i k=0 t k+1 t k θS0 (τ )γ(t i+1 -τ )dτ ≈ θi+1 S0 β 0 -θ0 S0 β i+1 + i k=0 θi+1-k S0 γ k h (7)
and

β k = β k-1 + hγ k .

Discretization of the micro problem of the F -phase

On a second hand, we wish to express the solution of the transient thermal problem on the F -cell with a recursion formula (in a back-in-time expression). To do so, we consider a spatial finite element discretization of the corresponding variational problem. It reads at time t i+1 :

w T (C ẇi+1 + Kw i+1 ) = -w T Cu θi+1 S0 (8)
C and K are the capacity and conductivity matrices using the material parameters ρ F c F and k F on the Ω F domain, w is the column vector of nodal values for the thermal field w, u is a unitary vector column of nodal values, so that u θi+1 S0 is a discretized uniform field with value θi+1 S0 . The thermal field w should be Ω-periodic and null on Γ, therefore the corresponding column vector can be built from a smaller size vector, w, getting rid of null values on the nodes on Γ and dependent values due to periodicity. Formally, we may use a boolean rectangular mapping matrix B (with BB T being the identity on a reduced subspace) to express these dependencies as w = B T w; so does the test functions: w = B T w . Therefore, the problem reads

(BCB T ) ˙ w i+1 + (BKB T ) w i+1 = -BCu θi+1 S0
Using now the time integration scheme, it transforms into:

(B 1 h CB T ) w i+1 = -BCu θi+1 S0 + 1 h BCw i
where C = C + hK is invertible, so one gets the expression of w i+1 , then of w i+1 = B T w i+1 :

w i+1 = Aw i -hAu θi+1 S0 with A = B T (B CB T ) -1 BC (9)
Using this recursion formula for w together with the time integration scheme allows to find the recursion formula for ẇ as well, leading after algebraic manipulations to:

ẇi+1 = A i (A -1) 1 h w 0 - i k=1 A k (A -1)u θi+1-k S0 -Au θi+1 S0 Noting that w 0 = θ 0 F 0 -θ 0 S0 = 0, one finally gets θi+1 F 0 = - i k=0 A k (A -1)u θi+1-k S0 (10)
Note that for a parabolic equation of a time-dependent problem such as this micro problem, spurious but bounded oscillations can appear for a too short time step due to the non satisfaction of the discrete maximum principle [START_REF] Ciarlet | Maximum principle and uniform convergence for the finite element method[END_REF][START_REF] Thomée | On the existence of maximum principles in parabolic finite element equations[END_REF]. The rule-of-thumb (obtained for a discretized 1D problem) is that [START_REF]Abaqus/Standard -User's Manual[END_REF][START_REF] Vermeer | An accuracy condition for consolidation by finite elements[END_REF] h ≥ (ρ F c F e 2 )/(6k F α) where e is the element size and α the parameter for the θ-method time integration scheme: here, α = 1 for implicit Euler method. Therefore, the bound can be stated as

h τ F ≥ π 2 6 e L F 2 = π 2 6 e l F 2 ε 2
Otherwise, a lumping of the capacity matrix (or a lower order estimates) allows to remove the bound [START_REF] Fujii | Some remarks on finite element analysis of time-dependent field problems, Theory and practice in finite element structural analysis[END_REF], for instance using a numerical integration with integration points located at nodes for the capacity model. In the following, we will always satisfy the previous criteria.

Identification of the two models

Identifying the two previous discretized versions, [START_REF] Auriault | Effective macroscopic description for heat conduction in periodic composites[END_REF] and ( 10), of the same term in ( 5) and ( 6):

τ =t τ =0 θS0 (τ )β(t -τ )dτ = -ρ F c F θF 0 F ( 11
)
leads to:

β 0 -hγ 0 = -ρ F c F (1 -A)u F hγ k = ρ F c F A k (A -1)u F for k = 1, .

.., i

There is still an indeterminacy for β 0 and γ 0 that can be overcome by ensuring the vanishing property of the memory function β. Indeed, using the incremental relationship β k+1 = β k + hγ k+1 , one gets

β i+1 = β 0 + i+1 k=1 hγ k = β 0 -ρ F c F (1 -A i+1 )Au F
If the spectral radius of the iterated operator A is ρ(A) < 1, then the limit of β i+1 for large values of i is β 0 -ρ F c F Au F . Nullifying this value for a vanishing memory function leads to choose

β 0 = ρ F c F Au F = 1 vol(Ω)
u T CAu so that for k = 0, . . . , i

β k = ρ F c F A k+1 u F = 1 vol(Ω) u T CA k+1 u (12)
A rough proof of the contraction property of the operator A is the following: within the subspace of Ω-periodic thermal fields that are null on Γ, consider any thermal field w = 0 and denote v = Aw with A = C -1

C and C = C + hK (C is symmetric definite positive, and K is symmetric positive, and definite on the considered subspace, provided that the measure of Γ is not null, i.e. that the phases are not separated). One has

(w -v) T C(w -v) = w T Cw + v T Cv -2v T Cw = w T Cw + v T Cv -2v T (C + hK)v = w T Cw -v T Cv -2hv T Kv
We can therefore conclude that v T Cv = w T Cw -(w -v) T C(wv) -2hv T Kv. Provided that v = w, one gets v T Cv < w T Cw. Assume now that w is an eigenvector of A associated to the eigenvalue λ = 1: v = λw and the previous inequality leads to λ 2 w T Cw < w T Cw, so λ 2 < 1. To complete the proof, consider the case where λ = 1 and v = w; the previous developments allow to conclude that v T Kv = 0 and therefore (on the dedicated subspace) v = 0, so w = 0 which is contradictory, so λ = 1. The final conclusion is therefore: ρ(A) < 1.

The simple overall pseudo-code to build the values of the memory function is summarized in the algorithm 1. Since this memory function is interpreted as a macroscopic homogenized material characteristics, it is determined with computations at the cell microstructure level, as the other homogenized material parameters are, such as macroscopic conductivity and macroscopic capacity.

Algorithm 1 Memory function determination

1: build unitary column vector u 2: n w ← u 3: compute r ← Cu 4: for k = 0, 1 . . . i do 5:

solve w T Cw = w T r to get w, Ω-periodic and null on Γ 6:

n old w ← n w , n w ← w , ρ k ← n w /n old w 7: compute r ← Cw 8:
compute β k ← 1 vol(Ω) u T r 9: end for Note that the characteristic time τ F 0 for the discretized-in-space problem is linked to the previous operators in the left hand side. Indeed, one can defined it as being the inverse of the minimal eigenvalue of (BCB T ) -1 (BKB T ). So, the spectral radius of A is ρ(A) = (1 + h/τ F 0 ) -1 < 1 and

τ F 0 = h ρ(A) -1 -1
More precisely, the restriction of A in the reduced subspace of admissible fields is

BAB T = 1 + h(BCB T ) -1 (BKB T )
Usually, the matrix A is not explicitly available, but its spectral radius can be obtained with the iterates of the power method. This is also what is done in the algorithm 1 for two successive vectors w k and w k+1 : when sufficiently iterated, the spectral radius is close to ρ k+1 = w k+1 / w k . The stabilization of this value therefore denotes that the evolution of the function β is found for "large" values of t/τ F therefore leading to an exponential decay for β(t) of the form exp(-t/τ F 0 ).

Memory effects on the different macroscopic fields

The previous developments consider θ S0 as a macroscopic temperature field. Indeed, it is independent of the fast coordinate. For the other phase, an average of the temperature is also a macroscopic field, i.e. θ F 0 F , as well as the average of the first temperature terms of both phases, i.e.

θ = θ S0 S + θ F 0 F = (1 -n)θ S0 + θ F 0 F
where n is the fraction of F -phase.

Considering the recurrence formula (10) (and its counterpart with no time derivative), one gets the linear combinations θi+1

F 0 F = i k=0 α k θi+1-k S0 and θ i+1 F 0 F = i k=0 α k θ i+1-k S0 with α k = -A k (A -1)u F .
In the special case where the microstructure is identically the same at any macroscopic point, and when ρ F c F is constant, the macroscopic problem [START_REF] Abdessamad | Memory effect in homogenization of a viscoelastic Kelvin-Voigt model with time-dependent coefficients[END_REF], can be written at the different time steps i + 1 -k and weighted with the coefficients α k to give leading when summing for k = 0, ..., i to:

α k C M S θi+1-k S0 = div x (k M S grad x (α k θ i+1-k S0 )) -ρ F c F α k θi+1-k F 0 F
C M S θi+1 F 0 F = div x (k M S grad x θ i+1 F 0 F ) - i k=0 ρ F c F α k θi+1-k F 0 F
whereas the first macroscopic problem ( 5) was

C M S θi+1 S0 = div x (k M S grad x θ i+1 S0 ) - i k=0 ρ F c F α k θi+1-k S0
i.e. with the same heredity term function in the right-hand-side. As a consequence, the macroscopic problem for the total average temperature is also of the same form:

C M S θi+1 = div x (k M S grad x θi+1 ) - i k=0 ρ F c F α k θi+1-k 4 Numerical examples 4.

Bilaminated composite

This 1D test case is defined in [START_REF] Auriault | Effective macroscopic description for heat conduction in periodic composites[END_REF] for which a Fourier analysis allows to have an explicit expression for the memory function 1 . It concerns a stack of 2 kinds of layers (S and F ) with a volumic fraction n of phase F , see figure 1. This example may serve as a validation test for the proposed approach to compute numerically the memory function, as well as to illustrate the dependence on the discretization parameters, i.e. the time step h and the finite element size e.

When the macroscopic temperature is considered constant, but a function of a parameter τ , a Fourier analysis gives the particular expression of the searched temperature for phase F as Φ(x, t; τ ) = φ(x, t)θ S0 (τ ) with

φ(x, t) = - ∞ q=0 u q (x)v q (t) u q (x) = 4 (2q + 1)π sin[(2q + 1)πx/L F ] v q (t) = exp[-(2q + 1) 2 t/τ F ]
To find the solution when the macroscopic temperature depends on time t, the use of Duhamel integral provides

w(x, t) = ∂ ∂t t 0 Φ(x, t -τ ; τ )dτ
With an initial condition θ S0 (x, 0) = 0, and v q (0) = 0, this leads to

w(x, t) = t 0 θS0 (τ )φ(x, t -τ )dτ ẇ(x, t) = θS0 (0)φ(x, t) + t 0 θS0 (τ )φ(x, t -τ )dτ θF 0 (x, t) = ẇ(x, t) + θS0 (t)
1 a typo in the original paper is corrected here, so intermediate steps are provided herein to obtain the analytical solution.

so that

-ρ F c F θF 0 (x, t) F = t 0 θS0 (τ ) ∞ q=0 ρ F c F u q (x) F v q (t -τ )dτ + + θS0 (0) ∞ q=0 ρ F c F u q (x) F v q (t) -ρ F c F n θS0 (t)
Hence

β(t) = ∞ q=0 ρ F c F u q (x) F v q (t) = ρ F c F n 8 π 2 ∞ q=0 1 (2q + 1) 2 exp[-(2q + 1) 2 t/τ F ] (13)
for which

β(0) = ρ F c F n 8 π 2 ∞ q=0 1 (2q + 1) 2 = ρ F c F n
If only one term of the development is kept for an approximate analytical solution, it reads:

β(t) ≈ ρ F c F n 8 π 2 exp[-t/τ F ] (14)
These solutions are depicted on Figure 2.

For validation purposes, a 1D discretization of Ω F = [0, L F = nL] is done with n e 2-node finite elements of uniform size e = nL/n e , and a time discretization of [0, τ F ] with n h time steps, also with a uniform time step h = τ F /n h . The error with respect to the exact solution is depicted on figure 3. It is concentrated at the beginning of the time evolution of β(t). This is due to the fact that the elementary solution that is to be captured by the numerical method is a somehow thermal boundary layer (boundary conditions with a different value than the initial condition), for which a fine discretization is useful to capture the good average. Such boundary layers can be captured using mesh local refinements or with enriched FE methods [START_REF] El-Zein | Exponential finite elements for diffusion-advection problems[END_REF][START_REF] Turner | Analytic enrichment of finite element formulations for capturing boundary layer behavior[END_REF]. Here, we focus on a semi-analytical approximation to get the first values of the memory effect function β(t) for the beginning of its evolution, i.e. for a "small" time t. Indeed, analytical solutions are available for an infinite half plane, i.e. locally to each point of the boundary Γ when the boundary layer is close to the boundary: in the direction of the normal to the surface, with a coordinate x, it is φ(x, t) ≈ erf(u); erf is the error function, and

u = x (4k F /ρ F c F )t = π 2 x/L F t/τ F
Using a similar development as before, this allows to get an approximate solution for the memory function:

β(t) ≈ ρ F c F 1 vol(Ω) Γ L 0 φ(x, t)dxdS = ρ F c F surf(Γ) vol(Ω) L 0 φ(x, t)dx
L is a normal distance to the boundary, allowing to have a significant integral on the whole domain Ω F . For the previous 1D problem, there is no surface integral but two boundaries, and L = L F /2 (each boundary is associated to half of the domain). For 3D problems, and for consistency reasons (exact integration of a constant field), one can choose

L = vol(Ω F )/surf(Γ) = n × vol(Ω)/surf(Γ) (15)
An analytic expression of the integral is available, so that

β(t) ≈ ρ F c F surf(Γ) L vol(Ω) erf 1 + √ π exp - 1 2 -1 (16)
L 4k F ρ F c F t = 2 π t/τ F L/L F
which is expected to be accurate for being small. This asymptotic solution is also depicted on figure 2 for the 1D test case, for which surf(Γ) L/vol(Ω) is replaced with n = L F /L. Indeed, we checked that the approximation ( 16) fits the exact solution at small time values (say for t/τ F ≤ 0.2), while the exponential decay approximation is correct for large time values (say for t/τ F ≥ 0.4).

A 3D connected phase with low conductivity

For 3D problems, approximate close forms of the solution -asymptotic solution ( 16) and the exponential decay exp (-t/τ F 0 )-are therefore driven by the correct choices of characteristic length L and characteristic time τ F 0 . Nevertheless, for a complex microstructure, their evaluation may not be trivial, so the previous transient numerical simulation on the microstructure is still of interest. The proposed test case is depicted on figure 4. The F -phase is arranged along a cubic connected lattice, and the S-phase is the matrix which surrounds. In the proposed case, with a scalar temperature field as unknown, only 1/8 of a cell can be simulated, discarding the periodicity conditions that are automatically satisfied due to the planar symmetries. The fraction of F -phase is n = 0.354.

Since the spectral radius of A is (1+h/τ F 0 ) -1 , the smaller h is, the slower the iterative power method converges (i.e. the simulation gives the initial evolution of β(t), not the asymptotic convergence rate when t is large). Therefore, an initial computation can be done with a somehow large time step h, not to provide the accurate memory function evolution, but to reach the exponential decay regime, to give the value of τ F 0 (depending on the spatial mesh size, but not on the time step). Then, a second simulation with a small time step will give the evolution of β(t) for small time t. Eventually, this last computation can be replaced with the approximation of the memory function for small time t (16); this requires a guess for the characteristic length L, that can be given with [START_REF] Ciarlet | Maximum principle and uniform convergence for the finite element method[END_REF].

In the following tests, we first compute a characteristic time value τ F 0 with a coarse time step and different mesh sizes to assess the precision of the estimate. Second, we run several simulations with different fine time steps and mesh sizes to study the convergence to the memory function, and compare the results with the previous approximation [START_REF] Thomée | On the existence of maximum principles in parabolic finite element equations[END_REF] for small time values. All pre and post-processings have been done with Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF] and Paraview [START_REF] Henderson | The ParaView Guide[END_REF] codes, and processing with the FE code Cast3M [START_REF]Cast3M web site[END_REF].

A first set of simulations is performed with 40 time steps only on a time interval [0, 10τ F ], with an estimated τ F obtained with a length L F taken as 0.4 times the cell edge length. Several meshes (with P1 4-node tetrahedral elements) are used: regular meshes with different refinement levels (meshes 1 to 3), and a mesh with a size gradient (mesh 4), the smallest elements being located on the interface Γ, figure 5. The obtained values as well as the meshes characteristics are reported in table 1, and time evolutions are plot on figure 6. The case of the overkilling mesh 3 is considered as the reference. Clearly, the mesh with a size gradient is not the most suited to estimate the exponential decay; this is due to the fact that at a long term evolution, the temperature field (i.e. the mode associated to the eigenvalue giving the spectral radius of A) is no more a boundary layer, but spreads on the whole cell, figure 7. A uniform mesh is therefore suited to capture it.

The second set of simulation is performed along a shorter time interval but with a much fine time discretization, also reported in table 1. The obtained memory functions are depicted on figure 8.

The result obtained with mesh 5 is considered as a reference. For this simulation, an adapted mesh such as mesh 4 is now suited to capture the solution, figure 7. The asymptotic solution for small time values is still a good approximation of the reference solution up to t/τ F ≤ 0. quantity to estimate, or to evaluate using the previous iterative power method on a uniform mesh an a coarse time step. Once obtained, it can be used to determine its associated characteristic length L F 0 = π τ F 0 ρ F c F /k F ; then one may choose τ F = τ F 0 and L F = L F 0 , which, together with the previously proposed estimation of L allows to reconstruct the asymptotic solution.

Applications

The applications of the proposed numerical modeling approach participate to the material design theme. Indeed, memory effects may be viewed as either a goal for design of functional materials (indeed, closedform expressions also allows optimization for goal-oriented design microstructures), or as a drawback that have to be checked for structural integrity.

For instance, thermal delayed effects obtained by highly heterogeneous materials can be a simpler alternative than those obtained by phase changes for applications in building thermal efficiency [START_REF] Sharma | Review on thermal energy storage with phase change materials and applications[END_REF][START_REF] Kuznik | A review on phase change materials integrated in building walls[END_REF]. If the characteristic time of memory function could be of the order of the day, the contrast in day and night temperatures could be counterbalanced.

Other diffusion-involved physics are also under concern: for instance the porous media flows [START_REF] Royer | Time analysis of the three characteristic behaviours of dual-porosity media. I: Fluid flow and solute transport[END_REF] that may exhibit such delayed effects at several scales, for double permeability media, to be taken into account for contaminant dispersion in soils, or for oil recovery [START_REF] Enamul Hossain | Inclusion of the memory function in describing the flow of shear-thinning fluids in porous media[END_REF], the blood perfusion in biological tissues or the compact bone porous media [START_REF] Rohan | Multiscale FE simulation of diffusion-deformation processes in homogenized dual-porous media[END_REF].

For the case of multifunctional designs, some microstructured composite materials may be used for other purpose-oriented applications, leading to material phase selections that exhibit a large contrast in thermal properties [START_REF] Ellis | Thermal conductivity and thermal expansion of graphite fiber/copper matrix composites[END_REF]. The thermomechanical behavior of these materials have to be checked and memory effects have to be taken into account for application with thermal transient loadings. Indeed, analyses for a phenomenological model (without memory effect) and an homogenized one may lead to quantitative different responses [START_REF] Auriault | Double conductivity media: a comparison between phenomenological and homogenization approaches[END_REF].

Conclusions

This article proposes a systematic numerical approach to derive the macroscopic memory effect function for composite material with two phases and a large contrast in their thermal properties. Using the periodic homogenization procedure, it derives the associated microscopic transient problem on a representative cell. The proposed numerical algorithm is related to the spectral properties of the incremental operator traducing the thermal transient behavior of the low thermal conduction phase.

Approximate simple forms of the memory function are provided, one for the short time memory effect, and one for the vanishing time effect. They do not require any 3D simulation, but they condense the microstructure complexity in two scalar parameters, namely a characteristic length and a characteristic time. Though estimates of these two parameters are provided, their validity on a large range of various complex microstructures is still to be validated. This validation is nevertheless exemplified in this article for two different cases: a simple 1D problem for which an exact analytical solution is available, and a 3D problem for which a 3D numerical simulation is performed.

Figure 1 :

 1 Figure 1: 1D test case of a bilaminated composite (infinite in vertical direction, periodic in horizontal direction).

Figure 2 :

 2 Figure 2: Analytical results for the 1D test case of a bilaminated composite.

Figure 3 :

 3 Figure 3: Error between the exact solution and the numerically obtained solution with different uniform discretizations.

Figure 4 :Figure 5 :

 45 Figure 4: Material microstructure for phase F . A group of 8 cells is depicted.

Figure 6 :

 6 Figure 6: Obtained evolutions of the memory function at large time, using coarse time steps. Influence of the spatial discretization.

Figure 7 :

 7 Figure 7: Left: mode for the long term temperature evolution, as computed on mesh 2. Right: transient solution at t/τ F = 2.67 10 -2 , as computed on mesh 4.

Figure 8 :

 8 Figure 8: Obtained evolutions of the memory function at small time. Influence of the spatial discretization and comparison with the approximation for small time values.

Table 1 :

 1 2. τ F 0 is still a Discretized problem characteristics and obtained exponential decay.

	Mesh	Nb of elem.	h/τF	τF 0/τF LF 0/ L
	1 (uniform)	230	0.25	0.905	2.55
	2 (uniform)	1894	0.25	0.968	2.64
	3 (uniform, ref.)	1.24 10 6	0.25	0.991	2.67
	4 (not uniform)	2123	0.25	0.935	2.59
	4 (not uniform)	2123	1/75	-	-
	5 (not uniform, ref.)	0.295 10 6	1/1500	-	-

* D. Dureisseix gratefully acknowledges support and partial funding from Carnot Institute 'Ingénierie@Lyon' (Carnot I@L), and B. Faverjon gratefully acknowledges the French Education Ministry, University of Lyon, CNRS, INSA of Lyon and LabEx iMUST for the CRCT and the out mobility grant.