© A. Al Takash et al., published by EDP Sciences 2020

1 Introduction

Polymers are widely used in industrial applications such as aerodynamics, internal-combustion 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 vibration-induced 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 space-time 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 non-incremental 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, non-convergence [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 space-time 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 one-dimensional 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 space-time problems where they used PGD as a space-time 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 fretting-fatigue. 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 Hyper-Reduction [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 non-incremental 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 well-known 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 space-time 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 multi-scale 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 space-time response denoted $T(\underset{\_}{x},t)$ can be written under a space-time 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 space-time problem into n spatial iterative problems to compute thespatial functions ($R_i(\underset{\_}{x})\text{for}i=\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⁴W.m⁻³; 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{for}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⁻³, k = 0.45 W.m⁻¹.C⁻¹, 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(1-A\times \exp (-\frac{5}{{\tau }_{\varphi }}t))\hfill \\ A=1-\exp (-0.7R_{\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.

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 space-time 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 space-time 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 & +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⁻⁴ $$E(n)(\%)=\frac{\Vert T_n-T_{n-1}\Vert }{\Vert T_{n-1}\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² 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 pre-calculated 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 (non-incremental 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 user-friendly 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⁻³. 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⁻¹.C⁻¹ (cf. Tab. 1). These time basis vectors are then normalized using the L² 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 time-saving 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 space-time 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².K⁻¹ and h_horizontal = 6.36 W.m².K⁻¹, 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⁻³, C_p = 0.075 Joule.kg⁻¹.C⁻¹ and k = 0.45 Wm⁻¹.C⁻¹, 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 Multi-times

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 non-homogeneous 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 re-written 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 & +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 non-homogeneous 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 multi-times 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-\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² 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=1-1\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, non-homogeneous 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 diffuso-thermo equations and finally by considering thermoviscoelasticity.

References

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 multi-times 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