Jose Aguado

Antonio Huerta

Francisco Chinesta

Elías Cueto

José V Aguado

Real-time monitoring of thermal processes by reduced order modeling

Keywords: Real-time, Heat Transfer, Monitoring, Model Reduction, Proper Generalized Decomposition, Harmonic Analysis

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

Introduction

Many thermal manufacturing processes require monitoring temperature and, moreover, being able to determine from these measurements if the process is running as designed. Typically these measurements are provided by some sensors (e.g. thermocouples) strategically placed. Consider, for instance, an industrial thermal process involving an external excitation, a heat source, moving on the surface of the considered part. This is typical of composite manufacturing processes. Today it is still a challenge to post-process these temperature measurements to monitor in real-time the correct evolution of this process, or to control the heat source, or to identify defects, material properties or power oscillations, etc. This paper proposes a methodology that circumvents some inherent difficulties and provides a fast strategy to determine temperature at a desired monitoring point. More precisely, the approach presented here allows to compute in real-time temperature at a point of interest given an arbitrary transient heat source traveling along a Neumann or Robin boundary.

The strategy proposed revolves on the reciprocity principle [START_REF] Love | A treatise on the Mathematical Theory of Elasticity[END_REF] extensively used in mechanics, dynamics, electromagnetic or wave scattering problems. Note that, in general, this principle is not applicable for the heat equation because of the lack of symmetry introduced by the first time derivative. This paper proposes to adopt the reciprocity principle also to the heat equation. For this, the heat equation is recast in the frequency domain. Details and proofs on this novel approach are presented here.

