Marwa Abida
email: marwa.abida@unicaen.fr

Florian Gehring

Jamel Mars

Alexandre Vivet

Fakhreddine Dammak

Mohamed Haddar

Hygro-mechanical coupling and multiscale swelling coefficients assessment of flax yarns and flax / epoxy composites

Keywords: A. Biocomposite, Natural fibres, B. Mechanical properties, D. Moisture

Introduction

The demand for bio-based materials in semi-structural and structural applications is constantly growing to conform to new environmental policies. Among them, plant fibre reinforced composites (PFRC) are attractive due to their relatively low densities and low environmental impact. Plant fibres appear to be suitable as reinforcing materials for polymers [START_REF] Bledzki | Composites reinforced with cellulose based fibre[END_REF][START_REF] Van De Velde | Wettability of natural fibres used as reinforcement for composites[END_REF][START_REF] Baley | Analysis of the flax fibres tensile behaviour and analysis of the tensile stiffness increase[END_REF][START_REF] Misnon | Potentiality of utilising natural textile materials for engineering composites applications[END_REF]. The durability of plant fibres when subjected to changing environments (relative humidity, temperature) is nevertheless a major issue.

The service life of PFRC could be strongly affected by exposure to moisture. Flax fibres and their associated products (yarns, fabrics and Flax Fibre Reinforced Composites (FFRC)) are characterised by a high affinity to water molecules. This hydrophilicity is mainly attributed to the flax fibre chemical composition and structure [START_REF] Hill | The Water vapor sorption behavior of natural fibers[END_REF][START_REF] Azwa | A review on the degradability of polymeric composites based on natural fibres[END_REF]. Flax fibres are considered as biocomposites constituted of cellulose 2 microfibrils held together by a matrix of hemicellulose and lignin [START_REF] Jayaraman | Manufacturing sisal-polypropylene composites with minimum fibre degradation[END_REF][START_REF] Baley | Influence of chemical treatments on surface properties and adhesion of flax fibre-polyester resin[END_REF][START_REF] John | Biofibres and biocomposites[END_REF]. It is the capacity of the molecular sites of the polymer to establish Hydrogen bonds with the water molecules that controls, in a way, the diffusion phenomenon [START_REF] Célino | The hygroscopic behavior of plant fibers: a review[END_REF]. Amorphous cellulose and hemicellulose are the principal contributors to the moisture absorption, due to their open structure containing hydroxyl (OH) and acetyl (C2H3O) groups [START_REF] Fotouh | Fatigue of natural fiber thermoplastic composites[END_REF]. In addition, water absorption in plant fibres could be due to capillary effects by the lumen for instance [START_REF] Bessadok | Effect of chemical treatments of Alfa (Stipa tenacissima) fibres on water-sorption properties[END_REF][START_REF] Morton | Physical Properties of Textile Fibres (Fourth Edition)[END_REF].

In literature, studies focus on influence of water absorption by immersion on the mechanical properties of PFRC [START_REF] Dhakal | Effect of water absorption on the mechanical properties of hemp fibre reinforced unsaturated polyester composites[END_REF][START_REF] Assarar | Influence of water ageing on mechanical properties and damage events of two reinforced composite materials: Flax-fibres and glass-fibres[END_REF][START_REF] Christian | Moisture diffusion and its impact on uniaxial tensile response of biobased composites[END_REF][START_REF] Alhuthali | Characterisation of the water absorption, mechanical and thermal properties of recycled cellulose fibre reinforced vinyl-ester eco-nanocomposites[END_REF][START_REF] Duigou | Long term immersion in natural seawater of Flax/PLA biocomposite[END_REF][START_REF] Mannberg | Service life assessment and moisture influence on bio-based composites[END_REF][START_REF] Kumar | A comparative analysis on tensile strength of dry and moisture absorbed hybrid kenaf/glass polymer composites[END_REF][START_REF] Sodoke | Wetting/drying cyclic effects on mechanical and physicochemical properties of quasi-isotopic flax/epoxy composites[END_REF]. Most of them indicated the decrease of PFRC performance with water absorption.

Nevertheless, these experimental conditions do not match the use conditions of those materials. Few studies focus on studying the effect of humid atmospheres on hygroscopic and mechanical behaviours of FFRC. Recently, Abida et al. [START_REF] Abida | Effect of hygroscopy on nonimpregnated quasi-unidirectional flax reinforcement behaviour[END_REF] showed a strong influence of water content on flax yarn fabric behaviour when exposed to a humid atmosphere. This behaviour modification can be explained by the physical and / or chemical modifications within the reinforcement [START_REF] Alix | Effect of chemical treatments on water sorption and mechanical properties of flax fibres[END_REF][START_REF] Placet | Influence of environmental relative humidity on the tensile and rotational behaviour of hemp fibres[END_REF][START_REF] Duigou | Hygroscopic expansion: A key point to describe natural fibre/polymer matrix interface bond strength[END_REF][START_REF] Péron | Study of hygroscopic stresses in asymmetric biocomposite laminates[END_REF]. The effect of water uptake on the surface properties of natural fibres such as flax and hemp is studied by Bismarck et al. [START_REF] Bismarck | Surface characterization of flax, hemp and cellulose fibers; Surface properties and the water uptake behavior[END_REF].

Water absorption by flax fibres causes a dimensional change. The ratio of length variation of a sample to its mass variation (due to water absorption) is defined as the hygroscopic expansion coefficient, also known as the swelling coefficient β . Le Duigou et al. [START_REF] Duigou | Hygroscopic expansion: A key point to describe natural fibre/polymer matrix interface bond strength[END_REF] showed that the residual stresses at the fibre / matrix interface are strongly influenced by the hygroscopic expansions and therefore this β coefficient.

Stuart et al. developed a modified Washburn equation (MWE) that takes into account swelling of flax and

viscose fibres [START_REF] Stuart | Modelling of wicking and moisture interactions of flax and viscose fibres[END_REF]. It was proven that wicking and swelling mechanisms of flax fibres with lumen and viscose fibres with no lumen are different. In order to model FFRC behaviour with consideration of hygro-mechanical coupling, this coefficient needs to be identified.

In a multi-scale experimental approach, the present study proposes to investigate the influence of water absorption on the hygroscopic and mechanical behaviours of flax yarns and composites reinforced with a quasi-unidirectional flax yarn fabric. The water absorption of flax yarn and flax UD fabric is studied, in order to analyse the reinforcement with the same microstructure as in composite. Firstly, a direct method for determining the hygroscopic expansion coefficients of flax yarns and composites reinforced with a quasi-unidirectional fabric of flax yarns is proposed. The method consists in sorption isotherm coupled with swelling measurements. The effect of water content on mechanical behaviour of flax/epoxy composite is then investigated.

Materials and methods

Materials

Quasi-unidirectional flax reinforcement

The studied material is a quasi-unidirectional flax fabric FUD 180 supplied by LINEO NV Company (Figure 1). This flax fabric is characterised by a basic weight of 180 g•m -2 . The abbreviation FUD corresponds to quasi-unidirectional reinforcement. The associated number is the theoretical basic weight of the fabric. On our sample batch, a basic weight of 207 ± 0.5 g•m -2 is measured on two different random samples of 1 m 2 surface each. This fabric is constituted of twisted yarns. As shown in Figure 1, few yarns in weft direction are used to facilitate the fabric handling. A density of 1.52 is measured on the quasi-UD flax fabric using a Helium pycnometer called Pycnomatic Evo from Thermo-Scientific. The value obtained is in accordance with those obtained by Amiri et al. [START_REF] Amiri | Standard density measurement method development for flax fiber[END_REF] and Legall et al. [START_REF] Gall | Recommended flax fibre density values for composite property predictions[END_REF]. The technical data of flax fabric are given in Table 1. The filling rate is the orthogonal projection of the fabric divided by the total area. It is determined theoretically using the flax density, the yarn metric number and the pick count. The balancing coefficient KG is defined as the warp pick count divided by the total pick count. This fabric is characterised by a pre-treatment of epoxy based sizing which justifies the choice of the matrix described in the next section.

Epoxy matrix

The matrix used for the elaboration of the composite laminates is a mixture of an epoxy resin developed by Hunstman and commercialised under the name Araldite LY1564 and a hardener developed by the same company and commercialised under the name Aradur 3487. As mentioned above, the resin was chosen for compatibility reasons with the flax fabric reinforcement. It is characterised by low toxicity and a good performance / cost ratio. The viscosity of the Araldite LY1564 resin at room temperature (T=25°C) is between 10 and 20 mPa.s and its density is between 1.1 and 1.2 g.cm -3 . The viscosity of the hardener Aradur 3487 at T=25 °C is between 30 and 70 mPa.s and its density is between 0.98 and 1 g.cm - adapted for contact moulding. Due to its low viscosity, impregnation and fibre / matrix adhesion is improved.

Specimens manufacturing

Flax yarns

The behaviour of elementary yarns constituting the flax reinforcement fabric is studied. Flax yarns are randomly extracted from the quasi-unidirectional flax fabric. The latter are handled very carefully in order to avoid their damage (untwisting of the yarn, partial rupture…). At least 30 specimens per water content are tested. An average length of 150 mm is chosen for flax yarn specimens. Yarn identification is performed by giving each sample a number so that it can be distinguished during conditioning and mechanical testing.

Flax fabric

Hygroscopic behaviour of quasi-unidirectional flax fabric reinforcement specimens is studied. Specimens of 80×25 mm 2 are cut randomly in the flax yarn fabric roll. Like flax yarns, flax fabric specimens are carefully handled in order to avoid their damage. At least 5 specimens are examined per relative humidity. The reproducibility of these tests is considered satisfactory.

Epoxy resin

