Introduction

When a material is subject to mechanical loading, heat sources appear inside the material. It has been widely studied for metallic materials machining, forming or laminating processes and fatigue characterization [START_REF] Boulanger | Calorimetric analysis of dissipative and thermoelastic effects associated with the fatigue behavior of steels[END_REF][START_REF] Amiri | Life prediction of metals undergoing fatigue load based on temperature evolution[END_REF]. Heat sources appearance occurs as well in composite materials [START_REF] Goidescu | Damage investigation in cfrp composites using full-field measurement techniques: Combination of digital image stereo-correlation, infrared thermography and x-ray tomography[END_REF][START_REF] Libonati | Damage assessment of composite materials by means of thermographic analyses[END_REF] under mechanical loading [START_REF] Rittel | On the conversion of plastic work to heat during high strain rate deformation of glassy polymers[END_REF], for instance in forming process [START_REF] Wu | Composite forming simulation of a threedimensional representative model with random fiber distribution[END_REF] or in fatigue tests [START_REF] Katunin | Evaluation of criticality of self-heating of polymer composites by estimating the heat dissipation rate[END_REF][START_REF] Islam | A thermographic and energy based approach to define high cycle fatigue strength of flax fiber reinforced thermoset composites[END_REF]. The interest for heat sources lies in their relationship to the thermomechanical properties of the material, to the fatigue life prediction, and to the damage assessment. Heat sources study is also important for investigation of new fabrication processes such as additive manufacturing (AM). Douellou et al. [START_REF] Douellou | Rapid characterization of the fatigue limit of additive-manufactured maraging steels using infrared measurements[END_REF] evaluates the mechanical dissipation in AM materials during fatigue loading. Additive manufacturing can also produce functional or multifunctional materials in which heat sources have a functional role. Neely et al. [START_REF] Neely | Soldered copper lap joints using reactive material architectures as a heat source[END_REF] describes reactive materials architectures that are designed for soldering in inaccessible areas. An energetic material is integrated in the material by additive manufacturing and is used for joining when ignited. The heat source must be prescribed for a use in welding at the appropriate fusion temperature.

Extension to the study of Functionally Graded Materials (FGM) properties might also be interesting as FGM are often designed to be used in severe temperature environment, under high thermal and mechanical loads [START_REF] Tian | Heat conduction investigation of the functionally graded materials plates with variable gradient parameters under exponential heat source load[END_REF].

For the design of functional materials as well as for structural parts, to assess damage or to contribute to the material characterization, the knowledge of the heat sources is then of main importance : the determination of the heat sources is needed to reach a target temperature at the surface in reactive materials joinings; it is also determined as an intrinsic data of the material to understand mechanisms under fatigue [START_REF] Munier | Determination of high cycle fatigue properties of a wide range of steel sheet grades from self-heating measurements[END_REF]. When the thermal effects reach the surface, thermal imaging is often used : the surface temperature gives information for heat sources reconstruction. The study of the thermal response of a component under a heat source is then very useful. We can also cite the application for non destructive damage assessment in laminated composites. Many techniques are used [START_REF] Cao | A novel damage characterization approach for laminated composites in the absence of material and structural information[END_REF] including thermography analysis.

This study is restricted to laminated composites. Laminated composite materials cover a very wide scope of application in industry, for instance in the transportation or energy fields as structural materials. For example in Carbon Fiber Reinforced Polymers [START_REF] Maier | Non-destructive thermography analysis of impact damage on large-scale cfrp automotive parts[END_REF], or in Glass Fiber Reinforced Polymers [START_REF] Colombo | Effect of delamination on the fatigue life of gfrp: A thermographic and numerical study[END_REF], the determination of heat sources is performed from the surface temperature field analysis, the temperature being obtained by IR thermography [START_REF] Movahedi-Rad | Fatigue damage in angle-ply gfrp laminates under tension-tension fatigue[END_REF] and the analysis by the solution of an inverse problem [START_REF] Ozisik | Inverse heat transfer: fundamentals and applications[END_REF]. Different methods have been developed, but in this framework, it appears that many computations have to be performed to identify the characteristics of the heat source in an iterative process. Thus, it calls for an efficient numerical tool to solve many direct problems involving heat diffusion phenomenon for different sets of parameters (position of the heat source, temporal evolution of the heat source ...). For this purpose, reduced order modeling can be advantageously used. In particular, the so-called Proper Generalized Decomposition (PGD) [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF] has shown interesting features to decrease significantly the computational costs. Indeed, it allows us to decrease the dimensionality of numerous physical problems by introducing the separated representations of solutions. It is firstly introduced under the name of ''radial approximation'' [START_REF] Ladevèze | Sur une famille d'algorithmes en mécanique des structures[END_REF] by P. Ladevèze within the framework of the LArge Time INcrement (LATIN) method. Concerning the thermal analysis, the PGD method has already been carried out for different issues: the moving thermal load on the surface is considered in [START_REF] Pruliere | On the solution of the heat equation in very thin tapes[END_REF] with a space/time separation. The non-linear transient thermal problem with a concentrated moving laser source is addressed in [START_REF] Favoretto | Reduced order modeling via pgd for highly transient thermal evolutions in additive manufacturing[END_REF]. It is applied in [START_REF] Ghnatios | Proper generalized decomposition based dynamic data-driven control of thermal processes[END_REF] for Dynamic Data-Driven Application Systems taking into account the prescribed temperature as parameters. The reciprocity principle for real-time monitoring of thermal processes (transient analysis) is developed in [START_REF] Aguado | Real-time monitoring of thermal processes by reduced-order modeling[END_REF]. The PGD is also used in the solution of an inverse heat conduction problem within the Bayesian framework [START_REF] Berger | Proper generalized decomposition model reduction in the bayesian framework for solving inverse heat transfer problems[END_REF]. The variables separation of the spatial coordinates x and z has already been used in [START_REF] Vidal | Thermo-mechanical analysis of laminated composite and sandwich beams based on a variables separation[END_REF] for laminates. For a review about the PGD and its fields of applications in our framework, the reader can refer to [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF][START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF][START_REF] Chinesta | First steps towards an advanced simulation of composites manufacturing by automated tape placement[END_REF][START_REF] Chinesta | 3 proper generalized decomposition[END_REF].

Note that an alternative way to decrease the computational cost consists in using the Proper orthogonal decomposition (see [START_REF] Białecki | Proper orthogonal decomposition and modal analysis for acceleration of transient fem thermal analysis[END_REF] for transient thermal problem).

This work aims at modeling composite beam structures regardless of the volume heat source location for unsteady thermal analysis with a low computational cost. For this purpose, the present approach is based on the separated representation where the temperature is written as a sum of products of unidimensional polynomials of 𝑥 coordinate, 𝑧 coordinate, the time 𝑡 and also the volume thermal heat position 𝑥 𝑠 . A piecewise fourth-order Lagrange polynomial of 𝑧 is chosen as it is particularly suitable to model composite structures. As far as the variation with respect to the axial coordinate is concerned, a 1D three-node beam Finite Element (FE) is employed. The functions of the volume heat source position are also piecewise quadratic. Finally, the deduced non-linear problem is solved using a classical fixed point method. Thus, four 1D linear problems are solved alternatively, in which the number of unknowns is much smaller than in a Layerwise approach. Once different sets of these four 1D functions are determined, the temperature fields can be deduced for any size and location of the heat sources by applying the superposition principle. Nevertheless, in our particular framework, the location of the heat source in the thickness direction z or its temporal evolution can change. Thus, a new strategy is developed to take into account these test cases. It relies on a socalled multi-resolution approach and comes from the preliminary stage of the Latin Method [START_REF] Bussy | The large time incremement method for numerical simulation of metal forming processes[END_REF][START_REF] Boisse | A new approach in non-linear mechanics : the large time increment method[END_REF]. Once the 1D functions are computed for a given configuration, they are used for subsequent computations with other heat source configurations. The new computations are limited to the corrections of some of the previously built 1D functions instead of performing the whole analysis. In this way, the computational cost will be reduced.

The outline of the article is the following. First, the unsteady heat conduction problem is described. Then, the formulation of the parametrized problem based on the variables separation is given. The particular assumption on the temperature field yields a non-linear problem which is solved by an iterative process. The FE discretization is also given. Then, different algorithms are shown to explain the multiresolution strategies. Finally, numerical results are provided to show the possibilities of the method and assess the accuracy of the involved results. Different heat source and composite structures are considered.

Unsteady heat conduction problem description

A laminated beam with 𝑁𝐶 layers is defined in a domain

 =  𝑥 ×  𝑦 ×  𝑧 = [0, 𝐿] × [-𝑏 2 ≤ 𝑦 ≤ 𝑏 2 ] × [-ℎ 2 ≤ 𝑧 ≤ ℎ 2 ]
expressed in a Cartesian coordinate (𝑥, 𝑦, 𝑧). The cross-section of the beam is rectangular. ℎ and 𝑏 are the height and the width, respectively. The central line of the beam is chosen as the 𝑥 axis. It is shown in Fig. 1.

Considering the geometry and the boundary conditions, a 2D problem will be considered hereafter.

Heat conduction problem

Thermal constitutive equation

The physical problem considered here involves the linear heat conduction equations. The constitutive equation is given by the Fourier law:

𝐐 = -𝝀 𝐠𝐫𝐚𝐝(𝜃) (1)
where 𝜃, 𝐐, 𝝀 and 𝐠𝐫𝐚𝐝 are the temperature, the heat flux, the thermal conductivity and the gradient operator, respectively. For the layer (𝑘), we have

𝝀 (𝑘) = [ 𝜆 (𝑘) 1 0 0 𝜆 (𝑘) 3 ]

Heat transfer equation

Considering a prescribed volume heat source 𝑟 𝑑 at the location 𝑀 𝑠 , the heat transfer equation can be written as

𝜌𝑐 𝑝 𝜕𝜃 𝜕𝑡 + 𝑑𝑖𝑣(𝐐) = 𝑟 𝑑 𝛿(𝑀 -𝑀 𝑠 ) (2)
where 𝜌 is the mass density, 𝑐 𝑝 the heat capacity per unit mass, 𝑡 the time, 𝑡 ∈ [0, 𝑡 𝑚𝑎𝑥 ]. 𝛿 is the Dirac delta function.

The weak form of the boundary value problem

In the presented problem, the prescribed surface heat surface and the convective heat transfer are not considered. For 𝛿𝜃 ∈ 𝛿𝛩 (𝛿𝛩 = {𝜃 ∈ 𝐻 1 ( × )∕𝜃 = 0 on 𝜕 𝜃 } with  = [0, 𝑡 𝑚𝑎𝑥 ]), the variational principle is given by: find 𝜃 ∈ 𝛩, with 𝜃(𝑡 = 0) = 𝜃 𝑖𝑛𝑖 such that:

∫  ∫  𝛿𝜃 𝜌𝑐 𝑝 𝜕𝜃 𝜕𝑡 𝑑𝑑𝑡 -∫  ∫  𝐠𝐫𝐚𝐝(𝛿𝜃) 𝑇 𝐐(𝜃)𝑑𝑑𝑡 = ∫  ∫  𝛿𝜃 𝑟 𝑑 𝛿(𝑀 -𝑀 𝑠 ) 𝑑𝑑𝑡 ∀𝛿𝜃 ∈ 𝛿𝛩 (3)
𝛩 is the space of admissible temperatures, i.e. 𝛩 = {𝜃 ∈ 𝐻 1 ( × )∕𝜃 = 𝜃 𝑑 on 𝜕 𝜃 }.

Application of the variable separation to the unsteady thermal analysis for a parametrized solution 𝒙 𝒔

The Proper Generalized Decomposition (PGD) was introduced in [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF] and is based on an a priori construction of separated variable representation of the solution. For the thermal analysis, the solution with separated coordinate variables has been developed in [START_REF] Pruliere | On the solution of the heat equation in very thin tapes[END_REF].

The following sections are dedicated to the introduction of this approach to build a parametric solution for unsteady analysis of a composite beam with heat sources at any locations on the 𝑥-axis. The position and the size in the 𝑧-direction and the temporal evolution of the heat source are fixed. The latter will be addressed in Section 4 considering a multi-resolution approach.

