Joint identification via deconvolution of the flux and energy relaxation kernels of the Gurtin-Pipkin model in thermodynamics with memory

In this paper we present a linear method for the identification of both the energy and flux relaxation kernels in the equation of thermodynamics with memory proposed by M.E. Gurtin and A.G. Pipkin. The method reduces the identification of the two kernels to the solution of two (linear) deconvolution problems. The energy relaxation kernel is reconstructed by means of energy measurements as the solution of a Volterra integral equation of the first kind which does not depend on the still unknown flux relaxation kernel. Then, flux measurements are used to identify the flux relaxation kernel.


Introduction
The linearized version of the model proposed in [9] by Gurtin and Pipkin to describe thermodynamical processes with memory depends on two memory kernels β(t), the energy relaxation kernel, and a(t), the flux relaxation kernel. These kernels are material properties of the body and have to be identified using suitable measurements. Several identification methods have been proposed in engineering and mathematical papers (see the references in section 2.1). In most of the cases, these methods assume β = 0. Here we are going to extend to the case β = 0 a linear algorithm first proposed in [16] for the identification of the kernel a(t) when β(t) = 0. The bonus of this algorithm is that it reduce the identification of the kernels to (linear) deconvolution problems.
In order to understand the identification algorithm, we shall shortly present the derivation of the Gurtin-Pipkin model in section 2. We shall see that the evolution in time of the temperature is described by the following equation Here θ = θ(x, t) is a function of the time t and of a variable x in a region Ω (∆ is the laplacian in the variable x).
Note that the variable x and also the time variable will not be explicitly indicated unless needed for clarity.
As explained in [4], the identification of the relaxation kernels can be obtained using a sample in the form of a bar. Hence, we shall consider the case Ω = (0, L) so that ∆θ = θ xx . But, the extension of the algorithm to general regions is important for applications to nondestructive testing. This extension is reserved to the future.
We state the assumptions which we use to justify the identification algorithm: Assumption 1 1. a(t) and β(t) are bounded and of class C 2 (0, +∞) (the derivatives have a continuous extension to t = 0).
2. a(0) > 0 (This condition is particularly important since it implies that signals propagate with finite velocity and for this reason eq. (1) is also called the hyperbolic heat equation).
In fact, the derivation of eq. (1) assumes also that a(t) and β(t) and their derivatives are integrable on [0, +∞) (and so a(t), β(t) tend to zero for t → +∞). Thermodynamics impose stronger conditions to the kernels (for example, it is proved in [6] that under natural conditions the kernels decreases and certain "positivity" conditions are proved in [1]). These properties are not used in the identification algorithm, but they can be used to improve the numerical implementation, as discussed in [16]. The properties of the kernels a(t) and β(t) (which are not used in the justification of our algorithm) imply that the system dissipates energy. The organization of the paper is as follows: in section 2 we describe informally the identification algorithm. In order to understand the rational under the algorithm, we describe also the derivation of eq. (1). Section 3 presents preliminary results on the solutions of eq. (1) (which of course depend on suitable initial and boundary conditions) while the algorithm is justified in section 4 while section 2.1 contains comments to the literature.

2
Informal description of the identification algorithm In order to understand the rationale behind the identification algorithm, we must shortly describe the derivation of eq. (1) in [9]. Note that eq. (1) is a linearized version of a "true" nonlinear model (see also [3] for a nice description of the nonlinear model) and so θ does not represent the absolute temperature but it is the perturbation of a stationary temperature of a nonlinear process. Hence, θ can be zero. The derivation of eq. (1) depends on the following two principles: 1. the (density of the) flux of heat q(x, t) at position x and time t depends on the past history of the gradient of the temperature θ(x, t) at the same position x, and it is given by Equation (2) replace the usual Fourier Law q(x, t) = −∇θ(x, t).
2. the internal energy is related to the temperature by the relation This relation replace the usual relation e(x, t) = b + cθ(x, t). In the following we assume that the scale has been chosen so to have c = 1.
Equation (1) follows when we balance the energy at every time t and position x, i.e. we impose d dt e(x, t) = −∇ · q(x, t) .
We note that eq. (1) can be written in several equivalent forms, in particular (we recall: c = 1). We did not yet specify the initial and boundary conditions. We fix an initial time t 0 , it is not restrictive to put t 0 = 0, and we assume that the measurements are taken for t > 0. So, in order to solve eq. (4) we need the initial condition θ(x, s) = ξ 0 (x, s) given for s < 0 , θ(x, 0) = ξ(x) .
The function ξ 0 (x, s) needs not be continuous so that in general ξ 0 (x, 0) does not even makes sense (and it is not equal to ξ(x)). Discontinuity at zero can be practically realized by putting the body in contact with a suitable source. This observation is used in the identification algorithm since we assume θ(x, s) = ξ 0 (x, s) = 0 s < 0 and ξ possibly different from zero. When θ(x, t) = ξ 0 (x, t) = 0 for t < 0, eq. (4) takes the form Due to the fact that a(t) is related to the flux, it is clear that its identification will require the measure of the flux (which can be measured only at the ends of the bar). The kernel β(t) is related to the relaxation of energy, and the quantity that can be directly measured is the temperature. So, in order to identify β(t) we measure the temperature. In the case of a bar, conceivably we can measure the temperature at each point.
The relaxation kernels are identified by taking the following measurements of the sample of the material: 1. first we impose boundary conditions θ(0, t) = θ(L, t) = 0 and an initial temperature θ(x, 0) = ξ(x). We measure (a) the total enery as a function of time. But, the quantity which can be measured is the temperature. So we estimate 2. then we repeat the same measurements but with ξ = 0 while θ(0, t) = f (t) and θ(L, t) = 0: (a) the total enery as a function of time. In fact we measure Remark 2 Note the use of a nonzero initial temperature. While there is no problem to apply a time varying temperature f (t) at one end of the bar, the initial condition is more difficult to realize. We shall see that we need a special initial condition which is easily realized.
The energy relaxation kernel β(t) can be computed from the measurements described in the items 1a and 2a even if a(t) is still unknown.
Once β(t) has been computed then a(t) is computed from the measurements described in the items 1b and 2b.