For pure resin plates manufacturing, the two constituents of the matrix (resin and hardener) are mixed in accordance with the mass proportions recommended by the manufacturer, the mixture is degassed for 30 min under vacuum in order to eliminate the air bubbles introduced during mixing. The resin is then casted and left to rest for 2 hours in order to allow the elimination of any air bubbles introduced during casting. The same thermal cycle as that of the composite laminates is respected, namely curing for 2 hours at 100 °C.

Composites

The composite laminates are manufactured following a two-step process. It consists of a hand lay-up of successive plies process (8 plies) followed by a curing by thermocompression. During the draping and impregnation phase (hand lay-up), the reinforcements are positioned so as to ensure the same orientation of the warp yarns of all the plies. This unidirectional stratification is respected during the fabrication of all plates required for this study. Prior to stratification, flax fabric bands of 330×250 mm² are firstly cut and then dried at 110 °C for 1 hour to evaporate water molecules and promote resin impregnation. At the end of the impregnation, a pressure of 7.1 ± 0.5 bar is applied on the composite laminates for 2 hours at room temperature to ensure the penetration and impregnation of the resin within the plies. The impregnated samples are then cured at 100 °C for 2 hours under pressure to respect the curing cycle of the resin (manufacturer recommendation).

The fibre volume fraction obtained is 55 % and the porosity volume fraction is between 2.5% and 3.5%.

Conditioning and monitoring protocol

Conditioning

After manufacturing, the specimens need to be dried in order to eliminate water molecules. Drying times up to 72 hours at 103 °C were previously investigated and no significant change in the lost mass was found after 24 hours which explains the drying time. Furthermore, all materials are dried at the same temperature (103 °C) and duration of 24 hours in order to have initially the same reference state. After drying, they are weighed and conditioned in an enclosed chamber at constant and regulated relative humidity (RH) according the standard ISO 483. The use of a saturated saline solution in a closed enclosure at constant temperature provides an atmosphere with constant RH. The saturated saline solutions are prepared from mineral salts. Four chambers with constant RH each are set up and kept in a room at 20 ±1 °C. Table 2 represents the aqueous saturated saline solutions used and the RH that each one of them provides. More details about the conditioning protocol are provided in [START_REF] Abida | Effect of hygroscopy on nonimpregnated quasi-unidirectional flax reinforcement behaviour[END_REF]. Some samples were immersed in water at T = 20 °C in order to reach higher water content at saturation. The latter conditioning is supposed to be an atmosphere with a relative humidity of 100 %. All samples were conditioned until water uptake stabilisation.

Gravimetric measurements

Flax yarn specimen weight after drying is around 10 mg depending on their length. Flax / epoxy composite specimens weigh around 10 g, and they are weighed on a 10 -4 g precision balance (Mettler Toledo). Water content in composite specimens is estimated assuming that all water molecules are evaporated and using the Equation (1).

Hygroscopic strain

An optical method which consists of using a profile projector equipped with a camera is referred to in order to measure the apparent diameter of the flax yarns. The measurements are carried out using a Baty Venture 2512CN machine with a precision of 1µm. The thickness of the projected profile is measured at locations previously identified to ensure that measurements are always performed at the same locations before and after conditioning. The cylindrical section of yarns was checked by performing measurements in two orthogonal directions at each point. Therefore, the profile thickness is assimilated to an apparent diameter, which makes it possible to assess the average swelling strain for each test specimen (by computing the average of the swelling ratio of each section). At least 10 measurements of the apparent diameter are made on each yarn.

For composite specimen, hygroscopic strain is measured in the three directions (

z y x ε ε ε , ,
).

Measurements in the out of plane direction ( z ε ) are conducted using a digital micrometre of 1 µm precision. For the other directions ( y x ε ε , ), a digital calliper of 10 µm precision is used for measurements.

Characterisation of mechanical properties

A servo-hydraulic machine INSTRON 8801 equipped with a force cell of 100 kN is used to perform monotonic tensile tests on composite specimens. The strain is measured with an axial extensometer with a gauge length of 50 mm. The test speed is set at 2 mm• min -1 so as to have a strain rate of the order of 10 - 4 s -1 , which is the standard strain rate for quasi-static tests. The warp yarns are aligned with the loading direction. In order to verify the repeatability and reproducibility, at least 3 samples of each configuration of the material are tested. Mechanical tests are carried out on dry samples as well on samples immersed in water or conditioned at various RH until water uptake stabilisation.

Results and discussions

Multiscale three-dimensional sorption

Contrary to 1D sorption where the edges of samples are sealed in order to force water molecules to diffuse in one specific direction (generally in the out of plane direction), in three-dimensional sorption no such preparation is needed. Thus, water molecules can penetrate and diffuse in the material through the 3 directions.

A detailed comparative study of 3D sorption behaviour of composites and its constituents at different scales (flax fabric, elementary yarns that constitute the fabric and epoxy resin) is presented in Table 3.

This study shows that water content within the materials is different from the relative humidity of the conditioning environment. A first study of the relationship between water content at saturation (WCS) and saturation time (ST) and RH of composite specimens shows their non-linear evolution with RH (Figure 2). Regardless of the relative humidity, the saturated water content of elementary flax yarns and flax fabrics is similar (Table 3). Therefore, the cavities between flax yarns within the fabric do not cause extra water absorption of the fabric due to capillarity effects. However, epoxy did not reach saturation after 45 weeks of conditioning. Furthermore, water content in epoxy resin remains low compared to the other studied materials. Thus, we can conclude that WCS of the composite is mainly governed by the flax absorption capacity. A classical law of mixtures based on the constituent moisture behaviour was applied to verify its reliability to obtain the experimental water content of the composite (Equation 2). According to the law of mixtures, the composite saturated water content is the average of the saturated water content of each constituent (flax fabric, epoxy and pores), weighted by its mass fraction. In this case, two hypotheses could be formulated. Either the pores are completely filled with water, which will lead to extra water absorption; or pores cannot be filled with water due to the presence of gas trapped in the voids. Under both assumptions, the computed composite water content at saturation is different from the experimentally measured one (Table 4). Hence, a classical law of mixtures based on the constituent moisture behaviour is not sufficient to explain the moisture absorption of the composite. The analysis should not only include the material behaviour (hygroscopic expansion of component) but also the structural effects (linked to the reinforcement microstructure). The reinforcement water absorption is not representative of the water absorption in the composite, due to extra phenomena in yarns and fabric (capillarity effects).

A comparison of saturation time of flax yarns and composites shows that the elementary yarns saturation time is lower than the composites, which could be attributed to the presence of epoxy that slows the diffusion of moisture in flax yarns.

υ ρ υ υ ρ ⋅ = ⋅ + ⋅ + (2)
8

Where wc is for water content, m υ is the mass fraction, pore v υ is the volume fraction of pores occupied by water molecules, composite ρ is the composite density. The composite water content is assessed when the pores are completely filled with water ( 3.5%

pore v υ =
) and when pores cannot be filled with water ( 0%

pore v υ = ).

Swelling coefficients

Flax yarns

As a reminder, the absorption of water by the flax fibres causes their swelling, which is characterised by the hygroscopic expansion coefficient. This paper focus on estimating the radial swelling coefficient of single flax yarns, noted r β . This parameter is very important to describe the hygro-mechanical coupling of composites as differential swellings of yarn and matrix could create interfacial stresses at the matrix / reinforcement interface and causing, thus, its damage. However, few studies in the literature are interested in the experimental assessment of this coefficient. Le Duigou et al. [START_REF] Duigou | Hygroscopic expansion: A key point to describe natural fibre/polymer matrix interface bond strength[END_REF] measured this coefficient on flax fibres using SEM observations and theoretical water content obtained from the sorption isotherms of Hill [START_REF] Hill | The Water vapor sorption behavior of natural fibers[END_REF]. However, studied composites are reinforced with a fabric constituted of flax yarns. Swelling mechanisms may be different (untwisting of yarns, etc.). Thus, determining the swelling coefficient of flax yarns is essential to understanding the swelling of FFRC composites. The evolution of the radial strain as a function of the water content is plotted in Figure 3. Each point of Figure 3 ). In conclusion, the free radial swelling of flax yarns appears to be similar to that of flax fibres.

Composite

Water absorption in composite material causes a variation in its volume. A study of the relationship between hygroscopic strains and water content according to the 3 directions is carried out. The swelling coefficients ' , L T T and β β β in the longitudinal, transverse and out of plane directions respectively are assessed (Figure 4 and Figure 5).

Even though the evolution of the transversal and longitudinal strains with water content may appear nonlinear, a linear regression is considered in order to assess the swelling coefficients in both directions. The swelling coefficient is the slope of the linear curve of the hygroscopic strain with water content. ). That is to say, moisture absorption of FFRC induces a shrinkage in the direction of the warp yarns. Further works are required to understand the origin of the flax/epoxy composite longitudinal shrinkage due the water absorption. However, this particular behaviour could be explained by the microstructure of the reinforcement. This behaviour could arise from structural effects due to the reinforcement microstructure (yarns or fabrics) or material (elementary flax fibres). On the one hand, water absorption by the flax yarns could lead to a longitudinal shrinkage of elementary yarns. As mentioned in section 2.1 the fabric is constituted of twisted yarns. Water absorption causes radial swelling of yarns. Consequently, the yarn twisting angle could increase due to a spiral spring-like contraction (structural effect). On the other hand, water absorption could lead to elementary flax fibre shrinkage (material behaviour). It could also be explained by the significant radial swelling at high moisture content which could lead to a longitudinal shrinkage of yarns and hence a shrinkage of the composite. Additional measures are needed to confirm these assumptions.

