Routes for Efficient Computational Homogenization of Nonlinear Materials Using the Proper Generalized Decompositions

Computational homogenization is nowadays one of the most active research topics in computational mechanics. Different strategies have been proposed, the main challenge being the computing cost induced by complex microstructures exhibiting nonlinear behaviors. Two quite tricky scenarios lie in (i) the necessity of applying the homogenization procedure for many microstructures (e.g. material microstructure evolving at the macroscopic level or stochastic microstructure); the second situation concerns the homogenization of nonlinear behaviors implying the necessity of solving microscopic problems for each macroscopic state (history independent nonlinear models) or for each macroscopic history (history dependant non-H

linear models). In this paper we present some preliminary results concerning the application of Proper Generalized Decompositions-PGD-for addressing the efficien solution of homogenization problems. This numerical technique could allow to compute the homogenized properties for any microstructure or for any macroscopic loading history by solving a single but highly multidimensional model. The PGD allows circumventing the so called curse of dimensionality that mesh based representations suffer. Even if this work only describes the firs steps in a very ambitious objective, many original ideas are launched that could be at the origin of impressive progresses.

Introduction
Homogenization approaches are now intensively used in numerous engineering applications. These approaches belong to the more general category of multiscale methods, in which it is usually distinguished two main families of approaches, namely hierarchical and concurrent, see for example [15,18,29].
In hierarchical approaches, the behavior at the highermacro-scale is obtained from the solution of a boundary value problem solved at the lower-micro-scale and posed on a representative volume element (RVE) of the material. This behavior can be determined under the form of macroscopic constitutive relations whose effective properties are identifie from the solution of the microscopic problem. Micromechanical approaches fall in this category [6,44]. The microscopic solution can also be incorporated in the macroscopic problem through a numerical scheme. It can be condensed at the macroscale as in the case of the variational multiscale method [20], or within a decomposition domain framework for structural mechanics applications [31]. The microscopic solution may also enrich the macroscale solution using strategies such as X-FEM and G-FEM, see e.g. [17] and [38] respectively. One of the major feature of hierarchical approaches is that the macroscopic and microscopic problems are not coupled. This means that the microscopic problems can be solved once for all, and consequently that only a finit number of microscopic solutions have to be known.
By contrast, the microscopic and macroscopic scales are solved simultaneously in concurrent approaches. This is required when the macroscopic behavior can not be obtained explicitly, which occurs as soon as microscopic nonlinear behaviors are involved. The area where the microscopic description is needed can be isolated but in this paper we will focus on problems where the whole domain is characterized by an underlying microstructure exhibiting a nonlinear behavior. These methods were reviewed in [24] in the continuum mechanics framework. The most popular approach is the FE 2 method initiated by Feyel in [16], which has also be named later on computational homogenization, see e.g. [19]. This method assumes that the microstructure is very small compared to the macroscopic domain. So at the macroscopic scale, it plays the role of a material point-an integration point if the finit element is used-at which a RVE is attached.
As noticed in [29], both the hierarchical and concurrent simulation approaches offer cost savings over large scale direct numerical simulation of the microstructure. However efforts are directed towards improving their numerical efficien y.
In hierarchical multiscale methods, one current challenge is the analysis of a representative volume element with a highly complex geometry. Powerful techniques such as computer tomography actually enable to obtain high resolution 3D images of material microstructures. Thus an image with 1024 × 1024 × 1024 voxels is now common, but its incorporation in a numerical model requires the use of an appropriate solution method.
In concurrent multiscale methods, such as the FE 2 approach, one main limitation is its computational cost, since several microscopic problems have to be solved at each integration point of the macroscopic domain. Therefore attempts have been made to build hierarchical approaches for nonlinear materials. One can thus defin approximate nonlinear constitutive relations, see e.g. [30]. Another possibility is to build a discrete material map (see Refs. [40,41,43] for the case of nonlinear elasticity). This approach has been used in these references in nonlinear elasticity. It consists in discretizing the macroscopic strain space and determine the microscopic solution for each node of this space. Even if the number of microscopic problems to be solved is large, this method is much more efficien than the FE 2 approach.
Recently, the authors have proposed a solution method based on the separated representation of the unknown fields coined Proper Generalized Decomposition-PGD-. This technique was successfully applied for solving multidimensional models encountered in the kinetic theory description of complex fluids involving steady state and transient linear and non-linear models [2][3][4]32]. Time-space separated representations were originally introduced many years ago by P. Ladeveze for addressing complex time-dependent nonlinear models within the LATIN framework (see [26] and the references therein). Some recent works proving the potentiality of this approach can be found in [27,33] among many others. In [2,3] we generalized the separated representation for treating general multidimensional models. Obviously, a firs possibility of separated representations consists of a separated representation of the physical space, by writing a generic function as a finit sums of functional products whose involved functions only depends on a single coordinate. This strategy was applied in [11] in high resolution linear homogenization. The computational complexity of separated representations scales linearly with the dimension of the space instead the exponential growing characteristic of mesh-based descriptions. Obviously, the use of the PGD allows for impressive computing cost savings, making possible the solution of models never until now solved (e.g. models define in spaces involving hundreds dimensions [3,12]).
Models usually encountered in computational mechanics that generally involve space and time as coordinates can be transformed, as we illustrate later, into highly multidimensional models. Obviously, the resulting models suffers of the curse of dimensionality if mesh-based representations are considered. However, the Proper Generalized Decomposition allows to circumvent efficientl such illness. For the sake of clarity, the use of a multidimensional modeling is illustrated and motivated in a quite simple physics. Imagine for example that you are interested in solving the heat equation without knowing the microscopic material thermal conductivity (because of its stochastic nature or simply because it has not been measured prior to solve the thermal model). You have three possibilities: (i) wait to know the conductivity before solving the heat equation (conservative solution !); (ii) solve the equation for many values of the conductivity (a sort of Monte Carlo) (brute force approach !); or (iii) solve the heat equation only once for any value of the conductivity (cleverest alternative !). Obviously the third alternative is the most exciting one. The computation of this general solution consists in introducing the conductivity as an extra coordinate, playing the same role than the standard space and time coordinates, even if it is not derivated. This procedure works very well, and can be extended for introducing many other extra-coordinates: source terms, initial conditions . . . [13]. This is the basic idea of spectral approaches in stochastic analysis, which are receiving a growing interest in computational science [34].
In the fiel of homogenization, two simpler scenarii can be considered: (i) defin a generic microstructure from a voxel description and assume the conductivity of each voxel as an extra-coordinate; and (ii) address non-linear homogenization by solving the model at the microscopic level for any macroscopic state or history. This paper has as main objective opening some routes for addressing both scenarios thanks to the use of Proper Generalized Decompositions. For the sake of simplicity, and because this work is a firs attempt in this ambitious objective, we will focus on the homogenization of linear and non-linear (time-independent and time-dependent) thermal models. We are aware that this contribution is opening many routes without closing any of them, but as just mentioned this is not more than a preliminary analysis of the potential application of PGD for some multidimensional descriptions related to computational homogenization. The application of these ideas to more realistic scenarios as well as to thermomechanical models is a work in progress.
This paper is be organized as follows: After illustrating the use of Proper Generalized Decomposition for addressing parametric models, it will be applied to the computation of effective conductivity of linear heterogeneous materials. The benefit of the PGD will be demonstrated for the analysis of random materials, for which numerous analyses are performed on mesoscale windows in order to defin the RVE, see e.g. [22,23]. Then, nonlinear thermal models will be addressed. In this case, the PGD enables to obtain an explicit representation of the nonlinear macroscopic behavior, for both history independent and history dependent microscopic behaviors. Finally, numerical examples are given in the last section to illustrate the performance of the approach.

