Issue 
Mechanics & Industry
Volume 21, Number 6, 2020



Article Number  606  
Number of page(s)  13  
DOI  https://doi.org/10.1051/meca/2020082  
Published online  26 October 2020 
Regular Article
A mixed PGDa priori time basis strategy for the simulation of cyclic transient thermal behavior
^{1}
Structure Dynamics Materials Group, Lebanese International University LIU,
Bekaa, Lebanon
^{2}
School of Engineering, International University of Beirut BIU,
Beirut,
Lebanon
^{3}
Institut Pprime, Département Physique et mécanique des matériaux, UPR CNRS 3346, ISAEENSMA,
86360
ChasseneuilduPoitou,
France
^{*} email: jeanclaude.grandidier@ensma.fr
Received:
15
February
2020
Accepted:
4
October
2020
The knowledge of the service life of polymers under cyclic loading, widely used in industrial applications, is required and usually based on the use of methods necessitating an accurate prediction of the stabilized cycle. This implies a large computation time using the Finite Element Method (FEM) since it requires a large number of cycles for polymers. To alleviate this difficulty, a model order reduction method can be used. In this paper, a mixed strategy is investigated. Through the Proper Generalized Decomposition Method (PGD) framework, this strategy combines the Fast Fourier Transform (FFT) to create a priori time basis and the FEM to compute the related spatial modes. The method is applied to 3D thermal problems under cyclic loadings. The robustness of the proposed strategy is discussed for various boundary conditions, multitimes, and different cyclic loadings. A large time saving is obtained proving the interest of this alternative strategy to deal with fatigue simulations.
Key words: Proper Generalized Decomposition / cyclic / heat problem / Model Order Reduction / Fast Fourier Transform / different time scales
© A. Al Takash et al., published by EDP Sciences 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Polymers are widely used in industrial applications such as aerodynamics, internalcombustion engine, turbines, biomechanics, etc [1,2]. Over the years, different strategies have been studied and proposed for the fatigue of polymers. These products are subjected to thermomechanical loading cycles. For example, in car industry, parts under the engine hood are subjected to vibrationinduced oscillations in a fluctuating temperature environment. Parking phases under the sun generating high temperature with constant loads (weight) alternate with vibratory phases at very high temperature generated by the engine. In such a context, structural parts must fulfill their functions throughout the life of an automobile. This implies that robust numerical predictions over several years must to be made, as it is not possible to perform tests over such periods. Consequently, numerical tools are a mandatory step for engineers. They are then confronted with numerical calculations of thermomechanical problems cycled over a large number of cycles with an additional difficulty in the case of polymers or metallic parts under high temperature [3] for example, which do not present real stabilized cycles, thermocreep being always active throughout the life of the material. Consequently, different numerical methods and strategies have been developed and used to predict the behavior of the material under cyclic loading. Models involving transient fields and different time scales must be solved in large time intervals using very fine time steps. The time step used is related to the smallest characteristic time scale: time associated with the cyclic loading, relaxation time, and the time associated with thermal diffusivity. To solve time problems, there are several numerical issues:

the incremental scheme with a step by step construction of the spacetime solution (of time) by a succession of spatial resolutions. Each spatial resolution can be solved by different methods such as the Finite Element Method (FEM) or the PGD method (Model Order Reduction Method). The time step must be smaller than the time evolutions of physical phenomena.

the nonincremental PGD where the space and time functions are iteratively found. The space functions can be decomposed or not under 1D functions. The time discretization must be smaller than the time evolutions of the physical phenomena but the time functions are calculated on the whole time domain.

