# PGD in thermal transient problems with a moving heat source: A sensitivity study on factors affecting accuracy and efficiency

Dominic Strobl, Jörg F Unger, Chady Ghnatios, Annika Robens-radermacher

► To cite this version:

Dominic Strobl, Jörg F Unger, Chady Ghnatios, Annika Robens-radermacher. PGD in thermal transient problems with a moving heat source: A sensitivity study on factors affecting accuracy and efficiency. Engineering Reports, 2024, 6 (11), 10.1002/eng2.12887. hal-05308695

HAL Id: hal-05308695

https://hal.science/hal-05308695v1

Submitted on 10 Oct 2025

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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# PGD in thermal transient problems with a moving heat source: A sensitivity study on factors affecting accuracy and efficiency

Dominic Strobl<fnref n="1" /> ID | Jörg F. Unger<fnref n="1" /> ID | Chady Ghnatios<fnref n="2" /> | Annika Robens-Radermacher<fnref n="1" />

## Correspondence

Dominic Strobl, Modellierung und Simulation, Bundesanstalt für Materialforschung und-prüfung, Unter den Eichen 87, Berlin 12205, Germany. Email: dominic.strobl@bam.de

## Abstract

Thermal transient problems, essential for modeling applications like welding and additive metal manufacturing, are characterized by a dynamic evolution of temperature. Accurately simulating these phenomena is often computationally expensive, thus limiting their applications, for example for model parameter estimation or online process control. Model order reduction, a solution to preserve the accuracy while reducing the computation time, is explored. This article addresses challenges in developing reduced order models using the proper generalized decomposition (PGD) for transient thermal problems with a specific treatment of the moving heat source within the reduced model. Factors affecting accuracy, convergence, and computational cost, such as discretization methods (finite element and finite difference), a dimensionless formulation, the size of the heat source, and the inclusion of material parameters as additional PGD variables are examined across progressively complex examples. The results demonstrate the influence of these factors on the PGD model's performance and emphasize the importance of their consideration when implementing such models. For thermal example, it is demonstrated that a PGD model with a finite difference discretization in time, a dimensionless representation, a mapping for a moving heat source, and a spatial domain non-separation yields the best approximation to the full order model.

## KEYWORDS

additive manufacturing, mapping for unseparable load, model order reduction (MOR), proper generalized decomposition (PGD), sensitivity analysis, thermal transient problem

## 1 | INTRODUCTION

Numerical simulations, especially finite element simulations, are widely used in the industry to achieve savings in the product time-to-market and costs.<fnref n="1" /> Such simulations often require a large computational effort which sums up even more when a model calibration or optimization process is required. In inverse problems, a large number of evaluations of the forward model is required – both in deterministic but even more pronounced in stochastic MCMC type of identification

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2024 The Authors. *Engineering Reports* published by John Wiley & Sons Ltd.

<footnote n="1">¹Modellierung und Simulation, Bundesanstalt für Materialforschung und-prüfung, Berlin, Germany</footnote>

<footnote n="2">²PIMM Laboratory, Arts et Métiers Institute of Technology, CNRS, Cnam, HESAM University, Paris, France</footnote>
procedures. Therefore, methods to improve the performance of the simulations are needed. Some proposed methods can improve the simulation performance based on mesh improvement techniques. These use, for example, adaptive mesh algorithms, such as a moving mesh algorithm used in Reference 2, octree-based voxelized conforming meshes in Reference 3, but also element-free algorithms, such as the element-free Galerkin method in References 4 and 5. Furthermore, semi-analytical approaches as in References 6 and 7, CPU-based implementations as in Reference 8, and data-driven approaches as in References 9 and 10 are used. Otherwise, reduced order modeling (ROM) is a promising approach to significantly improve the computational performance. This popular concept can decrease the computation time per model evaluation, potentially achieving a real-time simulation.¹¹ However, as usual, validation using experimental reference data is absolutely necessary.¹² Furthermore, providing real-time simulations using model order reduction methods enables an improved process control and builds the foundation for a digital twin of industrial manufacturing processes.¹³

The goal of reduced order modeling is to find a lower dimensional representation of the original problem, without compromising the physical behavior of the analyzed system. Several “a priori” model order reduction approaches have been developed and used, such as hyper-reduction in References 14 and 15. A widely used approach is the proper generalized decomposition (PGD),¹³,¹⁶,¹⁷ for example in data-driven applications,¹⁸,¹⁹ parameter estimations,²⁰,²¹ surgery simulations²²,²³ or contact problems.²⁴,²⁵ The main idea of the a-priori PGD approach is to approximate the solution field by a separated representation given by a number of modes for each model parameter which are computed based on the underlying PDE of the problem. The PGD model parameters could be the space, time, boundary conditions or any other model parameter, like a material parameter which may vary in a certain range. With the derived PGD solution, all possible solution fields in the defined parameter space can be approximated. The dimension of each PGD mode is defined by the corresponding PGD model parameter dimension. Except from the physical space, those are usually 1D-spaces allowing a separation of the complex multi-dimensional problem into a series of simple (1D) problems. For an efficient computation of the PGD modes, the separation of the whole problem (for example including external loads) with respect to the chosen PGD model parameters is necessary. A typical example, where the derived PGD weak form is hard to separate, is a moving load.²⁶ In that regard, novel approaches are derived overcoming this limitation, as Ghnatios et al.²⁷ with a mapping approach, Rubio et al.²¹ with a moving coordinate system, Huerta et al.²⁸ with a domain decomposition strategy and Favoretto et al.²⁹ with an asymptotic expansion.

The prediction of the temperature evaluation over time is extremely important for various technical problems for example in welding applications or metal-based additive manufacturing. For those cases, a fast-to-evaluate prediction model is required for an efficient process control and optimization.³⁰ Especially in welding problems, steep temperature gradients over time and space as well as moving heat sources occur, leading to challenges in the development of an appropriate reduced order model.

The aim of this article is to demonstrate and study those challenges by developing a PGD model for the thermal transient heat equation with a moving heat source. Several aspects are investigated and compared, namely the discretization method for the time integration to reduce temperature oscillations, a reformulation of the problem in a dimensionless form to achieve better convergence, the influence of different heat source sizes to get an understanding of occurring changes, the extension of the PGD variables by material parameters (Fourier number) required for calibration tasks and the separation of the spatial domain with its influence on the accuracy. In the case of a moving heat source, the mapping approach of Ghnatios et al.¹¹,²⁷,³¹ with a coordinate transformation is applied and the influence of the mapping on the accuracy is analyzed. These investigations are carried out using examples that gradually increase in complexity. Starting with a 1D thermal transient problem with a stationary heat source over a moving source up to a 3D problem. This 3D problem considers constant material parameters as well as movement of the heat source in one direction for simplicity. However, the methodology shown can be extended to temperature dependent material parameters as well as a free movement of the source in any direction. The findings are summarized in the conclusion section.

# 2 | PGD MODEL FOR THERMAL TRANSIENT MODEL

## 2.1 | Thermal transient model

The investigated thermal transient problem with both a standing and moving heat source can be described by the heat equation

$$ \rho c_p \frac{\partial T}{\partial t} - k \Delta T = q, \qquad (1) $$
where $T(\mathbf{x}, t)$ is the temperature field and $q(\mathbf{x}, t)$ the heat source with $\mathbf{x} = [x \ y \ z]^T \in \Omega$, $\Omega \in \mathbb{R}^3$ (here $\Omega = [0, L_x] \times [0, L_y] \times [0, L_z]$) and time $t \in [0, t_{\text{end}}]$. The density $\rho$, specific heat capacity $c_p$ and thermal conductivity $k$ are kept constant. Furthermore, the commonly used Goldak heat source model³² as depicted in Figure 1 with a movement in x-direction is applied to simulate the heat source

$$q(\mathbf{x}, t) = \frac{6\sqrt{3}\eta P}{(a_f + a_r)bc\pi^{\frac{3}{2}}} \exp^{-3((x-x_0+vt))/a_f} \exp^{-3((y-y_0)/b)^2} \exp^{-3((z-z_0)/c)^2}, \quad (2)$$

with the arc efficiency $\eta$, the thermal power $P$, the starting position $(x_0, y_0, z_0)$, the velocity $v$ and the Goldak spatial distribution parameters $a_f$, $a_r$, $b$, and $c$. In the following, $a_f$ and $a_r$ are set to be equal, such that the shape of the source is axisymmetric with respect to the x-axis. A property of the Goldak source is that the energy input drops to 5% intensity at its boundaries. Using $v = 0$, a standing heat source at location $(x_0, y_0, z_0)$ can be simulated. The initial condition is assumed to be the ambient temperature given by:

$$T(x, y, z, t = 0) = T_{\infty}. \quad (3)$$