The two transverse swelling behaviours are shown in Figure 5. In Figure 5 ). This difference could be attributed to the presence of the epoxy which has lower hygroscopic coefficient ( 6 . 0 = epoxy β ) [START_REF] Loos | Curing of Epoxy Matrix Composites[END_REF][START_REF] Agbossou | Thermal stresses and thermal expansion coefficients ofn-layered fiberreinforced composites[END_REF][START_REF] Lin | Investigation of the Moisture-Desorption Characteristics of Epoxy Resin[END_REF] which may block the swelling of the yarns. This constraining effect has been modelled by Joffre et al. [START_REF] Joffre | Swelling of cellulose fibres in composite materials: Constraint effects of the surrounding matrix[END_REF]. The swelling of flax yarns in the composite is therefore inferior to that of flax yarns that are free to deform. These differential swellings between yarns matrix are likely to induce interfacial stresses that can weaken the reinforcement / matrix interface in FFRC.

A ratio of 7 is measured between the weft and the out of plane swelling coefficients T β and T ′ β (Figure 5). A swelling coefficient in the weft direction of 019 . 0 13 . 0 ± = T β is assessed (Figure 5(b)). This dimensional stability could be explained by the presence of the few weft yarns that block the swelling in this direction. The swelling of the matrix and the radial swelling of the warp yarns should cause a higher swelling in the transverse direction. However, the assumed longitudinal behaviour of the flax yarns

( 0 < L β
) counterbalances the radial swelling and induces a dimensional stability of the composite in the direction of the weft yarns. This behaviour could also be the cause of internal stresses.

Measurements of the FFRC swelling coefficients, among the three directions, highlight the strong anisotropy of the hygroscopic behaviour of these composites.

Hygro-mechanical properties of composites

In the following, it is more suitable to reason in terms of water content and not in terms of conditioning atmosphere. The reason for this is that, contrary to relative humidity which is an atmospheric property, the water content is an intrinsic property of the material.

Tensile behaviour of composites

Tensile behaviour of flax / epoxy composites with different water content is investigated. As shown in Figure 6, water content dramatically affects the tensile behaviour of FFRC. As a matter of fact, for water

content % 9 . 4 < c w
, tensile curves are characterised by a linear elastic region followed by an apparent yield point. A change of slope involving both a rigidity decrease and the emergence of plasticity (in general, which includes viscoplasticity) is observed. This type of tensile curve was already observed by several authors [START_REF] Oksman | High Quality Flax Fibre Composites Manufactured by the Resin Transfer Moulding Process[END_REF][START_REF] Hughes | Deformation and fracture behaviour of flax fibre reinforced thermosetting polymer matrix composites[END_REF][START_REF] Lebrun | Tensile and impregnation behavior of unidirectional hemp/paper/epoxy and flax/paper/epoxy composites[END_REF][START_REF] Poilâne | Polymer reinforced by flax fibres as a viscoelastoplastic material[END_REF][START_REF] Yan | Flax fibre and its composites -A review[END_REF][START_REF] Rubio-López | Constitutive model to predict the viscoplastic behaviour of natural fibres based composites[END_REF][START_REF] Chilali | Effect of water ageing on the load-unload cyclic behaviour of flax fibre-reinforced thermoplastic and thermosetting composites[END_REF].

As water content increases, the behaviour progressively changes. The apparent yield point is still noticeable on tensile curves at higher water content ( % 9 . 4 ≥ c w

) but an inflection point appears after the yield stress evoking an apparent rigidification of flax / epoxy composites. This behaviour was only observed in some favourable experimental conditions by Richard et al. [START_REF] Richard | A viscoelastoplastic stiffening model for plant fibre unidirectional reinforced composite behaviour under monotonic and cyclic tensile loading[END_REF] at low loading rates and by Scida et al. [START_REF] Scida | Influence of hygrothermal ageing on the damage mechanisms of flax-fibre reinforced epoxy composite[END_REF] at high temperature. The explanation of this phenomenon is not well established and could be attributed to the reorientation of cellulose microfibrils in kink band areas, the spiral spring-like extension of cellulose microfibrils, and the shear-stress-induced crystallisation of the amorphous cellulose of the flax fibres [START_REF] Placet | Nonlinear tensile behaviour of elementary hemp fibres. Part I: Investigation of the possible origins using repeated progressive loading with in situ microscopic observations[END_REF].

Evolution of mechanical properties

In this section, the evolution of the mechanical properties of flax / epoxy composites is discussed.

Primarily, the evolution of Young's modulus is studied. Then, a closer analysis of tangent modulus is

presented. An investigation of the water content influence on the yield stress and on the appearance of the rigidification phenomenon is conducted in the last part of this section.

Young's modulus is computed in the first linear region of the tensile curve. A deformation between 10 %

( y σ × 1 . 0
) and 90% of the yield stress (

y σ × 9 . 0
) is set as the identification range. The yield stress is determined at the first slope change of tensile curves (Figure 6).

The evolution of Young's modulus with respect to water content appears to be non-linear and nonmonotonic (Figure 7). As water content increases, Young's modulus increases up to a maximal value then drops. In fact, Young's modulus is divided by 2 when c w increases from 0 to 10 %. However, the existence of an optimal water content ( 0

2 4 c w . % ≤ ≤
), for which the Young's modulus is maximal, is highlighted.

A more detailed examination of the FFRC tensile curves is then conducted. The analysis of the tangent modulus is presented in Figure 8. Tangent modulus is computed via the derivative of the stress with respect to the deformation. The evolution of tangent modulus with stress confirms the high influence of water content on FFRC's mechanical behaviour as well as the highly non-linear behaviour of flax / epoxy composites. In fact, for low water content ( ), tangent modulus decreases to a stress level (referred to as rigidification stress in the following) then increases.

The tangent modulus increases with stress at high water content highlights a rigidification phenomenon.

Using the tangent modulus evolution curves (Figure 8), a quantitative study of the moisture content effect on yield stress and rigidification stress is carried out. Yield stress decreases with c w in a non-linear way (Figure 9). A ratio of 5 when c w increases from 0 to 10% is measured. Furthermore, for high c w , yield stress is very low, meaning that composite behaviour becomes predominantly plastic, more specifically viscoelastic-viscoplastic [START_REF] Abida | Comportement hygroscopique et couplage hygromécanique dans les composites lin / époxy: approche expérimentale multi-échelle et modélisation[END_REF]. Water molecules could be considered as a plasticising agent for FFRC.

As mentioned above, for high c w ( 4 9 c w . % ≥ ), a rigidification phenomenon is observed for flax / epoxy composites. The evolution of the rigidification stress with water content is studied and given in Table 5. It is worthy of note that the rigidification stress is identified when tangent modulus is minimum. Water content increase lowers the rigidification stress. In fact, at a water content of 4.90 % the rigidification stress is 90 MPa and it drops sharply to 33 MPa when water content is 14.40 %. Rigidification phenomenon is still not well explained but is mostly linked to the rearrangement of microstructure of flax fibres. The rigidification stress drop indicates that the reorientation of cellulose microfibrils is enhanced with water absorption for high strains. The recent work conducted by Requile et al. on hemp / epoxy composites shows that, regardless of the conditioning atmosphere, the tangent modulus drops down to a minimum and then increases [START_REF] Réquilé | Deeper insights into the moisture-induced hygroscopic and mechanical properties of hemp reinforced biocomposites[END_REF]. However, the constant tangent modulus at low water content observed on FFRC is not noted. The influence of water content on behaviour of hemp / epoxy composites is different to the behaviour of flax / epoxy composites. Plant fibre reinforced composite behaviour depends on the kind of fibres used as reinforcement.

Conclusion

An experimental method to assess the radial hygro-expansion coefficients of flax yarns r In fact, flax yarn twisting angle could increase due to a spiral spring-like contraction with water content increase. The latter could also cause a shrinkage of elementary fibres. Another explanation could be the significant radial swelling of flax yarns that could induce a longitudinal shrinkage of flax yarns.

A good dimensional stability in the weft direction is measured ( ), the usual behaviour of FFRC is observed. Tensile curves exhibit a linear elastic region followed by an apparent yield point, which is still noticeable at higher water content ( 4.9 % c w ≥

). An inflection point appears after the yield stress evoking an apparent rigidification of flax / epoxy composites for high c w . This hygro-activated rigidification is related to the flax fibre microstructure rearrangement. The evolution of tangent modulus with respect to stress enhances the existence of the two types of behaviour (with or without rigidification). As water content increases, FFRC behaviour becomes predominantly plastic, more specifically viscoelastic-viscoplastic.

Figures caption

  weight of the dry sample before conditioning and the weight at saturation of sample after exposure respectively.