Another novelty is to apply the Proper Generalized Decomposition method (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. Part II: transient simulation using space-time separated representations[END_REF][START_REF] Chinesta | A short review on model order reduction based on Proper Generalized Decomposition[END_REF] to the frequency-domain heat equation. Combined with the reciprocity principle, this then allows real-time evaluation of the temperature at a specific point of a thermal system in the online stage. More specifically, the PGD computes for all range of frequencies a generalized transfer function solution of the frequency-domain heat equation.

Thus a strategy based on an offline and online phase is designed. In the offline phase the previously cited generalized transfer function is computed at the desired monitoring point. Then, given the external heat source excitation, a convolution is used to determine online and in real-time the temperature at the desired monitoring point for any instant. As it will be shown, the online approximation is so fast that it can be used for control purposes in real-time and on deployed devices. This opens vast possibilities related to real-time simulationbased monitoring and control [START_REF] Chinesta | PGD-based computational vademecum for efficient design, optimization and control[END_REF].

Representations in the frequency domain are appealing for analyzing responses of structures subjected to dynamic excitations. It is a powerful approach to study the response of structural systems when initial conditions can be neglected; that is, far enough from the initial transient response. This framework, which was integrated into the finite element framework from the very beginning, has been, and still is, extensively used. It is well described in most of the textbooks, some of them linked to many generations of scientists and engineers [START_REF] Bathe | Finite element procedures, Civil engineering series[END_REF][START_REF] Clough | Civil engineering series[END_REF][START_REF] Thomas | The finite element method: linear static and dynamic finite element analysis[END_REF][START_REF] Zienkiewicz | The finite element method[END_REF]. This approach is still today an active research area because several challenges are still open. For instance, separated representations in the frequency domain were considered in [START_REF] Barbarulo | Proper Generalized Decomposition applied to linear acoustic: A new tool for broad band calculation[END_REF][START_REF] Dorival | A substructured version of the variational theory of complex rays dedicated to the calculation of assemblies with dissipative joints in the medium-frequency range[END_REF][START_REF] Kovalevsky | The Fourier version of the Variational Theory of Complex Rays for medium-frequency acoustics[END_REF][START_REF] Kovalevsky | The variational theory of complex rays for three-dimensional Helmholtz problems[END_REF][START_REF] Ladevèze | The variational theory of complex rays for the calculation of medium-frequency vibrations[END_REF][START_REF] Ladevèze | Calculation of medium-frequency vibrations over a wide frequency range[END_REF][START_REF] Riou | The multiscale VTCR approach applied to acoustics problems[END_REF][START_REF] Rouch | The variational theory of complex rays: A predictive tool for medium-frequency vibrations[END_REF][START_REF] Rouch | Vibrational analysis of structures with stochastic interfaces in the medium-frequency range: Experimental validation on a touch screen[END_REF] for the so-called: variational theory of complex rays. Obviously, there have been many attempts considering such descriptions within the model reduction framework; the interested reader can refer to [START_REF] Hetmaniuk | Review and assessment of interpolatory model order reduction methods for frequency response structural dynamics and acoustics problems[END_REF] and the references therein. This however is not the aim of the present paper, thus, no exhaustive state-of-the-art on this topic is presented.

Finally, it is important to note that in spite of the large amount of scientific contributions using a frequency domain description in solid dynamics, this approach is not standard, to the author's knowledge, for thermal models subjected to dynamical forced thermal loads. Certainly, time domain approximations of thermal models are both efficient and robust and, moreover, model order reduction is successfully applied in this setting [START_REF] Bia Lecki | Proper orthogonal decomposition and modal analysis for acceleration of transient FEM thermal analysis[END_REF][START_REF] Chinesta | First steps towards an advanced simulation of composites manufacturing by automated tape placement[END_REF][START_REF] Ch | Proper Generalized Decomposition based dynamic data-driven control of thermal processes[END_REF][START_REF] Girault | Estimation of time-varying heat sources through inversion of a low order model built with the modal identification method from in-situ temperature measurements[END_REF][START_REF] Maday | The reduced basis element method: application to a thermal fin problem[END_REF][START_REF] Park | The use of the Karhunen-Loève decomposition for the modeling of distributed parameter systems[END_REF][START_REF] Prulière | On the solution of the heat equation in very thin tapes[END_REF][START_REF] Videcoq | Heat source identification and on-line temperature control by a branch eigenmodes reduced model[END_REF]. Whereas frequency domain approaches for thermal studies are scarce [START_REF] Feng | Real-time solution of heat conduction in a finite slab for inverse analysis[END_REF][START_REF] Monteyne | Heat transfer in a borehole heat exchanger: Frequency domain modeling[END_REF][START_REF] Pinnau | Model reduction techniques for frequency averaging in radiative heat transfer[END_REF][START_REF] Rylatt | Time and frequency domain investigation of the heat transfer to a synthetic air jet[END_REF][START_REF] Sierociuk | Modelling heat transfer in heterogeneous media using fractional calculus[END_REF].

2 Monitoring temperature at a surface point This section analyzes the representation of temperature at an arbitrary point of the boundary under an arbitrary transient external excitation, for instance, a traveling external heat source such as a moving laser. More specifically, the representation of the solution at the point of interest is first studied by means of a Green's function in the space-time domain. Then the applicability of the reciprocity principle for the transient heat equation with a forced excitation is discussed both in time and frequency domain. The use of reciprocity (only valid in the frequency domain) allows determining a representation of temperature at the desired point for any arbitrary external heat source.

Model problem in the space-time domain

Formally, the problem under consideration can be described as follows. Given a time interval I :=]0, T [ (where T can be chosen arbitrarily large) and a body Ω ⊂ R d , d ≤ 3, whose boundary ∂Ω is partitioned into Dirichlet, Γ D , and Robin/Neumann, Γ N , frontiers such that ∂Ω = Γ D ∪ Γ N and Γ D ∩ Γ N = ∅. The temperature evolution u(x, t), for x ∈ Ω and t ∈ I, is described by the transient heat equation:

         ρc p ∂ t u -∇ • K∇u = 0 in Ω × I, u = u D on Γ D × I, n • K∇u = -(u -u ext ) + q on Γ N × I, u = u 0 on Ω × {0}, (1)
where ρ is density (kg/m 3 ), c p is specific heat capacity (J/(kg K)), K is the thermal conductivity matrix (W/(m K)) and is the heat transfer coefficient (W/(m 2 K)), u ext is the external temperature (K), n is the exterior unit normal to Γ N (dimensionless) and q = q(x, t), for (x, t) ∈ Γ N × I, is the inflow forcing excitation (W/m 2 ). The international system of units of measurement is also employed for length (m) and time (s).

As noted in the introduction, q is typically the heat flux imposed by a laser. The objective in what follows is to determine a (fast) computable representation of the temperature at an arbitrary boundary point x 0 ∈ Γ N and at any instant t 0 , with 0 < t 0 < T .

Since the problem is linear, (1) can be further simplified. For instance, thermal diffusivity (thermal conductivity divided by density and specific heat capacity) can be considered the only material constant in the partial differential equation. Moreover, the increment of temperature with respect to the external one, i.e. (u -u ext ), can be defined as the unknown of the problem (that is, in practice, impose u ext = 0).

Finally, for the clarity of the presentation, the model problem studied is further simplified. However, these simplifications (considering unitary values of the coefficients) do not compromise the validity of the following developments.

More precisely, the following assumptions are used to define the model problem: homogeneous Dirichlet boundary conditions, canonical dimensionless form with an isotropic homogeneous material, and no convective heat exchanges. Note however, that these simplifications, are only done to simplify the presentation and do not hinder the application of the proposed methodology to real problems described by (1) as it will be shown in Section 5. Under these circumstances, problem (1) becomes

         ∂ t u -∇ 2 u = 0 in Ω × I, u = 0 on Γ D × I, n • ∇u = q on Γ N × I, u = u 0 on Ω × {0}.
(2)

Green's function and reciprocity in the space-time domain

The objective is to obtain the solution at an arbitrary point and time, (x 0 , t 0 ). Ideally the desired value u(x 0 , t 0 ) could be readily evaluated if the adjoint Green's function were known. The adjoint Green's function is the solution of

         ∂ t G + ∇ 2 G = 0 in Ω×]0, t 0 [, G = 0 on Γ D ×]0, t 0 [, n • ∇G = δ(x -x 0 )δ(t -t 0 ) on Γ N ×]0, t 0 [, G = 0 on Ω × [t 0 , T [, (3)
where the notation of G(x, t; x 0 , t 0 ) clearly identifies the parametric dependence on (x 0 , t 0 ). The representation of the desired temperature is then

u(x 0 , t 0 ) = Γ N t0 0 G(x, t; x 0 , t 0 ) q(x, t)dtdΓ + Ω u 0 (x) G(x, 0; x 0 , t 0 )dΩ.
See Appendix A for a detailed presentation. However, in general, the computation of the Green's function is by no means a trivial task, for instance when confronted to an arbitrary domain or inhomogeneous material properties. Consequently, this approach is not used in practice. An alternative, extensively employed in dynamics, is to use the reciprocity property [START_REF] Arias | Rayleigh wave correction for the BEM analysis of two-dimensional elastodynamic problems in a half-space[END_REF][START_REF] Love | A treatise on the Mathematical Theory of Elasticity[END_REF][START_REF] Vega | Dynamic response of underpasses for high-speed train lines[END_REF]. However, it is also well-known that it is not applicable to the heat equation because the operator is not self-adjoint. In order to recall this, the variational problem equivalent to (2) is presented: find u ∈ S such that

B(u, v) = L(q; v) ∀v ∈ V, (4)
with the appropriate spaces introducing the required regularity in space and time [START_REF] Parés | Bounds of functional outputs for parabolic problems. I. Exact bounds of the discontinuous Galerkin time discretization[END_REF][START_REF] Temam | Navier-Stokes equations. theory and numerical analysis[END_REF] V

:= v : v(•, t) ∈ H 1 (Ω), v(x, •) ∈ L 2 (I), v = 0 on Γ D × I ∩ v : v(•, t) ∈ H -1 (Ω), v(x, •) ∈ H 1 (I) , S := v : v ∈ V, v(x, 0) = u 0 , and
B(u, v) = Ω I v ∂ t u dtdΩ + Ω I ∇u • ∇v dtdΩ, L(q; v) = Γ N I qv dtdΓ.
Given two different exitations q 1 and q 2 , the corresponding solutions of Eq. ( 4) can be denoted u 1 and u 2 , respectively. Since both u 1 and u 2 belong to S the following expressions are also verified:

B(u 1 , u 2 ) = L(q 1 ; u 2 ) and B(u 2 , u 1 ) = L(q 2 ; u 1 ). (5)
However, the bilinear form is non-symmetric because of the time derivative, that is B(u, v) = B(v, u). Thus, subtracting both expressions in Eq. ( 5) attests that standard reciprocity is not satisfied because the left hand side terms do not cancel out. In conlusion, L(q 1 ; u 2 ) = L(q 2 ; u 1 ).

Therefore, in the space-time domain neither a Green's function approach or a reciprocity property can be used in practice to determine temperature at a given point and instance of the surface, say (x 0 , t 0 ) ∈ Γ N × I.

Space-frequency problem for an arbitrary excitation

Another alternative for studying this forced excitation problem is to consider harmonic analysis. In order to work in the frequency domain the Fourier transform and its inverse are used, namely

v(x, ω) = F[v] = +∞ -∞ v(x, t) e -iωt dt (6a) and v(x, t) = F -1 [v] = 1 2π +∞ -∞ v(x, ω) e iωt dω. (6b
     iωû -∇ 2 û = 0 in Ω, û = 0 on Γ D , n • ∇û = q on Γ N . (7)
Remark 2 (Long-term forced solution). Note that the harmonic solution is only concerned with the long-term forced solution and consequently it does not depend on the initial condition. Standard approaches should be used for evaluating the transient regime, which, in practice, decays rapidly to the obtained long-term solution.

The variational form associated to the strong form problem described in [START_REF] Bathe | Finite element procedures, Civil engineering series[END_REF] reads:

find û ∈ H 1 Γ D := v ∈ H 1 (Ω) : v = 0 on Γ D such that ∇û, ∇v + iω û, v = q, v ∀v ∈ H 1 Γ D , (8) where û
, v = Ω û v * dΩ , ∇û, ∇v = Ω ∇û • ∇v * dΩ and û, v = Γ N û v * dΓ (9)
denote, respectively, the L 2 scalar product of functions û and v in Ω and its traces over Γ N . Note also, that v * indicates the complex conjugate of v, since both û and v are, in general, in C.

It is important to observe that Eq. ( 8) is non-Hermitian but it is symmetric, the later property proves sufficient and also crucial for reciprocity.

Arbitrary excitation implies solving two problems with real excitation

Given any arbitrary excitation, q(x, t), it can always be decomposed in the sum of an even and odd function, namely q(x, t) = q e (x, t) + q o (x, t) = 1 2 q(x, t) + q(x, -t) + 1 2 q(x, t) -q(x, -t) .

Recalling the Fourier transform properties: F[q] = q = F[q e ] + F[q o ] = qe + iq o with qe (x, ω) ∈ R and qo (x, ω) ∈ R. Moreover, it is important to notice that, the decomposition of q(x, t) in the sum of an even and odd excitation produces two identical problems with a real excitation whose solutions are the real and the imaginary part of the solution of [START_REF] Bathe | Finite element procedures, Civil engineering series[END_REF]. More precisely, because of the linearity of [START_REF] Bathe | Finite element procedures, Civil engineering series[END_REF] and Remark 3, the solution of û of (7) can also be decomposed as û = ûe + iû o , which in turn are the solutions of

     iωû e -∇ 2 ûe = 0 in Ω, ûe = 0 on Γ D , n • ∇û e = qe on Γ N , and
     iωû o -∇ 2 ûo = 0 in Ω, ûo = 0 on Γ D , n • ∇û o = qo on Γ N . (10)
Then, in the time domain the original solution can be recovered by means of the inverse Fourier transform, namely

u = F -1 [û] = F -1 [û e ] + F -1 [iû o ] = u e + u o .
In summary, it is possible to find the solution of (7) -or (8)-for any arbitrary excitation q solving twice the same problem but with different real excitations, one corresponding to qe (x, ω) and the other to qo (x, ω).

Remark 3 (Imaginary excitation). Suppose u is solution of the following problem

     iωu -∇ 2 u = 0 in Ω, u = 0 on Γ D , n • ∇u = q on Γ N . Then, v = iu is solution of      iωv -∇ 2 v = 0 in Ω, v = 0 on Γ D , n • ∇v = iq on Γ N .
Hint: replace u = -iv in the first problem and obtain the second one. Obviously, this follows directly from linearity but it is explicitly recalled for didactic purposes.

Reciprocity in space-frequency holds for a real excitation

Given two real harmonic excitations q1 and q2 , the corresponding solutions of (8) are denoted by û1 and û2 , respectively. Then, since both solutions belong to H 1 Γ D , the following expressions hold:

∇û 1 , ∇û 2 + iω û1 , û2 = q1 , û2 , (11a)
∇û 2 , ∇û 1 + iω û2 , û1 = q2 , û1 . (11b)
Substracting (11b) from (11a) gives

∇û 1 , ∇û 2 -∇û 2 , ∇û 1 + iω û1 , û2 -û2 , û1 = q1 , û2 -q2 , û1 , (12)
which clearly shows that reciprocity is satisfied in the frequency domain, namely

q1 , û2 = q2 , û1 , (13)
if the following conditions hold

∇û 1 , ∇û 2 = ∇û 2 , ∇û 1 and û1 , û2 = û2 , û1 , (14)
which is precisely the case when both q1 and q2 are real. See Appendix B for a detailed proof of Eqs. [START_REF] Dorival | A substructured version of the variational theory of complex rays dedicated to the calculation of assemblies with dissipative joints in the medium-frequency range[END_REF].

Remark 4 (Reciprocity with convective heat flux). In the general case when convective heat fluxes are considered, see [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF], reciprocity also holds. Linearity is exploited solving for (u-u ext ) instead of the original temperature u. Accordingly, equations [START_REF] Chinesta | A short review on model order reduction based on Proper Generalized Decomposition[END_REF] are modified as follows:

∇û 1 , ∇û 2 + iω û1 , û2 + û1 , û2 = q1 , û2 , ∇û 2 , ∇û 1 + iω û2 , û1 + û2 , û1 = q2 , û1 .
Subtracting again these two previous equations the same equation ( 12) is obtained (recall the symmetry of the L 2 scalar products of functions). Thus reciprocity also holds for convective heat flux following the proof in Appendix B.

Using reciprocity to monitor temperature

Recall that the final objective is to monitor temperature at a given point for an arbitrary external excitation, i.e. evaluate u(x 0 , t 0 ) for (x 0 , t 0 ) ∈ Γ N × I.

For this purpose, it is necessary to determine û(x 0 , ω) for an arbitrary external excitation, q, in the space-frequency domain. The conclusion of Section 2.4 is that two problems, which are shown in strong form by [START_REF] Chinesta | PGD-based computational vademecum for efficient design, optimization and control[END_REF], with real excitations qe and qo , such that q = qe + iq o , must be solved to find û(x 0 , ω) = ûe (x 0 , ω) + iû o (x 0 , ω). Suppose, that for each frequency, ω ∈ R, one could also determine the corresponding solution ĥ(x, ω; x 0 ) under a Dirac flux imposed at the monitoring point x 0 , δ(x -x 0 ), namely

       iω ĥ -∇ 2 ĥ = 0 in Ω, ĥ = 0 on Γ D , n • ∇ ĥ = δ(x -x 0 ) on Γ N . (15)
Since all the excitations are real, reciprocity, see [START_REF] Dorival | A substructured version of the variational theory of complex rays dedicated to the calculation of assemblies with dissipative joints in the medium-frequency range[END_REF], holds and consequently δ(x -x 0 ), ûe = qe , ĥ and δ(x -x 0 ), ûo = qo , ĥ , that is, ûe (x 0 , ω) = ĥ(•, ω; x 0 ), qe (•, ω) and ûo (x 0 , ω) = ĥ(•, ω; x 0 ), qo (•, ω) ,

which implies that û(x 0 , ω) = Γ N ĥ(x, ω; x 0 ) q(x, ω) dΓ. (16)
See Remark 5 for a detailed deduction of [START_REF] Ch | Proper Generalized Decomposition based dynamic data-driven control of thermal processes[END_REF].

Once the temperature is monitored at the desired point in the space-frequency domain, the inverse Fourier transform is employed to obtain the desired final representation of temperature in the space-time domain. That is, the objective is to determine

u(x 0 , t) = F -1 û(x 0 , ω) = Γ N F -1 ĥ(x, ω; x 0 ) q(x, ω) dΓ.
This expression can be further simplified using the convolution theorem [START_REF] Duffy | Green's functions with applications[END_REF][START_REF] Kammler | A first course in Fourier analysis[END_REF], recall

F -1 ĥ(x, ω; x 0 ) q(x, ω) = t 0 h(x, τ ; x 0 ) q(x, t -τ )dτ, where h(x, τ ; x 0 ) = F -1 [ ĥ] = 1 2π +∞ -∞ ĥ(x, ω; x 0 ) e iωt dω. (17)
Thus, the representation of the temperature at the desired point

x 0 ∈ Γ N is u(x 0 , t) = Γ N t 0 1 2π +∞ -∞ ĥ(x, ω; x 0 ) e iωt dω q(x, t -τ ) dτ dΓ,
which can be written in a more compact form and for any instance

t 0 ∈ I =]0, T [ as u(x 0 , t 0 ) = t0 0 h(x, τ ; x 0 ), q(x, t 0 -τ ) dτ. (18)
Recall that q ∈ R. This is a compact and useful expression, it only requires knowledge of the external imposed excitation q up to the desired monitored instant t 0 (causality). Moreover, and this is a major point and advantage, [START_REF] Hetmaniuk | Review and assessment of interpolatory model order reduction methods for frequency response structural dynamics and acoustics problems[END_REF] does not reflect the decomposition of the excitation in even and odd contributions. Thus, if the transfer function, h(x, τ ; x 0 ), is known, (18) can be applied directly for any arbitrary excitation q(x, t).

However this expression presents a major drawback: the inverse Fourier transform of ĥ, solution to problem [START_REF] Feng | Real-time solution of heat conduction in a finite slab for inverse analysis[END_REF], must be known. This implies solving [START_REF] Feng | Real-time solution of heat conduction in a finite slab for inverse analysis[END_REF] for each frequency ω in the range needed by the arbitrary excitation. Thus, in general, the representation (18) cannot be used in practical problems.

The next section circumvents this drawback and proposes a methodology to obtain an expression for the generalized transfer function ĥ, solution of [START_REF] Feng | Real-time solution of heat conduction in a finite slab for inverse analysis[END_REF], for the all range of realistic frequencies. This expression can then be substituted in [START_REF] Hetmaniuk | Review and assessment of interpolatory model order reduction methods for frequency response structural dynamics and acoustics problems[END_REF] to determine the desired temperature. Moreover, (18) can be evaluated in real-time.

Remark 5 (Reconstruction of solution). To determine [START_REF] Ch | Proper Generalized Decomposition based dynamic data-driven control of thermal processes[END_REF] it is important to recall that qe and qo are real. Thus, the scalar products on Γ N , û(x 0 , ω) = ûe (x 0 , ω)+iû o (x 0 , ω) = ĥ(•, ω; x 0 ), qe (•, ω) +i ĥ(•, ω; x 0 ), qo (•, ω) , can be rewritten as

û(x 0 , ω) = Γ N ĥ(x, ω; x 0 ) qe (x, ω) dΓ + i Γ N ĥ(x, ω; x 0 ) qo (x, ω) dΓ = Γ N ĥ(x, ω; x 0 ) qe (x, ω) + iq o (x, ω) dΓ = Γ N ĥ(x, ω; x 0 ) q(x, ω) dΓ,
and do not present any complex conjugate. In fact, it is important to note that û(x 0 , ω) = Γ N ĥ(x, ω; x 0 ) q(x, ω) dΓ = ĥ(•, ω; x 0 ), q(•, ω) * = ĥ(•, ω; x 0 ), q(•, ω) .

Remark 6 (Inverse Fourier Transform of the generalized transfer function).

The inverse fourier transform of the generalized transfer function, ĥ(x, ω; x 0 ), see ( 17) is computed using the FFT algorithm. Only periodic signals with a finite number of harmonics can be exactly represented with the discrete Fourier transform, and thus non-periodic signals involving a continuous spectrum of frequencies can only be approximated. The range of frequencies included in the generalized transfer function must be chosen accordingly.

Remark 7 (Convective heat flux). In the general case described in (1), when convective fluxes are present, the computation of the generalized transfer function is modified. As noted in Remark 4, reciprocity also holds when convective fluxes are considered. Thus, the generalized transfer function problem originally described by [START_REF] Feng | Real-time solution of heat conduction in a finite slab for inverse analysis[END_REF] must have on the Neumann boundary the convective heat flux, namely

      
iω ĥ -∇ 2 ĥ = 0 in Ω, ĥ = 0 on Γ D , n • ∇ ĥ = -ĥ + δ(x -x 0 ) on Γ N .

The same development as before can be applied and, consequently, equations ( 16) and (18) also hold.

Computing the generalized transfer function

This section is aimed at computing a generalized transfer function, ĥ(x, ω; x 0 ), for a desired and predefined range of frequencies, I ω . This transfer function needs to be computed only once, and preferably offline. Since it is determined in the frequency domain, its inverse Fourier transform is later evaluated in order to use equation ( 18) as a simple and inexpensive post-process of any given excitation q(x, t).

A major contribution of the PGD approach is to view the frequency, ω, as a new coordinate [START_REF] Chinesta | PGD-based computational vademecum for efficient design, optimization and control[END_REF]. Thus, instead of solving problem [START_REF] Feng | Real-time solution of heat conduction in a finite slab for inverse analysis[END_REF] for each frequency, the objective is to solve, only once, a more general problem with ω as an extra coordinate, namely find ĥ(x, ω; x 0 ) satisfying

       iω ĥ -∇ 2 ĥ = 0 in Ω × I ω , ĥ = 0 on Γ D × I ω , n • ∇ ĥ = δ(x -x 0 ) on Γ N × I ω , (19)
where I ω is the predefined range of variation of ω. There is a price to pay, the increased dimensionality. However, as shown numerically in [START_REF] Chinesta | PGD-based computational vademecum for efficient design, optimization and control[END_REF] and references therein, the complexity of the PGD does not grow exponentially with the space dimension. On the contrary, numerical evidence shows that it scales linearly for a moderate number of dimensions. This is obvious for each resolution of the greedy algorithm. This is precisely what happens in this problem when ω is introduced as an extra coordinate. However, for large number of dimensions (> 100) or in non-symmetric problems the number of terms required for a given precision is sensible to the number of dimensions.

The weak problem equivalent to ( 19) is obtained using a weighted residual argument, namely, find ĥ for all v in the selected appropriate functional space such that

A( ĥ, v) = L(v) (20a)
with

A( ĥ, v) := Iω ∇ ĥ, ∇v dω + Iω iω ĥ, v dω (20b)
L(v) := Iω v * (x 0 , ω)dω. ( 20c
)
Note that formally, the required functional spaces must account for the singularity of the Dirac flux. Nevertheless, in practice, the Dirac delta is mollified and this allows to use the standard finite element functional setup.

The PGD approach assumes that the solution of (20), ĥ(x, ω; x 0 ), can be approximated by a rank-n separable function, ĥn (x, ω; x 0 ), namely, ĥ(x, ω; x 0 ) ≈ ĥn (x, ω;

x 0 ) = n s=1 X s (x) W s (ω) = ĥn-1 (x, ω; x 0 )+X n (x) W n (ω), (21) where, for s = 1, . . . , n, X s ∈ H 1 Γ D and W s ∈ L 2 (I ω ).
Note that all these functions give values in C.

A greedy algorithm [START_REF] Lee | A first course in combinatorial optimization[END_REF] is used to construct this approximation, that is, to determine the unknown functions X s and W s in Eq. [START_REF] Kovalevsky | The Fourier version of the Variational Theory of Complex Rays for medium-frequency acoustics[END_REF]. In other words, a progressive scheme is designed to compute each new term with all the information at hand. The sequence is stopped with an appropriate error estimator [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF][START_REF] Lavedèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF][START_REF] Moitinho De Almeida | A basis for bounding the errors of proper generalised decomposition solutions in solid mechanics[END_REF]. Note that each new term implies the computation of a product of unknown functions, X n and W n . Thus a nonlinear scheme must be designed. It is standard to use a fixed point alternating direction algorithm because it has proven robust in former works [START_REF] Chinesta | PGD-based computational vademecum for efficient design, optimization and control[END_REF][START_REF] Chinesta | A short review on model order reduction based on Proper Generalized Decomposition[END_REF]. To simplify notation each iterate approximating X n and W n is denoted by R and S. Hence, the nonlinear problem to solve for each new term of ĥn (x, ω; x 0 ) is obtained substituting in [START_REF] Kovalevsky | The Fourier version of the Variational Theory of Complex Rays for medium-frequency acoustics[END_REF] in (20a), in order to compute R and S (iterates of X n and W n ) such that

A(R S, v) = L(v) -A( ĥn-1 , v) (22)
with trial functions on the tangent manifold v = vR (x) S(ω) + R(x) vS (ω) ∀v R ∈ H 1 Γ D and ∀v S ∈ L 2 (I ω ). The alternating direction scheme, detailed below, consists in, for instance, updating the space function R from a given S assumed known, and then compute S from the just updated function R. This iteration continues until reaching convergence of both R and S. That is, the two stages for each iteration are:

1. Find R ∈ H 1 Γ D (S assumed known) such that A(R S, vR S) = L(v R S) -A( ĥn-1 , vR S) ∀v R ∈ H 1 Γ D . (23a) 2. Find S ∈ L 2 (I ω ) (R assumed known) such that A(R S, R vS ) = L(R vS ) -A( ĥn-1 , R vS ) ∀v S ∈ L 2 (I ω ). (23b)
Then at convergence, X n and W n are updated by R and S.

Updating the space function

For each new term in the series defined by [START_REF] Kovalevsky | The Fourier version of the Variational Theory of Complex Rays for medium-frequency acoustics[END_REF] and each iteration described by [START_REF] Ladevèze | The variational theory of complex rays for the calculation of medium-frequency vibrations[END_REF], equation (23a) must be solved. Taking advantage of the separated structure of the solution and also of A(•, •), see (20b), equation (23a) can be rewritten as, find R ∈ H 1 Γ D for all vR ∈ H 1 Γ D (S assumed known) such that

α S ∇R, ∇v R + iβ S R, vR = γ S v * R (x 0 ) - n-1 s=1 α S s ∇X s , ∇v R + iβ S s X s , vR , (24)
where the coefficients, which must be computed for each instance of S, are defined as

α S = S, S Iω , β S = ω S, S Iω , γ S = 1, S Iω , α S s = W s , S Iω , and β S s = ω W s , S Iω .
Bear in mind that in previous definitions •, • Iω denotes the L 2 scalar product of complex functions in I ω .

After the corresponding discretization of the spatial domain with a standard combination of piecewise linear shape functions, the system of linear equations induced by [START_REF] Ladevèze | Calculation of medium-frequency vibrations over a wide frequency range[END_REF] presents a conductivity and a mass matrix. These matrices are computed only once because they are constant for each iteration and for each term. Moreover, it is important to note that they are symmetric but non-Hermitian, which will preclude, for instance, Cholesky or conjugate gradient schemes.

Updating the frequency function

Similarly, the second stage, described by equation (23b), can also be rewritten using separated structure of the solution and also of A(•, •) as, find S ∈ L 2 (I ω ) for all vS ∈ L 2 (I ω ) (R assumed known from the previous stage) such that

α R S, vS Iω +iβ R ω S, vS Iω = γ R 1, vS Iω - n-1 s=1 α R s W s , vS Iω +iβ R s ω W s , vS Iω , (25)
where the coefficients, which must be computed for each instance of R, are defined as L 2 products over the spatial domain,

α R = ∇R, ∇R , β R = R, R , γ R = R * (x 0 ), α R s = ∇X s , ∇R , and β R s = X s , R .
The lack of derivatives with respect to ω in the problem ( 19) induces an point-wise algebraic equation for S. Piecewise discontinuous approximations of S will induce uncoupled scalar equations. Whereas continuous approximations over the one-dimensional range of frequencies lead to a symmetric but non-Hermitian matrix on the left-hand-side of [START_REF] Lavedèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF], as in the previous case.

Extension to multi-parametric and inverse problems

The approach presented here has a potentiality that exceeds real-time monitoring of temperature at a given location and can also be used for other thermal studies such as optimization, inverse analysis, non-destructive testing, etc. Here, two simple extensions are presented.

Multi-parametric models

As a simple demonstrator, thermal conductivity is chosen as an extra parameter. The underlying idea, already exploited in [START_REF] Chinesta | PGD-based computational vademecum for efficient design, optimization and control[END_REF][START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the Proper Generalized Decomposition[END_REF], is to solve multi-parametric models capitalizing the advantages of the PGD framework. A multi-parametric model is an extension of the procedure detailed in the previous section. Besides frequency as an extra-coordinate the generalized transfer function can encompass other parameters as extra coordinates. For instance, parameters characterizing the geometry, the constitutive behavior or the boundary conditions could be incorporated. The PGD methodology allows to compute efficiently a multiparametric solution defined in a high-dimensional space (spatial coordinates, frequency and other parameters). Multi-parametric models are of great interest in science and engineering because they make possible real-time simulation, optimization and inverse analysis, as illustrated in [START_REF] Ammar | Parametric solutions involving geometry: a step towards efficient shape optimization[END_REF][START_REF] Chinesta | PGD-based computational vademecum for efficient design, optimization and control[END_REF].

For instance, the thermal example that motivates this work can also be posed as an inverse analysis to find the actual conductivity of a certain material from the temperature provided by a thermocouple placed at location x 0 . Then, independently of the inverse identification method used, it is obvious that fast identification procedures can be envisaged if an approximation of temperature at the monitoring point x 0 for any instance t can be computed in real-time for any conductivity k, i.e. u(x 0 , t, k).

This approach can also be used for non-destructive testing. The monitored temperature being different from the computed one (obtained with the undamaged material parameters) triggers the inverse analysis to determine the "damaged" material parameter, which can be solved readily because a generalized solution for any material parameter is available.

The inverse problem is not solved in detail at this point because it is outside the scope of this work. Nevertheless, here, the generalized solution for any conductivity is provided. Once this solution is known, any standard inverse algorithm could be implemented. The key point is to determine the generalized transfer function, ĥ, for any value of the conductivity k ∈ I k , where I k is the desired range of conductivities. Problem ( 19) is now rewritten with the new parameter, conductivity, considered as an extra-coordinate. Thus the new problem consist in finding ĥ(x, ω, k; x 0 ) that satisfies

       iω ĥ -∇ • (k∇ ĥ) = 0 in Ω × I ω × I k , ĥ = 0 on Γ D × I ω × I k , n • k∇ ĥ = δ(x -x 0 ) on Γ N × I ω × I k .
The PGD approach must now determine an approximation ĥn (x, ω, k; x 0 ) to the solution of the previous problem, namely ĥ(x, ω, k; x 0 ) ≈ ĥn (x, ω, k;

x 0 ) = n s=1 X s (x) W s (ω) K s (k),
where extra separated functions must be determined; more precisely, those directly linked to conductivity, namely K s (k) for s = 1, . . . , n. The same greedy approach described earlier can be applied with now an extra stage in the nonlinear solve to determine each K s (k). Once the multi-parametric solution has been computed, it can be postprocessed in the same way explained in Section 2.6 in order to recover the solution at x 0 and any time t. The objective is that temperature at the monitoring point x 0 for any instance t can be computed in real-time for any conductivity k; that is, equation ( 18) is extended to approximate u(x 0 , t, k), namely

u(x 0 , t, k) = F -1 [û](x 0 , t, k) = t 0 h(x, τ, k; x 0 ), q(x, t 0 -τ ) dτ, (26)
where h = F -1 [ ĥ].

As noted earlier, the crucial point is to determine a reasonable approximation of the generalized transfer function ĥ(x, ω, k; x 0 ). The example presented in Section 5 is extended to demonstrate that such an approximation can be evaluated.

Inverse problem: an amplitude time-modulated calibration of excitation

This inverse problem considers that the amplitude of the power given by a laser varies with time because of uncontrolled power supply and has to be calibrated. Suppose an excitation defined by α(t)q(x, t) where, as assumed in previous sections, q(x, t) is given (and, thus, known) while its amplitude, which varies with time α(t) is not know. Temperature at the monitoring point x 0 for any instance t is obtained following the procedure described in Section 2.4 for this new excitation and equation [START_REF] Hetmaniuk | Review and assessment of interpolatory model order reduction methods for frequency response structural dynamics and acoustics problems[END_REF] becomes

u(x 0 , t) = t 0 α(t 0 -τ ) Γ N h(x, τ ; x 0 ) q(x, t -τ ) dΓdτ.
Then, given a discretization of the unknown function α(t), for instance

α(t) = n fit j=1 α i N i (t), (27)
where N i (t) are known interpolation functions, the coefficients α ∈ R n fit can be determined by a least-squares technique. This implies solving the normal equations Aα = b.

Given the series of instants {t 1 , t 2 , . . . , t m } (with m ≥ n fit ) at which temperature is going to be measured, the matrix of the normal equations is determined once and for all, during the offline phase, as

A = [a ij ] = m r=1 ψ i (t r )ψ j (t r ) with ψ i (t r ) = tr 0 N i (t r -τ ) Γ N h(x, τ ; x 0 ) q(x, t r -τ ) dΓdτ.
Then, the measured values of temperature, u meas (x 0 , t r ), at x 0 and for the series of instants {t 1 , t 2 , . . . , t m } allow to compute, in the online phase, the independent term boundary of a cylinder, is proposed. The outer boundary is also subjected to heat convection while the other boundaries are adiabatic. Figure 1 depicts the problem statement. The initial boundary value problem formerly described as

b = [b i ] = m r=1 ψ i (t r )u meas (x 0 , t r ) .
         ρc p ∂ t u -∇ • k∇u = 0 in Ω × I, n • k∇u = -(u -u ext ) + q on Γ Out × I, n • k∇u = 0 on ∂Ω/Γ Out × I, u = u 0 on Ω × {0}, (28)
where ρ = 1kg/m 3 is density, c p = 1J/(kg K) is specific heat capacity, k = 1 W/(m K) is isotropic thermal conductivity and = 1W/(m 2 K) is the heat transfer coefficient, u ext = 298K is the external temperature, Γ Out is the outer boundary where the laser impacts and with a radius of 1.0m,

q(ξ, t) = 500 exp -50(2ξ -πt) 2 W/m 2 (29)
is the inflow forcing excitation, and ξ is the local tangent coordinate along Γ Out . The temperature is measured (i.e. the point where temperature is monitored) at the middle point of the inner boundary. Finally, the thickness of the single-ply is 0.05m. The range of frequencies considered is f ∈ [-250, 250] Hz or, in terms of the angular frequency used in all previous sections, ω ∈ I ω := [-500π, 500π]. Such a large interval of frequencies has been chosen in order to be able to make a fair comparison with the reference solution. The time-step of a signal and the maximum frequency that can be computed from it are related by the Nyquist-Shannon theorem [START_REF] Shannon | Communication in the presence of noise[END_REF] as follows:

f max = 1/2∆t. ( 30
)
Figure 2: Single-ply composite cylinder: convergence of the generalized frequency transfer function.

The FE reference solution is obtained with a standard time-marching Crank-Nicolson scheme whose time-step is chosen for accuracy considerations. In practice, a time-step of 2ms is accurate enough, and thus the maximum frequency to be considered is 250 Hz. With that frequency, the time signal recovered after performing the inverse Fourier transform has the same time-step as the reference solution.

The range [-250, 250] Hz is clearly an overkill because the frequency range of the imposed heat flux q is in [-20, 20] Hz. This later range is determined because a Fourier transform of the heat flux seen by a point on Γ Out reveals that harmonics of frequency greater than 20 Hz transfer a negligible amount of energy to the system. Moreover, note that for a particular negative frequency, the transfer function must be the complex conjugate of its symmetric (positive) counterpart. This is also a consequence of the Nyquist-Shannon theorem. Here negative frequencies are also computed to numerically verify that the PGD method reproduces a symmetric real part and an anti-symmetric imaginary part.

As discussed in Section 3 and depicted in [START_REF] Kovalevsky | The Fourier version of the Variational Theory of Complex Rays for medium-frequency acoustics[END_REF], PGD is used to determine an approximation, say ĥn , of the transfer function, ĥ, solution of problem [START_REF] Thomas | The finite element method: linear static and dynamic finite element analysis[END_REF]. In fact, Figure 2 shows the relative residue of (20) as the number of modes n increases. The normalized residue is computed using the L 2 norm of the discrete residue normalized by the norm of the rigth-hand-side of (20a).

For demonstration purposes tolerances are taken small, beyond engineering accuracy. However, 19 terms induce negligible (below 10 -8 !) normalized relative residues. The average number of fixed-point iterations is 25 (same tolerance of 10 -8 ). Since 19 terms are necessary to reduce the residual norm below 10 -8 , the total amount of FE solves in the offline stage is around 475.

Note that spatial modes and frequency modes are localized. Figures 3 and4 show respectively the first three spatial and frequency modes. These modes are X 1 , W 1 , X 2 , W 2 , X 3 and W 3 . Since all of them are complex, the real part and the imaginary part are depicted. As expected, the real part of the frequency modes is symmetric, while the imaginary part is anti-symmetric.

Notice that a uniform spatial discretization is used with 770 bilinear quadrilateral elements of size 0.01m, this implies 930 nodes with scalar complex unknowns. A non-uniform discretization is used for frequency because it varies more rapidly near the origin. The mesh consists of 500 C 0 -continuous linear elements refined around the zero-frequency using a cubic polynomial ω 3 for the element length. Since frequency modes do not involve any derivative their corresponding algebraic equation can be solved point-wise. However, a FE discretization is introduced to approximate the frequency separated functions in a least-squares sense.

To further verify that the PGD approximation of the generalized transfer function ĥn is reasonable, it is evaluated at the extreme frequencies 0 and 250 Hz and then compared with a direct finite element (FE) resolution of ( 15) for those precise frequencies. Figure 5 depicts the difference between both approximations. The generalized solution gives approximations very close to those obtained with an FE computation, errors are always below 10 -5 .

Finally, once ĥn is determined and its inverse Fourier transform computed, h n = F -1 [ ĥn ], equation ( 18) is used to determine the temperature at the desired monitoring point x 0 for a given excitation. The evaluation of the PGD approximation and its inverse Fourier transform is performed only once for any excitation, it is the offline phase. The actual application of the convolution, see equation [START_REF] Hetmaniuk | Review and assessment of interpolatory model order reduction methods for frequency response structural dynamics and acoustics problems[END_REF], for any excitation in order to determine the temperature at the monitored point is the online or post-process stage. Here this is done for the imposed external flux q(ξ, t), see equation [START_REF] Moitinho De Almeida | A basis for bounding the errors of proper generalised decomposition solutions in solid mechanics[END_REF], in the time interval t ∈ [0, 1] seconds, which is the window of interest. The convolution integral is discretized in time using the time-step that comes naturally from the greatest frequency considered in the generalized transfer function. This follows the previous discussion on the use of the Nyquist-Shannon theorem [START_REF] Shannon | Communication in the presence of noise[END_REF].

The evolution with time of the temperature at x 0 is shown in Figure 6. The reference solution (solid red line) is computed with FE and a Crank-Nicolson time-marching scheme. This scheme uses a ratio ∆t/∆x 2 = 20 which has errors below 0.4 10 -4 compared with a reference solution using ∆t/∆x 2 = 1/2 to ensure an accurate transient response. Note that each time-step requires the resolution of a system of equations whose dimension is determined by the FE mesh used.

It is clear from this figure that the proposed method produces an accurate response. However, it is more important to note that the online phase for PGD-based approximation of temperature, which is determined at 500 instants (equispaced by 2ms), requires with MATLAB on a laptop 0.34s, which is almost a third of the physical time 1s. This confirms that given the generalized transfer function h n and an imposed heat flux q(ξ, t) along the outer boundary, the temperature at a point x 0 can be evaluated in real-time (actually faster than real-time!). Note that this 35 times faster than the full FE solution, which needs of around 12 seconds to be computed with a standard commercial code.

Multi-ply composite cylinder: response verification and imperfection influence

An imperfection is introduced in the previously studied problem. That is, the same geometry, see Figure 1, and the same equations, see [START_REF] Maday | The reduced basis element method: application to a thermal fin problem[END_REF], with the same parameters are considered. However, now a zero thickness imperfection of a prescribed length is introduced in the middle of the ply just between the inner boundary and the outer boundary and centered at the measuring point. This imperfection models a possible delamination of the composite and it is modeled as a perfect adiabatic boundary. Obviously this imperfection affects convergence but not drastically, see Figure 7 for a delamination length of 0.2m. Its influence is more clear when plotting modes or the generalized transfer function. In fact, the same comparison shown earlier (that is, the difference between a PGD approximation of the generalized transfer function ĥn and the direct FE) is shown in Figure 8 for the extreme frequencies 0 and 250 Hz. Recall that computations are still done for the whole range [-250, 250] Hz to further verify the symmetric nature of the solution. Again errors are always below 10 -5 . As discussed at the end of Section 2.3, the generalized transfer function problem (8) is non-Hermitian but it is symmetric. The later property proves sufficient for reciprocity, but it is well known that there is no proof of monotonic convergence for the PGD method when confronted to non-Hermitian operators. Thus, the convergence shown in Figure 7 shows a slight increase of the residual norm for the last computed terms.

The influence of the imperfection is even more clear when plotting temperature evolution at the measuring point. Figure 9 shows three cases with their corresponding (expensive) comparison with FE, the no defect case, same curves shown in Figure 6, and two defects, one of length 0.1m, the other with length 0.2m. Since the computer cost of the reference FE solution is similar to the one of the previous section, its overhead with respect to the prosed method is also of the same magnitude. The number of terms used in the PGD expansion is 21 and 25 for the short and long imperfection, respectively, recall that 19 were used with no imperfection. Thus, the imperfection does not increase the number of terms dramatically.

This methodology clearly shows that, in real-time, defects can be detected as the difference between measured and computed temperature values. Moreover, it further opens the possibility to determine the defect nature from its thermal signature.

Multi-parametric extension

As presented in Section 4.1 this PGD-based approach presented here has a potentiality that exceeds real-time monitoring of temperature at a given location. More precisely, thermal conductivity, k, can be chosen as an extra parameter. The generalized transfer function is now parametric in k and consequently, temperature at the monitoring point x 0 for any instance t can be computed in real-time for any conductivity k, i.e. u(x 0 , t, k), see equation [START_REF] Lee | A first course in combinatorial optimization[END_REF]. The single-ply example presented and discussed in Section 5.1 is further generalized for any conductivity k ∈ [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF][START_REF] Kammler | A first course in Fourier analysis[END_REF]W/(m K). In this section, the range of frequencies is up to 60 Hz instead of 250 Hz because there is no need to assess the computational performance of the method compared to the reference solution. In any case, the range of frequencies is taken such that the proper symmetries must be recovered in the computed solution. This is clearly seen in the modes shown in Figures 10 and11.

This problem however is much more challenging because there is an extra parameter. But more important, these difficulties are relevant because variations in thermal conductivity introduce major changes in the real behavior of the thermal field. This is clearly observed in the modes associated to conductivity in Figure 11. Note the large variations introduced close to the lower bound of the thermal range, recall k ∈ [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF][START_REF] Kammler | A first course in Fourier analysis[END_REF]W/(m K). That is, for low conductivities solutions must localize close to the heat source. This has a clear influence in the convergence process of PGD, see Figure 12. Although convergence to engineering precision (0.5 10 -2 ) is obtained with 14 modes, the rate of convergence is slower compared to space-frequency separated representations computed in the previous cases. Moreover, if further precision is required (beyond engineering accuracy), the algorithm fails to converge due to the non-Hermitian character of the operator. See [START_REF] Barbarulo | On a PGD model order reduction technique for midfrequency acoustic[END_REF] for further details and strategies to overcome this issue. Moreover, regarding the fixed-point convergence of the multi-parametric problem (i.e. convergence of each greedy algorithm), the average number of iterations per mode is now increased to 59. Consequently, a total of 826 FE solves are done during the offline phase to compute the necessary 14 terms of the PGD expansion. Note, that a brute force approach sampling (no functional approximation) of the generalized transfer function at the 121 frequency and 77 conductivity nodes would imply 9317 FE solves. Thus, the cost of the brute force computation of the transfer functions depends on the fidelity of the discretization of the frequency and conductivity spaces. Whereas, the cost of the computation of the PGD expansion (and the number of terms in the expansion) is largely independent of the fidelity of the discretization of the frequency and conductivity spaces, assuming the discretization is sufficiently fine.

Nevertheless, the PGD-based scheme converges globally although its local behavior, precisely near the lower bound for conductivity, shows less precision compared to larger values of k. This is better appreciated in the temperature evolution at the measuring point. Figure 13 shows this evolution for several values of the conductivity with their corresponding (expensive) comparison with FE. Results are in very good agreement with the reference FE solution and precision increases as k increases. Now, since the largest frequency considered for the evaluation of the generalized transfer function is 60Hz, the response at the measuring point x 0 is recovered with a time-step of 8ms, recall [START_REF] Monteyne | Heat transfer in a borehole heat exchanger: Frequency domain modeling[END_REF]. Consequently, the computational time needed for the online phase is faster than in previous examples, only 0.09 seconds. Whereas, the FE reference solution is computed with the same timemarching scheme described in previous sections (∆t = 2ms) because this implies that 20 ≤ k∆t/∆x 2 ≤ 400. Therefore, both solutions are not strictly comparable. The computational time associated to the FE solution remains unchanged (around 12 seconds) and the online stage of the proposed approach is now 133 times faster than the FE one. In any case, this approach clearly demonstrates its applicability to inverse problems where the measured temperature is used to determine thermal conductivity. Remember, that the temperature evolution at the measuring point is obtained extremely fast.

Amplitude identification

This example shows the applicability of the proposed approach for the model case of identification discussed in Section 4.2. First, some measured temperature is needed, that is, u meas (x 0 , t r ) for r = 1, . . . , m. In this case, the "measured" temperature is synthetically generated at m = 501 instances (∆t = 2ms) with an inflow forcing excitation equal to the one defined in (29) whose amplitude is modulated by [1 + cos(2πt)]/2. That is, a FE code with a Crank-Nicolson timemarching scheme is used to generate the temperature data at the monitoring point under an external heat source

q(ξ, t) = 250[1 + cos(2πt)] exp -50(2ξ -πt) 2 W/m 2 .
Second, this data is used to determine the laser input amplitude following the procedure described in Section 4.2. The amplitude of the excitation is assumed not known and it is approximated following a piecewise linear approximation with n fit = 50 (i.e. a uniform mesh of 51 nodes each 20ms), see [START_REF] Love | A treatise on the Mathematical Theory of Elasticity[END_REF]. Figure 14 shows the synthetically generated temperature at the monitoring point (left) and a comparison (right) between the approximated nodal values (blue markers) and the reference amplitude (solid red line). The coincidence is remarkable.

Conclusions

This paper proposes a novel approach to monitor in real-time thermal processes. This technique provides temperature evolution at a particular point in the domain, usually, where thermocouples will be placed. The online computation is drastically speed-up. This can be used in a large number of problems ranging from a simple surveillance of the process to a simulation-based control, as well as identification problems for defects, material properties or power oscillations.

This approach recourses to the reciprocity principle, which can be employed in thermal problems when working in the frequency domain where the equation is symmetric (but not Hermitian). Then a convolution of the heat source and a generalized transfer function is the only computational effort that the online stage requires. Note that the generalized transfer function is determined by means of the Proper Generalized Decomposition for all the range of frequencies excited by the heat source.

Numerical results show the potentiality and the accuracy of this approach. Examples of real-time monitoring, defect detection, multi-parametric evaluation and inverse calibration are shown.

A Green's function problem for a parabolic operator

Although what follows is well known it is not standard to write the Green Consequently, using ( 2) and (3) in the previous identity [START_REF] Parés | Bounds of functional outputs for parabolic problems. I. Exact bounds of the discontinuous Galerkin time discretization[END_REF], a representation for u(x 0 , t 0 ) with (x 0 , t 0 ) ∈ Γ N ×]0, T [ is obtained, namely u(x 0 , t 0 ) = Γ N t0 0 G(x, t; x 0 , t 0 ) q(x, t)dtdΓ + Ω u 0 (x) G(x, 0; x 0 , t 0 )dΩ, [START_REF] Park | The use of the Karhunen-Loève decomposition for the modeling of distributed parameter systems[END_REF] where now it can be clearly identified that the causality condition implies that the solution at time t 0 cannot depend on any of its values at later times. Recall T can be arbitrarily large.

B Reciprocity proof in the frequency domain for even and odd real excitations

This appendix is aimed to prove the reciprocity principle, recall [START_REF] Dorival | A substructured version of the variational theory of complex rays dedicated to the calculation of assemblies with dissipative joints in the medium-frequency range[END_REF] q1 , û2 = q2 , û1 ,

when both q1 and q2 are real. As discussed in Section 2.5, reciprocity, holds if the two conditions stated in equation ( 14) are verified, namely ∇û 1 , ∇û 2 = ∇û 2 , ∇û 1 and û1 , û2 = û2 , û1 .

Recalling that •, • is the L 2 scalar product of complex functions in Ω, see [START_REF] Chinesta | First steps towards an advanced simulation of composites manufacturing by automated tape placement[END_REF], these conditions are equivalent to Im ∇û 1 , ∇û 2 = 0 and Im û1 , û2 = 0.

Since û1 and û2 give values in C, one can define the real functions a 1 , b 1 , a 2 and b 2 such that û1 = a 1 + ib 1 and û2 = a 2 + ib 2 , and the previous expressions are equivalent to ∇a 1 , ∇b 2 = ∇b 1 , ∇a 2 (34a)

a 1 , b 2 = b 1 , a 2 . ( 34b
)
Recall also that û1 and û2 are the corresponding solutions of the weak problem [START_REF] Bia Lecki | Proper orthogonal decomposition and modal analysis for acceleration of transient FEM thermal analysis[END_REF] for the two excitations q1 and q2 and for any imposed frequency ω. By definition, a i and b i for i = 1, 2 also belong to space of trial and test functions, namely H 1 Γ D . Thus, equation ( 8) for the pair {û 1 , q1 } can be particularized as follows:

∇û 1 , ∇a 2 + iω û1 , a 2 = q1 , a 2 ,

∇û 1 , ∇b 2 + iω û1 , b 2 = q1 , b 2 , (35)
and likewise for {û 2 , q2 } ∇û 2 , ∇a 1 + iω û2 , a 1 = q2 , a 1 ,

∇û 2 , ∇b 1 + iω û2 , b 1 = q2 , b 1 . (36)
Recall the splitting of û1 into real and imaginary parts, û1 = a 1 + ib 1 , to also split equations [START_REF] Prulière | On the solution of the heat equation in very thin tapes[END_REF] into real and imaginary parts as

∇a 1 , ∇a 2 -ω b 1 , a 2 = q1 , a 2 , (37a)
∇b 1 , ∇a 2 + ω a 1 , a 2 = 0, (37b)

∇a 1 , ∇b 2 -ω b 1 , b 2 = q1 , b 2 , ( 37c
)
∇b 1 , ∇b 2 + ω a 1 , b 2 = 0, ( 37d
)
where the hypothesis that both q1 (x, ω) and q2 (x, ω) belong to R is used. Likewise, for û2 = a 2 + ib 2 each equation in ( 36) is split into real and imaginary parts as Recall now that the L 2 scalar product is symmetric for any pair of real functions, more specifically u, v = v, u for all u and v ∈ R. Then subtract (38b) from (37b), and (38d) from (37d) to obtain the desired conditions [START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the Proper Generalized Decomposition[END_REF]. Thus if [START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the Proper Generalized Decomposition[END_REF] are verified, then (14) also holds and reciprocity is demonstrated. Note that these results hold for any ω.

∇a 2 , ∇a 1 -ω b 2 , a 1 = q2 , a 1 (38a)
To further close these appendix, note that the other equations not used up to now also produce the same results. After subtracting (38a) from (37a) and (38c) from (37c), the following equations are obtained:

ω b 2 , a 1 -b 1 , a 2 = q1 , a 2 -q2 , a 1 , ∇a 1 , ∇b 2 -∇a 2 , ∇b 1 = q1 , b 2 -q2 , b 1 .
However, both left-hand-sides in the previous equations are zero because they correspond to [START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the Proper Generalized Decomposition[END_REF], which was just proven, and these equations become q1 , a 2 = q2 , a 1 and q1 , b 2 = q2 , b 1 , which corresponds to split (33) into real and imaginary parts using the definitions splitting of û1 = a 1 + ib 1 and û2 = a 2 + ib 2 .

Finally

  the normal equations Aα = b are solved.

Figure 1 :

 1 Figure 1: Single-ply composite cylinder: problem statement.

Figure 3 :

 3 Figure 3: Single-ply composite cylinder: first 3 PGD spatial modes real (left) and imaginary (right).

Figure 4 :

 4 Figure 4: Single-ply composite cylinder: first 3 PGD frequency modes.

Figure 5 :Figure 6 :

 56 Figure 5: Single-ply composite cylinder: difference between PGD and FE solutions, real (left) and imaginary (right) parts, for frequencies 0 Hz (top) and 250 Hz (bottom).

Figure 7 :

 7 Figure 7: Multi-ply composite cylinder: convergence of the generalized frequency transfer function.

Figure 8 :

 8 Figure 8: Multi-ply composite cylinder: difference between PGD and FE solutions, real (left) and imaginary (right) parts, for frequencies 0 Hz (top) and 250 Hz (bottom).

Figure 9 :

 9 Figure 9: Multi-ply composite cylinder: comparison of the temperature evolution at point x 0 for the proposed PGD-based solution (discontinuous) and the standard FE (solid) for different defect lengths.

Figure 10 :

 10 Figure 10: Single-ply composite cylinder: first 3 PGD spatial modes real (left) and imaginary (right).

Figure 11 :

 11 Figure 11: Single-ply composite cylinder: first 3 PGD frequency (left) and thermal conductivity (right) modes.

Figure 12 :

 12 Figure 12: Multi-parametric convergence of the generalized frequency transfer function.

Figure 13 :

 13 Figure 13: Multi-parametric comparison of the temperature evolution at point x 0 for the proposed PGD-based solution (discontinuous) and the standard FE (solid) for thermal conductivities.

Figure 14 :

 14 Figure 14: Amplitude identification: synthetically generated temperature measurements (left) used to calibrate the amplitude of the heat source (right), calibrated values (blue markers) and reference solution (solid red line).

0 G(∇ 2 u 0 G

 020 's function problem associated to the homogenous heat equation with homoge-neous Dirichlet and non-homogeneous Neumann boundary conditions. The adjoint Green's function problem associated to (2) is detailed in (3) where its last equation, namely G = 0 on Ω × [t 0 , T [, corresponds to the causality condition.Green's identity can be written asΩ t0 -∂ t u) -u(∇ 2 G + ∂ t G) dtdΩ = ∂Ω t0 (n • ∇)u -u(n • ∇)G dtdΓ + Ω uG| t=0 -uG| t=t0 dΩ.[START_REF] Parés | Bounds of functional outputs for parabolic problems. I. Exact bounds of the discontinuous Galerkin time discretization[END_REF]

∇b 2 ,

 2 ∇a 1 + ω a 2 , a 1 = 0 (38b) ∇a 2 , ∇b 1 -ω b 2 , b 1 = q2 , b 1 (38c) ∇b 2 , ∇b 1 + ω a 2 , b 1 = 0 (38d)

  ) For an odd function in time v o , vo is redefined as the imaginary part of F[v o ]. Thus, F[v o ] = iv o with vo ∈ R.

	etc). In what follows it is important to recall that, in general, v ∈ C; but, for an
	even function in time v e , i.e. v e (x, -t) = v e (x, t), F[v e ] = ve ∈ R; whereas, for
	an odd function in time v o , i.e. v o (x, -t) = -v o (x, t), F[v o ] is imaginary, i.e.
	F[v o ] ∈ iR. Applying the Fourier transform to problem (2) it becomes
	Remark 1 (Fourier transform properties). Fourier transforms have been largely
	studied and they hold a large number of properties (viz. linearity, translation,

Elías Cueto

Aragon Institute of Engineering Research (I3A), Universidad de Zaragoza, Maria de Luna 3, E-50018 Zaragoza, Spain. e-mail: ecueto@unizar.es web http://amb.unizar.es