The convective heat transfer at the boundary $\partial\Omega$ writes

$$-k\nabla T_s \cdot \vec{n} = h(T_s - T_{\infty}) \quad (4)$$

with the surface temperature $T_s$ and the heat transfer coefficient $h$. After applying these boundary conditions and integrating by part, the weak form of Equation (1) reads

$$\int_{\Omega} \rho c_p T^* \frac{\partial T}{\partial t} + k \nabla T^* \cdot \nabla T \, dV + \int_{\partial \Omega} h T^* (T - T_{\infty}) \, dA = \int_{\Omega} T^* q \, dV. \quad (5)$$

Often, Equation (5) is then reformulated again to improve the convergence properties of the numerical methods by considering the space, time, and temperature in dimensionless form. Hence, defining

$$\xi_x := \frac{x}{L_{\text{ref}}}, \quad \xi_y := \frac{y}{L_{\text{ref}}}, \quad \xi_z := \frac{z}{L_{\text{ref}}}, \quad \tau := \frac{t}{t_{\text{ref}}}, \quad \text{and} \quad \Theta := \frac{T}{T_{\text{ref}}} \quad (6)$$

with arbitrary $L_{\text{ref}}$, $t_{\text{ref}}$, $T_{\text{ref}} \in \mathbb{R}^+$. With a special choice of the reference temperature $T_{\text{ref}}$ as $T_{\text{ref}} = L_{\text{ref}}^2/k$ and the dimensionless spatial domain, Fourier and Biot number

$$\Omega = \{ \xi = [\xi_x \ \xi_y \ \xi_z]^T \in \mathbb{R}^3 | \mathbf{x} \in \Omega \}, \quad \text{Fo} = \frac{kt_{\text{ref}}}{\rho c_p L_{\text{ref}}^2}, \quad \text{Bi} = \frac{hL_{\text{ref}}}{k} \quad (7)$$

<figure><img src="image_2.png" /><figcaption>FIGURE 1 Goldak heat source model with a local coordinate system $\mathbf{x}' = [x' \ y' \ z']^T$ moving in the global coordinate system $\mathbf{x}$.³⁰</figcaption></figure>
the number of parameters can be reduced, and the dimensionless weak form can be obtained as

$$ \int_{\Omega} \mathrm{Fo}^{-1} T^* \frac{\partial \Theta}{\partial \tau} + \nabla T^* \cdot \nabla \Theta \, dV + \mathrm{Bi} \int_{\partial \Omega} T^* \Theta \, dA = \int_{\Omega} T^* q \, dV + \mathrm{Bi} \int_{\partial \Omega} T^* T_{\infty} \, dA. \quad (8) $$

## 2.2 | Reference solution

As a reference solution, the solution of the full order model is obtained by applying the finite element method solving the weak form given in Equations (5) and (8) coupled with a backward Euler approach for the time integration scheme

$$ \frac{\partial \Theta^{m+1}}{\partial \tau} = \frac{\Theta^{m+1} - \Theta^m}{\Delta \tau}. \quad (9) $$

Depending on the dimension of the problem, linear 1D elements or second order tetrahedrons are used to discretize the physical space. Since the full order model is used to quantify the accuracy of the reduced PGD model, the mesh influence on the reference solution has to be considered. Either the same mesh discretization is used in the full order model as in the PGD model thus the influence on the reduction compared to the full order model is analyzed or a convergence study is used to define an appropriate mesh independent reference solution. In the convergence analysis, the mesh is successively refined and in each refinement step the solution is compared with the previous solution. The convergence study is only carried out for the full order finite element solution and the resulting mesh is subsequently applied in the PGD model, since the same behavior of the PGD solution and the full order solution is assumed.

## 2.3 | PGD model

Deriving a PGD model to handle Equation (1) requires a separate representation of the temperature field T. Besides the standard space and time parameters, the source efficiency $\eta \in [\eta_{\min}, \eta_{\max}]$ is additionally chosen as a PGD parameter, where its domain boundaries $\eta_{\min}$ and $\eta_{\max}$ represent the possible inputs for the PGD model. This choice is motivated by the fact that those parameters are uncertain, and are usually identified using measurement data. Different examples for deriving a PGD model for transient thermal problems can be found in References 29 and 33. With the above-mentioned parameter selection, T is approximated as

$$ T \approx T^n(\mathbf{x}, t, \eta) = \sum_{i=1}^{n} F_1^i(\mathbf{x}) F_2^i(t) F_3^i(\eta) + G_x(\mathbf{x}) G_t(t) G_\eta(\eta), \quad (10) $$

with PGD modes $F_i^j$. The additional summand G, separated for each PGD variable, is the inhomogeneous initial condition, which holds

$$ G_x(\mathbf{x}) = G_\eta(\eta) \equiv 1, \quad \text{and} \quad G_t(t) \equiv T_\infty. \quad (11) $$

Inserting this approach into Equation (5) leads to a nonlinear, multidimensional problem

$$ \int_{\Omega} \int_{0}^{t_{\mathrm{end}}} \int_{\eta_{\min}}^{\eta_{\max}} \rho c_p T^* \frac{\partial T^n}{\partial t} + k \nabla T^* \cdot \nabla T^n \, dV \, dt \, d\eta = \int_{\Omega} \int_{0}^{t_{\mathrm{end}}} \int_{\eta_{\min}}^{\eta_{\max}} T^* q \, dV \, dt \, d\eta. \quad (12) $$

As long as the heat source q can be separated into the PGD parameters in the same way the temperature is sought, the integration over q in Equation (12) can be separated, too. For a stationary heat source with $v = 0$, the heat source is given by

$$ q(\mathbf{x}, t, \eta) = \frac{6\sqrt{3P}}{(a_f + a_t)bc\pi^{\frac{3}{2}}} \underbrace{\eta \exp^{-3((x-x_0)/a_t)^2} \exp^{-3((y-y_0)/b)^2} \exp^{-3((z-z_0)/c)^2}}_{q(\eta)}, \quad (13) $$
which can be easily separated. If $v \neq 0$, additional treatment is required (see Section 2.3.3). The PGD weak form Equation (12) is solved iteratively with an alternated direction fixed-point algorithm, as for example described in Reference 34, with

$$ 
\begin{aligned}
T^n(\mathbf{x}, t, \eta) &= \sum_{i=1}^{n-1} F_1^i(\mathbf{x}) F_2^i(t) F_3^i(\eta) + \mathbf{X}(\mathbf{x}) K(t) M(\eta) + G_x(\mathbf{x}) G_t(t) G_\eta(\eta) \\
&= T^{n-1} + \mathbf{XKM} + G_x G_t G_\eta,
\end{aligned} 
\quad (14) $$

where the new modes $\mathbf{X}$, $\mathbf{K}$, and $\mathbf{M}$ are computed in each enrichment step based on the old solution $T^{n-1}$. The numerical solution of the occurring separate PDE problems in the PGD solver can be handled individually using, for example finite elements or finite difference approaches.

### 2.3.1 | Time integration schemes

Especially in cases with steep temperature gradients, the discretization scheme of the time parameter influences the solution. Two different time integration schemes are discussed, in other words the problem regarding the time $t$ is solved by using the finite element (FE) method and by using the finite difference (FD) method.

In FD, the domain of interest is divided into a grid, often uniform, where discrete points represent locations in the domain. The derivatives in the PDEs are approximated by finite difference equations at these grid points. Different possibilities to approximate derivatives exist, since the time problem is one-directional, the first-order upwind scheme is applied

$$ \frac{\partial u}{\partial t} = \frac{u_i^{n+1} - u_i^n}{\Delta t}. \quad (15) $$

It is well-known that oscillations can be reduced with such an upwind approach. The FE method divides the domain into elements, and it uses a piecewise approximation with variable polynomial degree of the solution over these elements. The PDEs are integrated over these elements to obtain a system of equations. The FE method is very flexible in terms of handling irregular geometries or higher dimensional problems. Since the time problem is a 1D problem, standard 1D finite elements with linear ansatz functions are used. Nevertheless, finite element formulations are not one-directional like the upwind approach. Values of both sides of the 1D element will be used in the computation.

The two discretizations are compared in Section 4.1.1 to determine which method is more suitable for solving the temporal problem.

### 2.3.2 | Dimensionless formulation

As mentioned above, the dimensionless formulation is often used for improving the convergence properties, in the standard numerical approaches. The dimensionless PGD approach obtained by inserting the definition of the dimensionless variables Equation (6) into the PGD approach Equation (14) reads

$$ 
\begin{aligned}
\Theta &\approx \Theta^n(\xi, \tau, \eta) = \sum_{i=1}^{n-1} F_1^i(\xi) F_2^i(\tau) F_3^i(\eta) + X(\xi) K(\tau) M(\eta) + G_x(\xi) G_t(\tau) G_\eta(\eta) \\
&= \Theta^{n-1} + \mathbf{XKM} + G_x G_t G_\eta.
\end{aligned} 
\quad (16) $$