Illustrating the Solution of Multidimensional Parametric Models by Using the PGD
In what follows, the construction of the Proper Generalized Decomposition is illustrated by considering a quite simple parametric heat transfer equation: where (x, t, k) ∈ × I × and for the sake of simplicity the source term is assumed constant, i.e.f = cte. Because the conductivity is considered unknown, it is assumed as a new coordinate define in the interval . Thus, instead of solving the thermal model for different values of the conductivity parameter we prefer introducing it as a new coordinate. The price to be paid is the increase of the model dimensionality; however, as the complexity of PGD scales linearly with the space dimension the consideration of the conductivity as a new coordinate allows fast and cheap solutions.
The solution of (1) is searched under the form: In what follows we are assuming that the approximation at iteration n is already done: and that at present iteration we look for the next functional product X n+1 (x) · T n+1 (t) · K n+1 (k) that for alleviating the notation will be denoted by R(x) · S(t) · W (k). Prior to solve the resulting non linear model related to the calculation of these three functions, a model linearization is compulsory. The simplest choice consists in using an alternating directions fi ed point algorithm. It proceeds by assuming S(t) and W (k) given at the previous iteration of the non-linear solver and then computing R(x). From the just updated R(x) and W (k) we can update S(t), and finall from the just computed R(x) and S(t) we compute W (k). The procedure continues until reaching convergence. The converged functions R(x), S(t) and W (k) allow definin the searched functions: We are illustrating each one of the just referred steps.
Computing R(x) from S(t) and W (k) We consider the global weak form of (1): where the trial and test functions write respectively: and Introducing (5) and (6) where R n define the residual at iteration n that reads: Now, being known all the functions involving the time and the parametric coordinate, we can integrate (7) in their respective domains I × . Integrating in I × and taking into account the following notations Equation (7) reduces to: Equation (10) define an elliptic steady state boundary value problem that can be solved by using any discretization technique operating on the model weak form (finit elements, f nite volumes . . . ). Another possibility consists in coming back to the strong form of (10): that could be solved by using any collocation technique (fi nite differences, SPH . . . ).
Computing S(t) from R(x) and W (k) In the present case the test function writes: Now, the weak form reads (13) that integrating in the space × and taking into account the notation (9) results: Equation (14) represents the weak form of the ODE definin the time evolution of the fiel S that can be solved by using any stabilized discretization technique (SU, Discontinuous Galerkin, . . . ). The strong form of (14) reads: than can be solved by using backward finit differences, or higher order Runge-Kutta schemes, among many other possibilities.
Computing W (k) from R(x) and S(t) In the present case the test function writes: Now, the weak form reads that integrating in the space × I and taking into account the notation (9) results: Equation (18) does not involve any differential operator. The strong form of (18) reads: that represents an algebraic equation. Thus, the introduction of parameters as additional model coordinates has not a noticeable effect in the computational cost, because the original equation does not contain derivatives with respect to those parameters. Finally, note that other minimization strategies have been proposed, leading to more robust and faster convergence for building-up the PGD [13].