For

  corresponds to the averages of the radial swelling and the water content obtained on at least 30 yarns, and the bars represent the observed variability. The evolution of flax yarn radial swelling with water content appears to be linear. A high square of the linear correlation coefficient R is found (R 2 =0.993). The linear correlation coefficient is an indicator that allows evaluation of the quality of a simple linear regression. It measures the fit between the model and the experimental data or how well the regression equation is adapted to describe the distribution of points. ,32]. The swelling coefficient of flax yarn is similar to the one measured by Le Duigou et al. [25] on flax fibres ( 14 . 1 = r β

  The longitudinal hygroscopic deformations (in the direction of the warp flax yarns) as a function of the water content of the FFRC are plotted in the Figure4. The FFRC longitudinal swelling coefficient is

(

  (a), the out of plane hygroscopic strain referred to as transversal' ( T β ′ ) is computed. The hygroscopic strain in the weft direction called transversal strain ( T β ) is assessed in Figure 5(b).The swelling coefficient in the out of plane direction is found to be Figure5(a)). This value is lower than that of the radial swelling of flax yarns (

  modulus remains constant to a stress level (yield stress) after which, rigidity decreases until failure. For high water content (

)

  is developed. This work highlights the strong anisotropy of the studied composite. It shows that the radial hygro-expansion coefficient of flax yarns is 039 radial swelling is similar to radial swelling of flax fibres reveals that a composite hygro-expansion coefficient in the out of plane direction ( than the radial swelling coefficient of flax yarns. The differential swelling between the resin and the yarns blocks the swelling along the thickness. Furthermore, measurements show a negative swelling coefficient of FFRC in the warp direction ( longitudinal shrinkage when their water content increases. This behaviour could be attributed to the microstructure reinforcement. It could arise from structural effects due to the microstructure reinforcement (yarns or fabrics) or material (elementary flax fibres behaviour).

  swelling coefficients in the longitudinal and out of plane directions ( L β and T ′ β ) make it possible to explain the coefficient of hygro-expansion in the weft direction. Indeed, the negative longitudinal hygro-expansion coefficient of the weft yarns offsets the swelling of the warp yarns (radial) and the matrix. Tensile behaviour of flax / epoxy composites with different water content is investigated. The study points out the high non-linear behaviour of FFRC. At low water content ( 4.9 % c w <

Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 7 :Figure 8 :Figure 9 :

 2345789 Figure 1: quasi-unidirectional flax fabric (FUD180) Figure 2: evolution of water content at saturation (WCS) and saturation time (ST) with RH of composite samples Figure 3: radial swelling coefficient of flax yarns Figure 4: composite swelling coefficient in the warp direction Figure 5: assessment of the swelling coefficients of composites in the two transverse directions (a) in the out of plane direction (b) in the weft direction Figure 6: stress-strain curves of composites at different water content Figure 7: Young's modulus evolution with water content of flax / epoxy composites for fibre orientation θ=0°

Figure 2 :

 2 Figure 1: quasi-unidirectional flax fabric (FUD180)

Figure 3 :

 3 Figure 3: radial swelling coefficient of flax yarns

Figure 6 :

 6 Figure 6: stress-strain curves of composites at different water content

Figure 7 :

 7 Figure 7: Young's modulus evolution with water content of flax / epoxy composites for fibre orientation θ=0°

Figure 8 Figure 9 :

 89 Figure 8: smoothed tangent modulus

Table 1 :

 1 details on flax fabric used (FUD180)

	Warp Titration (mg/m)	Warp pick count (yarns/m)	Weft Titration (mg/m)	weft pick count (yarns/m)	Balancing coefficient KG	Filling rate (%)
	48	4068	31	304	0.93	94

Table 2 :

 2 relative humidities provided by each aqueous saturated saline solution

Aqueous saturated saline solution Relative humidity at 20°C

	Figures
	Potassium hydroxide (KOH)	8 %
	Potassium Carbonate (K2CO3)	45 %
	Sodium Chloride (NaCl)	75 %
	Potassium Sulphate (K2SO4)	97 %
	Figure 1:

. An optimum weight ratio of 34 : 100 is proposed by the manufacturer. Based on this ratio of the masses of the two constituents, the density of the matrix is around 1.14 g.cm -3 . This resin has been developed and

Funding

This research received no specific grant from any funding agency in the public, commercial, or not-forprofit sectors

Declaration of conflicting interest

The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.


Nicolas Barth

Daniel George
email: george@unistra.fr

Frédéric Bouyer

Aurélien Schwartz

Charles-Henri Lambert

Saïd Ahzi

Yves Rémond

An inverse method predicting thermal fluxes in nuclear waste glass canisters during vitrification and cooling

Keywords: Inverse method, finite element method (FEM), thermal stresses, nuclear waste can, vitrification, cooling

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

p. 2

Introduction

Disposing of ultimate nuclear wastes around the world remains challenging and has been considered with various possible strategies [START_REF] Birac | [END_REF], Wealer et al. (2019)). Lately in France, deep geological storage has received considerable attention, but the disposal of high-level nuclear waste (HLW) depends on available host bedrocks, also considering the type of stored waste products. The reference process currently in use for the management of such long-lived radionuclides is their vitrification as HLW packages (Bonniaud et al. (1980), Goel et al. (2019)) that could be stored under deep underground facility (see Acts 91-1381 of the 30 December 1991 and2006-739 of the 28 June 2006, as well as [START_REF] Andra | Andra research on the geological disposal of highlevel long-lived radioactive waste -results and perspectives[END_REF], [START_REF] Andra | Synthesis: Evaluation of the feasibility of a geological repository in an argillaceous formation -Meuse/Haute-Marne site[END_REF], [START_REF] Andra | Synthesis: Assets of granite formations for deep geological disposal[END_REF], Andra 305C/VA (2009), Wealer et al. (2019)). Other similar vitrification facilities exist around the world (Decamps and Dujacquier (1997), [START_REF] Gin | [END_REF], Gin et al. (2017), [START_REF] Vienna | Current understanding and remaining challenges in modeling long-term degradation of borosilicate nuclear waste glasses[END_REF]) leading to different processes developed to produce various kinds of HLW glasses (Tennant andMurphree (1982), Perez Jr. et al. (2001), Goel et al. (2019)). Different processes induce different thermal treatments that most significantly use different HLW glass block sizes. The resulting nuclear glasses have then specific thermal evolutions over several timescales, due to the accounting for the radionuclide decay rates. All such case-specific operational conditions make the following reported general approach adaptable to many kinds of vitrification, and very interesting to estimate the thermal response of HLW packages.

From the safety and environmental point of view, the manufacturing of these HLW packages induce internal cracking of the nuclear glass block during the vitrification and cooling processes, and their thermal treatments [START_REF] Peters | [END_REF], Faletti and Ethridge (1988), Kahl et al. (1991), Riley Jr. et al. (1999), Connelly et al. (2011)). Hence, numerous studies were carried out on small glass samples such as [START_REF] Perez | [END_REF], Doquet et al. (2013), Dubé et al. (2010), Mallet et al. (2013), Mallet et al. (2015), [START_REF] Ougier-Simonin | [END_REF]). Inside the glass matrix, the glassy phase transitions of the thermal expansion and the viscoelastic behaviors impact the thermal stress build-up for any of the relevant timescales. The behavior of the temperature gradients, thermal history and thermomechanical response of the HLW glass block can then be used to interpret trends of the vitrification process [START_REF] Peters | [END_REF], [START_REF] Barth | [END_REF], Barth et al. (2014)). Finally, the (thermo)mechanical response -using Continuum Damage Mechanics -can be calculated to predict and quantify relevant crack network developments of the glass block (Barth (2013), Barth et al. (2015), Dubé et al. (2010)).

p. 3

When disposed under long-term underground storage conditions, the HLW packages' exposure to water (leaching) also needs to be further investigated for confinement performance, which is particularly related to the internal cracking surface areas of the glass block (Seetharam et al. (2014), Gin (2014), Riley Jr. et al. (1999), [START_REF] Schlegel | [END_REF]). This alteration phenomenon was also studied using glass analogs (Dillmann et al. (2016), Frugier et al. (2008), Michelin et al. (2015), Parruzot (2014), Verney-Carron (2008)). Another challenge lies at the repository scale with HLW thermal power influencing the lithostatic stresses and/or hydrogeological interactions with the repository (Schneefuß et al. (1989), Tennant (1982), Kwon et al. (2013), Cho et al. (2014), Seetharam et al. (2014), Lee et al. (2014), Urpi et al. (2019)). It should be noted that thermally activated nuclear glass crack may heal by mechanical loadings (Doquet et al. (2015)), and that also other complementary effects on materials properties are linked to the radioactivity (Gin (2014), Gin et al. (2017)). At the larger scale, the repository configurations engineered with designed confinement structures, e.g. for a large number of HLW packages, need to be accounted for integrating all these important phenomena. By parts, depending on the multi-physics approach developed, this is abundantly reported in the literature so to evaluate ultimately the safety of underground repositories (Cha et al. (2017), [START_REF] Andra | Phenomenological evolution of a geological repository[END_REF], Hodgkinson (1980), Hodgkinson and Bourke (1980), Kar et al. (2008), Li et al. (2009), [START_REF] Sizgek | Three-dimensional thermal analysis of in-floor type nuclear waste repository for a ceramic waste form[END_REF], Li et al. (2011), [START_REF] Moncouyoux | Conference on the glass to contain high level radioactive wastes: definition, fabrication, storage and disposal [Conference sur le verre pour confiner les déchets de haute activité : Définition, fabrication, entreposage et stockage[END_REF], [START_REF] Avila | [END_REF]).

As a result, a complete dimensional and temporal multi-scale approach of thermal, mechanical and leaching phenomena would help in the prediction of the capacity of any HLW repositories to comply with given environmental standards. A combining effort towards such aim was proposed by the computational approaches of Bouyer and coworkers (Crevoisier et al. (2011), Repina et al. (2019), Repina et al. (2020)). However, some of the aforementioned thermal studies only deal with the HLW as a thermal source term, depending on the radionuclides mix. This is often the case as a point source or linked to the geometry approximation of the package [START_REF] Acar | Comparison of the oncethrough and closed nuclear fuel cycles with regard to waste disposal area required in a geological repository[END_REF], Finsterle et al. (2019)). In such configurations, the HLW packages with thermal power to dissipate due to the presence of radioactive waste are in thermal balance with their direct environment.

However, the thermal history of each HLW package has an effect on the ultimate internal structure of the HLW glass block. The track of this thermal history stems from the cast and cooling process conditions, to the temporary storing conditions, and finally to the ongoing repository characteristics (its design and the neighboring rock formation). Hence, to better predict the p. 4 ultimate internal structure of the HLW glass blocks, better tools and knowledge of the HLW packages' thermal history need to be exploited.