Inserting the approach into the dimensionless weak form Equation (8) the dimensionless PGD problem is obtained as

$$ 
\begin{aligned}
&\int_{\Omega} \int_{0}^{t} \int_{\eta_{\min}}^{\eta_{\max}} \mathrm{Fo}^{-1} T^{*} \frac{\partial \Theta^{n}}{\partial \tau} + \nabla T^{*} \cdot \nabla \Theta^{n} \, dV \, d\tau \, d\eta + \int_{\partial \Omega} \int_{0}^{t} \int_{\eta_{\min}}^{\eta_{\max}} \mathrm{Bi} \, T^{*} \Theta^{n} \, dA \, d\tau \, d\eta \\
&= \int_{\Omega} \int_{0}^{t} \int_{\eta_{\min}}^{\eta_{\max}} T^{*} q \, dV \, d\tau \, d\eta + \int_{\partial \Omega} \int_{0}^{t} \int_{\eta_{\min}}^{\eta_{\max}} \mathrm{Bi} \, T^{*} T_{\infty} \, dA \, d\tau \, d\eta.
\end{aligned} 
\quad (17) $$
The PGD modes in the dimensionless form can be computed in the same way as discussed for the standard form assuming a separable heat source $q$. The influence of the dimensionless formulation on the simulation in terms of accuracy and speed is analyzed in Section 4.1.2.

### 2.3.3 | Mapping

The development of a PGD model for the parametric problem Equation (1) with a moving heat source is more complex. This is because – to put the formulation in a similar form to Equation (16) – the problem has to be separable for each PGD variable. But with $v \neq 0$, the Goldak heat source Equation (2) is not factorable into space direction $x$ and time $t$, due to the squared exponent. Ghnatios et al.²⁷ derived a method to handle such a hardly separable parametric problem. Here, the main idea is to transform the hardly separable problem by a coordinate transformation to a separable one. For the sake of simplicity, the mapping is derived directly using the dimensionless formulation, in other words the dimensionless space and time domain $(\xi_x, \tau)$ is mapped into a $(s, r)$ domain in which the moving Goldak heat source function is separable into $r$ and $s$ demonstrated in Figure 2. The variable $r$ describes the heat source position using $\xi_{x_0} = x_0/L_{ref}$

$$r := \xi_{x_0} + \frac{t_{\text{ref}}}{L_{\text{ref}}} v \tau \quad \Leftrightarrow \quad \tau = \frac{r - \xi_{x_0}}{\frac{t_{\text{ref}}}{L_{\text{ref}}} v}, \quad r \in \left[ \hat{r}_{\text{on}} = \frac{r_{\text{on}}}{L_{\text{ref}}}, \hat{r}_{\text{off}} = \frac{r_{\text{off}}}{L_{\text{ref}}} \right]. \quad (18)$$

The activation and deactivation positions $r_{on}$ and $r_{off}$ for the load are the points where the heat source is switched on and off. These points are mathematically restricted by $r_{on} > a_f$ and $r_{off} < L - a_f$ to avoid division by zero in the mapping. This means that a heat source that is switched on or off at the edge of the spatial domain cannot be simulated without artificially expand the domain. Furthermore, the mapped spatial domain $s \in [s_{\min}, s_{\max}] = [0, 3]$ is split into three parts. These three parts comprise the domain behind the heat source $s \in [0, 1)$, the heat source itself $s \in [1, 2]$ and the domain in front of the head source $s \in (2, 3]$. The dimensionless mapping parameter $s$ is defined as

$$\begin{align}
s \in [0, 1] \quad : s &:= \frac{\xi_x}{h_1(r)} & \Leftrightarrow \xi_x &= sh_1(r), \nonumber \\
s \in [1, 2] \quad : s &:= \frac{\xi_x - h_1(r)}{h_g} + 1 & \Leftrightarrow \xi_x &= (s-1)h_g + h_1(r), \tag{19} \\
s \in (2, 3] \quad : s &:= \frac{\xi_x - h_g - h_1(r)}{h_2(r)} + 2 & \Leftrightarrow \xi_x &= (s-2)h_2(r) + h_g + h_1(r), \nonumber
\end{align}$$

with

$$\hat{a}_f = \frac{a_f}{L_{\text{ref}}}, \quad \hat{a}_r = \frac{a_r}{L_{\text{ref}}}, \quad \hat{b} = \frac{b}{L_{\text{ref}}}, \quad \hat{c} = \frac{c}{L_{\text{ref}}} \quad (20)$$

<figure><img src="image_3.png" /><figcaption>FIGURE 2 Mapping from the dimensionless space $T(\xi_x, \tau)$ (left) to the mapped space $T(s, r)$ (right) adapted from Reference.³⁰</figcaption></figure>
and

$$h_g = \hat{a}_f + \hat{a}_r, \quad h_1(r) = r - \hat{a}_r, \quad \text{and} \quad h_2(r) = \frac{L_x}{L_{\text{ref}}} - \hat{a}_f - r \qquad (21)$$

describing the load position. Following the mapping approach, the PGD parameters are now $\mathbf{s} = [s \ \xi_y \ \xi_z]^T$, $r$, and $\eta$ and the mapped PGD approach is given by

$$\begin{aligned} \Theta &\approx \Theta^n(\mathbf{s}, r, \eta) = \sum_{i=1}^{n} F_1^i(\mathbf{s}) F_2^i(r) F_3^i(\eta) + G_s(\mathbf{s}) G_r(r) G_\eta(\eta) \\ &= \sum_{i=1}^{n-1} F_1^i(\mathbf{s}) F_2^i(r) F_3^i(\eta) + \mathbf{S}(\mathbf{s}) R(r) M(\eta) + G_s(\mathbf{s}) G_r(r) G_\eta(\eta) \\ &= \Theta^{n-1} + \mathbf{S} \mathbf{R} \mathbf{M} + G_s G_r G_\eta. \end{aligned} \qquad (22)$$

Applying the coordinate transformation from the dimensionless space to the dimensionless mapped space given by

$$dV_{\xi} d\tau d\eta = \det(J) dV_s dr d\eta, \qquad (23)$$

with the Jacobian matrix $J$ and the mapped PGD approach Equation (22) to the dimensionless weak form Equation (17) yields the mapped dimensionless weak form

$$\begin{aligned} & \int_{\mathfrak{A}} \int_{\hat{r}_{\text{on}}}^{\hat{r}_{\text{off}}} \int_{\eta_{\text{min}}}^{\eta_{\text{max}}} \mathrm{Fo}^{-1} T^* \mathbf{B}_\tau \nabla_{sr} \Theta^n \det(\mathbf{J}) + \nabla_{\xi r} T^* \cdot \nabla_{\xi r} \Theta^n \det(\mathbf{J}) \, dV \, dr \, d\eta + \int_{\partial \mathfrak{A}} \int_{\hat{r}_{\text{on}}}^{\hat{r}_{\text{off}}} \int_{\eta_{\text{min}}}^{\eta_{\text{max}}} \mathrm{Bi} \, T^* \Theta^n \det(\mathbf{J}) \, dA \, dr \, d\eta \\ &= \int_{\mathfrak{A}} \int_{\hat{r}_{\text{on}}}^{\hat{r}_{\text{off}}} \int_{\eta_{\text{min}}}^{\eta_{\text{max}}} T^* q \det(\mathbf{J}) \, dV \, dr \, d\eta + \int_{\partial \mathfrak{A}} \int_{\hat{r}_{\text{on}}}^{\hat{r}_{\text{off}}} \int_{\eta_{\text{min}}}^{\eta_{\text{max}}} \mathrm{Bi} \, T^* T_\infty \det(\mathbf{J}) \, dA \, dr \, d\eta, \end{aligned} \qquad (24)$$

with the dimensionless mapped space $\mathfrak{A} := [s_{\text{min}}, s_{\text{max}}] \times [0, L_y/L_{\text{ref}}] \times [0, L_z/L_{\text{ref}}]$ and

$$\mathbf{B}_{\tau} := \left[ \frac{\partial s}{\partial \tau} \quad \frac{\partial r}{\partial \tau} \right], \quad \mathbf{B}_{\xi_x} := \left[ \frac{\partial s}{\partial \xi_x} \quad \frac{\partial r}{\partial \xi_x} \right], \quad \nabla_{sr} \Theta := \left[ \frac{\partial \Theta}{\partial s} \quad \frac{\partial \Theta}{\partial r} \right] \quad \text{and} \quad \nabla_{\xi r} \Theta := \begin{bmatrix} \mathbf{B}_{\xi_x} \nabla_{sr} \Theta \\ \frac{\partial \Theta}{\partial \xi_y} \\ \frac{\partial \Theta}{\partial \xi_z} \end{bmatrix}. \qquad (25)$$

