# Efficient solutions for transient thermal processes involving moving heat sources

Jose Aguado, Antonio Huerta, Francisco Chinesta, Elias Cueto, Adrien Leygue

► To cite this version:

Jose Aguado, Antonio Huerta, Francisco Chinesta, Elias Cueto, Adrien Leygue. Efficient solutions for transient thermal processes involving moving heat sources. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717102

HAL Id: hal-01717102

https://hal.science/hal-01717102

Submitted on 25 Feb 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Efficient solutions for transient thermal processes involving moving heat sources

Jose V. AGUADO<fnref n="1" />, Antonio HUERTA<fnref n="2" />, Francisco CHINESTA<fnref n="3" />, Elias CUETO<fnref n="4" />, Adrien LEYGUE<fnref n="5" />

**Résumé** — Many industrial processes involve transient heat transfer problems where the heat source moves over the domain. For instance, the Automated Tape Placement (ATP) in the context of composites manufacturing may be considered. In this process, the composite piece is built up by adding different layers successively, which are welded on the substrate by means of a local heat source, typically a laser, which moves all over the boundary. If one is interested in monitoring or controlling the process, some thermocouples need to be installed where knowing the temperature is of interest. However from a numerical point of view, having a moving heat source poses several problems. Classically, even when the problem remains linear, it would be needed to solve the transient problem by performing an appropriate time stepping, needing for the solution to solve many algebraic problems. In order to overcome these limitations, this work proposes solving the problem in the frequency domain because in this case it can be proved that the reciprocity principle is satisfied. Then, a unitary heat source may be applied where the temperature measurement is performed giving place to a transfer function which relates the measurement point and the boundary, where the heat flux applies. However this transfer function is only valid for a single frequency, and thus if the heat flux signal contains several frequencies, many transfer functions need to be computed. Instead, we prefer to compute a generalized transfer function by using the PGD method, i.e. the frequency is included as an extra-coordinate, like the physical space. Once this generalized transfer function has been computed, a simple and computationally cheap postprocessing suffices for obtaining the temperature response at the point of interest, for any heat source path on the boundary.

**Mots clés** — monitoring, processes, heat transfer, PGD, transient problem

# 1 Introduction

Consider a body $\Omega \subset \mathbb{R}^d$, $d \le 3$, whose boundary $\partial\Omega$ is partitioned in a Dirichlet, $\Gamma_D$, and a Neumann, $\Gamma_N$, part such that $\bar{\partial}\Omega = \bar{\Gamma}_D \cup \bar{\Gamma}_N$ and $\Gamma_D \cap \Gamma_N = \emptyset$. The time interval of interest is $I = [0, T]$. Under the assumption of a isotropic and homogeneous linear material, the temperature evolution $\hat{u}(x,t)$, $x \in \Omega$ and $t \in I$, is described by the transient heat equation :

$$ egin{cases} 
\displaystyle \frac{\partial \hat{u}}{\partial t} - \Delta \hat{u} = 0 & \text{in } \Omega \times I, \\ 
\hat{u} = 0 & \text{on } \Gamma_D \times I, \\ 
n \cdot \nabla \hat{u} = \hat{q} & \text{on } \Gamma_N \times I, \\ 
\hat{u} = \hat{u}_0 & \text{on } \Omega \times \{0\}, 
\end{cases} 
\quad (1) $$

where, without loss of generality, all material constants are suppressed for the sake of clarity, $n$ is the exterior unit normal to $\Gamma_N$ and $\hat{q} = \hat{q}(x,t)$ is a heat flux on the Neumann boundary, which could eventually include convective heat exchanges. Note that homogenous Dirichlet boundary conditions are imposed without loss of generality because the problem is linear.

Let us assume an imposed heat flux of frequency $\omega$, i.e. $\hat{q}(x,t) = \xi(x) \exp(i\omega t)$ for $x \in \Gamma_N$, and $\xi(x) \in \mathbb{R}$. Due to the linearity, we simply state that the response will have the same frequency but with a

<footnote n="1">¹ Ecole Centrale de Nantes, GeM, jose.aguado-lopez@ec-nantes.fr</footnote>

<footnote n="2">² Universitat Politecnica de Catalunya, LaCaN, antonio.huerta@upc.es</footnote>

<footnote n="3">³ Ecole Centrale de Nantes, GeM, francisco.chinesta@ec-nantes.fr</footnote>