Remark 1
The construction of each term in the sum (2) needs a certain number of iterations because of the non linearity of the problem related with the approximation (3). We denote by m i the number of iterations that were needed for computing the i-sum in (2). Let m = i=N i=1 be the total number of iterations involved in the construction of the separated approximation (2). It is easy to note that the solution procedure needs the solution of m 3D problems related to the construction of the space functions X i (x), i = 1, . . . , N, m 1D ordinary differential equations related to the construction of functions T i (t) and m linear systems related to the definitio of functions K i (k). In general m rarely exceeds ten. On the other hand the number of sums N needed to approximate the solution of a given problem depends on the solution regularity itself, but all the experiments carried out until now reveal that this number ranges from few tens to few hundreds. Thus, we can conclude that the complexity of the solution procedure is of some hundreds of 3D solutions (the cost related to the one dimensional problems being negligible with respect to the one related to the 3D problems). Now, if we assume a classical approach one should solve a 3D problem at each time step and for each value of the parameter k. In usual applications the complexity reaches millions of 3D solutions. In [14] we proved that the CPU time savings by applying the PGD can be of several orders of magnitude.

Formulation of the Homogenization Problem
Due to the microscopic heterogeneity, the macroscopic thermal modeling needs an homogenized thermal conductivity which depends on the microscopic details.
To compute this homogenized thermal conductivity in the linear case, one could isolate a representative volume element at position X, rve (X) (this representative volume element is supposed to be defined using for example methods which are proposed in [22,23,35]) and assume that the microstructure is perfectly define at such scale. Thus, the microscopic conductivity k(x) is known at each point in the microscopic domain x ∈ rve (X).
We can defin the macroscopic temperature gradient at position X, G(X), from: where We also assume the existence of a localization tensor L(x, X) such that Now, we consider the microscopic heat flu q according to the Fourier's law (22) and the macroscopic counterpart Q(X) that writes: from which the homogenized thermal conductivity can be define from As k(x) is perfectly known everywhere in the representative volume element, the definitio of the homogenized thermal conductivity tensor only requires the computation of the localization tensor L(x, X). Several approaches are proposed in the literature to defin this tensor, according to the choice of the boundary conditions. Our objective here is not to discuss this choice. The interested reader can fin some details in [21,22,37], or [36]. For sake of simplicity, we use essential boundary conditions corresponding to the assumption of uniform temperature gradient on the RVE boundary. Hence we are led to consider in the 3D case the solution of the three boundary value problems related to the steady state heat transfer model in the microscopic domain, define by: It is easy to prove that these three solutions verify where (·) T denotes the transpose. Thus, the localization tensor results finally The resulting non-concurrent homogenization procedure is illustrated in Fig. 1. As soon as tensor L(x, X) is known at each position X, the constitutive law relating the macroscopic temperature gradient and macroscopic heat flu becomes defined From the computational point of view the main difficult concerns the solution of the three boundary value problems using very fin meshes required to represent accurately all the microscopic details as well as the application of this procedure at different places for accounting for the spatial evolution of the microstructure of its stochastic nature [9].

Separated Representation of the Physical Space
As mentioned in introduction, recent advances in imaging techniques offer the opportunity to obtain high resolution images of material microstructures. This leads to very large finit element models which can reach one billion of degrees of freedom, see e.g. [5]. Nowadays, the performance of computers fails to solve, in a reasonable time and with classical solvers, problems of this size. In order to circumvent this difficult we propose to use a separated representation of the physical space.
Thus, the temperature field involved in the thermal models (25)- (27) are searched in a separated form. Therefore, at iteration n the solution writes and it is enriched until reaching convergence by computing the next term in the finit sums decomposition Applying the strategy described in the previous section we compute the function R(x) assuming known S(y) and W (z), then S(y) from the just computed R(x) and W (z) and finall W (z) from R(x) and S(x). After reaching convergence we can assign We assume that the microscopic conductivity can be written in a separated form, i.e. each component of the conductivity tensor k ij (x) can be expressed from: This representation can be performed by applying for example a multidimensional singular value decomposition (SVD) [25]. Following the procedure described in Sect. 2, the computation of R(x) implies the solution of the second order where the coefficient depend on the integrals of functions depending on y and z in their respective domains. Analogously, the calculation of functions S(y) and W (z) implies the solution of similar equations but implying different coefficient and source terms: and Thus, as proved in [11], models involving millions of nodes in each direction can be solved because the procedure based on separated representations avoids the definitio of a mesh. This enables degrees of resolution never until now reached, the main issue being the solution of the one dimensional boundary value problems. In general the discrete form of these problems leads to tridiagonal matrices that can be easily solved. In the next section we are describing a procedure, originally proposed in [7], specially adapted to such kind of solutions that avoids the solution of a large size linear system.