Where, due to the mapping the heat source is separable in $(\mathbf{s}, r, \eta)$

$$q(\mathbf{s}, r, \eta) = \frac{6\sqrt{3}P}{(\hat{a}_f + \hat{a}_r)\hat{b}\hat{c}\pi^{\frac{3}{2}}} \underbrace{\eta}_{q(r)} \underbrace{\exp^{-3((s-1)h_g - \hat{a}_r)/\hat{a}_t} \exp^{-3((\xi_y - \xi_{y0})/\hat{b})^2} \exp^{-3((\xi_z - \xi_{z0})/\hat{c})^2}}_{q(\mathbf{s})}, \qquad (26)$$

with $\xi_{y0} = y_0/L_{\text{ref}}$ and $\xi_{z0} = z_0/L_{\text{ref}}$. This mapping approach has another benefit when it comes to the discretization of the heat source. Due to the splitting of the spatial domain $\xi_x$ into three parts, the heat source itself is automatically refined adaptively.

A more detailed example on how to apply this mapping can be found in Reference 30. The effect of the mapping approach on the solution is discussed in Section 4.2.2.

### 2.3.4 | Material parameter as PGD variable

Due to the concept of PGD, other model parameters like material parameters can be easily added as PGD parameters allowing an efficient computation of the solution field with regard of its value.35 Having material parameters as PGD
variables enables, for example its identification based on measured data in a calibration process. Moreover, the resulting PGD model can directly handle different material parameters without resolving any equation systems. Using the dimensionless formulation, the original material parameters $ρ, c_p, k$ and $h$ are reduced to the dimensionless model parameters Fourier and Biot number. Here, the Fourier number $Fo$ and Biot number $Bi$ are added as PGD variables to the model to be able to adapt to any situation in combination with the heat source efficiency $η$. More specifically, scaling parameters $f$ for the Fourier number and $b$ for the Biot number are added, that is $Fo = f \cdot Fo_{ref}, f \in [f_{min}, f_{max}]$ and $Bi = b \cdot Bi_{ref}, b \in [b_{min}, b_{max}]$ with fixed values $Fo_{ref}, Bi_{ref} \in \mathbb{R}^+$ based on prior knowledge. Extending the PGD approach as in Equation (22)

$$ \Theta \approx \Theta^n(s, r, \eta, f, b) = \sum_{i=1}^{n} F_1^i(s) F_2^i(r) F_3^i(\eta) F_4^i(f) F_5^i(b) + G_s(s) G_r(r) G_\eta(\eta) G_f(f) G_b(b) \quad (27) $$

and inserting that into Equation (24) yields

$$ \begin{aligned} & \int_{\mathfrak{A}} \int_{\hat{r}_{on}}^{\hat{r}_{off}} \int_{\mathfrak{U}} \mathrm{Fo}^{-1} T^* \mathbf{B}_r \nabla_{sr} \Theta^n \det(\mathbf{J}) + \nabla_{\xi_r} T^* \cdot \nabla_{\xi_r} \Theta^n \det(\mathbf{J}) \, dV \, dr \, d\eta \, d\hat{f} \, dB + \\ & \qquad \int_{\partial \mathfrak{A}} \int_{\hat{r}_{on}}^{\hat{r}_{off}} \int_{\mathfrak{U}} \mathrm{Bi} \, T^* \Theta^n \det(\mathbf{J}) \, dA \, dr \, d\eta \, d\hat{f} \, dB \\ = & \int_{\mathfrak{A}} \int_{\hat{r}_{on}}^{\hat{r}_{off}} \int_{\mathfrak{U}} T^* q \det(\mathbf{J}) \, dV \, dr \, d\eta \, d\hat{f} \, dB + \\ & \qquad \int_{\partial \mathfrak{A}} \int_{\hat{r}_{on}}^{\hat{r}_{off}} \int_{\mathfrak{U}} \mathrm{Bi} \, T^* T_\infty \det(\mathbf{J}) \, dA \, dr \, d\eta \, d\hat{f} \, dB, \end{aligned} \quad (28) $$

which is simplified by defining $\mathfrak{U} := [\eta_{min}, \eta_{max}] \times [f_{min}, f_{max}] \times [b_{min}, b_{max}]$. Note, the extension of the PGD approach for material parameters like the Fourier number can also be done without applying the mapping approach using Equation (17). This extended model is then compared in Section 4.2.3 to reference solutions to analyze the accuracy of the PGD model for various material parameter values.

### 2.3.5 | Spatial domain splitting

The spatial domain plays a special role in PGD because it can be held together as a usual 3D space or decomposed into a 2D and a 1D space or even into three 1D spaces. Whether such a decomposition is possible depends mainly on the given geometry. Decomposing the spatial domain leads to tremendous savings in computation time, but potentially accompanied by a reduction in the accuracy of the model. Applying the PGD approach with a spatial domain divided into three 1D intervals

$$ \Theta \approx \Theta^n(s, \xi_y, \xi_z, r, \eta, f, b) = \sum_{i=1}^{n} F_1^i(s) F_2^i(\xi_y) F_3^i(\xi_z) F_4^i(r) F_5^i(\eta) F_6^i(f) F_7^i(b) + G_s(s) G_{\xi_y}(\xi_y) G_{\xi_z}(\xi_z) G_r(r) G_\eta(\eta) G_f(f) G_b(b) \quad (29) $$

to Equation (28) yields the dimensionless mapped weak form with decomposed spatial domain

$$ \begin{aligned} & \int_{\mathfrak{A}} \int_{\hat{r}_{on}}^{\hat{r}_{off}} \int_{\mathfrak{U}} \mathrm{Fo}^{-1} T^* \mathbf{B}_r \nabla_{sr} \Theta^n \det(\mathbf{J}) + \nabla_{\xi_r} T^* \cdot \nabla_{\xi_r} \Theta^n \det(\mathbf{J}) \, ds \, d\xi_y \, d\xi_z \, dr \, d\eta \, d\hat{f} \, dB \\ & \qquad + \int_{\partial \mathfrak{A}} \int_{\hat{r}_{on}}^{\hat{r}_{off}} \int_{\mathfrak{U}} \mathrm{Bi} \, T^* \Theta^n \det(\mathbf{J}) \, dA \, dr \, d\eta \, d\hat{f} \, dB \\ = & \int_{\mathfrak{A}} \int_{\hat{r}_{on}}^{\hat{r}_{off}} \int_{\mathfrak{U}} T^* q \det(\mathbf{J}) \, ds \, d\xi_y \, d\xi_z \, dr \, d\eta \, d\hat{f} \, dB + \\ & \qquad \int_{\partial \mathfrak{A}} \int_{\hat{r}_{on}}^{\hat{r}_{off}} \int_{\mathfrak{U}} \mathrm{Bi} \, T^* T_\infty \det(\mathbf{J}) \, dA \, dr \, d\eta \, d\hat{f} \, dB. \end{aligned} \quad (30) $$

Assuming the same separation as in Equation (29) for the heat source $q$, Equation (30) can be split into one-dimensional problems for each parameter $(s, \xi_y, \xi_z, r, \eta, f)$ in the fixed-point algorithm problem allowing an efficient way to compute the PDG modes involving only 1D integrals. As already discussed, the numerical solution of the single one-dimensional problems in the PGD solver can be handled individually using, for example finite elements or finite difference approaches. The influence of the splitting combined with different discretization methods on the accuracy of the solution is discussed in Section 4.3.
## 2.4 | Error measurement

Each model comparison of two models M1 and M2 is carried out by calculating the mean of the relative error for $\eta \in E = \{0.2, 0.4, 0.6, 0.8, 1\}$ over space $x$ at a given time point $t_{fixed}$

$$ \epsilon_x = \frac{1}{|E|} \sum_{\eta \in E} \frac{\|T_{M1}(x, t_{fixed}, \eta) - T_{M2}(x, t_{fixed}, \eta)\|_{L^2(x)}}{\|T_{M2}(x, t_{fixed}, \eta)\|_{L^2(x)}}, \quad (31) $$

or over time $t$ at a given space point $x_{fixed}$

$$ \epsilon_t = \frac{1}{|E|} \sum_{\eta \in E} \frac{\|T_{M1}(x_{fixed}, t, \eta) - T_{M2}(x_{fixed}, t, \eta)\|_{L^2(t)}}{\|T_{M2}(x_{fixed}, t, \eta)\|_{L^2(t)}}, \quad (32) $$

where $|E|$ represents the cardinality of $E$. If either model M1 or M2 is a PGD model, the error depends on the number of modes $n$ used to compute the PGD model solution. In the case of a moving heat source, the error is affected by the mapping approach resulting in the relative error over the mapped spatial domain $s$ at a given heat source position $r_{fixed}$