The aforementioned reported studies cover multi-physics approaches contributing by part to address the ultimate metric of the confinement of radionuclides. With the nature of the thermomechanical response of the HLW glass block, the crack network is the result of the thermal history of this glass structure. Mechanical loadings may further damage the thermally-induced crack network, under repository conditions or at the occasion of mishandlings. In this context, a predictive tool to determine the thermal history of each glass block is relevant over timescales corresponding to the cooling process, i.e. to the glass solidification and to the cooling end after tens or even hundreds of thousands of years.

To this aim, we present a numerical approach, which is applicable to any HLW vitrification process.

The purpose of the proposed tool is to more accurately predict the global thermal kinetics of the nuclear glass by imposing thermal fluxes corresponding to the process used over short and long periods of time as compared to only known external temperatures for a given process at the can surfaces. We model the thermal behavior using an inverse thermal method within the glass block during the complete vitrification process. This is carried out by considering the nuclear glass composite as a homogeneous media for the thermal conduction, during its cooling periods (short and long), even encompassing the integration of the cracking phenomena after solidification.

Materials and methods

In order to validate the predictive numerical model, full-scale HLW package vitrification experiments were conducted by casting about 400 kg of borosilicate glass within a container of stainless steel alloy 309S. In this experiment, the HLW package is equipped with 7 thermocouples.

An inactive glass is used to represent the nuclear glass. This glass without radionuclide, "SON68", is a borosilicate glass corresponding to the industrial "R7T7" nuclear glass. The experimental setup and the measured temperatures on the thermocouples as a function of time are presented in Fig.

1 (CEA (2005)). The heat pipe effect of the setup of 7 thermocouples inside the package is neglected, due to their small overall diameter. Standard Inconel® sheathed thermocouples were used, inside thermowells made of stainless steel of 5 mm diameter that are soldered to the HLW package envelope. The sampling rate for the thermocouple measurements is every 10s with a measurement numerical precision of the datalogger to the thousands of °C (precision of the thermocouples is of 0.1% of the measure). Initial data was removed when any thermocouple was p. 5 not relevant to the thermal response, such as thermocouple outside of glass, too far from the steel envelope, or in accidental contact with the centerline trickle for those located at the centerline of the canister. The main objective of the work was to know the overall temperature field of the nuclear glass, as close as possible to what it would be without the thermocouple experimental setup, and not the temperature at each single geometrical point. This was necessary in order to further calculate the solidification process and mechanical behavior of the glass after solidification.

Hence, seven thermocouples were enough to be relevant to the confinement efficacy of the glassy matrix. In addition, inserting many thermocouples within the glass block during the cooling and solidification process would impact greatly the solidification and cracking process that would be disturbed by the presence of the thermocouples. The hereby described vitrification process is the one based on the cold crucible process at Marcoule CEA Center (PEV facility) in France. While this process consists in two successive glass p. 6 castings in the container, as was previously reported for the case of a hot crucible process [START_REF] Barth | [END_REF]), the casting and thermal conditions are significantly different between these two experimental vitrification processes. For practical reasons, as it is not possible to realize an experiment with radioactive glass (R7T7 nuclear glass), in the numerical model, a thermal power source term "r" will be considered to homogeneously heat the glass volume due to radioactivity.

This source term is what drives the long-term thermal behavior, with a long-term delay of the cooling and solidification compared to the SON68 inactive glass.

Purpose of the inverse thermal method

In usual thermal models of HLW packages, Dirichlet-type thermal boundary conditions (BC) can be directly prescribed as external temperatures, using for instance collected data from thermocouples at different heights on the envelope of the HLW package. With appropriate locations to collect the temperature values, these may be enough to obtain the internal thermal state of a glass block for that experiment. However, this leads to an oversimplified model where the thermal fluxes are of the conduction type only. This approach was used in the case of the aforementioned variant vitrification process studies (hot crucible process) [START_REF] Barth | [END_REF], Barth et al. (2014)). Dirichlet-type BCs enable only the study of the experimentally driven scenarios and do not allow to properly investigate process parameters for more general cases.

The proposed inverse method consists in implementing Neumann-type BC, i.e. thermal fluxes, for the modeling of the thermal behavior. The objective is then to compute a simulated temperature field that is coherent with the thermocouple measurements. The approach leads to more complexity, since there is no unique way to define the expressions of the thermal fluxes over the relevant interfaces of the HLW package (such problem is mathematically ill-posed). In the proposed approach, many corrections to the conduction-convection-radiation model of the HLW package are grouped into a high number of fluxes superimposed to the thermal model properties and boundary conditions. This may create heat at some point of the transient analysis to correct the thermal response toward the observed temperature dynamics, which may then only be removed later by a flux of opposite sign. We made this choice to avoid multiplying the thermal sub-models to account for specific effects expected in the HLW package thermal balance as for example the ones provided in the Refs. [START_REF] Gogoi | Inverse analysis of an internal reforming solid oxide fuel cell system using simplex search method[END_REF], Das et al. (2017), Das (2014), Das et al. (2011), Das et al. (2008)). In these works, the thermal inverse analysis modeling approaches are conducted on many kinds of thermal parameters and BC. It can be noted that these approaches include the use of optimization techniques, such as nature-inspired evolutionary p. 7

algorithms [START_REF] Barth | [END_REF]), to assist with exploring the search space and finding good solutions. In the current work, we rely on fixed thermal model parameters and idealized convection and radiation fluxes, all from experiments, using the inverse analysis' Neumann-type BC to optimize the thermal response. It can be noted that the degree of ill-posedness is extremely high for the inverse problem under our pseudo-industrial vitrification conditions, and no regularization technique was implemented to stabilize the solution under noisy measurements, such as the one proposed by (Hon and Wei (2005)). It was applied only to one experiment, with all thermocouples used for identification of the thermal model parameters. The proposed inverse method is carried out manually following by part the transient thermal response and comparing it to the experimental results. This approach is compatible with further optimization scheme, more systematic or with a more complex objective function.

The applications of the proposed thermal modeling approach concerns the thermo-mechanical load of the HLW package, and the ill-posedness and potential divergence to the experimental results at the highest temperature have a rather limited impact since the glass block will likely be in the completely relaxed state, with little effect on the glass transition through the thermal history at high temperature.

Although the proposed inverse method is more complex than prescribing temperature in a deterministic numerical approach, the thermal modeling gains significant extrapolation capabilities for the vitrification process parameters. For example, a thermal forcing (or perturbation) in the thermal energy balance, due to the radioactive power source term, will directly shift the temperature field of the whole glass block with a physically relevant extrapolated behavior. This includes the shifting of the temperatures at the HLW package envelope, which would otherwise remain constant under Dirichlet-type BC. In this Dirichlet-type case, fixed temperatures at the envelope would lead to an overestimation of the thermal dissipation fluxes, then inappropriate for the prediction of the thermal gradients and thermal stresses occurring within the vitrification process.

Description of the proposed inverse method

We apply physically-based heat dissipation fluxes, i.e. natural convection and thermal radiation, whenever possible at the different external boundaries of the HLW package, accounting for both the casting, solidification, and cooling periods. Contrarily to other types of inverse methods, the liquid or solid glass thermal parameters are assumed determined by experiments, they are p. 8

homogeneous in space and functions of temperature. As such, they are not the calibration objectives of the inverse problem. Instead, the objective of the proposed inverse method is to determine additional, virtual, imposed thermal fluxes (if remaining degrees of freedom allows).

Such virtual fluxes are added at the external surface layers of the container, or the glass surface after castings. The aim of the inverse method, over the whole experimental process, is to make the temperatures converge to the measured values at the thermocouple locations.

The net total thermal fluxes at the HLW package external surface (convection and radiation fluxes, with the additional virtual fluxes) is hereafter called "homogenized". It will then encompass different complex thermal phenomena. These include, but not limited to: the heat exchanges between the lower and upper parts of the HLW package (due to external or internal thermal fluxes from the glass that has just been casted); a correction of the empirical relations used for convection as well as the thermal radiation; and less important process design characteristics unaccounted for with the model itself, such as the environing facility's sink temperatures map. The homogenized flux also depends on the internal cracking of the glass block during the thermal experiment, since the conduction is hindered by this phenomenon during the final period of the cooling. However, with the conduction not modeled to account explicitly for such behavior, this effect is directly integrated within the homogenized thermal flux, not the thermal conduction in glass. Altogether, the different thermal fluxes effectively implemented at different locations of the model of the HLW package are further described later in the current section.

Throughout the HLW vitrification and cooling process time, all thermal fluxes are fine-tuned to correctly address the main thermal phenomena we expected. The inverse method (IM) steps of the proposed algorithm are presented in the following resolution diagram, see Flowchart, where IM-i represents the inverse method step i.

At the end of the resolution diagram (IM-7), the output of the inverse method is assumed to accurately model the experiment without overexploitation of the objective experimental values, these being the thermocouple measurements. In the framework of the described optimization technique, overexploitation means the appearance of thermal variations (predicted fluxes) that become at some point less physically correct due to a too strict compliance to the objective measurements. The predictive approach reported here is assumed to produce accurate thermal results within a given range of the defined process parameters.

p. 9

Flowchart: Resolution diagram of the proposed inverse method (IM).

For the experimental vitrification process presented in Fig. 1, the thermal behavior was successfully simulated by the proposed inverse method during this HLW vitrification process with

No Yes

IM-1: We rely on experimental HLW vitrification temperature data (e.g. with inactive glass) using a scale 1:1 HLW package equipped with thermocouples at relevant locations. The collected data are used in IM-3 and IM-4. The thermal process is modeled in IM-2.

IM-2:

We imposed in a 3D FEM model of the HLW package the conduction properties in the glass and the stainless steel envelope from materials' property, with physically-based radiation-convection external BC of the envelope towards the thermal environment. The steel initial temperature is equal to the initial experimental data. Glass initial temperature, top glass surface BC; and fictive can fluxes were the variables of the inverse problem.

IM-3:

The thermal model of IM-2 is used to simulate the vitrification (direct method) using the 3 kinds of variables: i.e. prescribed initial glass temperature; glass free surface sink temperature; and fictive envelope volume fluxes. The latter two are time-dependent.

IM-4 -Test:

Is the direct approach of IM-3 yielding simulated temperatures close to the experimental temperatures? Evaluation scope: at the n th successive transient definition period of the inverse analysis latter two, time-dependent, parameters in IM-3.

Notes: An acceptable error depends on the foreseen application of the thermal analysis.

Moreover, a way to satisfy the test may be to shorten the current n th transient definition period of IM-3 parameters, e.g. as soon as the error hits a minimum in the simulation.

IM-6: Validation of the glass initial temperature

and n definition periods of the inverse method's variables of IM-3. New (n + 1) th definitions will be needed for the rest.

IM-5:

Update and/or adapt the three kinds of variables by going back to the model of IM-2, mainly through the last definition of the latter two kinds. Iterate until IM-6 is reached with all the updated homogenized thermal fluxes.

Note: This can also lead to going backwards before the n th definition sequence of IM-3, in case erroneous heat fluxes slowly emerge through the subsequent process times.

IM-7:

Ultimately, the end of the vitrification experiment is reached (e.g. at ambient temperature), the finely-tuned direct simulation of IM-2 then yield temperatures simulated at any point of the glass block, during the whole experimental time period. This temperature field correctly matches the temperature evolution at the thermocouples locations. Extrapolations are also possible with variant vitrification process parameters. Such is the case for the introduction of a source term accounting for the radioactive decay thermal power, or for any continuity of the simulation after the experiment's end.

p. 10 two castings. We used the following set of thermal BC implemented in the model at IM-2 step, see also Fig. 2: 1-24 layers, modeling the glass block, are filled at each corresponding process time to model the pouring down of liquid glass. During the pouring down of the glass, non-adiabatic thermal exchanges with the upper unfilled layers are assumed, if the layer fills within the same casting period. For the experiment under focus, there are 2 castings of 12 layers each. This modeling is a discretization resolution improvement compared to the pouring/filling technique described in [START_REF] Barth | [END_REF]). The initial temperature in each layer was prescribed at the temperature of the melted glass in the crucible, varied by a modulation using two additional optimization parameters of the inverse method: (a) a vertical temperature linear gradient, and (b) a radial one; as follows: a. The vertical gradient accounts for the cooling effects of the liquid glass trickle with the distance from the crucible, in the container filled with air, as well as the increase of the temperature of the glass free surface with pouring time (assumed a fixed glass temperature at the outlet of the crucible).

b. The radial gradient accounts for cooling effects of the liquid glass free surface, with the mass transfer from the center axis of the package to the cooled wall of the container. In fact, this radial gradient is implemented with two successive linear radial gradients from the centerline to the radius of the canister envelope; see Fig 2a andEq. (1).

The glass initial temperature at the beginning of the inverse analysis is applied as a uniform initial temperature of the glass ( g T0) equal to the glass melt in the crucible. Further in the inverse analysis, the temperature at the top of the glass block and at the centerline (R = 0 m) is kept close to this crucible temperature. To illustrate the prescribed initial temperatures later in the inverse method iterations, the linear gradients ultimately used at the end of the inverse analysis are provided as follows. The glass temperature is in K with the radius R in m, with +1.1 m in the z-axis in m being the top of the glass block, and with z = 0 m at HLW package bottom (steel envelope included):

g R R T z δ × -   = + × - -× × +     (1)
with δ = 0 if R < 0.19 m, or δ = 1 otherwise. This prescribed temperature is applied at the beginning of the filling of the glass casting's corresponding layer body which corresponds to p. 11 the estimated process time (the flow of glass from the crucible being tracked by a scale measure, then idealized linearly in this model). This corresponded to filling times at 0 s for the 1 st glass layer, 49 s for the 2 nd , 90 s for the 3 rd , 133 s for the 4 th , 190 s for the 5 th , 270 s for the 6 th , 390 s for the 7 th , 476 s for the 8 th , 590 s for the 9 th , 690 s for the 10 th , 750 s for the 11 th , and at 819 s for the last layer filling corresponding to the first casting. Then, the 2 nd casting filling of the relevant layers begins at 19150 s for the 13 th layer, the subsequent layers being filled at the following times : 19201; 19222; 19240; 19340; 19454; 19540; 19660; 19740; 19797; 19840; and 19990 s for the last layer up to 1.1 m height of the glass block.

2-The stainless steel container external surfaces are dissipating heat uniformly by natural convection and thermal radiation, with the temperature-dependent fluxes of Eqs. ( 3) and ( 6)

relying on sink temperatures "Tsink" that are the ones of the facility environment:

s Tsink = 313.15 K in our case for the steel, as per the ambient temperature at the end of the experimental vitrification process. The initial temperature of the steel canister is the one measured by the thermocouples closest to the container at the start of the experiment (in K)

which is approximated by the ambient temperature: and by Eq. ( 6). For the sake of simplicity, we identify both sink temperatures of the phenomena at the same value, Tsink, representing the facility thermal environment, or the internal cavity for the glass free surface. The natural convection is described by: ( )

sink S c T T h q - = (3)
with hc the film coefficient for the convection from the surface at the temperature TS, the convection being at ambient air temperature Tsink. For the glass, we define hc in W/K/m 2 as [START_REF] Cereser | Tuyauteries -transmission de la chaleur[END_REF]):

( )

5 0.0113 g c S sink h T T = + × - (4)
and for the steel, with the same units, cylindrical in shape [START_REF] Cereser | Tuyauteries -transmission de la chaleur[END_REF]):

0.25 1.3 0.43 s S sink c T T h   - = ×     (5)
For the sake of simplicity, we also used these same Eqs. ( 4) and ( 5) corresponding to convection in the air even for extrapolations to a radioactive glass casting process, which should have in practice environing neutral gaseous species. The thermal radiation is described by the Stefan-Boltzmann law:

( )

4 4 sink S T T a q - = (6)
with a the radiation constant (emissivity times the Stefan-Boltzmann constant, see Table 1) for the concerned surface towards a uniform sink temperature Tsink, assumed identical to the one for convection.

p. 14

Table 1: Thermomechanical properties used for the HLW package materials. The thermal properties are only temperature-dependent (no relaxation). Unless a formula is provided, the properties are constant to the first/last one specified here, and linearly interpolated in-between.

SON68 glass Stainless steel

Temperature

Phenomenology of the subsequent thermomechanical modeling

The HLW package thermomechanical state relies in this paper on the following multi-physics:

I. Thermal simulation: it can give the thermal history for any point within the HLW package.

II. Vitrification and thermal strain: to account for the time-dependent glassy transition temperature of the coefficient of thermal expansion.

III.

Vitrification and viscoelasticity: the elastic constants (shear modulus) glassy transition is also time-dependent, with possibly different relaxation dynamics compared with point II.

In the previous sections, we only reported the boundary conditions linked to the thermal modeling and simulation of the studied experimental HLW vitrification process (I). For this, we implement the whole approach within a 3D finite element method (FEM), using the ABAQUS code (Abaqus (2013)). All the defined thermal boundary conditions were scripted into ABAQUS FORTRAN 77 p. 15

subroutines named: HETVAL, DFLUX, and FILM. This way, it facilitates the change of values for the virtual thermal fluxes throughout the iterations of the proposed inverse method.

The geometry of the real HLW package is reported by (Moncouyoux et al. (1991), Verney-Carron ( 2008)). We designed the same geometric approximations for the thermomechanical FEM using a refinement of the mesh in the three cylindrical directions (see Fig. 3). The mesh is composed of 50763 total hexahedral elements (C3D8 for the thermal model in ABAQUS) with 50017 nodes, and same number of degrees of freedom for the thermal simulations. The experimental HLW vitrification process reported in Fig. 1 involved the inactive borosilicate glass SON68 onto which a heat source is added to account for the presence of radionuclides, simulating in this experiment the R7T7 nuclear glass. The composition of the SON68 nuclear glass is found in (Frugier et al. (2008)) as a weight percentage of each of the 26 constitutive oxides, or in (Jollivet et al. (2012)) in relation with the molecular mass of each of the constitutive oxides. The SON68 borosilicate glass has 69.1 g/mol, based on the molar mass percentage of the oxides, and an average 21.3 standard atomic mass. The glass is considered as a homogeneous material (with p. 16 homogenized thermomechanical properties) in the current continuum approaches, even if precipitates, cracks and bubbles exist within the experimental glass block.

For the thermal simulation, the Fourier's heat equation is used for conduction in the liquid or glassy SON68, as well as in the stainless steel container (liquid glass is assumed exempt of internal convection):

, , = ∇ • k∇T , , + (7)
where the thermal conduction properties are taken from Table 1 depending on the material location -either glass or steel -at the Cartesian coordinates (x,y,z); r is the internal heat source in the volume V of the body of interest. For the glass, this term is the aforementioned radioactivity source term (thermal power). For the stainless steel, the estimated homogenized fluxes in the inverse method are also led by r, which was appropriate since the stainless steel is not the main thermal resistance for conduction. The initial state for the temperature field over the body of interest was introduced with Eqs. ( 1) and ( 2). The interaction with the surface thermal BC q will be introduced within Eq. ( 8).