Remark 2
The main limitation in using this fully separated representation is precisely the necessity of build-up a separated approximation of the microscopic thermal conductivity (32) because complex microstructures imply thousands of terms in the sum (32) as well as the necessity of defin ing the conductivity in a grid of the representative volume element prior to apply a singular value decomposition-SVD-. The separated representation of the conductivity at the microscopic level is not mandatory, but without such a representation, integrals in the whole domain have to be computed, which can dramatically reduce the computational efficien y. In Sect. 3.3 we explore some possibilities for alleviating this drawback.

Efficien Solvers of One-Dimensional BVP
Consider the generic second order BVP implying the fiel with u(x = 0) = u l and u(x = L) = u r . The solution of this equation can be written as the addition of the general solution of the homogeneous equation and a particular solution of the complete equation, that is denoted by h h (x) and u c (x) respectively.
The solution u c (x) can be computed by integrating equation from u c (x = 0) = 0 and du c dx | x=0 = 0. The integration can be performed by using a standard backward finit difference scheme, avoiding the solution of any linear system. Now the solution of the homogeneous equation can be expressed as The coefficient α and β can be easily computed by enforcing the model boundary conditions: because u c (x = 0) = u h 1 (x = 0) = 0 and u h 2 (x = 0) = 1, and The advantage of such procedure lies in its simplicity to be implemented in any computing platform even when considering thousands of millions of discrete values [7,8,10].

Enriched Parametric Models
Many materials exhibit a random microstructure, needing the solutions of different realizations of the microstructure. Such situation is also encountered for the numerical determination of the RVE size, see e.g. [22,23,41]. Other times the microstructure is known everywhere, but as it could evolve in the macroscopic scale, one should solve a homogenization problem at each integration point (or nodal position) at the macroscopic scale. In this section we are describing an efficien way of circumventing this computational complexity.
For the sake of simplicity, a two dimensional thermal model define in ⊂ 2 is considered. It is characterized by a microstructure that evolves in the macroscopic scale (X). Thus, at each position X, a representative volume element rve (X) = [0, L] 2 can be extracted.
In the following, points in rve (X) are denoted by x. We assume that the microstructure in each rve (X) consists of M × M cells, each one characterized by a value of the thermal conductivity.
If h = L M , the cells boundaries are define by the seg- If we defin the characteristic function χ x k (x), k = 1, . . . , M by and analogously χ y k (y), k = 1, . . . , M by any cell C ij with i, j = 1, . . . , M, can be define by C ij = χ x i (x) × χ y j (y). Using this notation the thermal conductivity at the microscopic level can be expressed in the separated form: In what follows for the sake of simplicity we are assuming that the thermal conductivity is isotropic in each cell, assumption that allows writing Now, two possibilities can be considered for performing the homogenization: Now, by solving the resulting thermal model only once, but define in a space involving 2 + M × M dimensions, one could have direct access to the temperature at each point for any value of the thermal conductivities in each cell composing the microscopic representative element volume. Consequently, when we move from one position X to another one X , as soon as the different cell thermal conductivities are known at this new position X , the temperature fiel is determined from (46). This allows the computation of the homogenized thermal conductivity without the necessity of solving a thermal model at each location in the macroscopic domain. The price to be paid is the solution of a thermal model in a multidimensional space, but this solution, thanks to the proper generalized decomposition can be performed, only once and off-line.

Remark 3
As just argued in Sect. 2, the fact that the thermal model does not contain derivatives with respect to the thermal conductivity, implies that the introduction of the different cell conductivities as extra-coordinates has not a significan impact in the computational efforts needed to construct the separated representation. However, the number of terms involved in the separated representation, N , can increase significantl with the number of extra-coordinates introduced in the model. In [3] the solution of models containing hundreds of dimensions was performed in some minutes using a standard laptop, however, the solution of the models just described needed hundreds of terms and computations ranging from some hours to some days. In any case, as indicated previously, this multidimensional problem should be solved only once and off-line.

Advanced Non-concurrent Non-linear Homogenization
Until now only linear behaviors were addressed, however in practice the most serious difficultie appear as soon as nonlinear behaviors are considered. In this section, we focus on this issue, and more precisely in two scenarii: The firs one concerns non-linear behaviors independent on the thermal history, and the second one consider the more general case in which the thermal properties at a certain point and time depends not only on the present temperature, but also on its entire thermal history.