$$ \epsilon_s = \frac{1}{|E|} \sum_{\eta \in E} \frac{\|T_{M1}(s, r_{fixed}, \eta) - T_{M2}(s, r_{fixed}, \eta)\|_{L^2(s)}}{\|T_{M2}(s, r_{fixed}, \eta)\|_{L^2(s)}}, \quad (33) $$

and over the heat source position $r$ at a given mapped space point $s_{fixed}$

$$ \epsilon_r = \frac{1}{|E|} \sum_{\eta \in E} \frac{\|T_{M1}(s_{fixed}, r, \eta) - T_{M2}(s_{fixed}, r, \eta)\|_{L^2(r)}}{\|T_{M2}(s_{fixed}, r, \eta)\|_{L^2(r)}}, \quad (34) $$

These errors are necessary for an analysis in the mapped space.

In addition, the relative absolute error with a fixed position $x_{fixed}$, fixed time $t_{fixed}$ and fixed energy input $\eta_{fixed}$

$$ \epsilon = \frac{|T_{M1}(x_{fixed}, t_{fixed}, \eta_{fixed}) - T_{M2}(x_{fixed}, t_{fixed}, \eta_{fixed})|}{|T_{M2}(x_{fixed}, t_{fixed}, \eta_{fixed})|} \quad (35) $$

is introduced for a detailed analysis of the temperature on the plate.

The models $M_1$ and $M_2$ in the error measures $\epsilon_s$, $\epsilon_r$, and $\epsilon$ will later be extended with the material parameters Fourier and Biot number as variables, which are treated as fixed parameters in the error, denoted as $f_{fixed}$ and $b_{fixed}$.

## 3 | NUMERICAL IMPLEMENTATION

The whole code is implemented in python-version 3.8. The PGD module PGDrome developed by the authors and published as open source package 36 is used. The finite difference approach with the necessary matrices is implemented in this package. As the finite element tool, it uses FEniCS 37-version 2019.1. All full order simulations are done using FEniCS. The simulations are performed on an Ubuntu system with 8 cores and 7.7 GB of memory.

## 4 | NUMERICAL RESULTS AND DISCUSSION

In this section, the challenges described in Section 2 are discussed and examined using three examples of the introduced thermal transient problem with increasing complexity. Thereby, a plate—in 1D or 3D—is loaded by a standing or moving heat source. The parameters that are kept constant throughout all examples are listed in Table 1. To illustrate the solution with realistic material parameters, the material parameters of 1.4404 (AISI 316L) are chosen – an austenitic stainless steel often used in the chemical industry. The influence of these material parameters on temperature is neglected. A realistic average thermal power $P$ of the heat source for these material parameters is 5250 W as in Reference 30, but in 1D case
<table><caption>TABLE 1 Parametric problem's parameters.</caption><tr><td>Welding velocity v</td><td>0.01 m s<sup>-1</sup></td></tr><tr><td>Plate length L<sub>x</sub></td><td>0.1 m</td></tr><tr><td>Initial x-position x<sub>0</sub></td><td>0 m</td></tr><tr><td>Activation time t<sub>on</sub></td><td>2 s</td></tr><tr><td>Deactivation time t<sub>off</sub></td><td>8 s</td></tr><tr><td>Ambient temperature T<sub>∞</sub></td><td>25°C</td></tr><tr><td>Density ρ</td><td>7950 kg/m<sup>3</sup></td></tr><tr><td>Specific heat capacity c<sub>p</sub></td><td>475 J/(kgK)</td></tr><tr><td>Heat conductivity k</td><td>14.3 W/(mK)</td></tr><tr><td>Average thermal power P</td><td>2625 W</td></tr><tr><td>Activation position r<sub>on</sub></td><td>0.02 m</td></tr><tr><td>Deactivation position r<sub>off</sub></td><td>0.08 m</td></tr><tr><td>Last time point t<sub>end</sub></td><td>10 s</td></tr><tr><td>Reference length L<sub>ref</sub></td><td>0.1 m</td></tr><tr><td>Reference temperature T<sub>ref</sub></td><td>1372°C</td></tr><tr><td>Reference time t<sub>ref</sub></td><td>10 s</td></tr><tr><td>Minimal efficiency η<sub>min</sub></td><td>0</td></tr><tr><td>Maximal efficiency η<sub>max</sub></td><td>1.1</td></tr><tr><td>Number of elements in η-mesh</td><td>100</td></tr></table>

studies this value is reduced to 2625 W. This choice is due to the lack of a heat flux in y- and z-direction. Thus, this choice leads to an overall temperature in the 1D case similar to the full 3D case. Furthermore, the one-dimensional problems of the proper generalized decomposition method are solved here using first-order linear finite elements.

## 4.1 Case study 1: 1D PGD model with fixed heat source position

At first, the problem with a stationary heat source is analyzed. For the sake of simplicity, this parametric problem is considered in 1D and sketched in Figure 3. The x-mesh consists of 600 and the t-mesh of 100 linear elements. The temperature distribution of the reference solution over space at 8 s and over time at 0.05 m with heat source size $a_f = 0.005$ m and heat source efficiency $\eta = 1$ is depicted in Figure 4. Very steep temperature gradients in space can be observed with a peak temperature that propagates in time (with the heat source). The resolution of those steep gradients is challenging for numerical methods. Thus, the choice of the applied numerical method is of utmost importance and the effects of different time integration schemes are discussed first.

### 4.1.1 Time integration schemes

Two different PGD solutions, one per time integration scheme (FE and FD) as introduced in Section 2.3.1, are computed, and compared with a reference full order FE solution with backward Euler time integration, which uses the same x- and t-mesh as the reduced order model. Furthermore, the size of the Goldak heat source is kept constant at $a_f = 0.005$ m for this analysis and the computations are done with the non-dimensionless formulation. The heat source is positioned at $x = 0.05$ m. The relative error over space at fixed time $t = 8$ s and over time at fixed space $x = 0.05$ m is depicted in Figure 5 for both PGD solution approaches. On the left, the mean relative error over space and on the right, the mean relative error over time for increasing number of PGD modes is shown. In both cases, the PGD solution with FD in time (FE for all other PGD parameters) shows a better convergence behavior to the full order solution than the PGD solution
<figure><img src="image_4.png" /><figcaption>FIGURE 3 Boundary and geometry conditions of the problem with a fixed heat source.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4 Temperature distribution over space at 8 s (left) and over time at 0.05 m (right) with $a_f = 0.005$ m and $\eta = 1$.</figcaption></figure>

with FE in time. After a certain number of modes, the error of the latter case does not further decrease, since the PGD solution oscillates around the full order solution and therefore the convergence stagnates. Moreover, the relative absolute error in space at 8 s and in time at 0.05 m with $a_f = 0.005$ m and $\eta = 1$ is depicted in Figure 6 using 30 PGD modes. This shows that the solution using FE in time oscillates much more than the solution using FD in time, right before the steep temperature gradient occurs. Hence, the faster convergence of the PGD solution with FD in time can be seen. Thus, from now on, all PGD solutions are computed with the more accurate FD approach in time.

### 4.1.2 | Dimensionless formulation

Next, the effect of the dimensionless approach Section 2.3.2 is analyzed. Hence, a non-dimensionless and dimensionless PGD solution is computed and compared to the same reference non-dimensionless full order FE solution as above. Again,
<figure><img src="image_6.png" /><figcaption>FIGURE 5 Error over space $\epsilon_x$ at fixed time $t = 8$ s per mode of the PGD solution to the reference full order FE solution for different time integration schemes (left) and the corresponding error over time $\epsilon_t$ at fixed space $x = 0.05$ m (right).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6 Relative absolute error $\epsilon$ in space at 8 s (left) and in time at 0.05 m (right) with $a_f = 0.005$ m and $\eta = 1$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIGURE 7 Error over space $ε_x$ at 8 s (left) per mode of the PGD solution to the reference FE solution for the dimensionless formulation and non-dimensionless formulation and the corresponding error over time $ε_t$ at 0.05 m (middle), as well as the sum of used fixed-point iterations up to a given mode to reach the threshold of $10^{-9}$ (right).</figcaption></figure>

the size of the Goldak heat source is kept constant at $a_f = 0.005$ m and the convergence threshold of the fixed-point iteration per mode is $10^{-9}$. The comparison of both approaches over space at 8 s and over time at 0.05 m is depicted in Figure 7. The error plots show that the actual error of the models is the same per mode. This is expected since the fixed-point iteration in both PGD models converged and the dimensionless formulation does not change the PDE. But the difference is the number of fixed-point iterations needed to converge per mode. The third plot shows that the dimensionless solution requires less fixed-point iterations overall, thus converges faster and saves computation time. In addition, the dimensionless approach reduces the number of parameters when taking these into account, which in turn leads to savings in computation time. Therefore, from now on, all PGD solutions are computed using the dimensionless formulation.