an alternative method that would be applied to more specific situations is to predefine the time basis by solving representative problems. For example, with regard to cyclic phenomena, it is possible to use the time basis related to two known data: physical parameters and loads. Once the time basis has been constructed, the iterative technique will only concern a succession of resolution of spatial problems with the FEM. This way is investigated in this paper.
Let us detail these different issues. Using an incremental time scheme, researchers encounter the following difficulties: large computation time, stability requirements, storage memory, nonconvergence [4]. In this case, the standard incremental method becomes inefficient. Other numerical techniques must be investigated and developed.
To begin with, the Model Order Reduction (MOR) is used as an alternative method. MOR is one of the techniques which is classically divided into two schemes a priori and a posteriori. MOR has been developed to solve multiphysics, complex problems [5–7] and showed efficient results. PGD (a MOR method) [4,8] is an a priori method. PGD has been used by many researchers to deal with multiphysics, fatigue, etc... Nguyen et al. [9] used the PGD to solve coupled transient multiphysics problems in 2D with different characteristic times. Hammoud et al. [10] used PGD to solve the viscoelastic behavior of polymer under creep and cyclic loadings. It was shown that a combination of PGD and an adapted time discretization for the time functions could be efficiently used to predict the viscoelastic behavior. When the relaxation time is smaller or of the same order than the cycle time, the stabilized cycle is reached more or less quickly. But this was not the case when the relaxation time is larger than the cycle time. Ammar et al. [11] studied a transient simulation of mechanical behavior where the characteristic time of the mechanical response is less than the time of interest by using a very small time step. Also, they used a spacetime separated representation to simulate an integrodifferential model. In the case of viscoelastic behavior (time dependent behavior), the standard incremental simulations are inefficient. Not to mention, Beringhier et al. [12] solved a thermoviscoelastic coupled problem in the onedimensional case, where two different characteristic times were used. Different shape functions for space and time were used to adjust the discretization for different characteristic times whereas a coupled strategy is used. Then the displacement and temperature are simultaneously determined.
Bergheau et al. [13] discussed different strategies that can be used to solve coupled spacetime problems where they used PGD as a spacetime integrator for elastoplastic problems under cyclic loading. It has been shown that in the case of a cyclic loading, more modes are needed to reach convergence. From the author’s knowledge [14], an alternative method, the LArge Time INcrement (LATIN) is the first method applied to solve evolution problems under cyclic loading. Furthermore, Boisse et al. [14] used this method to solve elastoplasticity problems under different loading cases. Straightaway, Cognard et al. [15] used LATIN method in solving viscoplastic problems under cyclic load. The feasibility of the LATIN method is tested within a large number of cycles in the case of viscoplasticity. More recently, Comte et al. [16] presented an alternative numerical method to solve the structural evolution problem. The method issued from two methods, LATIN and wavelet transform.
So far, Montebello et al. [17] used Proper Orthogonal Decomposition (POD) to describe the phenomenon of frettingfatigue. As an illustration, POD generates a spatial basis from a set of time dependent fields [18]. POD has been used in materials science, thermal science [19,20] and other fields like fluid mechanics. In [21], a strategy based on the PGD method has been proposed to determine the most significant modes and a large time saving has been obtained of order 30 compared to Finite Element Method (FEM). A combination of the Fourier analysis and Model Reduction in specific POD has been used by Ichihashi et al. [22] to analyze the energy released from experimental data, where Fast Fourier Transform (FFT) generates the temporal characteristics and POD provides the significant information. Equally important, an A Priori HyperReduction [23] showed the ability to solve a nonlinear behavior involving internal variables to predict the fatigue life [24]. Not to mention, Ryckelynck et al. [25] developed an a posteriori error estimator of hyper reduced prediction for elastoviscoplastic problems.
In light of these observations, a new approach to solve physical problems under cyclic load is introduced. This paper aims to propose a numerical strategy to reduce computation time by avoiding the use of an incremental scheme. The objective of this paper is to demonstrate that it is possible to perform a nonincremental cycling calculation on a transient evolution of a 3D structure subject to a periodic loading. This work is a first validation step for the new proposed numerical strategy. Nevertheless, the studied configurations have been chosen to be general enough to prove the effectiveness of the basic idea of the method. This strategy is based on the knowledge of an a priori basis like in the wellknown POD method, but the basis is the time basis, and on the construction of the spatial basis within an iterative procedure like for the PGD method. The method can be seen as a LATIN method as the solution is sought under a spacetime representation. The originality of the approach is the construction of the a priori time basis. The way to construct the time basis arises from the observation that, in the context of alternating time problems, the solution consists of the combination of a transient part (deviation) and a fluctuation (cyclic) for the response of a 3D conduction problem under cyclic loading. With this in mind, for a thermal problem, the solution illustrates the combination of time scales present in the phenomenon, when two different scales control the phenomenon: the physical time τ_{ϕ} related to the material property and the cycle time τ_{c} associated with the loading. Based on this observation, the time basis is obtained as follows. For different physical times and cycle times, a solution is calculated by FEM. The different solutions are evaluated at a given spatial point and the FFT is applied to the different resulting time functions. The transformed functions are decomposed into peaks. Only the greatest peaks are retained. Finally, applying the Inverse Fast Fourier Transform (IFFT), a time basis is obtained and then an analytical form is given. The paper is organized as follows. Section 2 is dedicated to the presentation of the method called the mixed strategy and more particularly the construction of the time bases and spatial bases. The numerical results are presented in Section 3. In Section 4, a discussion of the obtained results is done. Finally, conclusions and perspectives are drawn in Section 5.
2 Presentation of the mixed strategy
2.1 Problem description
As presented in the introduction, to deal with a very large number of cycles in multiscale problems, the numerical difficulties come from the foundations of the Finite Element Method which solves the problem with an incremental scheme. Nevertheless the time response is controlled by the physical parameters of the material (such as density, specific heat and conductivity coefficient) and the loading frequency. The proposed numerical method first assumes that the spacetime response denoted $T(\underset{\_}{x},t)$ can be written under a spacetime decomposition written as follows: $$T(\underset{\_}{x},t)={\sum}_{i=1}^{n}{R}_{i}(\underset{\_}{x}){S}_{i}(t).$$(1)
Time basis vectors (S_{i}(t) for i = 1..n) are defined over the whole time domain. They are derived through the analysis of different time response of Problem (1) obtained by FEM calculations for a number of particular values of the data. These time basis vectors can be used for any other problem to be solved (see Problems (2) and (3)). Note that the time computation related to the time basis vectors creation can be considered as offline computation time. The decomposition finally allows to transform the incremental spacetime problem into n spatial iterative problems to compute thespatial functions (${R}_{i}(\underset{\_}{x})\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}i\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\mathrm{1..}n$). Spatial iterative problems then become independent of the number of cycles, which reduces the computation time if the number of time basis vectors is less than the number of FEM incremental steps. Note that online computation refers to the computation of spatial functions. The different studied problems are now presented.
Consider a cube of side L with k the conductivity coefficient, ρ the density and C_{p} the specific heat that are supposed constant in the domain. The temperature field $T(\underset{\_}{x},t)$ where $\underset{\_}{x}=(x,y,z)$, t ∈ [0, L_{t}] is governed by equation (2) $$\rho {C}_{p}\frac{\partial T}{\partial t}=k\Delta T+Q,$$(2)
where Q is the volumetric heat source.
The initial temperature is supposed to vanish: $$T(x,y,z,t=0)=0.$$(3)
Different boundary conditions (BC) and volumetric heat sources (Q) are considered, defining three particular problems:

Problem (1): null Dirichlet BC on each face of the cube and Q cyclic,

Problem (2): cyclic Robin BC on each face of the cube and Q = 0,