Remark 3
We observe that the integral L 0 θ(x, t) dx can be estimated since in practice the kernels can be determined from samples which have the form of a bar. It has an interest to see that identification of β(t) can be acieved also by measuring the temperature θ(L, t) when we impose the boundary condition θ x (L, t) = 0 (i.e. the right hand side is now insulated) and we measure This variant is examined in section 4.1.

References to previous results
Equation (1) has important applications in thermodynamics and viscoelasticity and in other applications (like nonfickian diffusion, which is often encountered in biological matherials). Hence the determination of the parameters in the equation has been widely studied in mathematics and in engineering journal. See the reference in [16]. The methods used in Engineering journals are mostly based on this idea: the kernels to be identified are assumed to belong to specific classes, for example Prony sums, and depend on few parameters. The equation is discretized and the parameters of the kernels are determined so that the discretized version of the solution best fit the experimental measure y(t). Usually y(t) is a measure of temperature or flux or, in viscoelasticity, it is a measure of the traction on the boundary. This method was first proposed to identify one memory kernel (see for example [4]) but of course it can be used as well to identify the two memory kernels in eq. (1). In this contest we mention in particular the paper [7] which is concerned wit three dimensioanl viscoelasticity.
In this case, different relaxation kernels correspond to the shear and bulk creep relaxation kernels.
The method proposed in [12] is as follows. It is noted that the observation y(t) is a linear functional of θ, which is unknown since eq. (1) cannot be solved, because the kernels are still unknown. But then, the equation of y and eq. (1) constitute a nonlinear system in the unknown (θ, N ). The solution of this system provides both θ and the unknown kernel N . Usually this system is solved by first computing the derivative of y so to get a nonlinear differential equations in the unknowns (θ, N ) in a suitable Hilbert spaces. The method, first proposed in the case that the equation depends on only one kernel was then extended to the identification of two kernels, for example in [5,10,11,13]. We note that this method requires the solution of highly nonlinear equation in a Hilbert space but it is noted in [8] (in the case of one unknown kernel, but the observation is easily extendable) that once the kernels and θ have been identified in a first interval [0, τ ] then, thanks to the Volterra structure of the problem, the identification up to T is reduced to a linear problem in a Hilbert space. In contrast with this, the method in this paper requires the solution of two Volterra integral equations of the first kind in R, since the quantity to be identified are the scalar functions β(t) and a(t).

Representation of the solutions
The solutions of eq. (6) have been studied by several authors. A nice reference is [2]. Here we need a Fourier type approach which is similar to the one used in [15]. We consider the operator A in L 2 (0, L) which is defined as follows: Let φ n (x) = γ 0 sin λ n x , γ 0 = 2 L , λ n = n π L .
The function φ n is an eigenvector of A whose eigenvalue is −λ 2 n and {φ n } is an orthonormal basis L 2 (0, L). We expand the solutions of eq. (6) with conditions θ(0) = ξ and The function θ n (t) solves We intoduce the functions z n (t) which solve We note that We compute explicitly the derivative under the integral sign and we see that The last integral can be manipulated as follows (when f (t) is differentiable): In the following, we shall use f ∈ C 1 with f (0) = 0 (consistent with θ(x, 0) = 0), i.e. Hence