The problem defined by Eq. (1)-Eq. ( 3) is considered as a parametrized problem where the heat source location 𝑥 𝑠 is defined in a bounded interval [0, 𝐿]. The thermal solution for a point 𝑀 of coordinate (𝑥, 𝑧) of the structure, at the time 𝑡, depends on the values of this parameter and is denoted 𝜃(𝑥, 𝑧, 𝑡, 𝑥 𝑠 ) (see Fig. 2).

The parametrized problem

The temperature solution 𝜃(𝑥, 𝑧, 𝑡, 𝑥 𝑠 ) is constructed as the sum of 𝑁 products of functions of only two spatial coordinates, time and heat source location (𝑁 ∈ N is the order of the representation) such that:

𝜃(𝑥, 𝑧, 𝑡, 𝑥 𝑠 ) = 𝑁 ∑ 𝑖=1 𝑇 𝑖 𝑥 (𝑥)𝑇 𝑖 𝑧 (𝑧)𝑇 𝑖 𝑡 (𝑡)𝑇 𝑖 𝑥 𝑠 (𝑥 𝑠 ) (4)
where 𝑇 𝑖 𝑥 (𝑥), 𝑇 𝑖 𝑧 (𝑧), 𝑇 𝑖 𝑡 (𝑡) and 𝑇 𝑖 𝑥 𝑠 (𝑥 𝑠 ) are defined in  𝑥 ,  𝑧 , [0, 𝑡 𝑚𝑎𝑥 ] and  𝑥 𝑠 , respectively. A classical quadratic FE approximation is used in  𝑥 and  𝑥 𝑠 and a fourth-order lagrangian expansion is used in  𝑧 .

The expression of the gradient of the temperature is reduced to:

𝐠𝐫𝐚𝐝(𝜃) = 𝑁 ∑ 𝑖=1 [ 𝑇 𝑖 𝑥,𝑥 (𝑥)𝑇 𝑖 𝑧 (𝑧)𝑇 𝑖 𝑡 (𝑡)𝑇 𝑖 𝑥 𝑠 (𝑥 𝑠 ) 𝑇 𝑖 𝑥 (𝑥)𝑇 𝑖 𝑧,𝑧 (𝑧)𝑇 𝑖 𝑡 (𝑡)𝑇 𝑖 𝑥 𝑠 (𝑥 𝑠 ) ] (5)
For a structure submitted to a localized volume heat source at 𝑥 = 𝑥 𝑠 , the parametrized problem can be formulated as:

find

𝜃 ∈ 𝛩 𝑒𝑥𝑡 with 𝜃(𝑡 = 0) = 𝜃 𝑖𝑛𝑖 (𝛩 𝑒𝑥𝑡 = {𝜃 ∈ 𝐻 1 ([0, 𝑡 𝑚𝑎𝑥 ] ×  𝑥 ×  𝑧 ×  𝑥 𝑠 )∕𝜃 = 𝜃 𝑑 on 𝜕 𝜃 }) such that: 𝑏( 𝜕𝜃 𝜕𝑡 , 𝛿𝜃) -𝑎(𝜃, 𝛿𝜃) -𝑑 𝑟 (𝛿𝜃) = 0, ∀𝛿𝜃 ∈ 𝛿𝛩 𝑒𝑥𝑡 (6)
where the bilinear forms 𝑎, 𝑏 and the linear forms 𝑑 𝑟 are defined by

𝑎(𝜃, 𝛿𝜃) = ∫ 𝑡 𝑚𝑎𝑥 0 ∫  𝑥 ∫  𝑧 ∫  𝑥 𝑠 𝐠𝐫𝐚𝐝(𝛿𝜃) 𝑇 𝐐(𝜃) 𝑑 𝑥 𝑠 𝑑𝑑𝑡 𝑏( 𝜕𝜃 𝜕𝑡 , 𝛿𝜃) = ∫ 𝑡 𝑚𝑎𝑥 0 ∫  𝑥 ∫  𝑧 ∫  𝑥 𝑠 𝛿𝜃 𝜌𝑐 𝑝 𝜕𝜃 𝜕𝑡 𝑑 𝑥 𝑠 𝑑𝑑𝑡 𝑑 𝑟 (𝛿𝜃) = ∫ 𝑡 𝑚𝑎𝑥 0 ∫  𝑥 ∫  𝑧 ∫  𝑥 𝑠 𝛿𝜃 𝑟 𝑑 𝛿(𝑀 -𝑀 𝑠 ) 𝑑 𝑥 𝑠 𝑑𝑑𝑡 (7)

The problem to be solved

The resolution of Eq. ( 6) is based on a greedy algorithm. If we assume that the first 𝑚 functions have been already computed, the trial function for the iteration 𝑚 + 1 is written as

𝜃 𝑚+1 (𝑥, 𝑧, 𝑡, 𝑥 𝑠 ) = 𝜃 𝑚 (𝑥, 𝑧, 𝑡, 𝑥 𝑠 ) + 𝑇 𝑥 (𝑥)𝑇 𝑧 (𝑧)𝑇 𝑡 (𝑡)𝑇 𝑥 𝑠 (𝑥 𝑠 ) (8)
where 𝑇 𝑥 , 𝑇 𝑧 , 𝑇 𝑡 and 𝑇 𝑥 𝑠 are the functions to be computed and 𝜃 𝑚 is the associated known sets at iteration 𝑚 defined by

𝜃 𝑚 (𝑥, 𝑧, 𝑡, 𝑥 𝑠 ) = 𝑚 ∑ 𝑖=1 𝑇 𝑖 𝑥 (𝑥)𝑇 𝑖 𝑧 (𝑧)𝑇 𝑖 𝑡 (𝑡)𝑇 𝑖 𝑥 𝑠 (𝑥 𝑠 ) (9)
The test function is [START_REF] Neely | Soldered copper lap joints using reactive material architectures as a heat source[END_REF] The test function defined by Eq. ( 10), the trial function defined by Eq. ( 8) are introduced into the weak form Eq. ( 6) to obtain the four following equations: [START_REF] Maier | Non-destructive thermography analysis of impact damage on large-scale cfrp automotive parts[END_REF] From Eq. [START_REF] Tian | Heat conduction investigation of the functionally graded materials plates with variable gradient parameters under exponential heat source load[END_REF] to Eq. ( 14), a coupled non-linear problem is derived. Thus, a non-linear resolution strategy has to be used. A classical fixed point method is carried out as shown in Algorithm 1.

Algorithm 1 for 𝑚 = 0 to 𝑁 𝑚𝑎𝑥 do Initialize T (0) 𝑥 , T (0) 𝑡 , T (0) 𝑥 𝑠 for 𝑘 =
= T (𝑘) 𝑥 , 𝑇 𝑚+1 𝑧 = T (𝑘) 𝑧 , 𝑇 𝑚+1 𝑡 = T (𝑘) 𝑡 , 𝑇 𝑚+1 𝑥 𝑠 = T (𝑘) 𝑥 𝑠 , Set 𝜃 𝑚+1 = 𝜃 𝑚 + 𝑇 𝑚+1 𝑥 𝑇 𝑚+1 𝑧 𝑇 𝑚+1 𝑡 𝑇 𝑚+1 𝑥 𝑠
Check for convergence end for

Discretization of the problem

A discretization of Eq. [START_REF] Tian | Heat conduction investigation of the functionally graded materials plates with variable gradient parameters under exponential heat source load[END_REF] to Eq. ( 14 [START_REF] Colombo | Effect of delamination on the fatigue life of gfrp: A thermographic and numerical study[END_REF] where

 𝑒 𝑥 𝑇 = [ 𝑇 𝑥,𝑥 𝑇 𝑥 ] ,  𝑒 𝑧 𝑇 = [ 𝑇 𝑧 𝑇 𝑧,𝑧 ] . The matrices [𝐍 𝑥 ], [𝐁 𝑥 ], [𝐍 𝑧 ], [𝐁 𝑧 ], [𝐍 𝑥 𝑠 ]
contain the interpolation functions, their derivatives and the jacobian components dependent on the chosen discrete representation.

Finite element problem to be solved on  𝑧

The functions 𝑇 (𝑘-1)

𝑥 , 𝑇 (𝑘-1) 𝑡 , 𝑇 (𝑘-1)
𝑥 𝑠 are assumed to be known. They will be denoted T𝑥 , T𝑡 , T𝑥 𝑠 , respectively and the function 𝑇 (𝑘) 𝑧 to be computed will be denoted 𝑇 𝑧 .

The variational problem defined on  𝑧 from Eq. ( 12) can be written under the following form: 15) in the variational Eq. ( 16) leads to the linear system