### 4.1.3 | Heat source size

Finally, the effect of different heat source sizes on the temperature field is studied using the findings of the previous analyses, in other words, using the dimensionless problem formulation and FD as time discretization whereas the other parameters are discretized via FE. For this, the property of a symmetric heat source is kept, that is $a_f$ and $a_r$ are always the same. Here, the size $a_f$ is chosen as $a_f \in \{0.002 \text{ m}, 0.005 \text{ m}, 0.009 \text{ m}\}$. The overall energy input into the system is kept the same for each $a_f$ by changing the thermal power $P$ respectively. Thus, the smaller the size of the heat source the higher the temperature rises, which leads to steeper gradients for smaller $a_f$. This behavior leads to slightly larger inaccuracies, which is reflected in the errors compared to the reference full order solution as depicted in Figure 8 over space at 8 s (left) and over time at 0.05 m (right). Thus, for smaller $a_f$ more modes are necessary to achieve the same accuracy. Overall, the PGD solution shows good approximation results with errors smaller than $10^{-3}$ for each investigated heat source size.

## 4.2 | Case study 2: 1D PGD model with moving heat source

In a second example, the problem with a moving heat source is analyzed. This parametric problem is also considered in 1D and sketched in Figure 9. In this whole subsection, the results of the previous Section 4.1 are taken into account.
<figure><img src="image_9.png" /><figcaption>FIGURE 8 Mean relative error over space $\epsilon_x$ at 8 s per mode of the PGD solution to the reference full order FE solution for different heat source sizes $a_f$ (left) and the corresponding error over time $\epsilon_t$ at 0.05 m (right).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9 Boundary and geometry conditions of the problem with a moving heat source.</figcaption></figure>

Thus, all the computations are performed with the finite difference upwind scheme in time considering the dimensionless formulation Equation (24) and using $a_f = 0.005$ m as the heat source size.

Applying the mapping approach to handle the moving heat source as discussed in Section 2.3.3 leads to a PGD solution which converges to a solution in the mapped space. But this does not need to be the same as the reference solution in the original space due to different meshes. Therefore, a mesh convergence study in both spaces is performed for the full order model.

### 4.2.1 | Mesh convergence

The analysis in the mapped space compares the full order FE solution in each refinement step to the previous solution. Here, the spaces $s$ and $r$ are analyzed separately. Thus, for the $s$-space analysis, the initial discretization has 6 elements in
<figure><img src="image_11.png" /><figcaption>FIGURE 10 Mesh convergence analysis in the mapped space of the *s*-mesh (left) and *r*-mesh (right) over heat source positions $\epsilon_r$ at *s* = 1.2. The full order FE solution in each refinement step is compared to the previous solution.</figcaption></figure>

the *s*-mesh, which are doubled in each refinement step, and 2560 elements in the *r*-mesh, which are fixed. In the *r*-space analysis the initial discretization has 10 elements in the *r*-mesh, which are doubled in each refinement step, and 480 elements in the *s*-mesh, which are fixed. For both spaces, eleven refinement steps are performed. The results of this study in the mapped space are depicted in Figure 10, which shows the errors over the heat source position *r* at *s* = 1.2 for the successively refined solutions to each other. It can be observed that in both spaces the solution converges nicely. The error in the *r*-mesh analysis reaches $10^{-4}$ and in the *s*-mesh analysis $10^{-7}$ after eleven refinement steps. Hence, a mesh with 6144 elements in the *s*-mesh and 2560 elements in the *r*-mesh is considered as converged and used later on.

In the mesh convergence study in the real space, both the *x*- and *t*-spaces are again considered separately. Moreover, the used *x*-space discretization is not equidistant. A refined grid around the current position of the heat source is used as automatically induced by the mapping approach. The initial discretizations are based on the converged solution in the mapped space. Starting from this, the meshes are both coarsened and refined three times to analyze which discretization should be applied in the real space. For each of these operations, the resulting full order FE solution is directly compared to the converged solution in the mapped space. This leads to an error representing how good the solutions in both the mapped and real space coincide with each other. The results of this study are depicted in Figure 11, which shows the errors over the time *t* at *s* = 1.2. It can be observed that the error in both cases is of magnitude $10^{-3}$ and that further refinements are not necessary. This though is because the refinement given by the mapping approach is applied in the real space. Besides that, using a finer temporal resolution mesh leads to problems in the solution, such that the resulting error does not decrease further. The full order FE solution in the real space is considered as converged with 6144 elements in the *x*-mesh using a refinement around the heat source and 2560 elements in the *t*-mesh.

Since the PGD solution approximates the full order solution in the mapped space, the error of the PGD solution to the full order solution in the real space is limited by the error of the full order solutions to each other. Hence, the magnitude of the error between the PGD solution and the full order FE solution in the real space reaches at most $10^{-3}$, when a sufficient number of modes is taken into account in the PGD solution.

### 4.2.2 | Mapping

With the information of an optimal mesh discretization, a PGD solution is computed with the PGD variables *s*, *r* and *η* (see Section 2.3.3). This model consists of 70 modes computed within a fixed-point algorithm using a convergence threshold of $10^{-10}$ and a maximum of 150 iterations per mode. The 1D problems regarding *s* and *η* are solved with first order finite elements, whereas the *r* problem is solved with the FD upwind scheme. For each number of modes up to 70, the resulting PGD solution is then compared to a reference full order FE solution in the mapped space. The relative errors over *s* at source position 0.064 m and over *r* at *s* = 1.2 are plotted in Figure 12. It can be observed that the error decreases using more modes in the typical non-monotonically behavior of the PGD approach. After reaching 70 modes, this error is around $10^{-2}$.
<figure><img src="image_12.png" /><figcaption>FIGURE 11 Mesh convergence analysis in the real space of the x-mesh (left) and t-mesh (right) over time $\epsilon_t$ at $s = 1.2$. The full order FE solution in each refinement step is compared to the converged solution in the mapped space.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 12 Error over $s \epsilon_s$ at source position 0.064 m (left) per mode of the PGD solution to the reference FE solution in the mapped space and the corresponding error over $r \epsilon_r$ at $s = 1.2$ (right).</figcaption></figure>

The difference in accuracy between stationary and moving cases in Figures 5 and 12 arises from the greatly increased complexity of the problem as soon as the heat source moves. Therefore, a lower accuracy than before was to be expected even with an increased total number of modes.

### 4.2.3 | Material parameter as PGD variable

Finally, the factor of the Fourier number $f$ is added as a PGD variable to the model as described in Section 2.3.4. The new $f$-space [$f_{min}, f_{max}$] = [0.5, 1.5] is discretized equidistantly using 100 linear finite elements. Once more, a
<figure><img src="image_14.png" /><figcaption>FIGURE 13 Error over $s \in \{0.6, 0.8, 1, 1.2, 1.4\}$ at source position 0.064 m (left) and over $r \in \{0.8, 1.0, 1.2, 1.4\}$ at $s = 1.2$ (right).</figcaption></figure>

PGD solution with the same settings as before is computed. This solution is evaluated at $\mathfrak{f} \in \{0.6, 0.8, 1, 1.2, 1.4\}$ and compared to reference full order FE solutions with the corresponding factor $\mathfrak{f}$ to show that the PGD model can adapt correctly to different material parameters. The results are depicted in Figure 13, which shows the error over $s$ at heat source position 0.064 m and over $r$ at $s = 1.2$. It can be observed that the error is decreasing slower as in Figure 12 but still reaching a suitable accuracy for, for example, a model calibration procedure. That is because the PGD model became more complex with the additional PGD variable $\mathfrak{f}$, hence needs more modes to achieve the same accuracy. When using only a few modes, the error is slightly increasing with decreasing $\mathfrak{f}$, since the maximal temperature increases and with it the maximal temperature gradient. Though, this effect can be compensated with more modes. Thus, the extended model is able to capture a change in the material parameters given by the Fourier number without the need of performing new computations. In that way, the derived PGD model can be used for efficient material parameter estimation with measured data.

## 4.3 | Case study 3: 3D PGD model with moving heat source