Identification of the kernels
Here we show that the proposed method enables us to identify the relaxation kernels by solving (linear) deconvolution problems. First we consider the identification of the energy relaxation kernel β(t).
Identification of the energy relaxation kernel. Note that Then it is easily computed that, when f (t) = 0: We use the special initial condition ξ(x) whose Fourier coefficients are ξ n = 1/λ n The function provided by this measurement is denote H(t):

Remark 4
We shall see that samples with this special initial condition are easily realized in practice without using the eigenvalues and eigenvectors of the problem.
Now we use formula (12) with ξ n = 0 to compute The left hand side is measured, H(t) has been already estimated and so the first and second terms on the right side are known. Hence, β(t) is computed from the deconvolution problem: Note that a(t), still unidentified, does not appear in the previous computations.
Identification of the flux relaxation kernel. It is not possible to measure the flux inside a body. So, in this case we are forced to measure the flux at L, when the boundary conditions are θ(L, t) = 0 .
and we measure the two fluxes produced either by ξ or by f . Let us first consider the flux when f = 0 and ξ(x) = 0. In this case we have This flux can be measured, provided that the initial condition ξ can be realized. As already stated and discussed below, we can realize that initial condition ξ whose Fourier coefficients are ξ n = 1/λ n . So, we measure the corresponding flux, that we denote K(t): Now we compute the flux throughot the right hand L due to the boundary temperature θ(0, t) = f (t) (f (t) as in eq. (11)) and ξ = 0. Using eq. (12) we see that is a known function since β(t) has already been identified. So we have In order to compute the flux throughout L we must compute the derivative respect to x. We note that d dx φ n (L) = (−1) n γ 0 λ n . So, the contribution of the second line is (15) and this is a known quantity. Now we consider the first addendum on the right side of eq. (14). We note that d dx (δ is the Dirac's delta). So, the flux at L is The relaxation kernel a(t) can be computed from this equality via deconvolution. We observe that the exchange of the series and the derivative in eq. (15) is easily justified since the sequence of continuous functions {z n (t)} is bounded on bounded intervals.
This ends the description of the identification procedure. It remains to understand whether it is possible to realize easily a source with the special temperature The positive answer was given in [17]. The temperature ξ(x) is the solution of the following problem: So, it is the stationary temperature achieved by a thermal body whose boundary temperatures are θ(0, t) = 0 , θ(L, t) = 1 γ 0 .
In fact, it is known that the heat equation η ′ = η xx with conditions η(0, t) = η(L, t) = 0 is exponentially stable, so that lim t→+∞ η(x, t) = 0 in L 2 (0, L) for every initial condition η(x, 0) = ξ(x) ∈ L 2 (0, L). Let θ solves The function solves the heat equation with the boundary condition put equal zero and the initial condition equal to ξ. Stability implies lim t→+∞ η(·, t) = 0 in L 2 (Ω) so that in order to realize a source with the special temperature ξ 0 needed in the identification process it is sufficient to apply the constant boundary temperatures at x = 0, x = L as described in (16), for a time large enough. Note that the same arguments apply not only to the solutions of the standard heat equation but also to the solutions of eq. (6), provided that the kernels are so chosen that the system is stable, as usually happens in practice.

Remark 5
The Fourier expansion has been used to justify the procedure, but the actual implementation of the procedure does not need any Fourier expansion, not even in order to realize the special initial condition needed in this procedure.

The variant in Remark 3
We show that the procedure in Remark 3 leads to the identification of β. We note that in this case the operator A has to be replaced by The eigenvectors of this operator are φ n (x) = γ 0 sin n + 1 2 with eigenvalues −λ 2 n where now λ n = n + 1 2 (π/L). Note that φ n (L) = γ 0 (−1) n .
Granted this new meanings of the symbols, the formulas for θ n (x, t), z n (t) are still given by eqs. (8) and (9). The manipulations in (10) still hold and θ(x, t) has the expansion θ(t) = ∞ n=1 θ n (t)φ n (x) in eq. (12) where now φ n is in (17). So, the algorithm can be justified as follows.
Convergence of this series is easily seen because {ξ n } ∈ l 2 and the same estimates as in [14] shows that z n (t) = cos λ n t + 1 λ n M n (t) where {M n (t)} is a sequence of continuous functions which is bounded on bounded intervals. We call H(t) the temperature θ(L, t) when So, β(t) is given by the following deconvolution problem (convergence of the numerical series is discussed below) The right hand side is known since H(s) is known from the previous measurement. Leibniz test shows that the numerical series in (18) converges but it converges slowly. In spite of this, the identification procedure can be performed since the sum of the series can be explicitly computed, This is easily seen by noting that for x ∈ (0, 1) we have The series converges also for x = 1 and so Abel Theorem shows that it converges to a function which is continuous also for x = 1. The sum of the series is obtained by computing both the sides for x = 1.