∫  𝑧 𝛿 𝑇 𝑧 𝝀 𝑥𝑡𝑥 𝑠 ( T𝑥 ,
𝜦 𝑥𝑡𝑥 𝑠 ( T𝑥 , T𝑡 , T𝑥 𝑠 ) 𝐪 𝑧 =  𝑥𝑡𝑥 𝑠 ( T𝑥 , T𝑡 , T𝑥 𝑠 , 𝜃 𝑚 ) (17)

T , T , T T , T , T

where 𝐪 𝑧 is the vector of the nodal temperatures associated with the finite element mesh in  𝑧 . As in the previous section, 𝜦 𝑥𝑡𝑥 Again, the functions 𝑇 (𝑘) 𝑧 , 𝑇 (𝑘-1)

𝑥

, 𝑇 (𝑘-1)

𝑥 𝑠
which are assumed to be known, will be denoted T𝑧 , T𝑥 , T𝑥 𝑠 , respectively and the function 𝑇 (𝑘) 𝑡 to be computed will be denoted 𝑇 𝑡 .

From Eq. ( 13), the following ODE can be deduced to compute 𝑇 𝑡 :

𝐴( T𝑥 ,
This equation can be solved using a fourth-order Runge-Kutta method. 𝐴, 𝐵, 𝐶 given in Appendix A.1.

Finite element problem to be solved on  𝑥

For the sake of simplicity, the functions 𝑇 (𝑘) 𝑧 , 𝑇 (𝑘) 𝑡 , 𝑇 (𝑘-1)

𝑥 𝑠

which are assumed to be known, will be denoted T𝑧 , T𝑡 , T𝑥 𝑠 , respectively and the function 𝑇 (𝑘) 𝑥 to be computed will be denoted 𝑇 𝑥 . The variational problem defined on  𝑥 from Eq. ( 11) can be written under the following form:

∫  𝑥 𝛿 𝑇

Explicit solution on  𝑥 𝑠

For this problem, the functions 𝑇 (𝑘) 𝑥 , 𝑇 (𝑘) 𝑧 , 𝑇 (𝑘) 𝑡 are assumed to be known and are denoted T𝑥 , T𝑧 and T𝑡 respectively. The function 𝑇 (𝑘) 𝑥 𝑠 to be computed will be denoted 𝑇 𝑥 𝑠 . The problem Eq. ( 14) can be written such that the unknown function 𝑇 𝑥 𝑠 could be computed in a straightforward manner as:

𝑇 𝑥 𝑠 (𝑥 𝑠 ) = 𝐹 𝑠 ( T𝑥 , T𝑧 , T𝑡 ) T𝑥 (𝑥 𝑠 ) - 𝑚 ∑ 𝑖=1 𝐹 𝑖 𝑚 ( T𝑥 , T𝑧 , T𝑡 , 𝜃 𝑚 )𝑇 𝑖 𝑥 𝑠 (𝑥 𝑠 ) 𝐺 𝑠 ( T𝑥 , T𝑧 , T𝑡 ) (25)
where

𝐺 𝑠 ( T𝑥 , T𝑧 , T𝑡 ) = ∫  𝑧 𝜌𝑐 𝑝 T 2 𝑧 𝑑𝑧 ∫ 𝑡 𝑚𝑎𝑥 0 𝜕 T𝑡 𝑑𝑡 T𝑡 𝑑𝑡 ∫  𝑥 T 2 𝑥 𝑑𝑥 + ∫ 𝑡 𝑚𝑎𝑥 0 T 2 𝑡 𝑑𝑡 ∫  𝑥  𝑥 ( T𝑥 ) 𝑇 [ ∫  𝑧 𝜮 𝐳 ( T𝑧 ) 𝑇 𝝀

Multi-resolution strategy

The aim of this section consists in describing different strategies to take into account the modification of parameters associated to the heat source or the geometry of the composite beam (number of layers, location of the heat source and time dependency of the loading). The main issue is to decrease the computational cost of numerous calculations for different configurations. For instance, it will be possible to change the location in the thickness direction and the temporal evolution of the heat source and the number of layers of the beam. For this purpose, a multi-resolution algorithm is given, using the advantages of the present variable separation method. The approach described here can be considered as an extension of the so-called preliminary stage introduced in [START_REF] Bussy | The large time incremement method for numerical simulation of metal forming processes[END_REF][START_REF] Boisse | A new approach in non-linear mechanics : the large time increment method[END_REF] to improve the convergence rate of the LATIN method. Once some couples (at least 2) are built, the process consists in updating: (i) only the 1D time functions 𝑇 𝑖 𝑡 ; (ii) both all the 1D zfunctions 𝑇 𝑗 𝑧 then the 1D time functions 𝑇 𝑖 𝑡 . The other basis functions are assumed to be known. Thus, assuming that 𝑁 1 couples are built, these problems to be solved can be written as The associated multi-resolution algorithms are given in Algorithm 2 and 3. The first step consists in building 𝑁 1 4-uplets for an initial set of parameters. Then, only 𝑇 𝑝 𝑧 or 𝑇 𝑝 𝑡 are updated, and the other functions are reused in the following calculations involving new sets of parameters. The associated problem is 1D. If the residual error remains high, only few new couples could be computed to achieve the convergence in the so-called enrichment stage.

Numerical results

In this section, the following FE discretization is chosen: (i) a classical 3-node quadratic FE for the thermal unknowns depending on the 𝑥-axis coordinate and the 𝑥 𝑠 location, (ii) a fourth-order layer-wise description for the transverse direction.

Unsteady diffusion tests are presented evaluating the efficiency and the properties of the algorithm and also validating our approach. To this purpose, a laminated composite beam is addressed. A preliminary test is first considered to illustrate the computation of the temperature with the variable separation. The explicit solution with respect to the location of the heat source allows us to build the temperature field for any size using the superposition principle (see Appendix B). The number of layers, the location of the heat source and the time dependency of the loading can change. The approach is assessed by comparing with 2D FEM solutions using a commercial code. A multiresolution approach is performed to take into account the modification of such parameters.

The main characteristics of the test cases are summarized as follows: 3). We have 𝜙 0 = 6000 W m -2 .

material properties: reference values: results are computed using the Ansys software. In the subsequent numerical results, the distribution of the temperature/heat flux along the thickness or the beam axis is such that the results path goes through the center of the heat source. [START_REF] Białecki | Proper orthogonal decomposition and modal analysis for acceleration of transient fem thermal analysis[END_REF] for a given value of 𝑧 = z

𝜆 1 1 = 𝜆 3 11 = 5 W m -1 K -1 ; 𝜆 1 3 = 𝜆 3 3 = 0.5 W m -1 K -1 𝜆 2 1 = 0.5 W m -1 K -1 ; 𝜆 2 3 = 5 W m -1 K -1 𝜌 = 10 kg m -3 ; 𝑐 𝑝 = 1 J kg -1 K -1
√ √ √ √ ∫  𝑥 ∫ 𝑡 𝑚𝑎𝑥 0 (𝑇 (𝑥, z, 𝑡) -𝑇 𝑟𝑒𝑓 (𝑥, z, 𝑡)) 2 𝑑𝑡𝑑 𝑥 ∫  𝑥 ∫ 𝑡 𝑚𝑎𝑥 0 𝑇 𝑟𝑒𝑓 (𝑥, z, 𝑡) 2 𝑑𝑡𝑑 𝑥

Preliminary test case

A first test case is carried out to assess the accuracy and the behavior of the method. A three-layer composite beam is submitted to a heat source in the layer 2 or 3. A ramped loading is considered as shown in Fig. 3(b). A nil temperature is imposed on two edges. Different sizes of heat source (in the axial direction) are addressed, but only the most severe case (𝐼 3 ) is presented below.

Firstly, a convergence study is carried out. Once the 4-uplets are built for a heat source in the layer 2 or 3 at any locations, the distribution of the temperature over the whole structure and for the whole time domain can be deduced for a given size of the heat source. A small domain 𝐼 3 is taken for the presented results. For the two types of error indicators (Eq. ( 29) and Eq. ( 30)), the convergence of the PGD method is rather monotonous (see Fig. 4).

In the following, 60 functions are built to compute the solution.

The distributions of temperature along the 𝑥-axis and the thickness are given in Fig. 5 and Fig. 6, respectively. It can be inferred from these figures that the accuracy of the results is very good despite the localized loading. The steep temperature gradient is well-captured (heat source in the middle layer, Fig. 5(a)). For further assessment, the distribution of the heat flux along the thickness is shown in Fig. 7. The results are very satisfactory. The use of a higher-order expansion through the thickness is justified as the variation of the heat flux along the 𝑧-axis is not linear. It can be also noticed that the distributions are rather different depending on the position of the heat source. The distribution of the temperature over the whole beam given in Fig. 9 and Fig. 10 allows also us to illustrate that. Finally, the variation of the temperature versus time is plotted in Fig. 8. It shows the capability of the present method to perform unsteady diffusion analysis.

Multi-resolution approach

In this section, we take advantages of both the explicit expression of the temperature with respect to the heat source location and also the variables separation to deduce an approach well-suited to perform numerous computations involving different configurations. Hereafter, different parameters can be changed :

• Modification of the temporal variation of the heat source (value of 𝑡 𝑟 ) • Modification of the position of the heat source (through the thickness) • Modification of the number of layers • Simultaneous modifications of the number of layers/heat source location/time dependence of the heat source First, an initial computation is carried out to build the functions (𝑇 𝑖 𝑥 , 𝑇 𝑖 𝑧 , 𝑇 𝑖 𝑥𝑠 , 𝑇 𝑖 𝑡 ) associated to a given configuration. It is described in each test case. For the present study, 60 functions are built. Then, these functions are used in the subsequent computations with a new set of parameters as described below. Thus, the new computations have a low computational cost as it is reduced to update only 1D functions, as 𝑇 𝑖 𝑧 and 𝑇 𝑖 𝑡 for our applications. The process is given in Algorithm 2 and 3. The results are illustrated in the following.

Modification of 𝑡 𝑟

In this example, the initial basis is built for 𝑡 𝑟 = 0.1 s and a nil temperature is imposed on all edges. The heat source is localized in layer 3. Due to the type of modification, only the time functions 𝑇 𝑡 are updated for the new value of 𝑡 𝑟 = 0.15 s. The distribution of temperature with respect to the time and along the thickness, and the transverse heat flux through the thickness are shown on Fig. 11. The agreement with the reference solution is very good.

Modification of the heat source location

The same test case as in Section 5.2.1 is considered. The initial basis involves a heat source in layer 2. Then, a new analysis is carried out with a heat source in layer 3. First, the 𝑇 𝑖 𝑧 functions are updated. From Fig. 12 and Fig. 13, it can be noticed that it is insufficient to obtain accurate results (dotted line). The value of the error indicator 𝐸𝑟𝑟𝑜𝑟 𝑥𝑧 is about 28%. Then, it is also needed to update the 𝑇 𝑖 𝑡 functions to deduce distribution of temperature and transverse heat flux in a very good agreement with the reference solution (full line). The error rate decreases to about 4%.

Modification of the number of layers

In this section, a modification of the number of layers is considered. Thus, the localization and the height of the heat source change. The initial basis is related to a three-layer beam with a heat source in layer 2. A nil temperature is prescribed on two edges. The new computation involves 6 layers with a heat source in layer 4. This new configuration does not induce any difficulties as only the z-problem and the integration over the thickness are involved owing to the variables separation. As previously, the z-functions are naturally updated, but a significant error rate remains (about 37%). Thus, the time functions are also updated. From Fig. 14 and Fig. 15, we can see that the throughthickness and axial distribution of the temperature is very close to the Ansys reference solution. The variation of the temperature with respect to time is also very accurate. Despite the significant difference on the distributions on thermal quantities between the initial configuration and the current one, the present method allows us to obtain very satisfactory results.

Modification of the number of layers/heat source location/time dependence of the heat source

To illustrate the wide range of validity of the multi-resolution approach, both geometry and heat source characteristics are changed. The initial configuration is described as : (i) geometry : 3 layers, (ii) heat source : located in layer 2 with a stepped time variation as in Fig. 3(a). The new configuration involves : (i) geometry : 6 layers, (ii) heat source : located in layer 4 with a ramped time variation as in Fig. 3(b), with 𝑡 𝑟 = 0.1 s.

The temperature is imposed on two edges. The strategy consists in updating both the z-functions, then the time functions. First of all, the first nine z-functions and time functions are shown in Fig. 17 and Fig. 16, respectively, for the initial configuration and the new one. We clearly see the influence of the new characteristics of the test case. The location and the size of the heat source has an important influence on the 𝑇 𝑖 𝑧 functions as it can be expected. A high variation occurs through the layer where the heat source is applied. Then, Fig. 18 and Fig. 19 show the accuracy of the strategy. Again, we notice that the computation of the z-functions and time functions is sufficient. Moreover, only two 1D problems have to be solved for a new configuration. Thus, the efficiency of the method is very good.

Remark.

For the test cases involved in the present section, the updating of the two functions 𝑇 𝑧 and 𝑇 𝑡 is sufficient to achieve an accurate solution. If needed, few new 4-uplets can be easily computed to improve the temperature field.

Comments about computational cost

The computational complexity of the present multi-resolution method is evaluated. It is compared with a classical approach in which a whole Layer-Wise computation is carried out for each new set of parameters. By assuming a direct band solver to solve these two types of methods, the estimation gives:

• Classical LayerWise Approach: ∼ 9 • PGD approach computation time for the 𝑇 𝑡 problem: ∼ 1 2 𝑁 3 4𝑢𝑝 𝑁 𝑡 • PGD approach computation time for the 𝑇 𝑧 problem:

∼ 1 2 (𝑁 𝑧 .𝐷 𝑧 ) 3 𝑁 3 4𝑢𝑝
where 𝐷 𝑧 is the order of expansion of the unknowns with respect to z, 𝑁 𝑧 the number of numerical layers, 𝑁𝐷 𝑥 the number of nodes in the x direction, 𝑁 𝑡 the number of time steps and 𝑁 4𝑢𝑝 is the number of 4-uplets built in the PGD process. A high number of sets of parameters is considered, thus, the initial computation cost of the first basis is neglected. This estimation is suitable when the number of nodes 𝑁𝐷 𝑥 and 𝑁 𝑧 are high. So, the most important gain of the PGD approach will be made when 9𝑁𝐷 𝑥 .𝑁 𝑡 ≫ 𝑁 3 4𝑢𝑝 and 9𝐷 3 𝑧 .𝑁 3 𝑧 .𝑁𝐷 𝑥 ≫ 𝑁 3 4𝑢𝑝 . The most important gain of the multiresolution approach will be made when the time and spatial discretization will be refined.

Conclusion

In the present study, a multi-resolution strategy based on a variable separation is performed to model laminated composite beams for unsteady diffusion problems. On the one hand, an explicit solution with respect to the heat source location is built. The solution for any size and location of the heat source can be easily computed with a low computational cost. On the other hand, the basis built for a given configuration can be reused to analyze new test cases. The computations involve only one or two 1D problems to be solved. It allows us to extend the application fields of the method (any temporal evolution of the heat source, location in the thickness of the structure, other composite laminates). The method is assessed by comparing with reference solutions, and the results are very satisfactory. Finally, this strategy avoids to perform numerous analyses for a fixed set of parameters in an iterative process framework (identification, inverse problems, reliability ...).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .

 1 Fig. 1. The laminated beam and coordinate system.

Fig. 2 .

 2 Fig. 2. The laminated beam and heat source.

geometry:= 4 .

 4 composite beam made of several layers with 𝐿 = 1 m. All layers have the same thickness. 𝑆 = 𝐿 ℎ boundary conditions: Temperature is imposed on two edges (𝜃 = 0 at 𝑧 = -ℎ∕2, 𝜃 = 0 at 𝑥 = 0) or on all edges. Different heat sources on 𝐼 1 = [0.3, 0.4], 𝐼 2 = [0.35, 0.4], 𝐼 3 = [0.38, 0.4] with a stepped or a ramped variation are considered (see Fig.

  mesh: A 1D mesh for the x-functions with 𝑁 𝑥 = 200 elements time problem : time domain [0, 𝑡 𝑚𝑎𝑥 ] with 𝑡 𝑚𝑎𝑥 = 0.4 s; The number of time steps : 𝑁 𝑡 = 350 numerical layers: 𝑁 𝑧 is the total number of numerical layers. We have 𝑁 𝑧 = 𝑁𝐶 (number of layers).

Fig. 3 .Fig. 4 .

 34 Fig. 3. Heat source intensity versus time.

Fig.

  Fig. Distribution of temperature along the beam length -𝑡 = 𝑡 𝑚𝑎𝑥 -NC = 3 layers -𝐼 3 = [0.38, 0.4] -heat source in layer 2/3.

Fig. 6 .

 6 Fig. 6. Distribution of temperature along the thickness at x𝑠 -𝑡 = 𝑡 𝑚𝑎𝑥 -NC = 3 layers -𝐼 3 = [0.38, 0.4] -heat source in layer 2/3.

Fig. 7 .

 7 Fig. 7. Distribution of heat flux along the thickness at x𝑠 -𝑡 = 𝑡 𝑚𝑎𝑥 -NC = 3 layers -𝐼 3 = [0.38, 0.4] -heat source in layer 2/3.

8 .

 8 Distribution of temperature with respect to the time -NC = 3 layers -𝐼 3 = [0.38, 0.4] -heat source in layer 2/3.

Fig. 9 .

 9 Fig. 9. Distribution of temperature over the beam at 𝑡 = 𝑡 𝑚𝑎𝑥 -NC = 3 layers -𝐼 3 = [0.38, 0.4] -heat source in layer 2.

Fig. 10 .

 10 Fig. 10. Distribution of temperature over the beam at 𝑡 = 𝑡 𝑚𝑎𝑥 -NC = 3 layers -𝐼 3 = [0.38, 0.4] -heat source in layer 3.

Fig. 11 .

 11 Fig. 11. Multi-resolution (only 𝑇 𝑡 updated) -ramped 𝑡 𝑟 = 0.15 s -3 layers -𝐼 3 = [0.38, 0.4] -heat source in layer 3.

Fig. 12 .

 12 Fig. 12. Multi-resolution (𝑇 𝑧 , then 𝑇 𝑡 updated) -initial heat source in layer 2-3 layers -𝐼 3 = [0.38, 0.4] -heat source in layer 3 -ramped.

Fig. 13 .

 13 Fig. 13. Multi-resolution (𝑇 𝑧 , then 𝑇 𝑡 updated) -initial heat source in layer 2-3 layers -𝐼 3 = [0.38, 0.4] -heat source in layer 3 -ramped.

Fig. 14 .

 14 Fig. 14. Multi-resolution (𝑇 𝑧 , then 𝑇 𝑡 updated) -initial configuration: 3 layers with a heat source in layer 2 -final configuration: 6 layers with a heat source in layer 4 -𝐼 1 = [0.3, 0.4] -ramped.

Fig. 15 .

 15 Fig. 15. Multi-resolution (𝑇 𝑧 , then 𝑇 𝑡 updated) -initial configuration: 3 layers with a heat source in layer 2 -final configuration: 6 layers with a heat source in layer 4 -𝐼 1 = [0.3, 0.4] -ramped.

Fig. 16 .

 16 Fig. 16. Multi-resolution (𝑇 𝑧 , then 𝑇 𝑡 updated) -initial configuration: 3 layers, constant heat source in layer 2 -final configuration: 6 layers, heat source in layer 4 with stepped time variation -𝐼 2 = [0.38, 0.4].

Fig. 17 .

 17 Fig. 17. Multi-resolution (𝑇 𝑧 , then 𝑇 𝑡 updated) -initial configuration: 3 layers, constant heat source in layer 2 -final configuration: 6 layers, heat source in layer 4 with stepped time variation -𝐼 2 = [0.38, 0.4].

Fig. 18 .

 18 Fig. 18. Multi-resolution (𝑇 𝑧 , then 𝑇 𝑡 updated) -initial configuration: 3 layers, constant heat source in layer 2 -final configuration: 6 layers, heat source in layer 4 with stepped time variation -𝐼 2 = [0.38, 0.4].

Fig. 19 .

 19 Fig. 19. Multi-resolution (𝑇 𝑧 , then 𝑇 𝑡 updated) -initial configuration: 3 layers, constant heat source in layer 2 -final configuration: 6 layers, heat source in layer 4 with stepped time variation -𝐼 2 = [0.38, 0.4].

Appendix A .

 . Matrices related to the discretized problemsA.1. Ordinary differential equation to be solved on [0, 𝑡 𝑚𝑎𝑥 ]𝐴( T𝑥 , T𝑧 , T𝑥 𝑠 ) =

  𝛿(𝑇 𝑥 𝑇 𝑧 𝑇 𝑡 𝑇𝑥 𝑠 ) = 𝛿𝑇 𝑥 .𝑇 𝑧 .𝑇 𝑡 .𝑇 𝑥 𝑠 + 𝑇 𝑥 .𝛿𝑇 𝑧 .𝑇 𝑡 .𝑇 𝑥 𝑠 + 𝑇 𝑥 .𝑇 𝑧 .𝛿𝑇 𝑡 .𝑇 𝑥 𝑠 + 𝑇 𝑥 .𝑇 𝑧 .𝑇 𝑡 .𝛿𝑇 𝑥 𝑠

  𝑏(𝑇 𝑥 .𝑇 𝑧 .𝑇 𝑥 𝑠 . , 𝑇 𝑧 .𝑇 𝑡 .𝑇 𝑥 𝑠 .𝛿𝑇 𝑥 ) -𝑎(𝑇 𝑥 .𝑇 𝑧 .𝑇 𝑡 .𝑇 𝑥 𝑠 , 𝑇 𝑡 .𝑇 𝑥 𝑠 .𝑇 𝑧 .𝛿𝑇 𝑥 ) = 𝑑 𝑟 (𝑇 𝑡 .𝑇 𝑥 𝑠 .𝑇 𝑧 .𝛿𝑇 𝑥 ) -𝑏( 𝜕𝜃 𝑚 𝜕𝑡 , 𝑇 𝑡 .𝑇 𝑥 𝑠 .𝑇 𝑧 .𝛿𝑇 𝑥 ) + 𝑎(𝜃 𝑚 , 𝑇 𝑡 .𝑇 𝑥 𝑠 .𝑇 𝑧 .𝛿𝑇 𝑥 ), ∀𝛿𝑇 𝑥 𝑏(𝑇 𝑥 .𝑇 𝑧 .𝑇 𝑥 𝑠 . 𝑇 𝑥 .𝑇 𝑡 .𝑇 𝑥 𝑠 .𝛿𝑇 𝑧 ) -𝑎(𝑇 𝑥 .𝑇 𝑧 .𝑇 𝑡 .𝑇 𝑥 𝑠 , 𝑇 𝑡 .𝑇 𝑥 𝑠 .𝑇 𝑥 .𝛿𝑇 𝑧 ) = 𝑑 𝑟 (𝑇 𝑡 .𝑇 𝑥 𝑠 .𝑇 𝑥 .𝛿𝑇 𝑧 ) -𝑏( 𝜕𝜃 𝑚 𝜕𝑡 , 𝑇 𝑡 .𝑇 𝑥 𝑠 .𝑇 𝑥 .𝛿𝑇 𝑧 ) + 𝑎(𝜃 𝑚 , 𝑇 𝑡 .𝑇 𝑥 𝑠 .𝑇 𝑥 .𝛿𝑇 𝑧 ), ∀𝛿𝑇 𝑧 𝑏(𝑇 𝑥 .𝑇 𝑧 .𝑇 𝑥 𝑠 . , 𝑇 𝑥 .𝑇 𝑧 .𝑇 𝑥 𝑠 .𝛿𝑇 𝑡 ) -𝑎(𝑇 𝑥 .𝑇 𝑧 .𝑇 𝑡 .𝑇 𝑥 𝑠 , 𝑇 𝑥 .𝑇 𝑧 .𝑇 𝑥 𝑠 .𝛿𝑇 𝑡 ) = 𝑑 𝑟 (𝑇 𝑥 .𝑇 𝑧 .𝑇 𝑥 𝑠 .𝛿𝑇 𝑡 ) -𝑏( 𝜕𝜃 𝑚 𝑏(𝑇 𝑥 .𝑇 𝑧 .𝑇 𝑥 𝑠 . , 𝑇 𝑡 .𝑇 𝑥 .𝑇 𝑧 .𝛿𝑇 𝑥 𝑠 ) -𝑎(𝑇 𝑥 .𝑇 𝑧 .𝑇 𝑡 .𝑇 𝑥 𝑠 , 𝑇 𝑡 .𝑇 𝑥 .𝑇 𝑧 .𝛿𝑇 𝑥 𝑠 ) = 𝑑 𝑟 (𝑇 𝑡 .𝑇 𝑥 .𝑇 𝑧 .𝛿𝑇 𝑥 𝑠 ) -𝑏( 𝜕𝜃 𝑚 𝜕𝑡 , 𝑇 𝑡 .𝑇 𝑥 .𝑇 𝑧 .𝛿𝑇 𝑥 𝑠 ) + 𝑎(𝜃 𝑚 , 𝑇 𝑡 .𝑇 𝑥 .𝑇 𝑧 .𝛿𝑇 𝑥 𝑠 ), ∀𝛿𝑇 𝑥 𝑠

	𝜕𝑇 𝑡
	𝜕𝑇 𝑡
	(11)
	𝜕𝑇 𝑡 𝜕𝑡 , (12)
	𝜕𝑇 𝑡

𝜕𝑡 𝜕𝑡 𝜕𝑡 , 𝑇 𝑥 .𝑇 𝑧 .𝑇 𝑥 𝑠 .𝛿𝑇 𝑡 ) + 𝑎(𝜃 𝑚 , 𝑇 𝑥 .𝑇 𝑧 .𝑇 𝑥 𝑠 .𝛿𝑇 𝑡 ), ∀𝛿𝑇 𝑡 (13) 𝜕𝑡

  1 to 𝑘 𝑚𝑎𝑥 do Compute T (𝑘)

	T (𝑘-1) 𝑡	𝑧 , T (𝑘-1) 𝑥 𝑠	from Eq. (12) (linear equation on  𝑧 ), T (𝑘-1) 𝑥 being known	,
	Compute T (𝑘) 𝑡	from Eq. (13), T (𝑘-1) 𝑥	, T (𝑘) 𝑧 , T (𝑘-1) 𝑥 𝑠	being known
	Compute T (𝑘) 𝑥 from Eq. (11) (linear equation on  𝑥 ), T (𝑘) 𝑧 , T (𝑘) 𝑡 ,
	T (𝑘-1) 𝑥 𝑠	being known
	Compute T (𝑘) 𝑥 𝑠 from Eq. (14) (linear equation on  𝑥 𝑠 ), T (𝑘) 𝑥 , T (𝑘) 𝑡 ,
	T (𝑘) 𝑧 being known
	Check for convergence
	end for
	Set 𝑇 𝑚+1 𝑥

  ) is introduced. Classical finite element approximation of the functions (𝑇 𝑥 , 𝑇 𝑧 , 𝑇 𝑥 𝑠 ) in  𝑥 ,  𝑧 and  𝑥 𝑠 respectively is used. The elementary vectors of degree of freedom (dof) associated with the finite element mesh in  𝑥 ,  𝑧 and  𝑥 𝑠 are denoted

	𝐪 𝑥 𝑒 , 𝐪 𝑧 𝑒 and 𝐪 𝑥 𝑠 𝑒 , respectively. The temperature fields and the associated gradient are determined from the values of 𝐪 𝑥 𝑒 , 𝐪 𝑧 𝑒 and 𝐪 𝑥 𝑠 𝑒 by
	𝑇 𝑥𝑒 = 𝐍 𝑥 𝐪 𝑥 𝑒 , 𝑇 𝑧𝑒 = 𝐍 𝑧 𝐪 𝑧 𝑒 , 𝑇 𝑥 𝑠 𝑒 = 𝐍 𝑥 𝑠 𝐪 𝑥 𝑠 𝑒 ,	 𝑒 𝑥 = 𝐁 𝑥 𝐪 𝑥 𝑒 ,  𝑒 𝑧 = 𝐁 𝑧 𝐪 𝑧 𝑒 ,

  T𝑡 , T𝑥 𝑠 ) 𝑧 𝑑𝑧 = ∫  𝑧

	𝛿 𝑇 𝑧	[ 𝐐 𝑟 𝑥𝑡𝑥 𝑠	( T𝑥 , T𝑡 , T𝑥 𝑠 )
	+ 𝐐 𝑚 𝑥𝑡𝑥 𝑠

( T𝑥 , T𝑡 , T𝑥 𝑠 , 𝜃 𝑚 ) ] 𝑑𝑧

[START_REF] Movahedi-Rad | Fatigue damage in angle-ply gfrp laminates under tension-tension fatigue[END_REF]

where

𝝀

𝑥𝑡𝑥 𝑠 , 𝐐 𝑟 𝑥𝑡𝑥 𝑠 and 𝐐 𝑚 𝑥𝑡𝑥 𝑠 are given in Appendix A.3. The introduction of the finite element approximation Eq. (

(

  𝑠 ( 𝑥 𝑡 𝑥 𝑠 ) and the equilibrium residual  𝑥𝑡𝑥 𝑠 ( 𝑥 𝑡 𝑥 𝑠 , 𝜃 𝑚 ) are obtained by summing the following elementary matrices and vectors: T𝑥 , T𝑡 , T𝑥 𝑠 ) = ∫ 𝑧 𝑒 𝐁 𝑇 𝑧 𝝀 𝑥𝑡𝑥 𝑠 ( T𝑥 , T𝑡 , T𝑥 𝑠 ) 𝐁 𝑧 𝑑𝑧

	𝜦 𝑒 𝑥𝑡𝑥 𝑠					(18)
	and
	 𝑒 𝑥𝑡𝑥 𝑠	( T𝑥 , T𝑡 , T𝑥 𝑠 , 𝜃 𝑚 ) = ∫ 𝑧 𝑒	𝐁 𝑇 𝑧	[ 𝐐 𝑟 𝑥𝑡𝑥 𝑠	( T𝑥 , T𝑡 , T𝑥 𝑠 )
		+ 𝐐 𝑚 𝑥𝑡𝑥 𝑠	( T𝑥 , T𝑡 , T𝑥 𝑠 , 𝜃 𝑚 ) ]	𝑑𝑧	(19)
	3.5. Ordinary differential equation to be solved on [0, 𝑡 𝑚𝑎𝑥 ]

  The introduction of the finite element approximation Eq. (15) in the variational Eq. (21) leads to the linear system𝜦 𝑧𝑡𝑥 𝑠 ( T𝑧 , T𝑡 , T𝑥 𝑠 ) 𝐪 𝑥 =  𝑧𝑡𝑥 𝑠 ( T𝑧 , T𝑡 , T𝑥 𝑠 , 𝜃 𝑚 )(22)where 𝐪 𝑥 is the vector of the nodal temperatures associated with the finite element mesh in  𝑥 , 𝜦 𝑧𝑡𝑥 𝑠 ( T𝑧 , T𝑡 , T𝑥 𝑠 ) the matrix obtained by summing the elements' conductivity and capacity matrices involved in 𝜦 𝑒 𝑧𝑡𝑥 𝑠 ( T𝑧 , T𝑡 , T𝑥 𝑠 ).  𝑧𝑡𝑥 𝑠 ( T𝑧 , T𝑡 , T𝑥 𝑠 , 𝜃 𝑚 ) is the equilibrium residual obtained by summing the elements' residual load vectors  𝑒

						 𝑥	𝛿 𝑇 𝑥	[ 𝐐 𝑟 𝑧𝑡𝑥 𝑠	( T𝑧 , T𝑡 , T𝑥 𝑠 )
						+ 𝐐 𝑚 𝑧𝑡𝑥 𝑠	( T𝑧 , T𝑡 , T𝑥 𝑠 , 𝜃 𝑚 ) ]	𝑑𝑥	(21)
	where 𝝀 𝑧𝑡𝑥 𝑠 , 𝐐 𝑟 𝑧𝑡𝑥 𝑠	and 𝐐 𝑚 𝑧𝑡𝑥 𝑠	are given in Appendix A.2.
	T𝑥 𝑠 , 𝜃 𝑚 )				𝑧𝑡𝑥 𝑠	( T𝑧 , T𝑡 ,
	𝜦 𝑒 𝑧𝑡𝑥 𝑠	( T𝑧 , T𝑡 , T𝑥 𝑠 ) = ∫ 𝐿 𝑒	𝐁 𝑇 𝑥 𝝀 𝑧𝑡𝑥 𝑠 ( T𝑧 , T𝑡 , T𝑥 𝑠 ) 𝐁 𝑥 𝑑𝑥	(23)
	and
	 𝑒 𝑧𝑡𝑥 𝑠	( T𝑧 , T𝑡 , T𝑥 𝑠 , 𝜃 𝑚 ) = ∫ 𝐿 𝑒	𝐁 𝑇 𝑥	[ 𝐐 𝑟 𝑧𝑡𝑥 𝑠	( T𝑧 , T𝑡 , T𝑥 𝑠 )
				+ 𝐐 𝑚 𝑧𝑡𝑥 𝑠	( T𝑧 , T𝑡 , T𝑥 𝑠 , 𝜃 𝑚 ) ]	𝑑𝑥	(24)

𝑥 𝝀 𝑧𝑡𝑥 𝑠 ( T𝑧 , T𝑡 , T𝑥 𝑠 ) 𝑥 𝑑𝑥 = ∫

(

  ∫  𝑥 𝑠 T 2 𝐶 𝑟 (𝑡, T𝑥 , T𝑧 , T𝑥 𝑠 ) = 𝑟 𝑑 (𝑡) ∫  𝑠𝑜𝑢𝑟𝑐𝑒 𝑧 T𝑧 𝑑𝑧 ∫  𝑥 T𝑥 𝑠 T𝑥 𝑑𝑥 𝐶 𝑚 (𝑡, T𝑥 , T𝑧 , T𝑥 𝑠 , 𝜃 𝑚 ) = A.2. FE problem to be solved on  𝑥 The matrices involved in Eq. (21) are given as: 𝝀 𝑧𝑡𝑥 𝑠 ( T𝑧 , T𝑡 , T𝑥 𝑠 ) = ∫ T𝑧 , T𝑡 , T𝑥 𝑠 , 𝜃 𝑚 ) = ∫ A.3. FE problem to be solved on  𝑧The matrices involved in Eq. (16) are given as:𝝀 𝑥𝑡𝑥 𝑠 ( T𝑥 , T𝑡 , T𝑥 𝑠 ) = ∫ 𝑧 𝜮 𝐱 ( T𝑥 ) 𝑇 𝝀 𝜮 𝐱 ( T𝑥 )𝑑𝑥

	𝐐 𝑟 𝑧𝑡𝑥 𝑠	( T𝑧 , T𝑡 , T𝑥 𝑠 ) = ∫  𝑠𝑜𝑢𝑟𝑐𝑒 𝑧	T𝑧 𝑑𝑧 ∫	0	𝑡 𝑚𝑎𝑥	𝑟 𝑑 (𝑡) T𝑡 𝑑𝑡	[	0 T𝑥 𝑠	]	(A.5)
	𝐐 𝑚 𝑧𝑡𝑥 𝑠			0	𝑡 𝑚𝑎𝑥	∫  𝑧 ∫  𝑥 𝑠	𝐠𝐫𝐚𝐝(𝛿𝜃) 𝑇 𝝀𝐠𝐫𝐚𝐝(𝜃 𝑚 )
		+ 𝛿𝜃𝜌𝑐 𝑝	𝜕𝜃 𝑚 𝜕𝑡	𝑑 𝑥 𝑠 𝑑 𝑧 𝑑𝑡	(A.6)
		0	𝑡 𝑚𝑎𝑥	T 2 𝑡 𝑑𝑡
		∫  𝑥 𝑠 𝑑𝑥 𝑠 ∫  + ∫ T 2 𝑥 𝑠 𝑡 𝑚𝑎𝑥 0 T𝑡 𝜕𝑡 [ 𝜌𝑐 𝑝 0 ] 𝑑𝑡 ∫  𝑥 𝑠 𝜕 T𝑡	T 2 𝑥 𝑠	𝑑𝑥 𝑠 ∫  𝑥	T 2 𝑥 𝑑𝑥	(A.7)
				0				0
	𝐐 𝑟 𝑥𝑡𝑥 𝑠
											𝐵( T𝑥 , T𝑧 , T𝑥 𝑠 ) = ∫  𝑥 𝑠	T 2 𝑥 𝑠	𝑑𝑥 𝑠 ∫  𝑥	 𝑥 ( T𝑥 ) 𝑇
											[
											×	∫  𝑧	𝜮 𝑚
											∑ 𝑖=1 ∫  𝑥 ∫  𝑥 𝑠 T𝑥 𝑇 𝑖 𝑥 𝑑𝑥 T𝑥 𝑠 𝑇 𝑖 𝑥 𝑠 𝜕𝑇 𝑖 𝑑𝑥 𝑠 ∫  𝑧 𝑡 𝜕𝑡 𝑚	𝜌𝑐 𝑝 T𝑧 𝑇 𝑖 𝑧 𝑑𝑧
											∑ 𝑖=1 [	∫  𝑥 𝑠	T𝑥 𝑠 𝑇 𝑖 𝑥 𝑠	𝑑𝑥 𝑠 ∫  𝑥	 𝑥 ( T𝑥 ) 𝑇 ]
											∫  𝑧	𝜮 𝐳 ( T𝑧 ) 𝑇 𝝀 𝜮 𝐳 (𝑇 𝑖 𝑧 )𝑑𝑧	 𝑥 (𝑇 𝑖 𝑥 ) 𝑑𝑥𝑇 𝑖 𝑡 (𝑡)
											(A.3)
											[	0	𝑡 𝑚𝑎𝑥 𝑑𝑥 𝑡 𝑚𝑎𝑥 T 2 𝑡 𝑑𝑡 ∫  𝑥 𝑠 T 2 𝑥 𝑠 0 T𝑡 𝜕 T𝑡 𝜕𝑡 T 2 𝑥 𝑠 𝑑𝑡 ∫  𝑥 𝑠 0 0 ]	𝑑𝑥 𝑠 ∫  𝑧	𝜌𝑐 𝑝	T 2 𝑧 𝑑𝑧
											0 1
		𝑥 𝑠	𝑑𝑥 𝑠 ∫  𝑧	𝜌𝑐 𝑝	T 2 𝑧 𝑑𝑧 ∫  𝑥	T 2 𝑥 𝑑𝑥	(A.1)	(A.4)

𝐳 ( T𝑧 ) 𝑇 𝝀𝜮 𝐳 ( T𝑧 )𝑑𝑧 ]  𝑥 ( T𝑥 ) 𝑑𝑥 (A.2) 𝑠 ∫  𝑧 𝜮 𝐳 ( T𝑧 ) 𝑇 𝝀 𝜮 𝐳 ( T𝑧 )𝑑𝑧 + ∫

(

  T𝑥 , T𝑡 , T𝑥 𝑠 ) = ∫  𝑥 T𝑥 T𝑥 𝑠 𝑑𝑥 ∫

	0	𝑡 𝑚𝑎𝑥	𝑟 𝑑 (𝑡) T𝑡 𝑑𝑡	[	1 0	]	(A.8)
	𝐐 𝑚 𝑥𝑡𝑥 𝑠

(

  T𝑥 , T𝑡 , T𝑥 𝑠 , 𝜃 𝑚 ) = ∫

	0	𝑡 𝑚𝑎𝑥	∫  𝑥 ∫  𝑥 𝑠	𝐠𝐫𝐚𝐝(𝛿𝜃) 𝑇 𝝀𝐠𝐫𝐚𝐝(𝜃 𝑚 )
	+ 𝛿𝜃𝜌𝑐 𝑝	𝜕𝜃 𝑚 𝜕𝑡	𝑑 𝑥 𝑠 𝑑 𝑥 𝑑𝑡	(A.9)
	Appendix B.

Computation of the temperature with the superposi- tion principle

  The field 𝜃(𝑥, 𝑧, 𝑡, 𝑥 𝑠 ) being computed, it is possible to deduce easily the temperature for a heat source applied on a domain 𝐼 𝑘 such that 𝑇 (𝑥, 𝑧, 𝑡) = ∫ 𝐼 𝑘 𝜃(𝑥, 𝑧, 𝑡, 𝑥 𝑠 ) 𝑑𝑥 𝑠 =

	𝑁
	∑ 𝑖=1	𝑇 𝑖 𝑥 (𝑥)𝑇 𝑖 𝑧 (𝑧)𝑇 𝑖 𝑡 (𝑡) ∫ 𝐼 𝑘	𝑇 𝑖 𝑥 𝑠

(𝑥 𝑠 ) 𝑑𝑥 𝑠 (B.1)


Joel Minguet Lopez
email: joel.minguetlopez@cea.fr

Quentin Rafhay

Manon Dampfhoffer

Lucas Reganaz

Niccolo Castellani

Valentina Meli

Simon Martin

Laurent Grenouillet

Gabriele Navarro

Thomas Magis

Catherine Carabasse

Tifenn Hirtzlin

Elisa Vianello

Damien Deleruyelle

Jean-Michel Portal

Gabriel Molas

François Andrieu

S1R optimization for high-frequency inference on binarized spiking neural networks

Keywords: 1S1R, crossbar, ovonic threshold switch (OTS), resistive random-access memory (RRAM), spiking neural networks

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

Introduction

Deep learning accelerators have attained impressive performance on various cognitive tasks, such as image or audio recognition. However, its enormous electrical consumption a synapse an input spike, which flows through the dendrite until reaching the soma. This input spike is integrated into the soma of the neuron and contributes to increasing its membrane potential. When the neuron membrane potential reaches a certain threshold, it generates a spike that is transmitted to other neurons (post-neurons) through the axon. As a consequence, the neuron returns to its resting state (Figure 1A). The emerging bio-inspired spiking neural networks (SNNs) provide the opportunity to emulate these characteristics in artificial systems. Indeed, while classical ANNs rely on layer-by-layer multiply-and-accumulate operations between input activations and weights, SNNs computation is mostly based on the accumulation of spikes weighted by synaptic connections.

Therefore, only Accumulate operations are required in SNNs. In addition, SNNs can be processed in an event-based manner to benefit from their natural sparsity. Indeed, synaptic weight readings are triggered by the spikes received by a neuron. One should notice that the synaptic weights are not read in the absence of input spikes. Therefore, a low spiking activity on the SNN results in very low energy consumption compared to an ANN where all synaptic weights are read at each inference.

Various NVM devices have been considered for NN synaptic weight hardware implementation in crossbar arrays. [START_REF] Chen | [END_REF][4][5][6][7][8] Nevertheless, at the moment, there is no ideal candidate due to NVM technologies imperfection, which can strongly limit the overall accelerator inference accuracy. [9] Accordingly, considerable efforts have been invested in the precedent years to develop innovative techniques to deal with memory nonidealities. [10] Among others, resistive random-access memory (RRAM) devices are one of the most promising technologies proposed for NN hardware implementation. [10][11][12] However, the RRAM intrinsic operating variability remains challenging. In this context, innovative bit-error correcting codes and adaptive programming schemes have been developed. [START_REF] Jain | ISSCC[END_REF][START_REF] Chou | VLSI Circuits[END_REF] In addition, innovative and aggressive programming patterns have been considered to reduce device variability. Nevertheless, an important degradation of the device operating endurance capabilities appears as a counter-part. [START_REF] Nail | [END_REF]16] In this work, we focus on low-precision Binarized Spiking Neural Networks (BSNNs), which benefit from very high energy efficiency while promising satisfactory tolerance to RRAM nonidealities. [17] In this type of network, the synaptic weights are implemented by binarized values (-1 and +1) after the training process. Moreover, the standard 1T1R (1 Transistor 1 Resistor) architectures used to implement synaptic weights with memory devices are replaced by a denser 1S1R (1 Selector 1 Resistor) stack, where the memory (1R) is co-integrated in series with a back-end selector (1S). About one order of magnitude improvement in terms of memory density can be achieved with 1S1R with respect to conventional 1T1R architectures. [18] However, RRAM co-integration with a back-end selector implies several challenges. First, the process integration complexity increases. Second, achieving high memory capacity while preserving reduced electrical consumption becomes challenging due to a trade-off between programming voltages and leakage currents. Third, to achieve high precision on device resistive states implies a degradation of its programming endurance capabilities due to a trade-off between memory window and endurance capabilities. To deal with those challenges, a specific design of both memory stack and applied programming conditions for the application of interest is required. [19][20][21][START_REF] Lopez | presented at[END_REF] In this context, the co-integration of an HfO 2 -based OxRAM device with an ovonic threshold switch (OTS) selector on single memristor crossbar arrays to implement the BSNN synaptic weights (Figure 1B)) is proposed. OTS co-integration with HfO 2 -based OxRAM has been satisfactorily demonstrated in the previous years. [18,[20][21][START_REF] Lopez | presented at[END_REF][START_REF] Lopez | [END_REF] Through stack design and applied programming conditions adaptation, the 1S1R dynamic switching capabilities have been elucidated, the 1S1R binarized window margin optimized, and its programming endurance capabilities enlarged. [START_REF] Lopez | presented at[END_REF] In particular, the 1S1R pertinence for standard low-precision synaptic weight encoding during network training on-chip has been demonstrated. [18] However, beyond 1S1R-based crossbar arrays pertinence for low-precision neural network training, it remains essential to elucidate the 1S1R ability to perform reliable high-frequency inference on-chip. Inference in Adv. Electron. Mater. 2022, 2200323 Figure 1. A) Biological neuron dynamics. An input spike flows through the neuronal dendrite until it is integrated into the soma and contributes to increase the membrane potential. When the neuron membrane potential reaches a certain threshold, an output spike is generated and transmitted to other neurons. Then, the neuron returns to its resting state. B) A 1S1R-based single memristor crossbar array is proposed in this work to implement the binarized spiking neural network synaptic weights.