For r, the radioactive decay for the simulated R7T7 glass follows rates as function of the radionuclides mix, taken close to one of the scenarios of an internal CEA report by Godon et al. (2004) [unpublished], also available in (Gin et al. (2017)). It can be noted that we used values initially higher than a similar power term introduced by (Advocat et al. (2008)). In our study, for one HLW package, the power of the whole R7T7 glass block of 400 kg yields the following: (i) 2500

W when manufactured, (ii) 644 W at 50 years, (iii) 260 W at 100 years, (iv) 26 W at 1000 years, (v) 0.04 W at 10 4 -10 5 years (and neglected values later).

The heat Eq. ( 7) is applied to the glass block assuming an oversimplified morphology of the glass block, as a homogeneous body. The reality is that the glassy matrix contains precipitates or bubbles which may increase in volume and move during the cooling periods, also with the possibility of cracks nucleation at a more advanced solidified state. The viscosity of the glass is assumed too high for any convection to occur in liquid glass. Studies of semi-transparent materials such as glass also mention the existence of internal thermal radiation (Siegel (1998)), which is not implemented here for the borosilicate glass under focus that is very opaque [START_REF] Barth | [END_REF]).

For this latter phenomenon and the small defects of the glass, we assumed their effects to be already included in the thermal conductivity parameters applied to Eq. ( 7) that were determined from experiments using this complex glass. For the sake of simplicity about the conduction p. 17 parameters of Eq. ( 7), we also neglected that the glassy temperature should be determined after the thermal history. Then, in our case, we fixed the glass transition temperature for the thermal properties (for instance cp varies noticeably between the liquid and glassy states, here at 530 °C).

Moreover, the initial temperature state is also challenging since the glass flow may become anisotropic from the centerline of the axis of the can. Even though the heat Eq. ( 7) may appear too simply implemented, thus neglecting the complex process phenomenology, and in particular glass bubbles/cavitation and cracking. Nevertheless, we assumed that the proposed inverse analysis will significantly correct the thermal response of the glass block with a homogenized thermal conductivity parameter including these different effects.

The FEM is based on the weak form of the heat Eq. ( 7), where the BC from the inverse method are induced through r and q, is then expressed as (Abaqus ( 2013), [START_REF] Barth | [END_REF]):

δ d δ d δ d δ d V V V S U T V T k T V r T V q T S ρ + ∇ ⋅ ∇ = - ∫ ∫ ∫ ∫ & (8)
with the time rate of the specific internal energy, δT the variational field satisfying the essential boundary conditions in this discretization and q the outgoing heat flux from the surface S (outgoing by convention, opposite to the convention (Abaqus (2013))). Further details of the thermal modeling approach can be found in [START_REF] Barth | [END_REF]). It can be noted that the implicit time scheme is used (ABAQUS/Standard), so that the forward simulation time steps are varying and will be much larger than the experimental sampling rate (10 s). Consequently, comparisons with experiments will induce small differences when linearly interpolating simulated values and comparing them to experiments. It is then recommended to compare the values at each simulation time steps only, with a weight inversely proportional to the simulation step.

The thermomechanical modeling of the HLW package is uncoupled from the thermal modeling. It follows the viscous constitutive model described by (Barth et al. (2014)) and reported for different vitrification processes and for a different design where the stainless steel container was not considered in the reported mechanical response. The thermomechanical constitutive relation was implemented in the subroutine called UMAT for the used FEM code. The mechanical BC were kept similar to the ones in the approach by (Barth et al. (2014)), although less constraining in the current work with a designed canister that is fully 3D without any symmetry (for the purpose of applying non-axisymmetric BC in future works).

We report in Table 1 the main required properties for the thermal and thermomechanical modeling of the HLW package materials. The other parameters related to viscoelasticity p. 18

relaxations can be found in our previous work (Barth et al. (2014)). The material properties are extracted from the following references: (Moncouyoux et al. (1991)); Internal CEA note (Langlois (1998) [unpublished]); (Barthelemy (2003)); Internal CEA report (Godon et al. (2004) [unpublished]); [START_REF] Barth | [END_REF], Barth et al. (2014), Dubé et al. (2010), [START_REF] Matteï | Rayonnement thermique des matériaux opaques[END_REF]).

Results and discussion

Inverse method reference simulation

In the current study of the experimental vitrification process, we are matching the simulation temperatures in the proposed inverse analysis to two thermocouples locations corresponding to:

(i) the thermocouples closest to the stainless steel before the glass casting is filling the volume below them; then (ii) when the glass is above the height of the thermocouple, the temperature collected at the middle of the casting heights at a radius R = 16 cm from the centerline. The rest of the container thermal activity is predicted with thermal fluxes from the inverse method using the glass thermal parameters obtained from the literature (see Table 1), and by respecting as much as possible the other known parameters and temperatures of the vitrification process experiment.

The number of iterations required between step IM4 and IM5 of the inverse analysis is very small (typically of the order of 4 or 5) in order to get the temperature field corresponding to the experimental measurements. The obtained thermal results from the inverse method are illustrated for one of the casting in Fig. 4, where the whole casting period is translated into nominal operating times and casting rates.

During the glass pouring phase into the canister, the numerical simulation yields large temperature differences compared to the experimental measurements, which is comprehensible since the pouring process, and the glass free surface evolution, are not directly taken into account in our simulations (each layer of the numerical model is filled instantaneously). This implies several variabilities such as: (i) a glass flow with rupture of the axial symmetry, thus the temperature values may not rise regularly; (ii) locations of thermocouple assumed to be at the interface between the two castings, but not knowing precisely if this was experimentally achieved, and (iii) the glass free surface initial convexity with its central flow, which explains how several thermocouples of the same height, but different radii, are not observed wet by the glass free surface at the same time of the experiment, see Fig. 4. Nevertheless, even with such complex effects not specifically accounted for in the numerical model, they are integrated into the inverse p. 19

analysis by the homogenized thermal fluxes (but may not be studied separately). A rapid convergence of the numerical model towards the measured experimental values is observed. At the pouring phase, the allowed tolerance is relatively loose since the aim of these results is not to model perfectly the thermal history at high temperatures (i.e. for a relaxed liquid state of the glass) during the casting time, but rather to approximate the thermal evolution during the cooling and solidification of the glass. The end of the cooling process is hence presented in Fig. 5. 1b. Here, the process period only covers the simulated time of the second casting.

In the following, the different castings (Casting 1, Casting 2) are denoted by C#1 and C#2 respectively. We also denote by R the radius from the glass (cylindrical) axis to the considered point (See Fig. 1). Some fluctuations of experimental temperature values at the core thermocouple of Casting #2 (see Fig. 5b), are observed and can be due to cavitation effects/bubbles. From the numerical point of view, the temperatures of the thermocouples at the core of the castings glass blocks were difficult to exactly accommodate in the optimization procedure. This could be improved by taking into account a glassy transition for the glass heat conduction properties, and by investigating more appropriately the upper and lower interfaces/surfaces temperatures during the vitrification process. Nevertheless, after the complete filling of the container, the comparison p. 20 of temperatures between thermocouples and simulated ones is very good. We then assume that the thermal kinetics was correctly modeled for the given process with our homogenized fluxes. To evaluate the quality of the inverse analysis process, a comparison can be made between the numerically predicted temperatures and the objective experimental values. This is presented in Fig. 6 by illustrating the difference between the experimental and numerical temperatures of the thermocouple located at R = 16 cm from the container center line. We observe that despite large variations observed at punctual times where the melted glass is poured into the canister (as described above), the correlation between the numerical predictions and experimental temperatures for both the first and second casting is very good. Differences around 20 degrees are observed, i.e. ±3.5 % in absolute values in K, which shows the accuracy of the proposed inverse method.

Figure 6: Matching of the inverse analysis results to the thermocouple measurements -Inverse method quality check at the position of the objective temperatures (located at the middle of the heights of the two castings, for R = 16 cm). Note that the sharp peaks here during the two castings (can fill up), at 0-1 h and 5-6 h, are outside the boundaries of the y-axis (e.g., see Fig. 4 for the second casting). We note the robustness of the convergence algorithm with time.

Process parameters extrapolation

Once the homogenized thermal fluxes were estimated and the inverse thermal method validated with the thermocouples temperatures, it is then possible to simulate the vitrification process of both the inactive SON68 borosilicate glass (the experiment used to identify the fluxes), and the long-term extrapolation to R7T7 nuclear glass through the use of the additional thermal power p. 22 source r (see Fig. 7). The nuclear glass power source is only neglected after 100000 years, not represented here but relevant for the mechanical response of the glass block. In Fig. 7, we represent only the first decades of the solidification and cooling, for which the glass block thermal response is more sensitive to the thermal source term behavior than its natural cooling thermal BC (unchanged over the period). At one arbitrary time of the process, at 16.5 h past the end of the castings (i.e. 2×10 -3 year in Fig. 7), the entire temperature field is also represented in Fig. 8 for this HLW package. For the two simulations, the vitrification process was slightly adapted to normalize the process parameters as for: (i) the standby period between the two castings, and (ii) the final process times, which were out of the scope of the experimental period, so to fit identical thermal environment conditions (at the end of the vitrification process with Tsink = 315 K).

Thermomechanical results

From our previous works (Barth et al. (2014)) a damage mechanics approach coupled to the viscoelastic constitutive modeling was used (see also (Barth (2013), Barth et al. (2015), Barth et al. (2018))). We showed that the cooling rates were of the first importance on the thermal stress build-up inside the glass block, and this remains true within the current thermomechanical study with a more realistic HLW package design. Moreover, for the vitrification process studied in (Barth et al. (2014)), the inverse method was not applied: the experimental temperatures for SON68

were directly applied to the edges of the HLW package -by Dirichlet boundary conditions. This did not allow us at that time to extrapolate the case of R7T7 glass. In the current work, this is in contrary accounted for, and the thermal stress build-up in the glass block was calculated for these two kinds of HLW processes using the inverse thermal method (Neumann boundary conditions).