Problem (3): cyclic Dirichlet BC on each face of the cube and Q cyclic with different cycle times.
Q cyclic depends only on time in a triangular form with $R=\frac{{Q}_{\text{min}}}{{Q}_{\text{max}}}=0$, Q_{max} = 10.10^{4}W.m^{−3}; different values of the cycle time (τ_{c}) are considered as shown in Figure 2. The cyclic Robin BC is given by ϕ = h(T − T_{∞}) where ϕ is the heat flux, h is the convection coefficient and T_{∞} is cyclic (triangular form with $R=\frac{{T}_{\text{min}}}{{T}_{\text{max}}}=0$, T_{max} = 50°C, τ_{c} = 20s). Two different constant h values are considered depending on the faces, respectively. We will consider a value for the top and bottom of the cube and an other value for all the other faces.
The cyclic Dirichlet BC is similar to T_{∞} with τ_{c} = 50s.
A schematic representation of Q and T_{∞} with respect to time is illustrated in Figure 1. Note that Q and T_{∞} are supposed uniform in space.
Problem (1) is used to generate the analytical expression of the time basis and to validate the method. These analytical expressions are then considered in Problem (2) and Problem (3) leading to different temperature gradients to discuss the domain of validity of the mixed strategy. Problem (3) involves two different cycle times. For all the problems, we assume the geometry is a cube with L = 50mm. Note that the point (x, y, z) = (25, 25, 25) is the centerof the cube.
Let us detail the time basis generation (S_{i}(t) for i = 1..n) and the spatial modes computation (${R}_{i}(\underset{\_}{x})\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}i=\mathrm{1..}n$) before presenting the numerical results.
Fig. 1 A cyclic on time in a triangular form: A = Q or A = T_{∞}. 
2.2 Time basis generation
To generate the time basis, Problem (1) is first considered with τ_{c} = 10 s, L_{t} = 200 s (leading to 20 cycles) and supposing the thermal problem has two characteristic times: the cycle time of the volumetric heat source Q (denoted by τ_{c}) and the time associated with the physical properties $\left(\frac{\rho {C}_{p}}{k}{L}_{c}^{2}\right)$ (denoted by τ_{ϕ}) where L_{c} represents a characteristic length. The evolution of the temperature is controlled by the ratio ${R}_{\tau}=\frac{{\tau}_{\varphi}}{{\tau}_{c}}$ between the physical time and the cycle time, the ratio $\frac{{L}_{c}^{2}}{k}$ being assumed given.
Remark. L_{c} could be equal to L or another characteristic length in a simulation with spatial gradients due to gradients in material properties for example. L_{c} will be here equal to one half of L due to cube symmetries.
Figure 2 shows the different behaviors of the temperature associated with different values of R_{τ} :

When the physical time is smaller than the cycle time (R_{τ} << 1), the evolution of the temperature is very fast, and the stabilized cycle is quickly reached.

By way of contrast, when the physical time and cycle time are of the same order (R_{τ} = 10 for the illustrative example), the temperature reaches the stabilization after 4 cycles.