SNNs requires a huge amount of repeated and frequent reading operation of 1S1R devices to sense the network synaptic weights. While reading is expected to be nondestructive, it can induce erratic switching of the selector and memory devices and affect network performance. This reliability issue has never been addressed in the literature, and the link between device physics and circuit accuracy has not been studied and clarified so far.

In this work, we first present an experimental statistical study on OTS switching probability when repeated sub-threshold voltages are applied and propose to reduce its variability through device lifetime by optimizing the applied voltage on the devices. To gain insight into the OTS switching operation microscopic mechanisms, a Monte Carlo statistical model based on Bernoulli's conduction point is implemented. The overall 1S1R reading reliability being directly linked to the OTS switching variability, and general guidelines on 1S1R reading conditions (applied voltages, reading frequency) optimization for low reading bit error rate (BER) are provided based on this experimental and theoretical statistical study. To evaluate the 1S1R pertinence for BSNN inference hardware implementation, training simulations are performed on one hidden layer fully connected BSNN for an image classification task on the MNIST dataset. By introducing errors in the synaptic weights during training, an optimized BSNN tolerance to synaptic errors for the devices of interest is demonstrated. Based on this analysis, guidelines for an optimized system footprint, a reduced electrical consumption, maximized inference frequency, and maximized BSNN accuracy for the MNIST task are provided.