History Independent Non-linear Behaviors
Using the notation introduced previously we assume that in a certain representative volume element rve (X) located at position X, the temperature gradient G(X) is known (from the previous homogenized macroscopic thermal model solution). We are interested in knowing the associated heat flu Fig. 2 Concurrent homogenization of non-linear history independent thermal models Q(X) at that position. In the present case, the non-linearity is coming from the fact that the thermal conductivity depends on the temperature itself, i.e. k(T ). Usual concurrent homogenization proceeds as illustrated in Fig. 2. The macroscopic temperature and the macroscopic temperature gradient are used for computing at the microscopic level the averaged heat flu that corresponds to the macroscopic heat flux If one solves the thermal model in rve (X) it is easy to prove that ∇T (x) = G.
The main difficult lies in the fact that as soon as the gradient of temperature changes at position X a new thermal model must be solved in the associated representative volume element rve (X). In the linear case, only a 3D thermal model for each space coordinate must be solved in each representative volume element rve (X). But now, as just argued, the complexity increases because the dependence of the homogenized thermal properties on the thermal state and not only on the microstructure.
For simplifying the solution process, the thermal model in rve (X) could be solved for any macroscopic gradient G. Thus, as soon as the macroscopic gradient is known, the macroscopic heat flu Q could be obtained without the necessity of solving a new thermal model. The main idea lies in introducing the components of G (denoted by G 1 , G 2 , G 3 ) as extra-coordinates in the thermal model.
The main issue is how to transfer G from the boundary conditions to the partial differential equation itself. For this purpose, a new temperature fiel θ(x) is define as: The thermal model (47) becomes: by applying the proper generalized decomposition as described in Sect. 2. In mechanics, (for example non linear elasticity), this procedure suffice for determining the macroscopic stress from the macroscopic deformation tensor E. This approach is sufficien because prescribing a displacement fiel u = E · x on the boundary of rve ensures that the resulting microscopic displacement fiel verifie = E. Obviously, a displacement given by E · x + U 0 (with U 0 a uniform translation) could be prescribed on the boundary without modifying the solution because = E is still satisfie as the elastic constants depend on the strain but not on the displacement. In mechanical models only the gradient of the displacement plays a role. Thus, if a new displacement fiel is define as: The introduction of the resulting expression of u(x) into the non-linear elastic problem, would make the translation fiel U 0 disappears because (i) it is constant and then it vanishes when derivatives apply, and (ii) the elastic coefficient only depends on the gradient of the displacement (strain) and not on the displacement itself.
Other models however exhibit such a dependence on the fiel itself. For example in thermal models, the Fourier's law relates the temperature gradient with the heat flu trough the material conductivity that usually depends on the temperature itself. The same situation can be found in the Fick's diffusion law relating the concentration gradient and the chemical species flu es through the so-called diffusion coefficien that is usually assumed depending on the concentration itself. In all these cases the procedure just described for computing the macroscopic flu from the macroscopic gradient (in our case Q from G) fails: Prescribing the temperature G · x + T 0 on the boundary of rve ensures condition ∇T = G. However, the validity of the computed Q cannot be guaranteed because the computed temperature fiel depends on T 0 , due to the dependence of the thermal conductivity on the temperature itself.
This problem of transferring the macroscopic temperature to the microscopic problem is addressed in [36]. In this reference, the microscopic temperature over the RVE is linked to the macroscopic temperature through an equation which expresses the consistency of the stored heat at the macro and the micro level.
So we can see that as in [36], the data of the macroscopic problem are composed with the macroscopic temperature and the gradient of temperature. Therefore, the temperature fiel is searched under the following form, after applying the change of variable for transferring the boundary condition to the equation itself:

History Independent Non-linear Behaviors: A Simplifie Procedure
We have just proven that the accurate homogenization of non-linear behaviors, in absence of internal variables (and then in absence of history dependencies), needs to transfer two key information from the macroscopic scale to the microscopic one: (i) the gradient of the macroscopic temperature and (ii) the macroscopic temperature itself. The resulting microscopic problem is non-linear (the thermal properties depend on the microscopic temperature) and as just described, one equation links the microscopic temperature to the macroscopic one. Even if the procedure is conceptually simple, its numerical implementation deserves certain difficultie since the solution is nonlinear with respect to the gradient of the macro-scopic temperature. Thus, many authors proposed a simplifie modeling in which the thermal properties at the microscopic level (that depend on the existing local temperature) are frozen to values corresponding to the macroscopic temperature. That is, the thermal conductivity at any point in the microscopic RVE is fi ed according to the macroscopic temperature, and then it does not evolve with the local microscopic temperature. This approach has been rigorously justifie from the asymptotic expansion method for materials with periodic microstructure in [42]. In a few words, in this paper, the authors have shown that for a given 0-th order temperature-i.e. macroscopic temperature-the effective conductivity of the non-linear temperature dependent problem is equal to that of the linear problem with thermal properties fi ed at the macroscopic temperature. Therefore this paper provides a justificatio of an approach which was widely applied in the context of computational homogenization of non-linear thermal models, see for example [1,28,39]. Now, with the thermal properties frozen at the microscopic scale, the non-linearity of the microscopic model disappears. Thus, the microscopic solution only needs the prescription of the macroscopic temperature gradient G(X) as soon as the thermal conductivity is assumed known everywhere in the microscopic representative volume.
The simplest possibility for computing a general solution is solving the microscopic linear models define in Sect. 3.1. (two in the 2D case and 3 in the 3D case) for the material thermal conductivity related to any macroscopic temperature (X). The solution of these microscopic problems (T 1 (x), T 2 (x) and T 3 (x) in Sect. 3.1.) for any macroscopic temperature suggests the separated representation: These solutions allow computing the homogenized conductivity tensor for any value of the macroscopic temperature . Through the dependence of the conductivity tensor to the macroscopic temperature, this homogenization problem displays similarities with homogenization in nonlinear elasticity, for which non-concurrent multiscale approaches have been recently proposed [40,41,43]. In these references, the macroscopic nonlinear behavior is computed after the discretization of the macroscopic strain space. Then, in a second step, a continuous representation is determined by an interpolation technique. It can be seen that our approach provides automatically this interpolation and therefore enables to compute in an elegant and efficien way the homogenized behavior of non linear models.