For larger physical times (R_{τ} >> 1), the stabilization is not reached at the end of the simulation.
Different physical times are used, illustrating different responses. It is important to analyze the effect of physical time on the frequency domain.
Preliminary, Problem (1) is solved with the FEM using ABAQUS software (with 960 degrees of freedom and2, 000 time steps). In the first place, a fixed spatial point has been chosen (the center point in this case). Nevertheless, the effect of the chosen spatial point has to be discussed. Notably, an efficient basis can be generated by considering different points.
A FFT is held on the time spectrum for different problems with different physical and cycle times. Let us note that to evaluate theeffect of the physical time, the specific heat was changed whereas the conductivity and the density were kept constant. Table 1 gives the value of the physical time for different specific heat values with ρ = 950 kg.m^{−3}, k = 0.45 W.m^{−1}.C^{−1}, and L_{c} =25 mm (L_{c} being equal to one half the size of the cube). Small specific heat leads to small physical time, the physical time and the specific heat being proportional.
To highlight the effect of the spatial points, a study with a large physical time τ_{ϕ} = 1, 000 s is done. The evolution of the frequency spectrum is checked for 4 different points chosen from the cube.
Figure 3 depicts the amplitude of the FFT for four different points. First, the observation shows that the spectrum is a combination of transient and Gaussian peaks, the most compelling evidence is that the transient part is related to the physical characteristic time and the Gaussian peaks refer to the effect of cyclic load. To emphasize, since in this case, the applied load is the heat source, thus the maximum energy will be concentrated at the center of the cube which explains the highest amplitude at the center ((x, y, z) = (25, 25, 25)). Let us note that this variation shows that a single point is capable of capturing the information for different points since they are varied by the magnitude only. Accordingly, to build the basis from the first peak which is associated with the effect of the physical time, the first peak is converted for all different characteristic times, considered in these simulations, on the time domain using the IFFT. An analytical expression of the related time basis is generated within a fitting step. To clarify the time basis creation, an example, where R_{τ} = 10, τ_{c} = 10 s is detailed. First, the solution to this particular problem is solved with FEM (or another method like PGD). The FFT of the solution is shown in Figure 4a. Different peaks are seen where each peak is related to a time basis vector.
For peak 1 related to the average of the evolution and the effect of the transient (deviation), a time basis vector is generated from the IFFT of the peak and fitted within the function as shown in Figure 4b. The same procedure has been developed for different R_{τ} so that we would be able to construct a priori time basis vectors. The results obtained are compatible with the fact that as the physical time increases the evolution of the temperature decreases (see Fig. 3) and the variation of the amplitude of the peaks becomes more exponential in the frequency domain as illustrated in Figure 5. Thus, the transformation from frequency (at 0 Hz) to time domain shows the effect of the physical time.
The following empirical equation is developed by curve fitting for the first time basis vector: $$\begin{array}{l}{S}_{1}^{\varphi}(t)=a(1A\times \mathrm{exp}(\frac{5}{{\tau}_{\varphi}}t))\hfill \\ A=1\mathrm{exp}(0.7{R}_{\tau}),\hfill \end{array}$$(4)
where the coefficient a refers to the amplitude of the time basis which can be captured within the spatial modes, set to be 1, and ${S}_{1}^{\varphi}(t)$ corresponds to the time basis related to the physical time τ_{ϕ} and the cycle time τ_{c} through the ratio R_{τ}.
The other peaks refer to the effect of the cycle and physical times. To generate a time basis for these peaks, FFT has been held for different characteristic times, then converted using IFFT in time domain where an empirical equation has been generated by curve fitting as shown in equation (5): $${S}_{i+2}(t)=b\times cos((2i+1)w\times t+\theta ({R}_{\tau})),$$(5)
where i = 0...n, n is the number of time basis vectors, w = 2 × π × f, and $f=\frac{1}{{\tau}_{c}}$. To calculate θ(R_{τ}), an empirical equation has been developed: $$\theta ({R}_{\tau})=\{\begin{array}{ll}\pi ta{n}^{1}(1.04\times {R}_{\tau}),\hfill & \text{}\text{}0\le {R}_{\tau}\le 9.41\hfill \\ 1.007\pi ta{n}^{1}(1.329\times {R}_{\tau}),\hfill & \text{}\text{}9.41\le {R}_{\tau}\le 100,\hfill \end{array}$$
To sum up, to create the dictionary of time basis vectors, the first time basis vector is related to the effect of the physical time (at 0 Hz) and the other time basis vectors are related to the frequency of the load applied knowing that the physical time affects the evolution of the magnitude of the other peaks as illustrated in Figure 5.
Fig. 2 Evolution of the temperature at the center point of the cube. 
The physical time τ_{ϕ} with respect to the specific heat C_{p}.
Fig. 3 The variation of the amplitude of the peaks in the frequency domain for 4 different points. 
Fig. 4 Time basis. (a) The evolution of the temperature in the time and frequency domains at the center point (x, y, z) = (25, 25, 25) and a zoom on the second peak. (b) Generation of the time basis for the first peak. 
Fig. 5 Variation of the amplitude of the peaks in the frequency domain. 
2.3 Spatial modes computation
The time basis vectors are denoted by S_{i}(t), i = 1..n and generated solving Problem (1). They are known in an analytical form parameterized by τ_{c} and R_{τ}. Knowing S_{i} (t), i = 1..n, we look forward to the associated spatial modes ${R}_{i}(\underset{\_}{x}),i=\mathrm{1..}n$, $\underset{\_}{x}=(x,y,z)$ such that the spacetime solution is written in the separated form given by equation (1). The spatial modes are obtained via the following procedure which is based on the PGD algorithm [4,8] and more particularly an alternative direction iteration scheme.
Assume that the first k spatial modes are known and compute the next spatial mode k + 1 that is ${R}_{k+1}(\underset{\_}{x})$, which is the solution of the Galerkin variational form of equation (2) with the real and virtual fields written as follows: $$T(\underset{\_}{x},t)={\sum}_{i=1}^{k}{R}_{i}(\underset{\_}{x}){S}_{i}(t)+{R}_{k+1}(\underset{\_}{x}){S}_{k+1}(t),$$(6) $${T}^{*}(\underset{\_}{x},t)={R}_{k+1}^{*}(\underset{\_}{x}){S}_{k+1}(t).$$(7)
The field ${R}_{k+1}^{*}$ is assumed kinematically admissible.
After integrating by parts and taking into account the boundary conditions (null Dirichlet BC), this leads to the following equation: $$\begin{array}{lll}\hfill & \hfill & \rho {C}_{p}{\displaystyle {\int}_{{\Omega}_{\underset{\_}{x}}\times {\Omega}_{t}}{T}^{*}}\frac{\partial T}{\partial t}d\Omega dt+k{\displaystyle {\int}_{{\Omega}_{\underset{\_}{x}}\times {\Omega}_{t}}\overrightarrow{\nabla}}T\overrightarrow{\nabla}{T}^{*}d\Omega dt\hfill \\ \hfill & \hfill & ={\displaystyle {\int}_{{\Omega}_{\underset{\_}{x}}\times {\Omega}_{t}}Q}{T}^{*}d\Omega dt,\hfill \end{array}$$(8)
where ${\Omega}_{\underset{\_}{x}}$ is the spatial domain (the cube in our simulations) and Ω_{t} is the time domain(the time interval [0, L_{t}]). For the purpose of simplicity, Q is written under a spacetime separated representation as follows: $Q={q}_{x}(\underset{\_}{x}){q}_{t}(t)$ where q_{t} is cyclic.
To find ${R}_{k+1}(\underset{\_}{x})$, it is enough to solve a spatial problem. Using the FEM, it leads to the following discrete form: $$\begin{array}{ll}\hfill & \left(\rho {C}_{p}{\chi}_{t}^{k+1,k+1}[M]+k{\gamma}_{t}^{k+1,k+1}[K]\right)\left\{{R}_{k+1}\right\}\hfill \\ \hfill & ={\sum}_{i=1}^{k}\left(\rho {C}_{p}{\chi}_{t}^{i,k+1}[M]\left\{{R}_{i}\right\}+k{\gamma}_{t}^{i,k+1}[K]\left\{{R}_{i}\right\}\right)\hfill \\ \hfill & \text{\hspace{1em}}+{q}_{t}^{k+1}[M]{q}_{x},\hfill \end{array}$$(9)
with ${\gamma}_{t}^{i,k}={\displaystyle {\int}_{{\Omega}_{t}}{S}_{i}}(t){S}_{k}(t)$, ${\chi}_{t}^{i,k}={\displaystyle {\int}_{{\Omega}_{t}}({S}_{k}(}t)\frac{\partial {S}_{i}(t)}{\partial t})dt$, ${q}_{t}^{k}={\displaystyle {\int}_{{\Omega}_{t}}{S}_{k}}(t){q}_{t}(t)dt$, [K] is referred to the conductance matrix, [M] the mass matrix, R_{k+1} are the nodal values of the spatial function and q_{x} are the nodal values of the spatial function associated with the volume heat source, Q.
Let us detail the expressions of the matrices. [M] and [K] are obtained by assembling the elementary matrices [M^{e}] and [K^{e}], themselves constructed from the classical finite element interpolation functions written under a vector form $\{{N}^{e}(\underset{\_}{x})\}$. As ${R}_{k+1}^{e}(\underset{\_}{x})=\{{N}^{e}(\underset{\_}{x})\}{R}_{k+1}^{e}$, $[{M}^{e}]={\displaystyle {\int}_{{\Omega}_{\underset{\_}{x}}}\{{N}^{e}(}\underset{\_}{x}){\}}^{t}\{{N}^{e}(\underset{\_}{x})\}\text{d}\underset{\_}{\text{x}}$ and $[{K}^{e}]={\displaystyle {\int}_{{\Omega}_{\underset{\_}{x}}}{\{\nabla {N}^{e}(\underset{\_}{x})\}}^{t}}\{\nabla {N}^{e}(\underset{\_}{x})\}\text{d}\underset{\_}{\text{x}}$.
To choose the number of modes to be used, the following criteria is used with a stopping value of 10^{−4} $$E(n)\text{\hspace{0.17em}}(\%)=\frac{\Vert {T}_{n}{T}_{n1}\Vert}{\Vert {T}_{n1}\Vert}\times 100,$$(10)
where T_{k} is is the solution computed using the mixed strategy with k modes (n = k in Eq. (1)) and ‖ ‖ is the L^{2} norm i.e. ${\Vert T(\underset{\_}{x},t)\Vert}^{2}={\displaystyle {\int}_{{\Omega}_{t}}{\displaystyle {\int}_{{\Omega}_{\underset{\_}{x}}}T}}{(\underset{\_}{x},t)}^{2}\text{d}\underset{\_}{\text{x}}\text{dt}$.
Remark. The equations (8) and (9) to compute spatial modes depend on Q and BC and are different for Problem (2) and Problem (3). They will be given in equations (12) and (14) for Problem (2) and (3) respectively.
Figure 6 illustrates the algorithm of the proposed approach. Note that the time basis generation is considered as an offline step. Then, the computation time of the mixed strategy takes only into account the spatial modes computation, the online step.
Remark. The FEM is used here to compute spatial modes, but the PGD method with a spatial (x × y × z) decomposition could also be used.
Let us note that to evaluate the accuracy of this approach, the time relative error at a given spatial point denoted by ${\underset{\_}{x}}_{g}$ is computed using the following equation for all cases: $$tRE(\%)=\frac{\Vert {T}_{ref}({\underset{\_}{x}}_{g},t){T}_{method}({\underset{\_}{x}}_{g},t)\Vert}{\Vert {T}_{ref}({\underset{\_}{x}}_{g},t)\Vert}\times 100,$$(11)
where T_{ref} is the FEM solution, T_{method} is the solution computed using the mixed strategy, and ${\Vert T({\underset{\_}{x}}_{g},t)\Vert}^{2}={\displaystyle {\int}_{{\Omega}_{t}}T}{({\underset{\_}{x}}_{g},t)}^{2}\text{dt}$.
Remark. At the end of the iterative process, the following question arises: Are they enough precalculated time basis vectors? The number of significant time functions given by the Fourier transform is a finite number n. To verify that this number n is sufficiently large, one solution will be to use a new time function obtained by another method (nonincremental PGD method for example) and then to see if the result is modified. In this paper, as the n time functions were sufficient in the studied simulations, this technique has not been discussed.
Remark. It is important to note that this article presents a preliminary work on a new numerical strategy for dealing withfatigue simulations. To compute spatial modes, we use the 3D finite element method implemented on Matlab^{®} software which induces some restrictions on the number of degrees of freedom that can be processed (memory storage). This is the reason why it is not possible to process an industrial case. The choice of the Matlab^{®} software is due to our previous research on the PGD method. Moreover, this platform is userfriendly to manage the implementation of academic problems. To move to industrialization it will be necessary to build a spatial solver to solve (9) with specific research to be carried out on robustness, speed and parallelism, which is no longer a mechanical problem.
Fig. 6 Presentation of the mixed strategy algorithm. 
3 Numerical results
3.1 Validation of the strategy
Throughout this section, the thermal model with a null initial and boundary conditions with a cyclic heat source is considered. The cyclic heat source, Q, has a triangular form with $R=\frac{{Q}_{\text{min}}}{{Q}_{\text{max}}}=0$, τ_{c} = 10 s and a maximum amplitude equal to 100, 000 W.m^{−3}. First, a priori time basis vectors issued from the analytical expressions are constructed. Then, the time basis is used a priori to solve a spatial problem henceforth, and the solution is built from a priori time basis and spatial modes. A large physical time τ_{ϕ} equal to 1, 000 s is used by considering the specific heat C_{p} equal to 757 Joule.kg^{−1}.C^{−1} (cf. Tab. 1). These time basis vectors are then normalized using the L^{2} norm. They are called time modes. The number of modes are chosen based on the stopping criteria (Eq. (10)). Figure 7 depicts the convergence of equation (11) in function of the number of modes. Thus, 2 time modes will be used, the second peak being the most significant. The mixed strategy generates an accurate result. To check the percentage error, different points were considered where the time relative error is computed using equation (11) and presented in Table 2. The results show that the maximum time relative error is 2.8%.
The same procedure was repeated for a new physical time τ_{ϕ} = 100 s, 4 time modes were required. Accurate results are also obtained between the mixed strategy and the FEM. A time relative error less than 2.6% is obtained.
To sum up, accurate results were obtained with high timesaving of order 54 (CPU time for the approach = 1.932 s and CPU time for FEM = 104.6s). It can be explained by the fact that just 2 linear equations for τ_{ϕ} = 1, 000 s are solved with our approach (number of the a priori time basis used) instead of 2,000 linear equations with FEM (fixed time discretization with Δt = 0.1 and the total time 200 s). Let us recallthat the CPU of the mixed strategy is the online computation time, the generation of the analytical expression of the time basis isnot included. The finite element CPU time, like that of the PGD, includes the complete calculation of the spacetime solution.
The approach is able to capture a large physical time problem.
Fig. 7 Convergence with respect to the number of modes for τ_{ϕ} = 1, 000 s. 
3.2 Efficiency of the time basis for various conditions
3.2.1 Different spatial gradients
In this section, Problem (2) is considered. Two different values of convection coefficient, h_{vertical} = 4.6 W.m^{2}.K^{−1} and h_{horizontal} = 6.36 W.m^{2}.K^{−1}, are used to illustrate the effect of the density variation and geometry. These values correspond to a study in [26]. Let us note that at the top of cube (z direction), h_{horizontal} is considered and at the surfaces of the cube (x, y directions) h_{vertical} is considered. The physical properties are as follows: ρ = 950kg.m^{−3}, C_{p} = 0.075 Joule.kg^{−1}.C^{−1} and k = 0.45 Wm^{−1}.C^{−1}, accordingly, the physical time is τ_{ϕ} = 0.01 s. The time basis vectors are updated with the values of the physical and the cycle times through the analytical expressions.
Algorithm to compute spatial modes for Problem (2):
Taking into account Robin conditions, the variational formulation reads as: $$\begin{array}{lll}\hfill & \hfill & {\displaystyle {\int}_{{\Omega}_{t}}{\displaystyle {\int}_{{\Omega}_{\underset{\_}{x}}}{T}^{*}}}\rho {C}_{p}\frac{\partial T}{\partial t}\text{d}\underset{\_}{\text{x}}\text{dt}+{\displaystyle {\int}_{{\Omega}_{t}}{\displaystyle {\int}_{{\Omega}_{\underset{\_}{x}}}k}}\nabla T\nabla {T}^{*}\text{d}\underset{\_}{\text{x}}\text{dt}\hfill \\ \hfill & \hfill & +{\displaystyle {\int}_{{\Omega}_{t}}{\displaystyle {\int}_{\partial {\Omega}_{b}}h}}T{T}^{*}\text{dSdt}={\displaystyle {\int}_{{\Omega}_{t}}{\displaystyle {\int}_{\partial {\Omega}_{b}}h}}{T}_{\infty}(t){T}^{*}\text{dSdt}.\hfill \end{array}$$(12)
As for equation (8), an equation analogous to equation (9) can be derived for equation (12).
Numerical result.
The method converges with 10 modes. A good agreement, as shown in Figure 8, and a large time saving of order 50 are obtained with a maximum of time relative error less than 2% in the whole domain as shown in Table 3.
Fig. 8 Solution from the mixed strategy compared to the Finite Element solution for Robin boundary conditions at (x, y, z)= (25, 25, 25). 
3.2.2 Multitimes
In this case,the possibility of introducing two different cycle times is discussed, that is Problem (3) is considered with τ_{c} = 20s for Q and τ_{c} = 50s for the cyclic temperature on the boundary (denoted by T_{b}(t)).
Straightaway, a priori time basis is constructed using the analytical expressions derived in Section 2 from Problem (1)by updating the physical time and cycle time. In this case, the bases will be related to two different cycle times: 20s and 50s. The position of these peaks can be generated using FFT of the load applied, and later the IFFT by modifying the phase angle to fit with the physical time. As an illustration, Figure 9 pinpoints the location of the peaks. ${S}_{1}^{b}$, ${S}_{2}^{b}$... ${S}_{n}^{b}$ and ${S}_{1}^{q}$, ${S}_{2}^{q}$...${S}_{n}^{q}$ denote the time bases issued from the Dirichlet boundary conditions and the heat source respectively, where n refers to the number of basis vectors. To alleviate the notations, n is used for the two different kinds of time bases but n can be different for ${S}_{i}^{b}$ and ${S}_{i}^{q}$.
The first investigation of these bases show that the major information are contained in the first and second peaks. The results show that for loads with different cycle times and amplitudes, different time bases are generated by changing the cycle time of the analytical expressions. However, these bases are essential to construct the solution.
Algorithm to compute spatial modes for Problem (3).
Let us specify, the nonhomogeneous Dirichlet conditions are applied along the boundary ∂Ω: $T(\underset{\_}{x},t)={T}_{b}(t)$. In order to apply the separated form, we consider a function ${T}^{D}(\underset{\_}{x},t)$ that satisfies the Dirichlet conditions as presented previously in Section 3.2.2. The separated representation of the boundary condition is under the form: $${T}^{D}(\underset{\_}{x},t)=G(\underset{\_}{x}){T}_{b}(t),$$(13)
where $G(\underset{\_}{x})=1$ at $\underset{\_}{x}$ ∈ ∂Ω and $G(\underset{\_}{x})=0$ at $\underset{\_}{x}$ ∈ Ω\∂Ω.
Besides, the heat source is written under a separated form $Q={q}_{x}(\underset{\_}{x}){q}_{t}(t)$ where q_{t} is cyclic.
To find the spatial modes, the discrete form of equation (2) (Eq. (9)) is rewritten as follows: $$\begin{array}{ll}\hfill & \left(\rho {C}_{p}{\chi}_{t}^{k+1,k+1}[M]+k{\gamma}_{t}^{k+1,k+1}[K]\right)\left\{{R}_{k+1}\right\}\hfill \\ \hfill & ={\sum}_{i=1}^{k}\left(\rho {C}_{p}{\chi}_{t}^{i,k+1}[M]\left\{{R}_{i}\right\}+k{\gamma}_{t}^{i,k+1}[K]\left\{{R}_{i}\right\}\right)\hfill \\ \hfill & \text{\hspace{1em}\hspace{1em}}+{q}_{t}^{k+1}[M]{q}_{x}\left(\rho {C}_{p}{\chi}_{t,b}^{k+1}[M]k{\gamma}_{t,b}^{k+1}[K]\right)\left\{G\right\},\hfill \end{array}$$(14)
with ${\gamma}_{t,b}^{k}={\displaystyle {\int}_{{\Omega}_{t}}{S}_{k}}(t){T}_{b}(t)$, ${\chi}_{t,b}^{k}={\displaystyle {\int}_{{\Omega}_{t}}({S}_{k}(}t)\frac{\partial {T}_{b}(t)}{\partial t})dt$, G are the nodal values of the spatial function associated with ${T}^{D}(\underset{\_}{x},t)$. Two last terms are added to equation (9) to take into account the nonhomogeneous Dirichlet condition.
Numerical result.
The number of modes needed is 22. The time bases are alternatively based on the location of the peaks that can be determined from the decomposition of the load. The results are shown in Figure 10. There are accurate results compared tothe FEM that leads us to conclude that the mixed strategy method is able to solve problems under different cycle times. Table 4 depicts the relative error for different spatial points, where the maximum relative error is equal to 1.8%.
Remark. Since in this case, we have different cycle times due to the loads applied (through heat source and boundary conditions), the bases are related to both loads applied leading to two different values of R_{τ}.
Effect of the order.
Considering the model encountering two different loads with different cycle times (τ_{c} = 50 s, τ_{c} = 20 s) and τ_{ϕ} = 10 s, the bases are established a priori using the analytical expression. To check the effect of the order of the time basis, 3 different cases are discussed:

case 1: the time bases vectors are chosen alternatively according to their related frequency.

case 2: the time basis vectors issued from τ_{c} = 50 s are the first followed by the time basis vectors issued from τ_{c} = 20 s,

case 3: the time basis vectors issued from τ_{c} = 20 s are the first followed by the time basis vectors issued from τ_{c} = 50 s.
The results are shown in Figure 11. The evolution of the ratio shows that within case 2 and case 3, the number of modes is 9, and 10, respectively. However, for case 1, 20 modes are needed. But regarding the relative error compared to the FEM, it is noticed that for case 2 and case 3 the relative error is higher than 18%. However, for case 1 the relative error is 1.8%, thus the order of the basis is important: the time bases have to be ordered with increasing frequency. For this case, the time saving is 40 (where CPU using FEM is 80s, and using the approach is 2s) and an accurate prediction is obtained as illustrated in Figure 10.
Fig. 9 Bases issued from different cycle times imposed by a volumetric heat source (right) and a Dirichlet boundary condition (left). 
Fig. 10 The evolution of temperature using the approach (using the order of bases from case 1) compared to FEM solution at (x, y, z) = (25, 25, 25). 
Fig. 11 Convergence with respect to the number of modes for multitimes case: effect of the time basis vectors order. 
4 Discussion
4.1 Larger time domain
Throughout this section, the model encounters a load with a larger time domain with the following properties: τ_{ϕ} = 100 s, τ_{c} = 20 s leading to R_{τ} = 5 and the total time of the time domain is L_{t} = 2, 000 s (10 times larger than the fitting domain). The bases are used directly from the analytical expressions generated in the time basis generation section by taking into account the effect of cycle time and physical time. Let us note that we here assume that the analytical expressions are valid for large time domains, larger than the fitting time domain.
For example, the analytical expression of the first time basis vector leads to: $${S}_{1}^{100}=1\mathrm{exp}(\frac{5}{100}t),$$
and the second time basis vector: $${S}_{2}=cos(\frac{2\pi}{20}\times t+1.76).$$
The results obtained with our approach are compared to the FEM solution. The first set of analyses highlighted the impact of the physical time, where the stabilized cycle is obtained after 5 cycles. Interestingly, only 4 modes are needed to obtain an accurate solution with an error less than 3.6%. The error iscalculated in different spatial points (near boundaries and at the center) using L^{2} norm and are reported for some spatial points as shown in Table 5.
This is an important issue that a priori bases generated from a short time L_{t} = 200 s are used for larger time domain L_{t} = 2, 000 s. Let us note that this is valid only if the model used is valid for a larger domain. Thus, this method represents a viable technique to solve a physical problem with a large time domain. A large time saving of order 80 has been obtained in this case.
4.2 Different cyclic loadings
This sectionaims to check the effect of introducing a load with many frequencies. To do this, an Heaviside load is studied. This case calls into the following questions:

Can this approach take into account the effect of this type of load (many frequencies)?

What is the time basis to be used?
The model now encounters an Heaviside load with a physical time equal to 100 s, as shown in Figure 12. This load has different plateaus at T =0°C where the FFT of this load leads to generate only one fundamental frequency related to the cycle time τ_{c} = 20 s and all the location of the other peaks is related to this frequency. The time basis vectors need to be updated by adding $\frac{\pi}{2}$ to the analytical expression previously defined in Section 2 from Problem (1) and to modify the phase angle by doing the direct fitting with FFT. Thus, we can use directly the time basis vectors. As an example, the first two time basis vectors used are: $${S}_{1}=11\times exp(\frac{5}{100}t),$$ $${S}_{2}=1\times cos(\frac{2\pi}{20}t+(1.4+\frac{\pi}{2})).$$
The evolution of temperature using the approach is compared to a reference solution leading to a relative error less than 5%.
Fig. 12 Heaviside load (a), comparison between the FEM solution and the approach solution using an Heaviside load (b). 
5 Conclusions
Solving physical problems under cyclic loadings may lead to a large computation time when incremental methods are used. The present study focuses on finding an innovative way to enhance the computation time. The evidence from this study intimates that due to the presence of different time scales in different problems, the solution will be affected by the relationship between these characteristic times (physical time and cycle time) which leads to thinking about the correlation and the effect of each time scale. Therefore, this study is the first step towards enhancing and developing a priori time basis. A new approach based on building a priori time basis was presented. It leads to the generation of a dictionary which will be used for a certain family of problems. To sum up, our study provides that an a priori time basis can be generated by only using FFT analysis and the knowledge of time scales. This approach has been investigated with different problems such as homogeneous boundary conditions, nonhomogeneous boundary conditions, Robin boundary conditions, and different types of applied loads. To illustrate, we succeeded in finding a priori time basis when two different loads with different cycle times and amplitudes are considered. Accurate results compared to FEM have been obtained with a large time saving of around 50 times. Note that with classical PGD, we have obtained a time saving of order 30 compared to FEM. This means that the mixed strategy was able to reduce effectively the computation time compared to PGD and FEM respectively. The results are then encouraging. The method can be useful, for example, in optimization problems where the best solution for a criterion must be found. We consider a given set of parameters and the corresponding solution to the problem. The proposed method then makes it possible to calculate very quickly the solutions of the problem for parameter sets close to the given parameters.
Future studies will be firstly performed on nonlinear problems with a nonlinearity due to the use of nonlinear material parameters and can be enhanced to other non linearities such as coupling terms between diffusothermo equations and finally by considering thermoviscoelasticity.
References
 T.L. Anderson, NASA STI/Recon Technical Report A 92, 809 (1991) [Google Scholar]
 T.L. Anderson, Fracture mechanics: fundamentals and applications, CRC Press, 2017 [CrossRef] [Google Scholar]
 R. Ahmed, P.R. Barrett, M. Menon, T. Hassan, Int. J. Solids Struct. 126–127, 90 (2017) [Google Scholar]
 F. Chinesta, R. Keunings, A. Leygue, The proper generalized decomposition for advanced numerical simulations: a primer, Springer International Publishing, 2014 [CrossRef] [Google Scholar]
 F. Chinesta, P. Ladeveze, E. Cueto, Arch. Computat. Methods Eng. 18, 395 (2011) [CrossRef] [Google Scholar]
 P. Benner, S. Gugercin, K. Willcox, SIAM Rev. 57, 483 (2015) [CrossRef] [Google Scholar]
 E. Stein, R. de Borst, T.J.R. Hughes, Encyclopedia of computational mechanics, Vol. 1, Wiley, 2004 [CrossRef] [Google Scholar]
 F. Chinesta, A. Ammar, A. Leygue, R. Keunings, J. NonNewtonian Fluid Mech. 166, 578 (2011) [CrossRef] [Google Scholar]
 T.L. Nyuyen, Ph.D. thesis, ENSMA (2012) [Google Scholar]
 M. Hammoud, M. Beringhier, J.C. Grandidier, Comptes Rendus Mécanique 342, 671 (2014) [CrossRef] [Google Scholar]
 A. Ammar, A. Zghal, F. Morel, F. Chinesta, Comptes Rendus Mécanique 343, 247 (2015) [CrossRef] [Google Scholar]
 M. Beringhier, M. Gueguen, J.C. Grandidier, Arch. Comput. Methods Eng. 17, 393 (2010) [Google Scholar]
 J.M. Bergheau, S. Zuchiatti, J.C. Roux, É. Feulvarch, S. Tissot, G. Perrin, Comptes Rendus Mécanique 344, 759 (2016) [CrossRef] [Google Scholar]
 P. Boisse, P. Bussy, P. Ladeveze, Int. J. Numer. Methods Eng. 29, 647 (1990) [Google Scholar]
 J.Y. Cognard, P. Ladevèze, Int. J. Plast. 9, 141 (1993) [Google Scholar]
 F. Comte, H. Maitournam, P. Burry, T.M.L. Nguyen, Comptes Rendus Mécanique 334, 317 (2006) [CrossRef] [Google Scholar]
 C. Montebello, Ph.D. thesis, Université ParisSaclay (2015) [Google Scholar]
 A. Chatterjee, Current Sci. 78, 808 (2000) [Google Scholar]
 Y.C. Liang, H.P. Lee, S.P. Lim, W.Z. Lin, K.H. Lee, C.G. Wu, J. Sound Vibr. 252, 527 (2002) [CrossRef] [Google Scholar]
 M.O. Efe, H. Ozbay, Proper orthogonal decomposition for reduced order modeling: 2D heat flow, in: Proceedingsof 2003 IEEE Conference on Control Applications, CCA 2003 (2003), Vol. 2, pp. 1273–1277 [Google Scholar]
 A. Al Takash, M. Beringhier, M. Hammoud, J.C. Grandidier, Comput. Phys. 375, 950 (2018) [CrossRef] [Google Scholar]
 F. Ichihashi, S.M. Jeng, K. Cohen, Proper orthogonal decomposition and Fourier analysis on the energy release rate dynamics of a gas turbine combustor, in 48th AIAA Aerospace Science Meeting (2010), Vol. AIAA201022 [Google Scholar]
 D. Ryckelynck, F. Chinesta, E. Cueto, A. Ammar, Arch. Comput. Methods Eng. 13, 91 (2006) [Google Scholar]
 D. Ryckelynck, D. Missoum Benziane, Comput. Methods Appl. Mech. Eng. 199, 1134 (2010) [Google Scholar]
 D. Ryckelynck, L. Gallimard, S. Jules, Adv. Model. Simul. Eng. Sci. 2, 19 p. (2015) [Google Scholar]
 N.T. Nguyen, X. Huang, Biomed. Microdevices 8, 133 (2006) [CrossRef] [PubMed] [Google Scholar]
Cite this article as: A. Al Takash, M. Beringhier, M. Hammoud, J.C. Grandidier, A mixed PGDa priori time basis strategy for the simulation of cyclic transient thermal behavior, Mechanics & Industry 21, 606 (2020)
All Tables
All Figures
Fig. 1 A cyclic on time in a triangular form: A = Q or A = T_{∞}. 

In the text 
Fig. 2 Evolution of the temperature at the center point of the cube. 

In the text 
Fig. 3 The variation of the amplitude of the peaks in the frequency domain for 4 different points. 

In the text 
Fig. 4 Time basis. (a) The evolution of the temperature in the time and frequency domains at the center point (x, y, z) = (25, 25, 25) and a zoom on the second peak. (b) Generation of the time basis for the first peak. 

In the text 
Fig. 5 Variation of the amplitude of the peaks in the frequency domain. 

In the text 
Fig. 6 Presentation of the mixed strategy algorithm. 

In the text 
Fig. 7 Convergence with respect to the number of modes for τ_{ϕ} = 1, 000 s. 

In the text 
Fig. 8 Solution from the mixed strategy compared to the Finite Element solution for Robin boundary conditions at (x, y, z)= (25, 25, 25). 

In the text 
Fig. 9 Bases issued from different cycle times imposed by a volumetric heat source (right) and a Dirichlet boundary condition (left). 

In the text 
Fig. 10 The evolution of temperature using the approach (using the order of bases from case 1) compared to FEM solution at (x, y, z) = (25, 25, 25). 

In the text 
Fig. 11 Convergence with respect to the number of modes for multitimes case: effect of the time basis vectors order. 

In the text 
Fig. 12 Heaviside load (a), comparison between the FEM solution and the approach solution using an Heaviside load (b). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.