Technological Details

OxRAM is co-integrated with an OTS back-end selector in a 4kb 1S1R array configuration, where a transistor is used to limit the current on the devices. The OTS is composed of a 10 nm-thick Ge-Se-Sb-N (GSSN) alloy, sandwiched between two Carbon electrodes, and is used as a selector device. Then, a 10 nm-thick HfO 2 -based OxRAM deposited by Atomic Layer Deposition (ALD) is used as the memory device. This HfO 2 is deposited on a TiN inter-layer that separates the selector and the memory elements. A 10 nm Physical Vapor Deposition (PVD) Ti top layer acts as an oxygen scavenging layer for the memory. A TiN additional layer completes the top electrode. After the etching process, the memory dot is capped with a SiN layer for passivation. Afterward, the top contact and top metal line end the overall integration process. The 1S1R memory array is thus integrated into the Back-End-Of-Line (BEOL) of a 130 nm CMOS process between the fourth and fifth metal layers. Figure 2A provides an SEM image of the 1S1R devices integrated into the Back-End-Of-Line (BEOL), as well as a TEM cross-section for the stack of interest, which demonstrates a consistent co-integration of all the deposited layers.

Results and Discussion

1S1R Programming and Reading Operation

Figure 2B provides typical 1S1R current-voltage characteristics after the device firing operation, which is required for both OTS and OxRAM initialization. The 1S1R switching voltages are strongly impacted by the resistive state of the OxRAM. If the OxRAM is in the Low Resistive State (LRS) (resp. High Resistive State (HRS)), V th-LRS (resp. V th-HRS > V th-LRS ) is required for the 1S1R switching. The 1S1R read window margin is defined by (V th-HRS -V th-LRS ). Figure 2C provides experimental programming endurance characteristics for the stack of interest for eight 1S1R devices. No error is observed until 10 4 programming cycles, demonstrating the ability to encode binarized weights on the memory. Since only inference is targeted on chip, weights have only to be programmed once per application, and thus 10 4 cycles are sufficient. Nevertheless, 1S1R stack adjustments and optimized programming patterns as well as OTS material engineering are paths of improvement to enhance endurance for more cycling demanding applications. [START_REF] Lopez | presented at[END_REF][START_REF] Lopez | [END_REF] In this context, the OxRAM resistive state can be read by applying a certain reading voltage (V read ) caught between V th-HRS and V th-LRS . If the OxRAM is in LRS, the OTS switching occurs, and I LRS is read. If the OxRAM is in HRS, the OTS device remains at the OFF state, and I HRS is read. For a 200 ns long standard trapezoidal reading pulse, Figure 2D shows experimental 1S1R reading disturb characteristics for a given relaxation time t relax ≈100 µs between the applied pulses. No error is observed up to 10 9 cycles. However, provided the same reading pattern and V read , Figure 2E presents experimental 1S1R HRS reading disturb characteristics when t relax ≈1 µs. Through the cycles, the read HRS current progressively increases until an erratic OTS switching occurs after 10 7 cycles, which corresponds to the first failing cycle noted Cycle 1st-fail . I HRS is read on the subsequent cycles, which suggests that the OTS can nonetheless relax back to the OFF state and that this phenomenon is reversible. The apparition of reading fails on 1S1R devices being linked to the OTS switching reliability at V read , and the main goal of the following sections is to explore the OTS dynamic switching behavior. The reading pulses characteristics provided in Figure 2D,E are used for the following sections.