Remark 4
The gradient of the temperature field given by (53) allows computing the homogenized macroscopic conductivity tensor within a non-concurrent framework. However, because of the non-linearity of the macroscopic homogenized thermal model (i.e. the macroscopic homogenized thermal conductivity depends on the macroscopic temperature) we should determine an expression of the macroscopic tangent conductivity. This can be easily obtained from the derivatives of the different components of the homogenized conductivity tensor with respect to the macroscopic temperature. From (53) we can actually notice that for this purpose it suffice to take the derivative of the spatial gradient of the microscopic temperatures ∇T 1 , ∇T 2 and ∇T 3 involved in the definitio of K(X) with respect to the macroscopic temperature . Obviously, this procedure is very easy to use, whereas in [40,43] the tangent operators are computed once the interpolation of the solutions obtained in the discretized macroscopic space is built.

History Dependent Non-linear Behaviors
In the previous section we assumed that the thermal model coefficien depends on the instantaneous temperature, but for the sake of generality we extend the previous procedure (which consists in introducing all the thermal loads as extracoordinates) to models in which the thermal parameters also depend on the thermal history. This scenario is the equivalent of mechanical models involving internal variables, as elastoplasticity for example.
Usual concurrent homogenization proceeds as illustrated in Fig. 3. The macroscopic temperature and the macroscopic temperature gradient histories are used for computing at the microscopic level the averaged heat flu at present time that corresponds to the macroscopic heat flux In what follows, we assume, without looking for any physical interpretation, that the thermal conductivity at a certain location x ∈ rve (X) depends on the thermal history at that position, i.e. k = k(T (τ )), 0 ≤ τ ≤ t . Now, as in the previous case, we should solve the thermal model define in rve (X), but in this case the thermal model must be integrated in the whole time interval 0 ≤ t ≤ t max : We would like to introduce the components of G(t) and T 0 (t) as extra-coordinates in the partial differential equation describing the thermal model. However, as these potential candidates become coordinates, which still depend on the time, they cannot be at present assumed as coordinates in the model (the model coordinates must be independent).
Let G i (t), i = 1, 2, 3 be the components of vector G(t). We could approximate the evolution of each one of these components using a set of discrete values. The simplest choice consists of using a simple linear finit element interpolation: and N j (t) the standard linear finit element shape functions.
We proceed in a similar way for approximating T 0 (t) Now, the discrete values G i j , i = 1, 2, 3 and j = 1, . . . , M; and (T 0 ) j , j = 1, . . . , M, all them being independent, can be considered as extra-coordinates in the transient thermal model. By solving it only once, and possibly off-line, we have direct access to the solution of the thermal model for any prescribed macroscopic history of the thermal load G(t) and T 0 (t). The price to be paid is the solution of a transient model in a space of dimension 3 + 1 + 4 × M.
As proposed earlier, the solution of the resulting multidimensional model is searched in a separated form (to allow a linearly scaling complexity with respect to dimensionality). Thus, the solution results in: y, z, t, G 1  1 , . . . , G 1  M , . . . , G 3  1 , . . . , Remark 5 As in the previous case, prior to introduce these extra-coordinates in the transient thermal model, they should be transferred from the boundary condition into the PDE itself. This transfer is performed as in the former case by applying the change of variable θ(

History Dependent Non-linear Behaviors: A Simplifie Procedure
Even in the case of history dependent models, the complexity of the computation can be drastically reduced when considering the simplifying hypothesis proposed in Sect. 4.2 (thermal properties frozen to the values dictated by the macroscopic scale). As just described, the proper solution of this kind of models needs the introduction of many time evolutions: Those related to each component of the macroscopic temperature gradient and the one related to the time evolution of T 0 that is prescribed on the domain boundary to enforce that the microscopic averaged temperature equals the macroscopic one at any time. Thus, the use of M values for describing the time evolution of these fields implies the introduction, in the general 3D case, of 3 × M + M extra-coordinates as just described in the previous section. However, if the thermal properties are assumed evolving as dictated by the macroscopic temperature evolution, the non-linearity at the microscopic level disappears, and it suffice to compute the homogenized conductivity tensor as described in Sect. 4.1. for any evolution of the macroscopic temperature, i.e.: y, z, t, 1 , . . . , M ) that reduces the dimensionality of the problem with respect to problem (57) by removing 3 × M coordinates.