Lastly, the 3D parametric problem is considered, which is sketched in Figure 14. The parameters for this analysis are listed in Table 2. Two different PGD models are considered. In the first PGD model, the spatial domain is divided into three 1D intervals, where the number of elements of the meshes is chosen as $s : 800, \xi_y : 200, \xi_z : 20$. In this case, the problems regarding space and time are solved with FD and 100 modes are computed. A contour plot of the corresponding PGD solution at 8 s with $\eta = 0.8, \mathfrak{f} = 1, \mathbf{b} = 1$ and 100 modes, mapped back to the physical space $(x, y, z)$, is shown in Figure 15. In the second PGD model, the spatial domain is kept together as a 3D domain, where the number of elements of the mesh is chosen as $s : 90, \xi_y : 20, \xi_z : 4$. In this case, the problem regarding space is solved with second order FE and 40 modes are computed. These models are compared to a reference full order FE solution, where the number of elements of the space mesh is chosen as $s : 90, \xi_y : 20, \xi_z : 4$, with quadratic elements, $\eta = 0.8, \mathfrak{f} = 1$ and $\mathbf{b} = 1$. The relative error of the PGD solutions to the reference solution mapped back to the original space up to 20% error is depicted in Figure 16
<figure><img src="image_15.png" /><figcaption>FIGURE 14 Derived parametric problem for a moving heat source in the 3D-case.</figcaption></figure>

<table><caption>TABLE 2 Parameters for the 3D parametric problem.</caption><tr><td>Plate length L<sub>x</sub></td><td>0.1 m</td></tr><tr><td>Plate width L<sub>y</sub></td><td>0.05 m</td></tr><tr><td>Plate height L<sub>z</sub></td><td>0.008 m</td></tr><tr><td>Initial x-position x<sub>0</sub></td><td>0 m</td></tr><tr><td>Initial y-position y<sub>0</sub></td><td>0.025 m</td></tr><tr><td>Initial z-position z<sub>0</sub></td><td>0.008 m</td></tr><tr><td>Goldak front length parameter a<sub>f</sub></td><td>0.005 m</td></tr><tr><td>Goldak rear length parameter a<sub>r</sub></td><td>0.005 m</td></tr><tr><td>Goldak width parameter b</td><td>0.005 m</td></tr><tr><td>Goldak height parameter c</td><td>0.002 m</td></tr><tr><td>Average thermal power P</td><td>5250 W</td></tr><tr><td>Number of elements in r-mesh</td><td>600</td></tr><tr><td>Number of elements in η-mesh</td><td>20</td></tr><tr><td>Minimal Fourier number factor f<sub>min</sub></td><td>0.5</td></tr><tr><td>Maximal Fourier number factor f<sub>max</sub></td><td>1.5</td></tr><tr><td>Number of elements in f-mesh</td><td>20</td></tr><tr><td>Minimal Biot number factor b<sub>min</sub></td><td>0.5</td></tr><tr><td>Maximal Biot number factor b<sub>max</sub></td><td>1.5</td></tr><tr><td>Number of elements in b-mesh</td><td>20</td></tr></table>

on the left, for the first PGD model with the divided domain, and in Figure 16 on the right, for the second PGD model with the 3D domain, as a contour plot.

In Figure 16 on the left, it can be observed that the error in the zone with the highest temperature is the lowest and rises in y-direction in the heat-affected zone. Furthermore, an error peak is located at the beginning of the plate in x-direction. Both of these areas with a high error are heated up only by heat conduction, which seems to need many modes for an accurate approximation. In comparison to this, the error depicted in Figure 16 on the right is significantly lower. Only at the beginning of the plate in x-direction a slightly different error occurs. This error discrepancy is expected, since the splitting of the 3D spatial domain into single 1D domains is another approximation to the real solution, which mainly increases computation time and the required number of modes in the PGD solution. Whether a division of the spatial domain should be used or not depends heavily on the problem, as the error discrepancy may not be as pronounced in general. In this case, the computation of both PGD models, each of which is an equivalent of 8000 3D finite element solutions with different parameters, requires approximately 5 h, whereas the full order model requires approximately 1 h
<figure><img src="image_16.png" /><figcaption>FIGURE 15 Contour plot of the 3D PGD solution at $t = 8$ s with $\eta = 0.8$, $\mathfrak{f} = 1$, $\mathfrak{b} = 1$ and 100 modes, where the spatial domain is divided in three 1D intervals.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 16 Contour plot of the relative error $\epsilon$ at $t = 8$ s between the PGD model with a divided spatial domain with 100 modes and a reference full order model (left) and the corresponding relative error between the PGD model with a 3D spatial domain with 40 modes and a reference full order model (right) with $\eta = 0.8$, $\mathfrak{f} = 1$ and $\mathfrak{b} = 1$ for up to 20% error.</figcaption></figure>

for a single parameter set. However, any change in the Fourier number, for example by changing the material parameters, in the Biot number, for example by changing the heat transfer coefficient for the convection boundary condition, or in the energy input leads to a completely new simulation with the full order model. For the PGD model, changing such a value only requires milliseconds to compute the new solution. Therefore, after computing the full order FE model 5 times, the PGD model is computationally more efficient. This speed up is a great benefit in parameter estimation problems, since in these problems, the model is evaluated a lot of times for different parameter settings.

5 | CONCLUSIONS

In the present article, various challenges in developing a reduced order model using the proper generalized decomposition for a temperature field with steep gradients described by a transient thermal problem are discussed. The investigations and studies illustrated with examples of increasing different complexity demonstrate the significant influence of various
aspects—from the formulation of the problem over the discretization to the choice of the considered parameters—on the convergence and accuracy, as well as the computational speed-up of the derived PGD model. The following conclusions can be drawn from this work for the analyzed thermal problems:

* A more accurate PGD solution is reached by using the finite difference method in the temporal discretization. This effect is influenced by the steepness of the gradients in time. The benefits are particularly clear for steep gradients.

* A dimensionless formulation of the thermal problem saves computational time, since fewer fixed-point iterations are required for the same accuracy and the number of variables is reduced.

* When dealing with a moving heat source, the presented mapping approach leads to an accurate and efficient PGD approximation of the temperature field. Additionally, it provides an integrated adaptive mesh refinement for the heat source location.

* The PGD model is efficiently extended by additional variables, for example material and process parameters. Therefore, the great benefits of the PGD model, like real-time computation speed, are capable of forming the basis of parameter estimation and model calibration procedures as well as in future digital twins. Thus, these simulations lead to great savings in computation time and energy cost for the industry.

* Splitting the 3D spatial domain into several one-dimensional domains is always a compromise between a loss in accuracy and an immense gain in computational speed, which should be discussed on a case-by-case basis. In the selected problem, not splitting the spatial domains leads to a PGD model with smaller approximation error, even when using a lower number of PGD modes, and a similar overall computational cost, as compared to splitting the domain.

Overall, a good agreement to reference solutions from full order models can be reached by the derived PGD models with enormous reductions in computational time. Nevertheless, this discussion underlines the importance to reconsider the multiple influences of numerical as well as problem dependent factors on a problem specific level. For applications in welding or WAAM, the methodology needs to be further extended to allow non-linear model parameters, a layer build-up, and flexible movement of the heat source along complex manufacturing paths. In conclusion, a PGD model is a reduced order model involving approximations, but depending on the use case, a PGD model can be very powerful regarding its superior computation speed. Thus, when setting up such a model, especially for the thermal problem, this article presents useful tips and tricks to improve the implementation, the convergence, and the overall performance of the model.

## AUTHOR CONTRIBUTIONS

**Dominic Strobl**: Investigation; methodology; software; writing – original draft; writing – review and editing. **Jörg F. Unger**: Conceptualization; project administration; funding acquisition; writing – review and editing. **Chady Ghnatios**: Methodology; supervision; writing – review and editing. **Annika Robens-Radermacher**: Methodology; conceptualization; project administration; software; writing – review and editing.

## CONFLICT OF INTEREST STATEMENT

The authors declare no potential conflict of interests.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study and this study's code is openly available in Zenodo at http://doi.org/10.5281/zenodo.10102489. Python-version 3.8 is available at http://www.python.org, the finite element tool FEniCS³⁷-version 2019.1 is available at https://fenicsproject.org and the PGD module PGDrome³⁶ is available at https://github.com/BAMresearch/PGDrome.

## ACKNOWLEDGMENT

Open Access funding enabled and organized by Projekt DEAL.

## ORCID

Dominic Strobl https://orcid.org/0000-0001-8873-7792

Jörg F. Unger https://orcid.org/0000-0003-0035-0951
REFERENCES

1. Aromaa S, Viitaniemi J, Leino SP. *Virtual Prototyping in Human-Machine Interaction Design*. VTT Technical Research Centre of Finland; 2014:29.

2. Hu Z, Liu Z. Heat Conduction Simulation of 2D Moving Heat Source Problems Using a Moving Mesh Method. *Adv Theor Math Phys.* 2020;2020:1-16.

3. Moreira CA, Caicedo MA, Cervera M, Chiumenti M, Baiges J. An accurate, adaptive and scalable parallel finite element framework for the part-scale thermo-mechanical analysis in metal additive manufacturing processes. *Comput Mech.* 2023;29. doi:10.1007/s00466-023-02397-6