2D Monte Carlo Model for OTS Switching Operation

To study the OTS reading reliability from a phenomenological point of view, a novel 2D OTS reading endurance model is presented in this section. The main goal here is to statistically quantify the OTS switching statistics leading to the apparition of reading errors during device lifetime, which hence suggests a Monte Carlo simulation approach.

Conceptually, we focus on an OTS filamentary field-driven switching theory, where the OTS switching operation is described by the nucleation of metastable domains within the chalcogenide layer. [START_REF] Karpov | [END_REF] In particular, after the initialization of the device, the application of a certain electrical field at the OTS terminals leads to the random nucleation of conductive GSSN dots in the chalcogenide layer. At a microscopic scale, it can be foreseen that these local conductive dots correspond to metastable metavalent bonding formations, resulting from a change in the local order and GSSN bond alignment. [25] The alignment of these conductive dots implies the appearance of a conductive filament and leads to OTS switching to the ON state. This approach is schematically summarized in Figure 3A.

Figure 3B schematically illustrates the OTS 2D Monte Carlo reading endurance model deployment, which remains compatible with the conceptual switching approach introduced previously. A homogeneous OTS bulk is considered at the beginning of each reading disturb simulation. Through the application of a virtual pulse, conductive dots can randomly appear in the GSSN matrix for every reading cycle. Between two successive cycles, a certain relaxation time is allowed for the device. On the one hand, the nucleation of a conductive defect on the constriction layer is driven by a Bernouilli probability (P n ) for each bin of the GSSN matrix. On the other hand, the disappearance of a conductive defect in one bin is assumed to follow an exponential probabilistic distribution characterized by a time τ d . Thus, conductive dots on the constriction zone are considered reversible and can thus switch back to a nonconductive state at every reading cycle. The vertical alignment of a column of conductive bins forms a conductive path, which is the condition chosen for reading endurance failure. To emulate the read cycling process, each position on the constriction zone is scanned at each V th-HRS > V th-LRS ) is required for 1S1R switching operation. In this context, the application of a reading voltage V read between V th-HRS and V th-LRS allows to identify the OxRAM resistive state. C) 1S1R programming endurance capabilities. Up to 10 4 programming cycles are demonstrated without errors. D) Experimental 1S1R reading disturb characteristics at V read = 2.25V. A 100 µs relaxation time between pulses is allowed to the device. No error is observed up to 10 9 cycles. E) Experimental 1S1R high resistive state current reading disturb characteristics at V read = 2.25 V. A short 1µs relaxation time between pulses is allowed to the device, which increases the probability of OTS erratic switching with respect to longer relaxation times. In this case, the OTS switches after 10 7 reading cycles, leading to the appearance of a reading fail.