<footnote n="4">⁴ Universidad de Zaragoza, I3A, ecueto@unizar.es</footnote>

<footnote n="5">⁵ Ecole Centrale de Nantes, GeM, adrien.leygue@ec-nantes.fr</footnote>
phase lag, namely

$$ \hat{u}(x,t) = u(x; \omega)e^{i\omega t}. \qquad (2) $$

Note, that now, in the frequency domain formulation the unknown $u(x; \omega) \in \mathbb{C}$. Moreover, to emphasize the parametric dependence of the solution $u$ on $\omega$, this is explicitly indicated in the expression $u(x; \omega)$. Introducing (2) into (1) the frequency domain problem is defined as

$$ \left\{ \begin{array}{ll} i\omega u - \Delta u = 0 & \text{in } \Omega, \\ u = 0 & \text{on } \Gamma_D, \\ n \cdot \nabla u = \xi & \text{on } \Gamma_N. \end{array} \right. \qquad (3) $$

It can be proved, although it is not shown here, that (3) satisfies the reciprocity principle, which roughly can be stated as the temperature measured at a point A when an external heat flux is applied at B (i.e. heat source) is the same that if the heat flux was applied in A and the temperature measured in B. Then, when a problem involving a moving load has to be addressed we would prefer a model related to a unitary heat source applied where the temperature measurement is performed. The result of such computation gives place to a transfer function, denoted as $h(x; \omega)$, which relates the measurement point and the rest of the boundary (except the Dirichlet boundary), where the heat flux applies :

$$ \left\{ \begin{array}{ll} i\omega h - \Delta h = 0 & \text{in } \Omega, \\ h = 0 & \text{on } \Gamma_D, \\ n \cdot \nabla h = 1 & \text{on } \Gamma_M \subset \Gamma_N, \end{array} \right. \qquad (4) $$

where $\Gamma_M$ represents the boundary where the temperature is measured. We enfatize the parametric dependence of $h$ on $\omega$ by denoting $h(x; \omega)$. Now, the variational form associated to Eq. (4) reads : find $u \in \mathcal{H}_{\Gamma_D}^1$ such that

$$ (\nabla h, \nabla v) + i\omega(h, v) = \langle 1, v \rangle \quad \forall v \in \mathcal{H}_{\Gamma_D}^1, \qquad (5) $$

where both $h, v \in \mathbb{C}$,

$$ \mathcal{H}_{\Gamma_D}^1 = \{ v \in \mathcal{H}^1(\Omega) : v = 0 \text{ on } \Gamma_D \}, $$

and as a matter of notation

$$ (z, w) = \int_{\Omega} z \cdot \bar{w} d\Omega, $$

$$ \langle z, w \rangle = \int_{\Gamma_N} z \cdot \bar{w} d\Gamma, $$

being $\bar{w}$ the complex conjugate of $w$.

## 2 PGD frequency-space separated representation

One of the key points of this work is to compute $h(x, \omega)$, $x \in \Omega$ and $\omega \in \mathcal{K}$ by using the PGD framework [1, 2, 3], in which $\omega$ will be considered as an extra-coordinate [4], like the physical space $x$, instead of computing $h(x)$ for each $\omega$. We define :

$$ h(x, \omega) \approx \sum_{i=1}^{N} X_i(x) \cdot W_i(\omega) , \quad x \in \Omega, \omega \in \mathcal{K} \qquad (6) $$

Eq. (6) is built-up by performing a Greedy algorithm, a sort of enrichment procedure where new terms are progressively added to the series until the convergence quantified by using an appropriate error estimator is reached [5, 6].

Thus, suppose that we have already computed $n$ terms of the solution, and we seek a new one for enriching the solution. The $n+1$ approximation reads :

$$ h^{n+1} = h^n + S(x) \cdot R(\omega) = \sum_{i=1}^{n} X_i(x) \cdot W_i(\omega) + S(x) \cdot R(\omega). \qquad (7) $$
In what follows, the dependence of $X_i$, $S$ on $x$ and $W_i$, $R$ on $\omega$ will be frequently omitted to alleviate the notation. The variational form of the problem to be solved reads : find $h \in \mathcal{U}$ such that

$$ \int_{\mathcal{K}} (\nabla h, \nabla v) \, d\omega + \int_{\mathcal{K}} i\omega(h, v) \, d\omega = \int_{\mathcal{K}} \langle 1, v \rangle \, d\omega \quad \forall v \in \mathcal{U}, \qquad (8) $$

being

$$ \mathcal{U} = \{h : h(\cdot, \omega) \in \mathcal{H}^1(\Omega), h(x, \cdot) \in \mathcal{L}^2(\mathcal{K}), h = 0 \text{ on } \Gamma_D \times \mathcal{K}\}. $$

Given the separated nature of the solution, as a test function we consider

$$ v = v_S(x) \cdot R(\omega) + S(x) \cdot v_R(\omega). $$

By introducing Eq. (7) into (8)

$$ \int_{\mathcal{K}} \int_{\Omega} R \nabla S \cdot \nabla \bar{v} \, d\boldsymbol{x} \, d\omega + \int_{\mathcal{K}} i\omega \int_{\Omega} R S \cdot \bar{v} \, d\boldsymbol{x} \, d\omega = \\ = \int_{\mathcal{K}} \int_{\Gamma_N} 1 \cdot \bar{v} \, d\boldsymbol{x} \, d\omega - \sum_{i=1}^{n} \int_{\mathcal{K}} \int_{\Omega} W_i \nabla X_i \cdot \nabla \bar{v} \, d\boldsymbol{x} \, d\omega - \sum_{i=1}^{n} \int_{\mathcal{K}} i\omega \int_{\Omega} W_i X_i \cdot \bar{v} \, d\boldsymbol{x} \, d\omega \quad \forall v \in \mathcal{U}. \quad (9) $$

As stated before, now we perform a fixed point algorithm by first computing $S$ from the assumption of function $R$ is known, and vice versa, i.e.

1. Compute $S$ from $R$ known. Then, $v = v_S(x) \cdot R(\omega)$ and $S$ can be computed from Eq. (9) because all integrals defined on $\mathcal{K}$ can be computed.

2. Compute $R$ from $S$ known. Then, $v = S(x) \cdot v_R(\omega)$ and $R$ can be computed from Eq. (9) because all integrals defined on $\Omega$ can be computed.

Details are omitted for the sake of brevity, and we refer to the bibliography for information [1, 2, 3].

### 3 Coming back to the transient solution

In Section 2 a solution valid not only for a single frequency but for all of them in a given range has been computed. Here we assume that we need to know the temperature in a certain location $\Gamma_M$ on the boundary, and there is an excitation moving on the domain boundary. As we have previously computed a transfer function $h(x, \omega)$ by applying a unit flux in $\Gamma_M$, then by reciprocity :

$$ 0 = \langle 1, u \rangle - \langle q, h \rangle \qquad (10) $$

from where $u(x, \omega)$ can be computed. Then the transient response can be recovered by performing the Inverse Fourier transform, i.e.

$$ \hat{u}(x, t) = \int_{-\infty}^{+\infty} u(x, \omega) e^{i\omega t} \, d\omega \qquad (11) $$

### 4 Numerical example

Aiming to demonstrate the ability of the proposed method for solving transient models, here we solve an example defined in a 2D geometry which involves a flux moving over the upper boundary of a plate domain.

Fig. 1 shows schematically the proposed model to be solved. A punctual heat flux moves at a constant velocity, $v = 1.0 \, m/s$, over the upper boundary. The heat flux vanishes on the other domain boundaries, except on both the left and right edges where a homogeneous Dirichlet condition has been considered. A measurement is performed at the middle point of the bottom boundary. The domain dimensions are 1.0 m long and 1 m high. The heat flux intensity is considered to be unitary and, for the sake of simplicity, constant through the time. Regarding the material properties, we consider unit values. As it has been explained throughout this paper, a unitary heat flux can be applied where the measurement is going to be performed. Once the generalized transfer function, $h(x, \omega)$, has been computed via the PGD method, thanks to the reciprocity principle, we are able to reconstruct the solution at the measurement point, $u(x, \omega)$, and then $u(x, t)$, by performing a simple postprocessing. Thus we divide the study in two parts :
<figure><img src="image_1.png" /><figcaption>Fig. 1 – Sketch of the numerical example</figcaption></figure>

- The off-line computation of the generalized transfer function, $h(x, \omega)$, related to a unitary flux applied at the measurement point.

- The postprocessing step, i.e. the exploitation of the generalized transfer function to reconstruct $\hat{u}(x, t)$.

A symmetric range on frequencies spanning from $f \in [-60, 60]$ Hz is considered. The PGD solution in the form shown in Eq. (6) was achieved in only 9 terms, for a residual error of 0.01%. The first three space terms (modes) are shown in Fig. 2. The real part is placed in the left column while the imaginary part appears in the second one.

Fig. 3 shows the first three frequency modes. Each plot contains two curves, being the blue one the real part and the red one the imaginary part of each frequency mode. As expected, the real part is symmetric with respect the zero frequency, while the imaginary part is anti-symmetric. It can be noticed by looking at line markers that frequency mesh is not uniform, it is finner around the zero frequency because solution is quite sharp around this location while it becomes smoother as it gets away from the zero frequency.

After reconstructing the solution for certain frequencies (here we consider -60,0,60 Hz), we can compare the result against direct FEM computation. It is worthy to remind that once the generalized transfer function has been computed, extracting the solution for a particular frequency implies only a particularization, whereas a direct FEM computation requires solving a problem for each value of the frequency. Fig. 4 depicts the error map for both the real and imaginary parts, obtaining a good agreement.

Once the generalized transfer function has been verified, we use it for computing the transient response at the point $M$ when the load moves over the upper boundary, in the time interval $t \in (0, 1]$ of interest, as shown in Fig. 5. The reference solution has been computed by finite elements, using an explicit integration method. Notice that for the reference solution it is needed to solve at least as many problems as load positions are considered on the boundary, even if in the linear case this is not a major difficulty. Here we have considered 21 positions of the flux. Fig. 5 shows that both the PGD and the reference solutions are very close one another.

# 5 Conclusions

In this work we propose the use of an advanced harmonic technique for solving the linear transient heat equation. The technique here proposed does not present major advantages for fixed thermal sources, but when these sources move, the combination of this procedure and the reciprocity principle that only works in the frequency domain, allows very fast computations. Thus, two are the main ingredients : the reciprocity principle and the PGD method. However, let us suppose that the heat source presents a certain evolution along the time, which in a Fourier sense means that an arbitrary set of frequencies are being excited, and consequently the model has to be solved for each involved frequency. And here it is where the second ingredient shows up, because the PGD method allows an efficient computation
<figure><img src="image_2.png" /><figcaption>Fig. 2 – First three PGD space modes : real (left column) and imaginary (right column) parts</figcaption></figure>

of parametric solutions. This parametric solution can be seen as a sort of Computational Vademecum because it contains the whole necessary information, and thus, it can be computed once only and then the only work to do is a postprocessing which involves computationally cheap operations.

Références

[1] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings. A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of compex fluids. Part II : Tran-
<figure><img src="image_3.png" /><figcaption>Fig. 3 – First three PGD frequency modes : real (blue) and imaginary (red) parts</figcaption></figure>

sient simulation using space-time separated representation, Journal of Non-Newtonian Fluid Mechanics,
144, 98-121, 2007.

[2] F. Chinesta, P. Lavedeze, E. Cueto. *A short review in Model Order Reduction based on Proper Generalized Decomposition*, Archives of Computational Methods in Engineering, 18, 395-404, 2011.

[3] F. Chinesta, A. Leygue, F. Bordeu, J.V. Aguado, E. Cueto, D. Gonzalez, I. Alfaro, A. Ammar and A. Huerta : *PGD-based computational vademecum for efficient design, optimization and control*, Archives of Computational Methods in Engineering, 2012 (in press).

[4] A. Barbarulo, P. Ladeveze, H. Riou and L. Kovalesky : *Proper generalized decomposition applied to linear acoustic : a new tool for broad band calculation*, Journal of Sound and Vibration, 2012.

[5] A. Ammar, F. Chinesta, P. Diez and A. Huerta : *An error estimator for separated representations of highly multidimensional models*, Computer methods in Applied Mechanics and Engineering, 199, 1872-1880, 2010.

[6] P. Lavedèze, L. Chamoin : *On the verification of model reduction methods based on the proper generalized decomposition*, Computer methods in Applied Mechanics and Engineering, 200, 2032-2047, 2011.

[7] O.C. Zienkiewicz and R.L. Taylor : *The finite element method*, McGraw Hill, Vol. I., 1999.
<figure><img src="image_4.png" /><figcaption>Fig. 4 – PGD vs Direct FEM calculation : error map $|h - h_{ref}|$. Real part (left column) and imaginary part (right column) for frequencies -60 Hz (first row), 0 Hz (second row), 60 Hz (third row)</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 5 – Temperature evolution at point *M* during the heat flux path on the upper boundary</figcaption></figure>