Homogenization of Linear Models
In [11] the use of separated representations for computing the homogenized behavior of linear thermal models was already addressed. In this paper we retained the use of such technique (see Sect In what follows, we focus on the homogenization of materials whose microstructure can be accurately described by a microscopic representative volume element consisting of a certain number of cells, according to the description addressed in Sect. 3.4. For the sake of simplicity in what follows, and without loss of generality, only the 2D case is considered. More precisely, we assume a 2D microscopic RVE vre = (0, 1) × (0, 1) composed of 5 × 5 cells whose conductivity can take any value in the interval [k min , k max ], with k min = 1 and k max = 2 in the following.
According to the procedure described in Sect. 3.4, the linear thermal model is solved in rve for two different boundary conditions: T (x ∈ ∂ rve ) = x in the firs case and T (x ∈ ∂ rve ) = y in the second case. The resulting solutions are denoted T 1 (x) and T 2 (x) respectively. As described in Sect. 3.2 these solutions suffic for definin the homogenized conductivity tensor. On the other hand, in order to compute solutions whose validity does not depend on the particular value of the conductivity of each cell, both solutions are computed by considering the thermal conductivities of all the cells as extra-coordinates, i.e. both solutions are searched under the form: for j = 1, 2. Thus, instead of solving two thermal models for each different microstructure, the solution of two multidimensional models (define in a space of dimension 2 + 5 × 5 = 27) suffice to give access to the temperature fiel for any possible microstructure define on this grid. Moreover, thanks to the PGD, the solution of these multidimensional models does not introduce major difficulties With these solutions computed, and as soon as the microstructure is given, (value of the conductivity in each cell), the homogenized conductivity tensor can be computed.
In order to validate the solutions, a particular microstructure is define by taking random values of the components of the conductivity tensor k ij that we denote by k g ik (see Fig. 4). Then, the solution is obtained by particularizing the general PGD solution Temperatures T 1 (x, y) and T 2 (x, y) obtained using the PGD are compared with finit element solutions define on the same microstructure. The PGD results were obtained by considering a 101 nodes F.E. discretization for functions of x and y and a three nodes discretization for the functions depending on the different conductivities. The classical F.E. solution was obtained by mean of a 100 × 100 quadrilateral mesh (which leads to the same geometrical accuracy as the PGD). Figure 5 depicts both solutions: the one related to the PGD on the left and the one associated with the application of the finit element on the right. The comparison of both solutions illustrates the performance of the PGD for computing the solution (the norm of the difference of both solutions is lower than 1 percent). It is important to recall that such solution define in a space of dimension 27 cannot be computed by using a standard mesh based discretization technique. It was also verifie that an increasing number of terms N in the sum (59) makes the norm of the difference between the PGD and the FEM solutions decrease. Finally Fig. 6 depicts the functions X i (x), Y i (y), i 11 (k 1,1 ) and i 55 (k 5,5 ) involved in the decomposition (59) for N = 83.

Homogenization of Nonlinear History Independent Models
In this section we focus on the homogenization of a simple non-linear thermal model. The representative volume element rve = (0, 1) × (0, 1) is composed of two materials of different conductivities k 1 and k 2 . The zone ω = (x l , x r ) × (y d , y t ) ⊂ involves a material with conductivity k 2 whereas the material occupying the complementary region − ω involves a material of conductivity k 1 (see Fig. 7). The region of conductivity k 2 can be expressed by the following separate form: where both characteristic functions are define by: and χ y (y) = 1 if y d < x < y t 0 otherwise (63) In the computations, we considered the following values for the geometrical parameters: x l = 0.4, x r = 0.6, y d = 0.4 and y t = 0.6. The microscopic conductivities are supposed to be linearly temperature dependent. As justifie in Sect. 4.1 , their values are obtained considering that the microscopic temperature equals the macroscopic temperature i.e. k 1 = 10 + 9 × and k 2 = 90 + , with ∈ [0, 20]. One can notice that this choice for the conductivities implies that when = 10 both conductivities take the same value and the microstructure becomes homogeneous.
With both conductivities define and frozen during the microscopic step, the solution of the microscopic thermal model can be computed as previously for two different boundary conditions T (x ∈ ∂ rve ) = x and T (x ∈ ∂ rve ) = y, leading to the homogenized conductivity tensor. If one proceeds in this manner, as soon as the macroscopic temperature evolves in time, the homogenized conductivity tensor should be recomputed for the thermal properties related to the new macroscopic temperature. As discussed in Sect. 4.1 this procedure is very expensive from a computational point of view. An alternative lies in the solution of those microscopic models for any value of the macroscopic temperature, i.e. the calculation of T 1 (x, y, ) and T 2 (x, y, ) allowing to solve two higher dimensional problems from which the temperature field can be com-puted. Then, the homogenized conductivity tensor can be obtained for any macroscopic temperature.
In order to illustrate the feasibility of this approach, both solutions T 1 (x, y, ) and T 2 (x, y, ) are computed. Then, the resulting temperature fiel for a given macroscopic tem-perature = g (obtained by particularizing the general solution, T 1 (x, y, = g )) is compared with the one computed by applying the finit element method (with the thermal properties define at temperature g ). These results are depicted in Fig. 8 (left PGD, right FE) proving the ability of the multidimensional model to represent accurately the solution of any particular scenario, since the norm of the difference between the PGD and FE solutions is very small. Finally Fig. 9 shows the most significan functions X i (x), Y i (y) and θ i ( ) (i = 1, . . . , 100) related to the proper generalized decomposition. In this simulation, 1001 nodes were considered for the discretization of functions of x and y and 200 for those depending on the macroscopic temperature .
From these results we can conclude on the ability of the PGD to capture the main features of the homogenized thermal model related to non-linear behaviors. By solving a slightly higher dimensional model involving an extracoordinate (the macroscopic temperature at the location of the RVE rve ), the necessity of solving a microscopic model for each value of the macroscopic temperature can be avoided.

Homogenization of Non-linear History Dependent Models
In this case we consider the scenario described in the previous example, but also assume that the thermal properties depend on the thermal history. Thus, a transient thermal model should be solved. Moreover, as the thermal properties at each time depend on the thermal history, a microscopic thermal model should be solved at each time step. In addition, the thermal properties in the RVE are assumed fi ed (equal to (τ ), 0 ≤ τ ≤ t), following the simplifie model described in Sect. 4.4. To alleviate this solution one could try to solve the microscopic thermal model for any thermal history, from which the thermal properties in rve could be computed and frozen. Following the description given in Sect. 4.4, the approximation of (t) in the whole time interval was interpolated from a polynomial of degree 4 that can be expressed from 5 nodal values of , denoted by i , i = 1, . . . , M where M was set to M = 5 in our simulations.
Finally, the solutions for the two different boundary conditions are searched in the separated form T 1 (x, y, 1 , . . . , M ) and T 2 (x, y, 1 , . . . , M ) according to (58). In this case, two problems define in a space involving x, y, t and the f ve temperatures i , i = 1, . . . , 5 (8 dimensions) need to be solved in order to compute the homogenized conductivity tensor.
Functions depending on the space (x and y) where discretized by employing 101 nodes, whereas functions of time were integrated using also a coarse description involving 101 nodes. The length of the space and time interval were fi ed arbitrarily to one. Finally the f ve extra-coordinates related to the intermediate macroscopic temperatures i , i = 1, . . . , 5 were discretized by employing 10 nodes uniformly distributed in the interval of variation of those temperatures ([0, 10] here).
After computing both solutions T 1 (x, y, 1 , . . . , M ) and T 2 (x, y, 1 , . . . , M ) for any macroscopic thermal history described by the fourth order polynomial just defined these solutions are particularized for a thermal history given by a linear thermal evolution from zero to 10 and then it decreases again linearly until reaching a null temperature at the fina time. This problem, as soon as the thermal history is fi ed, was also solved by mean of the finit element method. To validate the proposed approach we depicts in Fig. 10 the solutions computed at the time t = 0.25, t = 0.5, t = 0.75 and t = 1 obtained by particularizing the solution T 1 (x, y, 1 , . . . , M ) for the macroscopic thermal history just defined i.e. T 1 (x, y, 1 = 0, 2 = 5, 3 = 10, 4 = 5, 5 = 0). The corresponding solutions, computed by applying the finit element solution and enforcing T (x ∈ ∂ rve , t) = x are depicted in Fig. 11. Both solutions agree in minute detail.
These results prove the ability of the proposed technique for solving non-linear homogenization problems related to models involving internal variables implying a history dependance. Obviously a deeper analysis should be carried out for validating this approach in the case of real dependences of thermal parameters on the thermal history as well as for analyzing the sensibility of the results to the regularity in that thermal dependence.

Conclusions
This paper enlarges the domain of applicability of proper generalized decomposition to computational homogenization. In particular it proves that in the linear case it could be possible to defin the homogenized thermal conductivity tensor for a generic microstructure composed of cells by assuming the conductivity of those cells as extra-coordinates. Even if the solution of the resulting multi-dimensional model could be expensive, it is important to recall that this solution should be carried out only once. Using a laptop and a handmade code in Matlab we solved models involving a 2D microstructure composed of 10 × 10 cells. This model was define in a space involving 2 + 10 × 10 = 102 dimensions, and its solutions were performed without major difficulties The computational resources available at present allows expecting new developments in this direction.
In the non-linear case, we proposed a fully non-concurrent homogenization procedure by solving the microstructure for a generic macroscopic solicitation (macroscopic All these developments were possible thanks to the introduction of a certain number of extra-coordinates in the models. The redoubtable curse of dimensionality related to multidimensional models was circumvented by employing separated representations, whose numerical complexity scales linearly with the dimension of the space in which the model is defined instead the exponential growing characteristic of mesh based representations. All the procedures were described by assuming a fully separated representation involving the separation of each one of the space coordinates x, y and z. However, the separated representation of the physical space is not mandatory, and consequently all the techniques proposed throughout this paper remain applicable if one considers space functions X i (x) instead of X i (x) · Y i (y) · Z i (z). The only difference when considering the former grouped physical space representation is that the computation of the space functions requires the solution of 3D problems instead of the solution of three one-dimensional problems involved in fully separated space representations.
The results presented in this work could be extended to mechanical models. This extension constitutes a work in progress. The numerical examples addressed here are quite simple, but the objective was more proving the potential of the Proper Generalized Decomposition in the fiel of computational homogenization, rather than addressing models of practical interest.