reading cycle, and its resistive state is updated based on P n and τ d simulation parameters. The first cycle at which a failure is detected corresponds to Cycle 1st-fail , whereas the number of failures over a given amount of cycles is used to estimate the switching probability.

OTS Switching Reliability Optimization

Figure 4A presents the Cycle 1st-fail dependence with relaxation times and applied reading voltages on OTS devices. The boxplots represent the (25%; 75%) percentiles. Increasing V read closer to the OTS intrinsic switching voltage lowers the median of Cycle 1st-fail . In addition, as the relaxation time between reading cycles increases, Cycle 1st-fail median is shifted to higher values for a given V read . This suggests a more robust OTS OFF state when a longer relaxation time between applied pulses is achieved. Moreover, Figure 4B The OTS switching operation is driven by an applied electrical field at the device terminals, which induces the nucleation of randomly spaced conductive dots on the bulk of the OTS. [START_REF] Karpov | [END_REF] If the conductive defects are aligned, a conductive filament appears on the device and leads to the OTS switching. B) OTS Monte Carlo reading endurance model deployment. The nucleation of a conductive defect on the constriction layer is driven by a certain probability (P n ). The disappearance of a conductive dot on the constriction zone is driven by a certain recovery time (τ d ). The alignment of several defects implies the apparition of a conductive filament on the constriction zone and leads to a reading fail. Monte Carlo model for OTS switching operation is used to provide physical insights on Cycle 1st-fail distributions. To capture the experimental data, the following dual behavior of the switching probability P n along cycling is considered:

• Under mild reading conditions (either low reading voltage V read or long rest time between subsequent reading cycles), the OTS is assumed to have sufficient time to relax between reading cycles so that the switching probability remains constant and the OTS state is unaltered. This is supported by the flat behavior of the OTS sub-threshold current, as provided in Figure 2D. • Under aggressive reading conditions (either high reading voltage V read or short rest time between subsequent reading cycles), gradual aging of the OTS is assumed and leads to a gradual increase of the switching probability. This is supported by the drift of the OTS sub-threshold current, as observed in Figure 2E.

Once the statistic of first switching is clarified, the next step consists in extending the study to subsequent OTS switching events through the device lifetime. Thus, OTS cycle-to-cycle switching probability has to be extracted. In this aim, the OTS current is measured ten times per decade when sub-threshold voltages are applied. The experience is repeated multiple times, and OTS switches are detected. Then, the OTS switching probability is calculated by quantifying the percentiles of switching events among all the experiments at a fixed reading cycle (Figure 5A)). The main goal here is to quantify the evolution of the OTS switching probability in the sub-threshold regime regarding the applied reading conditions (reading voltages and reading frequency). Figure 5B presents the evolution with the cycling of OTS switching probability at V read = 2.35V, as a function of the relaxation time between reading operations. ≈110 tests are considered per condition. Again, shorter OTS relaxation time leads to a larger OTS switching probability.

Adv. Electron. Mater. 2022, 2200323 Figure 5. A) OTS switching probability description over cycling. Multiple pulses are applied in sub-threshold regime. By extracting the switching events at every cycle, the OTS switching probability is calculated. B) Evolution of the OTS switching probability with cycling at V read = 2.35 V, as a function of the relaxation time between reading operations. More frequent OTS switching is evidenced for shorter relaxation times between reading pulses. C) Experimental and simulated OTS switching probability after 10 6 reading cycles (P SwitchOTS @ 10 6 cycle ) evolution with the median Cycle 1st-fail . Cycle 1st-fail appears to be a direct image of the device P OTS @ 10 6 cycle . D) Experimental and simulated P SwitchOTS @ 10 6 cycle evolution with the reading frequency. Three different reading voltages (V read = 2.15V; 2.25V; 2.35V) are considered. The higher the reading frequency, the larger the resulting switching probability on the devices for a certain V read .

Moreover, the experimental and simulated OTS switching probabilities after 10 6 reading cycles (P SwitchOTS @ 10 6 cycle ) as function as the median Cycle 1st-fail is presented in Figure 5C. Cycle 1st-fail appears to be a direct image of P SwitchOTS @ 10 6 cycle . In addition, the experimental and simulated P SwitchOTS @ 10 6 cycle evolution with the reading frequency is illustrated in Figure 5D, together with its dependence with V read . In agreement with precedent results, the higher the reading frequency, the larger the resulting OTS switching probability for a certain reading voltage. Thus, the lower V read , the farther to the OTS intrinsic switching voltages, and so the minor the resulting switching probability for a given reading frequency. In this context, a similar drift in the probability P n is required for the simulation to reproduce the experimental results. All in all, successful agreement between experimental and simulated data is evidenced.

1S1R Reading Reliability Optimization

Based on the previous experimental and simulated OTS read operation reliability analysis, OTS switching statistics, under applied pulses in the sub-threshold regime, have been elucidated. In this section, we propose to estimate the overall 1S1R reading BER, which is directly impacted by OTS erratic switching events during the reading operation. Figure 6A illustrates the influence of the OTS reading voltages and reading frequency on the OTS switching probability after 10 6 reading cycles. The symbol size represents the resulting switching probability. This figure generalizes the trends described in the previous section and shows that a trade-off between reading frequency and voltage amplitude is required to prevent OTS erratic switching. During 1S1R reading operation, the OTS behavior depends on the state of the OxRAM memory device: when the OxRAM is in LRS (resp. HRS), the OTS is expected to switch (resp. to remain in the OFF state) at V read . On the one hand, 1S1R HRS BER is governed by the OTS probability to switch from OFFto-ON state at V read , which corresponds to the OTS switching probability studied in the previous sections. On the other hand, LRS BER represents the OTS probability to remain in the OFF state at V read . Practically, LRS BER is obtained by flipping HRS BER (P OTS-ON = 1-P OTS-OFF ) and shifting it by the 1S1R window margin. Indeed, this shift corresponds to the read window margin (WM) and is equal to the voltage that drops on the memory in the HRS. The higher the memory resistance, the larger the threshold voltage shift (WM = V th-HRS -V th-LRS ). It should be noticed that, for a first-order approximation, only the mean value for both OxRAM LRS and HRS resistive states is considered. The variability in the 1S1R memory window, which could be at the origin of a dissymmetrical behavior between LRS and HRS, is not taken into account here. In this context, high reading voltage increases the OTS switching probability and so degrades HRS BER while improving the LRS BER. Figure 6B,C) identify the optimal V read allowing minimizing the overall 1S1R reading BER for both 10 and 4 MHz reading frequencies. A standard 1S1R read window margin of 700 mV is considered for the analysis. [18] The lower the 1S1R reading frequency, the lower the optimal reading BER. Altogether, a 1S1R reading BER of ≈10 -1 (resp. ≈4.5 × 10 -2 ) is estimated after 10 6 reading operations at 10 MHz (resp. 4 MHz) reading frequency. Improved performances could be expected for a larger window margin or lower reading frequency.

Binarized Spiking Neural Network Figures of Merit

To explore the benefit of a BSNN inference hardware implementation with 1S1R-based crossbar arrays, training simulations on a fully-connected BSNN with one hidden layer for an image classification task on the MNIST dataset are performed in this section. Figure 7A presents the considered BSNN topology, where the amount of neurons on the hidden layer is a variable of the study. Neurons of SNNs integrate spike inputs over time and thus are simulated with temporal dynamics. However, in our case, no significant network performance improvement is observed using several simulation timesteps (Figure 7B). Therefore, aiming to keep energy efficiency as high as possible, we chose to simulate our BSNN with a unique timestep. To perform the classification with our BSNN, the MNIST images are converted into spikes by using a pixel intensity threshold. If the pixel intensity is above (resp. below) the threshold, the input is set to +1 (a spike) (resp. 0). In this context, input sparsity is engineered by adapting the global pixel intensity threshold on the input images. The location of the pixel in the image is not considered for sparsity engineering. The BSNN is trained with standard deep learning techniques on the MNIST training dataset and evaluated on the validation dataset. The resulting maximal attainable BSNN accuracy evolution with memory BER is provided in Figure 7C. The number of neurons on the hidden layer is demonstrated to have an important influence on the network performance. [18] The higher the number of neurons, the bigger the number of learnable parameters of the network (binarized weights), and so the higher the accuracy, up to some extent. Indeed, by increasing too much the number of neurons on the hidden layer, the network may be subject to the overfitting phenomenon. Unfortunately, the reading BER on the devices of interest is demonstrated to degrade the network performance. In this context, we propose to improve the BSNN tolerance to parasitic bit errors by adapting the network training strategy. [START_REF] Hirtzlin | presented at AICAS[END_REF] To this aim, the BSNNs of interest are retrained, including bit errors in the weights during the training process, which allows anticipating the apparition of reading errors during the network inference. In particular, at each training iteration, the weights are randomly switched to their opposite with a probability equal to the target BER. Figure 7D presents the resulting maximal attainable BSNN accuracy evolution with memory BER during testing for different target BER during training. The network performance improves for the BER introduced during network training. However, the network becomes less performant for smaller BERs than the one used for training. This trend is generalized in Figure 7E, which presents the optimized inference accuracy evolution with BER with a training strategy adapted to each case. The 1S1R experimental optimized BER characteristics in this work are now demonstrated to be perfectly tolerated by the networks and not to induce any important accuracy degradation.