4. Álvarez-Hostos JC, Ullah Z, Storti BA, Tourn BA, Zambrano-Carrillo JA. An overset improved element-free Galerkin-finite element method for the solution of transient heat conduction problems with concentrated moving heat sources. *Comput Methods Appl Mech Eng.* 2024;418:116574. doi:10.1016/j.cma.2023.116574

5. Álvarez-Hostos JC, Ullah Z, Storti B, Tourn BA, Fachinotti VD. Solving heat conduction problems with a moving heat source in arc welding processes via an overlapping nodes scheme based on the improved element-free Galerkin method. *Int J Heat Mass Transf.* 2022;192:122940. doi:10.1016/j.ijheatmasstransfer.2022.122940

6. Pittner A, Karkhin V, Rethmeier M. Reconstruction of 3D transient temperature field for fusion welding processes on basis of discrete experimental data. *Weld World*. 2015;59:497-512. doi:10.1007/s40194-015-0225-4

7. Yang Y, Knol MF, van Keulen F, Ayas C. A semi-analytical thermal modelling approach for selective laser melting. *Addit Manuf.* 2018;21:284-297. doi:10.1016/j.addma.2018.03.002

8. Proell SD, Munch P, Wall WA, Meier C. A highly efficient computational framework for fast scan-resolved simulations of metal additive manufacturing processes on the scale of real parts. arXiv preprint, arXiv:2302.05164 2013. doi:10.48550/arXiv:2302.05164

9. Du Y, Mukherjee T, DebRoy T. Physics-informed machine learning and mechanistic modeling of additive manufacturing to reduce defects. *Appl Mater Today*. 2021;24:101123. doi:10.1016/j.apmt.2021.101123

10. Zhu Q, Liu Z, Yan J. Machine learning for metal additive manufacturing: predicting temperature and melt pool fluid dynamics using physics-informed neural networks. *Comput Mech.* 2021;67:619-635. doi:10.1007/s00466-020-01952-9

11. Ghnatios C, Cueto E, Falco A, Duval JL, Chinesta F. Spurious-free interpolations for non-intrusive PGD-based parametric solutions: application to composites forming processes. *Int J Mater Form.* 2021;14(1):83-95. doi:10.1007/s12289-020-01561-0

12. Sudnik W, Radaj D, Erofeev W. Computerized simulation of laser beam welding, modelling and verification. *J Phys D Appl Phys.* 1996;29(11):2811-2817. doi:10.1088/0022-3727/29/11/013

13. Chinesta F, Cueto E, Abisset-Chavanne E, Duval JL, Khaldi FE. Virtual, digital and hybrid twins: a new paradigm in data-based engineering and engineered data. *Arch Comput Method E.* 2020;27(1):105-134. doi:10.1007/s11831-018-9301-4

14. Cosimo A, Cardona A, Idelsohn S. Improving the k-compressibility of Hyper Reduced Order Models with moving sources: applications to welding and phase change problems. *Comput Methods Appl Mech Eng.* 2014;274:237-263. doi:10.1016/j.cma.2014.02.011

15. Cosimo A, Cardona A, Idelsohn S. Global-Local ROM for the solution of parabolic problems with highly concentrated moving sources. *Comput Methods Appl Mech Eng.* 2017;326:739-756. doi:10.1016/j.cma.2017.08.036

16. Chinesta F, Leygue A, Bordeu F, et al. PGD-based computational vademecum for efficient design, optimization and control. *Arch Comput Method E.* 2013;20:31-59. doi:10.1007/s11831-013-9080-x

17. Chinesta F, Ammar A, Cueto E. Recent Advances and New Challenges in the Use of the Proper Generalized Decomposition for Solving Multidimensional Models. *Arch Comput Method E.* 2010;17(4):327-350. doi:10.1007/s11831-010-9049-y

18. González D, Masson F, Poulhaon F, Leygue A, Cueto E, Chinesta F. Proper Generalized Decomposition based dynamic data driven inverse identification. *Math Comput Simul.* 2012;82(9):1677-1695. doi:10.1016/j.matcom.2012.04.001

19. Ghnatios C, Barasinski A. A nonparametric probabilistic method to enhance PGD solutions with data-driven approach, application to the automated tape placement process. *Adv Model Simul Eng Sci.* 2021;8:20. doi:10.1186/s40323-021-00205-5

20. Coelho Lima I, Robens-Radermacher A, Titscher T, Kadoke D, Koutsourelakis PS, Unger JF. Bayesian inference for random field parameters with a goal-oriented quality control of the PGD forward model's accuracy. *Comput Mech.* 2022;70(6):1189-1210. doi:10.1007/s00466-022-02214-6

21. Rubio PB, Louf F, Chamoin L. Fast model updating coupling Bayesian inference and PGD model reduction. *Comput Mech.* 2018; 62:1485-1509. doi:10.1007/s00466-018-1575-8

22. Niroomandi S, González D, Alfaro I, et al. Real-time simulation of biological soft tissues: a PGD approach. *Int J Numer Meth Bio.* 2013;29(5):586-600. doi:10.1002/cnm.2544

23. Mena A, Bel D, Alfaro I, González D, Cueto E, Chinesta F. Towards a pancreatic surgery simulator based on model order reduction. *Adv Model Simul Eng Sci.* 2015;2:31. doi:10.1186/s40323-015-0049-1

24. Claus S, Kerfriden P. A stable and optimally convergent LaTIn-CutFEM algorithm for multiple unilateral contact problems. *Int J Numer Methods Eng.* 2018;113(6):938-966. doi:10.1002/nme.5694

25. Giacoma A, Dureisseix D, Gravouil A. An efficient quasi-optimal space-time PGD application to frictional contact mechanics. *Adv Model Simul Eng Sci.* 2016;3:12. doi:10.1186/s40323-016-0067-7

26. Benner P, Grivet-Talocia S, Quarteroni A, Rozza G, Schilders W, Silveira LM. *Model Order Reduction Volume 2: Snapshot-Based Methods and Algorithms*. Berlin, Boston; 2021.

27. Ghnatios C, Abisset E, Ammar A, Cueto E, Duval JL, Chinesta F. Advanced separated spatial representations for hardly separable domains. *Comput Method Appl M.* 2019;354:802-819. doi:10.1016/j.cma.2019.05.047
28. Huerta A, Nadal E, Chinesta F. Proper generalized decomposition solutions within a domain decomposition strategy. *Int J Numer Methods Eng.* 2018;113(13):1972-1994. doi:10.1002/nme.5729

29. Favoretto B, Hillerin dCA, Bettinotti O, Oancea V, Barbarulo A. Reduced order modeling via PGD for highly transient thermal evolutions in additive manufacturing. *Comput Method Appl M.* 2019;349:405-430. doi:10.1016/j.cma.2019.02.033

30. Strobl D, Unger JF, Ghnatios C, et al. Efficient bead-on-plate weld model for parameter estimation towards effective wire arc additive manufacturing simulation. *Weld World*. 2024. doi:10.1007/s40194-024-01700-0

31. Ammar A, Ghnatios C, Delplace F, et al. On the effective conductivity and the apparent viscosity of a thin rough polymer interface using PGD-based separated representations. *Int J Numer Methods Eng.* 2020;121:5256-5274. doi:10.1002/nme.6448

32. Goldak J, Chakravarti A, Bibby M. A new finite element model for welding heat sources. *Metall Trans B.* 1984;15(2):299-305. doi:10.1007/BF02667333

33. Chinesta F, Ammar A, Leygue A, Keunings R. An overview of the proper generalized decomposition with applications in computational rheology. *J Non-Newtonian Fluid Mech.* 2011;166(11):578-592. doi:10.1016/j.jnnfm.2010.12.012

34. Chinesta F, Ladevze P. *Separated Representations and PGD-Based Model Reduction: Fundamentals and Applications*. 1st ed. Springer; 2014.

35. Ghnatios C, Masson F, Huerta A, Leygue A, Cueto E, Chinesta F. Proper generalized decomposition based dynamic data-driven control of thermal processes. *Comput Methods Appl Mech Eng.* 2012;213-216:29-41. doi:10.1016/j.cma.2011.11.018

36. Robens-Radermacher A, Strobl D. PGDrome. *Zenodo*. 2023. doi:10.5281/zenodo.10075334

37. Alnaes MS, Blechta J, Hake J, et al. The FEniCS project version 1.5. *Arch Numer Software*. 2015;3:9-23. doi:10.11588/ans.2015.100.20553

**How to cite this article:** Strobl D, Unger JF, Ghnatios C, Robens-Radermacher A. PGD in thermal transient problems with a moving heat source: A sensitivity study on factors affecting accuracy and efficiency. *Engineering Reports*. 2024;6(11):e12887. doi: 10.1002/eng2.12887