Table 2 presents the results of the developed thermal stresses following the proposed approach.

The thermal stresses are the ones at the end of the vitrification and cooling process, for the two above cases (the inactive SON68 glass with thermal experimental data and the extrapolated thermal process for R7T7).

p. 24

For prescribed temperature BC related to the hot crucible vitrification process, without the stainless steel container (see [START_REF] Barth | [END_REF])), the viscoelastic stress build-up was reported at 193 MPa (tension) (Barth et al. (2014)). Using the current work's more realistic geometrical approximation (with canister), the hot crucible vitrification process with container yielded maximum stresses of 200 MPa (tension), which shows the good correlation between the two validated cases with and without stainless steel canister, its mechanical role being then minor.

Table 2: Developed stresses within glass with the thermomechanical approach. The results correspond to the maximum stress value along the centerline of the cylindrical glass block, at the end of the cooling process (the glass block is then uniformly at ambient temperature). The current work's geometry refers to Fig. 3, whereas the one denoted "Cylindrical glass block (no steel)" refers to the simplified geometric approximations of (Barth et al. (2014)). In (Barth et al. (2014)), the vitrification process is the hot crucible process with prescribed temperature BC at the envelope, whereas the one presented in the current paper is the cold crucible experimental process with the proposed inverse method. The cold crucible vitrification process was significantly different because the thermal variations imposed to the glass block were different (although the thermal stress build-up follows the same global trends as in the work of (Barth et al. (2014))). For this work model with canister, the maximum stress of inactive SON68 borosilicate glass, at the centerline, was 152 MPa (tension). For p. 25 the extrapolation of the R7T7 nuclear glass, with an initial thermal power of 2500 W/HLW package (and extending the cooling time for the R7T7 nuclear glass to reach the same ambient temperature at the end of the cooling period), the resulting thermal stress build-up was rather similar, with a value of 145 MPa (tension) at the core of the glass block at ambient temperature.

We observe for thermally active R7T7 that, even if the glassy transition is lower and induces more thermal strain compared to the inactive SON68, the stress build-up predictions are then lower than in the case of the experimental vitrification process with inactive glass.

The proposed thermal inverse method enabled not only to determine precisely the developed stresses within the container with good accuracy compared to known experimental data (fixed case), but also to predict these stresses for other cases without an exact knowledge of experimental data, i.e. more interestingly the case of R7T7 nuclear glass.

Conclusion

A thermal inverse method was presented and applied to the thermal fluxes on the envelope and inside a glass block of a high-level nuclear waste package. The inverse method allows us to simulate the thermal behavior and thermal fluxes of the whole glass block, based on an experimental vitrification process using a set-up of thermocouples in the container.

The vitrification experiment was carried out in the experimental facility of the Marcoule CEA center, France, using inactive borosilicate glass within a scale 1:1 container. The interest of the proposed inverse method is that the behavior of an industrial nuclear glass package (radioactive R7T7 glass) can be numerically extrapolated from this experiment to cooling periods on the longterm timescale. This is up to 1 million years to obtain negligible heat source term from radioactive decay, and thus completing the cooling and thermal loading case of the glass block.

We conducted two thermomechanical analyses on the long-term timescale, one on the experimental package of the inactive borosilicate glass, and one on the thermally activated nuclear waste package. The time-dependent viscoelastic stress build-up inside the glass block for these two cases was rather similar, with a smaller stress build-up for the activated nuclear glass (calculated at the same uniform ambient temperature after sufficient time). With such thermal stress build-up analyses, it is possible to predict the driving force behind damage accumulation under long-term storage conditions of the high-level nuclear waste packages.

  Figure 1: (a) General experimental vitrification illustration (CEA (2005)), and thermocouple locations (stars) reported in a schematic of the HLW package for the current experiment (glass block height is 1.1 m). The thermocouples in the experiment were separated by 120° which is not represented by the axisymmetric schematic here. (b) Experimental results (temperatures) for the thermocouples within the 400 kg glass block, without radionuclides, as a function of experimental time.

  , when the pouring down finishes, the glass free surface is associated with the same natural convection and radiation as for the external steel surfaces, Eqs. (3) and (6), except that the g Tsink becomes a parameter to be optimized by the inverse method (each casting glass free surface has one such sink temperature parameter).4-The following additional virtual thermal fluxes are implemented to conduct the rest of the degrees of freedom of the inverse method. On the long term, such fluxes become 0 and the homogenized thermal fluxes at the border of the canister are only those of physically-based convection and radiation. The virtual thermal fluxes were more conveniently applied in the used finite element code as volume heat fluxes of the stainless-steel envelope, instead of the aforementioned surface fluxes. They are each optimized throughout the inverse method as a function of the intended role for each of them, based on the location in the package: a. a uniform axisymmetric flux throughout the height corresponding to each casting b. an axisymmetric flux for the already filled glass block, including the current casting, with capability to have a gradient of this flux value along the height p. 12 c. an axisymmetric flux for the unfilled canister, with the same capability to implement a gradient of the flux value along the height d. a uniform flux for the bottom of the canister We dissociated the contributions 4a to 4c to better mimic the heating of the internal cavity of the canister by convective air flow inside the canister, not modeled explicitly in the proposed approach, only through the inverse analysis. This led for instance to the sharp temperature increase of the thermocouple near the envelope in the middle of casting 2, following the process time of the 1 st glass cast, see Fig. 1b. To conduct the inverse analysis, it is possible to merge all 4a to 4c flux values in a unique axisymmetric volume flux, function of the height and time only. With the technical limitations of the commercial finite element code we use (the ones of its subroutines), only a uniform volume flux was implemented for the stainless steel next to each layer, thus leading to 24 resulting virtual fluxes for 4a to 4c.

  Figure 2: (a) Glass block casting temperature g T0 with the different linear gradients to model the pouring of the liquid glass during the castings for the inverse method (after Eq. (1), and elements along this longitudinal section are shown in the zoom, see also Fig. 3); the represented temperature map has been iterated to match the thermocouple measurements; (b) other thermal fluxes implemented for this vitrification process thermal model.

Figure 3 :

 3 Figure3: FEM design -Full-3D geometric approximations for the thermomechanical FEM approach (about 51000 hexahedral elements). The 24 successive casting layers are represented with different colors, here with different volumes to accommodate further studies where the process parameter of the number of castings is changed (between 1 to 8 castings to fill the entire container). On the right-hand side are illustrated different layers: at the surface, mid-height of a casting, interface between the two castings, and partially at the bottom of the canister. The total height of the glass block is 1.1 m.

Figure 4 :

 4 Figure 4: Temperature convergence in the vicinity of the 2 nd glass casting sensors -Thermal of the proposed inverse method applied to the experimental vitrification process from Fig.

Figure 5 :

 5 Figure 5: Thermal results of the proposed inverse method applied to the vitrification process. The illustrated process period is the same between first and second casting here. (a) Results for the 1 st Casting thermocouples only; (b) results for the rest of the thermocouples in the 2 nd casting, and continuation of Fig. 4 toward the final cooling period of the vitrification process.

Figure 7 :

 7 Figure 7: Long-term global cooling of the two vitrification cases -Extrapolations of the vitrification process time, extended, and of the type of nuclear class, with the thermal contribution r for the radionuclide decay for R7T7. The thermal environment for the whole period, 100 years represented here, is kept the same as the one at the end of the vitrification process (Tsink = 315 K, thermal radiation in vitrification facility and air natural convection).

Figure 8 :

 8 Figure 8: Temperature map 16.5 hours after the glass castings -Extrapolation of the vitrification process for R7T7nuclear glass, see also Fig. 7. Here the time is fixed at 16.5 h past the end of the castings (i.e. 2×10 -3 year). At the end of this extrapolated vitrification process, the only dissipating thermal fluxes are with a sink temperature of Tsink = 315 K, namely the thermal radiation to the vitrification facility and air natural convection from the sides. (Inverse method's virtual fluxes are then negligible.)

  to replace the R7T7 glass g T0 K Initial glass temperature, used for the thermal conduction s T0 313.15 K Initial stainless steel temperature, used for the thermal pTsink inside the HLW canister (convection and radiation only active beyond glass casting) s Tsink 313.15 K Stainless steel envelope's Tsink toward the ambient air and facility from the HLW package sides. U & W/kg Time rate of the specific internal energy V m 3 Volume of the body of interest (x,y,z) 3D position vector in the body of interest

Time since the end of Casting #2 [years, up to 100 years]

		913	R7T7 glass, initially 2500 W/package, glass block core temperature	640	°C
		863	R7T7 glass, initially 2500 W/package, mid-height surface	590	°C
		813	temperature					540	°C
	Temperature [K]	513 563 613 663 713 763	SON68 glass, experiment, glass block core temperature SON68 glass, experiment, mid-height surface temperature	240 290 340 390 440 490	°C °C °C °C °C °C
		463						190	°C
		413						140	°C
		363						90	°C
		313						40	°C
		1E-3	1E-2	1E-1		1E+0	1E+1	1E+2
		1 day	1 week	1 month	months	1 year	10 years	100 years
				p. 23

Acknowledgements

The authors wish to acknowledge the CEA (LCLT) -French Alternative Energies and Atomic Energy Commission -for its financial support as well as providing experimental data (under VESTALE research program of CEA, internal program DEN/2815-C15190). The Andra -French National Radioactive Waste Management Agency -and AREVA NC are also gratefully acknowledged for supporting this study. The funding sources were involved for providing the study designs (not released), the collection of data, and the release of the paper content, but had no involvement in the analysis and interpretation of data, the writing of the paper, or the decision to submit the article for publication.