Therefore, the ability to apply the optimal V read on the devices in a Crossbar environment remains a key condition not to degrade the overall network accuracy. Based on 28 nm technological node resistive rules in metal lines, and while the optimal V read is applied on the word-line extremity, Figure 8A provides the applied reading voltages on the 1S1R devices as a function of the device position on the word-line (column index n). The higher the number of devices per wordline on the Crossbar, the stronger the applied reading voltages degradation within the Crossbar due to IR voltage drop phenomenon. In this context, Figure 8B) provides the 1S1R reading BER evolution with the device position on the word-Adv. Electron. Mater. 2022, 2200323 line. Both 10 MHz and 4 MHz reading frequencies are considered for the analysis. In order to preserve acceptable reading BER on the devices (Figure 7B), the maximal readable number of devices per word-line (n limit ) can be identified. Therefore, a trade-off exists between the 1S1R reading frequency and the maximal readable Crossbar size due to IR voltage drop phenomenon on metal lines. In particular, the lower the 1S1R reading frequency, the lower the optimal reading BER, and so the larger the readable Crossbar arrays due to better tolerance to IR drop phenomenon. First, an adapted computation partitioning into smaller Crossbar arrays can be considered to prevent network performance degradation due to IR voltage drop on large Crossbar structures. Second, adapted network training strategies, including Crossbar IR drop phenomenon, allow preventing network accuracy degradation during inference. [START_REF] Fouda | [END_REF][START_REF] Liu | presented at ICCAD[END_REF][START_REF] Lepri | [END_REF] Third, including the IR voltage drop constraint when designing the circuit architecture, enhances its robustness to the phenomenon. [START_REF] Lepri | [END_REF][30][START_REF] Lepri | presented at IRPS[END_REF] In this context, Figure 9A quantifies the trade-off between crossbar area and system accuracy. A 1S1R reading frequency of 4 MHz is considered for the analysis. The provided area is calculated by adding the respective area contributions from both input-hidden weights and hidden-output weights. First, assuming one driver transistor height per bit-line and wordline, the overall periphery area is calculated for 28 nm highvoltage CMOS for the network topologies of interest. The area contribution of the circuit dedicated to BSNN neuron implementation is not considered for the analysis. Second, assuming a CD min metal width and space between metal lines, the 1S1R crossbar area is quantified for the various network topologies. Increasing linearly the number of neurons on the hidden layer leads to a quadratic increase in the number of synaptic weights and so of the overall crossbar area. Third, the equivalent 1T1R array area is estimated for the various network topologies.

In this context, the peripherals' footprint is demonstrated to be negligible compared to the actual crossbar area. In addition, one decade of area reduction is demonstrated for 1S1R-based Crossbar in comparison with 1T1R arrays. Moreover, network topologies with a large number of neurons on the hidden layer present degraded accuracy due to IR voltage drop phenomenon in the arrays. Altogether, the amount of neurons on the hidden layer of 1024 is observed to optimize the area-accuracy tradeoff, promising high tolerance to BER while strongly reducing the overall system footprint. Given this topology, we propose to Adv. Electron. Mater. 2022, 2200323 Figure 8. A) IR voltage drop on metal lines influence on the applied reading voltages on 1S1R devices as a function of their position on a Crossbar word-line (n value). 28 nm technological node resistive rules in metal lines are used for the analysis. The via resistivity is supposed to be negligible compared to metal line resistivity. The optimal V read voltages identified in the previous section are applied in the word-line extremity. B) 1S1R reading BER evolution with the device position on the Crossbar word-line. Both 10 MHz and 4 MHz reading frequencies are considered for the analysis. The lower the 1S1R reading frequency, the lower the optimal reading BER and so the larger the readable Crossbar arrays due to better tolerance to IR drop phenomenon. Crossbar array sizes required to accommodate the input-hidden synaptic weights are indicated. The higher the amount of neurons on the hidden layer, the higher the required amount of columns on the Crossbar, and so the more the IR voltage drop phenomenon is important within the array. The Crossbar arrays required for hidden-output weights accommodation are not susceptible to lead to a significant IR voltage drop phenomenon, given a very low amount of columns on the array that corresponds to the amount of neurons on the output layer. Figure 9. A) Evaluation of the trade-off between crossbar area and system accuracy for the 1S1R BER of interest, considering a 1S1R reading frequency of 4 MHz. Assuming one driver transistor height per bit-line and word-line, the overall periphery area is calculated for 28 nm high-voltage CMOS. The 1S1R crossbar area is provided (CD min metal width and space between metal lines). The equivalent 1T1R area is estimated, demonstrating on order of magnitude improvement on system area for 1S1R crossbar architecture. An amount of neurons on the hidden layer of 1024 is observed to optimize the area-accuracy trade-off, showing high tolerance to BER while strongly reducing the overall system footprint. B) Optimal BSNN activations sparsity identification, preventing any network accuracy degradation on the MNIST recognition task. C) BSNN electrical consumption per read bit estimation. One order of magnitude improvement in comparison with standard BNN is estimated.

optimize the overall BSNN electrical consumption by tuning the network activations sparsity. In particular, the activations sparsity is defined by the average number of spikes received by a synapse per image classification. Therefore, the sparsity corresponds here to the average number of weights that are read per inference. Thus, the lower the sparsity, the lower the number of weight readings per inference, and so the lower the overall electrical consumption. In our fully connected network, there are two contributions to the overall sparsity: the inputhidden synapses contribution and the hidden-output synapses contribution. Since the number of input-hidden weights is almost two orders of magnitude higher than the number of hidden-output weights in the given topology, the contribution of the hidden-output synapses to the overall sparsity is negligible. Particularly, the sparsity in the input-hidden synapses is the number of spikes fired by neurons in the input layer (previous layer in the network). In this context, the sparsity of our network is determined by the sparsity of the input layer, meaning the sparsity of the input image. To engineer the latter, the threshold that is used for the image binarization is modified. Figure 9B identifies the minimal BSNN sparsity, preventing the network performance degradation. The electrical consumption per reading bit evolution with the network sparsity is presented in Figure 9C). Approximately 1.4pJ electrical consumption per reading bit is estimated for this work. In comparison with standard Binarized Neural Networks (BNN), where all the weights are read at each inference, about one order of magnitude of energy consumption improvement is demonstrated for BSNN for an equivalent network accuracy on the task of interest. [18] Finally, general guidelines for BSNN figures of merit (area, operating frequency, and accuracy) optimization are provided in Figure 10. Assuming an adapted computation partitioning into smaller Crossbar arrays, IR voltage drop issues are not considered here. When the reading frequency is high, the erratic OTS switching probability (Section 3.

Conclusion

1S1R capabilities for synaptic weight storage for Binarized Spiking Neural Network high-frequency inference hardware implementation are demonstrated. By crossing statistical experimental data on memory arrays with Monte Carlo simulations, the OTS switching dynamics are elucidated. Stochastic formation of local conductive dots in the selector leading to the formation of a conductive path allows catching the OTS switching probability distribution when repeated sub-threshold pulses are applied. Based on this analysis, 1S1R reading conditions are optimized for low reading BER during high-frequency inference. Focusing on the MNIST handwritten digit recognition task, general guidelines for system footprint and electrical consumption reduction and inference frequency and accuracy maximization network are provided. Overall, the 1S1R array of interest is demonstrated to perfectly sustain a 1 MHz inference with 97% accuracy, with an estimated circuit area lower than 0.01 mm 2 . This opens the path to 1S1R exploitation for realtime image inference tasks.

  Adv. Electron. Mater. 2022, 2200323

Figure 2 .

 2 Figure 2. A) Technological details. SEM top view of the 1S1R devices integrated into the BEOL, together with a TEM cross-section of the 1S1R devices studied in this work. B) 1S1R typical current-voltage characteristics after the firing process. The 1S1R switching voltages are function of the OxRAM resistive state. If OxRAM is in low resistive state (LRS) (rep. high resistive state (HRS)), V th-LRS (resp. V th-HRS > V th-LRS ) is required for 1S1R switching operation. In this context, the application of a reading voltage V read between V th-HRS and V th-LRS allows to identify the OxRAM resistive state. C) 1S1R programming endurance capabilities. Up to 10 4 programming cycles are demonstrated without errors. D) Experimental 1S1R reading disturb characteristics at V read = 2.25V. A 100 µs relaxation time between pulses is allowed to the device. No error is observed up to 10 9 cycles. E) Experimental 1S1R high resistive state current reading disturb characteristics at V read = 2.25 V. A short 1µs relaxation time between pulses is allowed to the device, which increases the probability of OTS erratic switching with respect to longer relaxation times. In this case, the OTS switches after 10 7 reading cycles, leading to the appearance of a reading fail.

  Figure4Apresents the Cycle 1st-fail dependence with relaxation times and applied reading voltages on OTS devices. The boxplots represent the (25%; 75%) percentiles. Increasing V read closer to the OTS intrinsic switching voltage lowers the median of Cycle 1st-fail . In addition, as the relaxation time between reading cycles increases, Cycle 1st-fail median is shifted to higher values for a given V read . This suggests a more robust OTS OFF state when a longer relaxation time between applied pulses is achieved. Moreover, Figure4Bpresents experimental OTS Cycle 1st-fail distributions for two distinct applied voltages for a corresponding t relax of ≈1 µs. One unique device is considered per distribution. Each point on the distribution corresponds to a unique sub-threshold reading disturb test until 10 8 reading cycles. In this context, the

Figure 3 .

 3 Figure 3. A) Conceptual description of the OTS switching operation considered in this work.The OTS switching operation is driven by an applied electrical field at the device terminals, which induces the nucleation of randomly spaced conductive dots on the bulk of the OTS.[START_REF] Karpov | [END_REF] If the conductive defects are aligned, a conductive filament appears on the device and leads to the OTS switching. B) OTS Monte Carlo reading endurance model deployment. The nucleation of a conductive defect on the constriction layer is driven by a certain probability (P n ). The disappearance of a conductive dot on the constriction zone is driven by a certain recovery time (τ d ). The alignment of several defects implies the apparition of a conductive filament on the constriction zone and leads to a reading fail.

Figure 4 .

 4 Figure 4. A) First measured cycle when the OTS switches under repeated applied pulses under threshold (called Cycle 1st-fail ) as function as applied voltage. Two various relaxation times between reading cycles are considered. The higher V read , the lower the median and maximal Cycle 1st-fail . For a given V read , the larger the relaxation time between reading pulses and the higher appears Cycle 1st-fail . B) OTS Cycle 1st-fail experimental and simulated distributions for V read = 2.25 V and V read = 2.4 V. A fixed relaxation time of 1 µs between reading cycles is considered. High dispersion Cycle 1st-fail experimental distributions are satisfactorily captured by adjusting the OTS global probability to switch on the Monte Carlo model.

  Adv. Electron. Mater. 2022, 2200323

Figure 6 .

 6 Figure 6. A) OTS reading voltages and reading frequency influence on the resulting OTS switching probability after 10 6 sub-threshold reading cycles. The BER is represented by the symbol size. B) Identification of the optimal reading voltages for 1S1R reading BER minimization for the devices of interest. The reading frequency is fixed at 10 MHz. A 1S1R read window margin of 700 mV is considered. C) Identification of the optimal reading voltages for 1S1R reading BER minimization for the devices of interest. The reading frequency is fixed at 4 MHz. A 1S1R read window margin of 700 mV is considered.

Figure 7 .

 7 Figure 7. A) Binarized Spiking Neural Network (BSNN) considered in this work. Fully connected neural network with one hidden layer of X neurons for MNIST handwritten recognition. X values belong to [512, 1024, 2048, 4096]. B) BSNN network accuracy evolution with the number of inference timesteps, focusing on MNIST handwritten recognition task. No significant network performance improvement is observed using several simulation timesteps. C) Maximal attainable BSNN accuracy for a traditional training method as a function of the memory BER for the topologies of interest, considering one inference timestep. The 1S1R BERs are demonstrated to induce an important degradation in accuracy. D) An adapted BSNN training strategy is considered, where the memory BERs are included during the training process. 1024 neurons on the hidden layer are considered for this analysis. Again, one inference timestep is considered. The network accuracy strongly improves for the bit error of interest. However, the network becomes less performant for smaller BERs than the one used for training. E) Inference accuracy evolution with BER. Each point is obtained with a training strategy adapted to the considered BER. Again, one inference timestep is considered. No considerable accuracy degradation is observed for the 1S1R devices of interest (BER ≈10 -1 at 10 MHz read frequency and BER ≈4.5 × 10 -2 at 4 MHz read frequency).

  3) increases due to the creation of a conductive path in the OTS, leading to degraded 1S1R BER (Section 3.4) what is detrimental for the neural network accuracy (Section 3.5). Adapted training and a larger circuit area (by increasing the number of neurons in the hidden layer) make the circuit more robust to failure.

Adv. Electron. Mater. 2022, 2200323

© 2022 Wiley-VCH GmbH

© 2022 Wiley-VCH GmbH 2200323 (8 of 11) www.advelectronicmat.de

Acknowledgements

This work was partially funded by the European project ANDANTE and StorAIge, the French IPCEI program, as well as MIAI @ Grenoble Alpes (ANR-19-P3IA-0003) program.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflict of Interest

The authors declare no conflict of interest.

Adv. Electron. Mater. 2022, 2200323