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

Hajer Lamari, Amine Ammar, Patrice Cartraud, Grégory Legrain, Francisco Chinesta, Frédéric Jacquemin

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

## ► To cite this version:

Hajer Lamari, Amine Ammar, Patrice Cartraud, Grégory Legrain, Francisco Chinesta, et al. Routes for Efficient Computational Homogenization of Nonlinear Materials Using the Proper Generalized Decompositions. Archives of Computational Methods in Engineering, 2010, 17 (4), pp.373-391. 10.1007/s11831-010-9051-4. hal-01528268

HAL Id: hal-01528268

https://hal.science/hal-01528268

Submitted on 28 May 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Routes for Efficient Computational Homogenization of Nonlinear Materials Using the Proper Generalized Decompositions

H. Lamari · A. Ammar · P. Cartraud · G. Legrain · F. Chinesta · F. Jacquemin

H. Lamari · P. Cartraud · G. Legrain
GeM, UMR CNRS-Centrale de Nantes-Université de Nantes,
1 rue de la Noe, BP 92101, 44321 Nantes cedex 3, France
e-mail: Hajer.Lamari@ec-nantes.fr

P. Cartraud
e-mail: Patrice.Cartraud@ec-nantes.fr

G. Legrain
e-mail: Gregory.Legrain@ec-nantes.fr

A. Ammar
Arts et Métiers ParisTech, 2 boulevard du Ronceray, BP 93525,
49035 Angers cedex 01, France
e-mail: Amine.AMMAR@ensam.eu

F. Chinesta
EADS Corporate Foundation International Chair, GeM,
UMR CNRS-Centrale de Nantes-Université de Nantes,
1 rue de la Noe, BP 92101, 44321 Nantes cedex 3, France
e-mail: Francisco.Chinesta@ec-nantes.fr

F. Jacquemin
GeM, UMR CNRS-Centrale de Nantes-Université de Nantes,
58, rue Michel Ange, BP 420, 44606 Saint Nazaire cedex, France
e-mail: Frederic.Jacquemin@univ-nantes.fr

**Abstract** 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-

linear models). In this paper we present some preliminary results concerning the application of Proper Generalized Decompositions—PGD—for addressing the efficiency 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 first steps in a very ambitious objective, many original ideas are launched that could be at the origin of impressive progresses.

## 1 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 higher—macro—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 identified 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 solu-
tion using strategies such as X-FEM and G-FEM, see e.g.
[17] and [38] respectively. One of the major feature of hier-
archical 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 re-
quired when the macroscopic behavior can not be obtained
explicitly, which occurs as soon as microscopic nonlinear
behaviors are involved. The area where the microscopic de-
scription 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 be-
havior. These methods were reviewed in [24] in the contin-
uum mechanics framework. The most popular approach is
the FE² 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 at-
tached.

As noticed in [29], both the hierarchical and concurrent
simulation approaches offer cost savings over large scale di-
rect numerical simulation of the microstructure. However
efforts are directed towards improving their numerical ef-
ficien 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 reso-
lution 3D images of material microstructures. Thus an im-
age 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² ap-
proach, one main limitation is its computational cost, since
several microscopic problems have to be solved at each in-
tegration point of the macroscopic domain. Therefore at-
tempts have been made to build hierarchical approaches for
nonlinear materials. One can thus defin approximate non-
linear constitutive relations, see e.g. [30]. Another possibil-
ity 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² 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 multidimen-
sional models encountered in the kinetic theory description
of complex fluids involving steady state and transient lin-
ear and non-linear models [2–4, 32]. Time-space separated
representations were originally introduced many years ago
by P. Ladeveze for addressing complex time-dependent non-
linear models within the LATIN framework (see [26] and
the references therein). Some recent works proving the po-
tentiality of this approach can be found in [27, 33] among
many others. In [2, 3] we generalized the separated represen-
tation for treating general multidimensional models. Obvi-
ously, a firs possibility of separated representations consists
of a separated representation of the physical space, by writ-
ing a generic function as a finit sums of functional prod-
ucts whose involved functions only depends on a single co-
ordinate. This strategy was applied in [11] in high resolu-
tion linear homogenization. The computational complexity
of separated representations scales linearly with the dimen-
sion of the space instead the exponential growing character-
istic 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 mechan-
ics that generally involve space and time as coordinates can
be transformed, as we illustrate later, into highly multididi-
mensional models. Obviously, the resulting models suffers
of the curse of dimensionality if mesh-based representations
are considered. However, the Proper Generalized Decompo-
sition 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 equa-
tion without knowing the microscopic material thermal con-
ductivity (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 conductiv-
ity 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 coordi-
nate, 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 computa-
tional science [34].
In the field of homogenization, two simpler scenarios can be considered: (i) define 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 first 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 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 define 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.

## 2 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:

$$ \frac{\partial u}{\partial t} - k \Delta u - f = 0 \qquad (1) $$

where $(\mathbf{x}, t, k) \in \Omega \times I \times \mathcal{I}$ and for the sake of simplicity the source term is assumed constant, i.e. $f = \text{cte}$. Because the conductivity is considered unknown, it is assumed as a new coordinate defined in the interval $\mathcal{I}$. 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:

$$ u(\mathbf{x}, t, k) \approx \sum_{i=1}^{i=N} X_i(\mathbf{x}) \cdot T_i(t) \cdot K_i(k) \qquad (2) $$

In what follows we are assuming that the approximation at iteration $n$ is already done:

$$ u^n(\mathbf{x}, t, k) = \sum_{i=1}^{i=n} X_i(\mathbf{x}) \cdot T_i(t) \cdot K_i(k) \qquad (3) $$

and that at present iteration we look for the next functional product $X_{n+1}(\mathbf{x}) \cdot T_{n+1}(t) \cdot K_{n+1}(k)$ that for alleviating the notation will be denoted by $R(\mathbf{x}) \cdot S(t) \cdot 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 fixed 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(\mathbf{x})$. From the just updated $R(\mathbf{x})$ and $W(k)$ we can update $S(t)$, and finally from the just computed $R(\mathbf{x})$ and $S(t)$ we compute $W(k)$. The procedure continues until reaching convergence. The converged functions $R(\mathbf{x})$, $S(t)$ and $W(k)$ allow defining the searched functions: $X_{n+1}(\mathbf{x}) = R(\mathbf{x})$, $T_{n+1}(t) = S(t)$ and $K_{n+1}(k) = W(k)$. We are illustrating each one of the just referred steps.

Computing $R(\mathbf{x})$ from $S(t)$ and $W(k)$ We consider the global weak form of (1):

$$ \int_{\Omega \times I \times \mathcal{I}} u^* \left( \frac{\partial u}{\partial t} - k \Delta u - f \right) d\mathbf{x} dt dk = 0 \qquad (4) $$

where the trial and test functions write respectively:

$$ u(\mathbf{x}, t, k) = \sum_{i=1}^{i=n} X_i(\mathbf{x}) \cdot T_i(t) \cdot K_i(k) + R(\mathbf{x}) \cdot S(t) \cdot W(k) \qquad (5) $$

and

$$ u^*(\mathbf{x}, t, k) = R^*(\mathbf{x}) \cdot S(t) \cdot W(k) \qquad (6) $$

Introducing (5) and (6) into (4) it results

$$ \int_{\Omega \times I \times \mathcal{I}} R^* \cdot S \cdot W \cdot \left( R \cdot \frac{\partial S}{\partial t} \cdot W - k \cdot \Delta R \cdot S \cdot W \right) d\mathbf{x} dt dk = - \int_{\Omega \times I \times \mathcal{I}} R^* \cdot S \cdot W \cdot \mathcal{R}^n d\mathbf{x} dt dk \qquad (7) $$

where $\mathcal{R}^n$ define the residual at iteration $n$ that reads:

$$ \mathcal{R}^n = \sum_{i=1}^{i=n} X_i \cdot \frac{\partial T_i}{\partial t} \cdot K_i - \sum_{i=1}^{i=n} k \cdot \Delta X_i \cdot T_i \cdot K_i - f \qquad (8) $$
Now, being known all the functions involving the time
and the parametric coordinate, we can integrate (7) in their
respective domains $I \times \mathfrak{S}$. Integrating in $I \times \mathfrak{S}$ and taking
into account the following notations

$$
\begin{align*}
w_1 &= \int_{\mathfrak{S}} W^2 dk \qquad s_1 = \int_I S^2 dt \qquad r_1 = \int_{\Omega} R^2 d\mathbf{x} \\
w_2 &= \int_{\mathfrak{S}} kW^2 dk \qquad s_2 = \int_I S \cdot \frac{dS}{dt} dt \\
r_2 &= \int_{\Omega} R \cdot \Delta R d\mathbf{x} \qquad w_3 = \int_{\mathfrak{S}} W dk \qquad s_3 = \int_I S dt \\
r_3 &= \int_{\Omega} R d\mathbf{x} \qquad w_4^i = \int_{\mathfrak{S}} W \cdot K_i dk \tag{9} \\
s_4^i &= \int_I S \cdot \frac{dT_i}{dt} dt \qquad r_4^i = \int_{\Omega} R \cdot \Delta X_i d\mathbf{x} \\
w_5^i &= \int_{\mathfrak{S}} kW \cdot K_i dk \qquad s_5^i = \int_I S \cdot T_i dt \\
r_5^i &= \int_{\Omega} R \cdot X_i d\mathbf{x}
\end{align*}
$$

Equation (7) reduces to:

$$
\begin{align}
& \int_{\Omega} R^* \cdot (w_1 \cdot s_2 \cdot R - w_2 \cdot s_1 \cdot \Delta R) \, d\mathbf{x} \nonumber \\
&= - \int_{\Omega} R^* \cdot \left( \sum_{i=1}^{i=n} w_4^i \cdot s_4^i \cdot X_i - \sum_{i=1}^{i=n} w_5^i \cdot s_5^i \cdot \Delta X_i \right. \nonumber \\
& \qquad \left. - w_3 \cdot s_3 \cdot f \right) \, d\mathbf{x} \tag{10}
\end{align}
$$

Equation (10) define an elliptic steady state boundary value
problem that can be solved by using any discretization tech-
nique operating on the model weak form (finit elements,
f finite volumes ...). Another possibility consists in coming
back to the strong form of (10):

$$
w_1 \cdot s_2 \cdot R - w_2 \cdot s_1 \cdot \Delta R \\
= - \left( \sum_{i=1}^{i=n} w_4^i \cdot s_4^i \cdot X_i - \sum_{i=1}^{i=n} w_5^i \cdot s_5^i \cdot \Delta X_i - w_3 \cdot s_3 \cdot f \right) \tag{11}
$$

that could be solved by using any collocation technique (fi-
nite differences, SPH ...).

Computing $S(t)$ from $R(\mathbf{x})$ and $W(k)$ In the present case the test function writes:

$$
u^*(\mathbf{x}, t, k) = S^*(t) \cdot R(\mathbf{x}) \cdot W(k) \quad (12)
$$

Now, the weak form reads

$$
\int_{\Omega \times I \times \mathfrak{S}} S^* \cdot R \cdot W \cdot \left( R \cdot \frac{\partial S}{\partial t} \cdot W - k \cdot \Delta R \cdot S \cdot W \right) d\mathbf{x} dt dk \\
= - \int_{\Omega \times I \times \mathfrak{S}} S^* \cdot R \cdot W \cdot \mathcal{R}^n d\mathbf{x} dt dk \tag{13}
$$

that integrating in the space $\Omega \times \mathfrak{S}$ and taking into account
the notation (9) results:

$$
\begin{align}
& \int_I S^* \cdot \left( w_1 \cdot r_1 \cdot \frac{dS}{dt} - w_2 \cdot r_2 \cdot S \right) dt \nonumber \\
&= - \int_I S^* \cdot \left( \sum_{i=1}^{i=n} w_4^i \cdot r_5^i \cdot \frac{dT_i}{dt} - \sum_{i=1}^{i=n} w_5^i \cdot r_4^i \cdot T_i \right. \nonumber \\
& \qquad \left. - w_3 \cdot r_3 \cdot f \right) dt \tag{14}
\end{align}
$$

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:

$$
\begin{align}
w_1 \cdot r_1 \cdot \frac{dS}{dt} - w_2 \cdot r_2 \cdot S \nonumber \\
= - \left( \sum_{i=1}^{i=n} w_4^i \cdot r_5^i \cdot \frac{dT_i}{dt} - \sum_{i=1}^{i=n} w_5^i \cdot r_4^i \cdot T_i - w_3 \cdot r_3 \cdot f \right) \tag{15}
\end{align}
$$

than can be solved by using backward finit differences, or
higher order Runge-Kutta schemes, among many other pos-
sibilities.

Computing $W(k)$ from $R(\mathbf{x})$ and $S(t)$ In the present case the test function writes:

$$
u^*(\mathbf{x}, t, k) = W^*(k) \cdot R(\mathbf{x}) \cdot S(t) \quad (16)
$$

Now, the weak form reads

$$
\int_{\Omega \times I \times \mathfrak{S}} W^* \cdot R \cdot S \cdot \left( R \cdot \frac{\partial S}{\partial t} \cdot W - k \cdot \Delta R \cdot S \cdot W \right) d\mathbf{x} dt dk \\
= - \int_{\Omega \times I \times \mathfrak{S}} W^* \cdot R \cdot S \cdot \mathcal{R}^n d\mathbf{x} dt dk \quad (17)
$$

that integrating in the space $\Omega \times I$ and taking into account
the notation (9) results:

$$
\begin{align}
& \int_{\mathfrak{S}} W^* \cdot (r_1 \cdot s_2 \cdot W - r_2 \cdot s_1 \cdot k \cdot W) \, dk \nonumber \\
&= - \int_{\mathfrak{S}} W^* \cdot \left( \sum_{i=1}^{i=n} r_5^i \cdot s_4^i \cdot K_i - \sum_{i=1}^{i=n} r_4^i \cdot s_5^i \cdot k \cdot K_i \right. \nonumber \\
& \qquad \left. - r_3 \cdot s_3 \cdot f \right) \, dk \tag{18}
\end{align}
$$
Equation (18) does not involve any differential operator. The strong form of (18) reads:

$$ (r_1 \cdot s_2 - r_2 \cdot s_1 \cdot k) \cdot W = - \left( \sum_{i=1}^{i=n} (r_5^i \cdot s_4^i - r_4^i \cdot s_5^i \cdot k) \cdot K_i - r_3 \cdot s_3 \cdot f \right) \quad (19) $$

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 = \sum_{i=1}^{i=N}$ 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, \dots, 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.

# 3 Homogenization of Linear Heterogeneous Models

## 3.1 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 $\mathbf{X}$, $\Omega_{rve}(\mathbf{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 $\mathbf{k}(\mathbf{x})$ is known at each point in the microscopic domain $\mathbf{x} \in \Omega_{rve}(\mathbf{X})$.

We can defin the macroscopic temperature gradient at position $\mathbf{X}$, $\mathbf{G}(\mathbf{X})$, from:

$$ \mathbf{G}(\mathbf{X}) = \langle \mathbf{g} \rangle = \frac{1}{|\Omega_{rve}(\mathbf{X})|} \int_{\Omega_{rve}(\mathbf{X})} \mathbf{g}(\mathbf{x}) d\Omega \quad (20) $$

where $\mathbf{g}(\mathbf{x}) = \nabla T(\mathbf{x})$.

We also assume the existence of a localization tensor $\mathbf{L}(\mathbf{x}, \mathbf{X})$ such that

$$ \mathbf{g}(\mathbf{x}) = \mathbf{L}(\mathbf{x}, \mathbf{X}) \cdot \mathbf{G}(\mathbf{X}) \quad (21) $$

Now, we consider the microscopic heat flu q according to the Fourier's law

$$ \mathbf{q}(\mathbf{x}) = -\mathbf{k}(\mathbf{x}) \cdot \mathbf{g}(\mathbf{x}) \quad (22) $$

and the macroscopic counterpart $\mathbf{Q}(\mathbf{X})$ that writes:

$$ \mathbf{Q}(\mathbf{X}) = \langle \mathbf{q}(\mathbf{x}) \rangle = \langle -\mathbf{k}(\mathbf{x}) \cdot \mathbf{g}(\mathbf{x}) \rangle = \langle -\mathbf{k}(\mathbf{x}) \cdot \mathbf{L}(\mathbf{x}, \mathbf{X}) \rangle \cdot \mathbf{G}(\mathbf{X}) \quad (23) $$

from which the homogenized thermal conductivity can be defined from

$$ \mathbf{K}(\mathbf{X}) = \langle -\mathbf{k}(\mathbf{x}) \cdot \mathbf{L}(\mathbf{x}, \mathbf{X}) \rangle \quad (24) $$

As $\mathbf{k}(\mathbf{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 $\mathbf{L}(\mathbf{x}, \mathbf{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:

$$ \left\{ \begin{array}{l} \nabla \cdot (\mathbf{k}(\mathbf{x}) \cdot \nabla T^1(\mathbf{x})) = 0 \\ T^1(\mathbf{x} \in \partial \Omega_{rve}) = x \end{array} \right. \quad (25) $$

$$ \left\{ \begin{array}{l} \nabla \cdot (\mathbf{k}(\mathbf{x}) \cdot \nabla T^2(\mathbf{x})) = 0 \\ T^2(\mathbf{x} \in \partial \Omega_{rve}) = y \end{array} \right. \quad (26) $$
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Homogenization procedure of linear heterogeneous models</figcaption></figure>

and

$$ \left\{ \begin{array}{l} \nabla \cdot (\mathbf{k}(\mathbf{x}) \cdot \nabla T^3(\mathbf{x})) = 0 \\ T^3(\mathbf{x} \in \partial \Omega_{\text{rve}}) = z \end{array} \right. \qquad (27) $$

It is easy to prove that these three solutions verify

$$ \left\{ \begin{array}{l} \mathbf{G}^1 = \langle \nabla T^1(\mathbf{x}) \rangle^T = (1, 0, 0) \\ \mathbf{G}^2 = \langle \nabla T^2(\mathbf{x}) \rangle^T = (0, 1, 0) \\ \mathbf{G}^3 = \langle \nabla T^3(\mathbf{x}) \rangle^T = (0, 0, 1) \end{array} \right. \qquad (28) $$

where $(\cdot)^T$ denotes the transpose. Thus, the localization tensor results finally

$$ \mathbf{L}(\mathbf{x}, \mathbf{X}) = (\nabla T^1(\mathbf{x}) \quad \nabla T^2(\mathbf{x}) \quad \nabla T^3(\mathbf{x})) \qquad (29) $$

The resulting non-concurrent homogenization procedure is illustrated in Fig. 1. As soon as tensor $\mathbf{L}(\mathbf{x}, \mathbf{X})$ is known at each position $\mathbf{X}$, the constitutive law relating the macroscopic temperature gradient and macroscopic heat flux becomes defined

From the computational point of view the main difficult concerns the solution of the three boundary value problems using very fine 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].

## 3.2 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 finite 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

$$ T^n(x, y, z) = \sum_{i=1}^{i=n} X_i(x) \cdot Y_i(y) \cdot Z_i(z) \qquad (30) $$

and it is enriched until reaching convergence by computing the next term in the finite sums decomposition

$$ T^{n+1}(x, y, z) = T^n(x, y, z) + R(x) \cdot S(y) \cdot W(z) \qquad (31) $$

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 final $W(z)$ from $R(x)$ and $S(x)$. After reaching convergence we can assign $R(x) \to X_{n+1}(x)$, $S(y) \to Y_{n+1}(y)$ and $W(z) \to Z_{n+1}(z)$.

We assume that the microscopic conductivity can be written in a separated form, i.e. each component of the conductivity tensor $\mathbf{k}_{ij}(\mathbf{x})$ can be expressed from:

$$ \mathbf{k}_{ij}(\mathbf{x}) \approx \sum_{p=1}^{p=M} (\mathbf{k}_{ij}^1)_p(x) \cdot (\mathbf{k}_{ij}^2)_p(y) \cdot (\mathbf{k}_{ij}^3)_p(z) \qquad (32) $$

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
BVP

$$ \alpha_x \frac{d^2 R(x)}{dx^2} + \beta_x \frac{d R(x)}{dx} + \gamma_x R(x) = f_x(x) \quad (33) $$

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:

$$ \alpha_y \frac{d^2 S(y)}{dy^2} + \beta_y \frac{d S(y)}{dy} + \gamma_y S(y) = f_y(y) \quad (34) $$

and

$$ \alpha_z \frac{d^2 W(z)}{dz^2} + \beta_z \frac{d W(z)}{dz} + \gamma_z W(z) = f_z(z) \quad (35) $$

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 definition 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 defining 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 efficiency. In Sect. 3.3 we explore some possibilities for alleviating this drawback.

### 3.3 Efficient Solvers of One-Dimensional BVP

Consider the generic second order BVP implying the first order $u(x), x \in [0, L]$

$$ a \frac{d^2 u}{dx^2} + b \frac{du}{dx} + cu = f(x) \quad (36) $$

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

$$ a \frac{d^2 u^c}{dx^2} + b \frac{du^c}{dx} + cu^c = f(x) \quad (37) $$

from $u^c(x=0) = 0$ and $\frac{du^c}{dx}|_{x=0} = 0$. The integration can be performed by using a standard backward finite difference scheme, avoiding the solution of any linear system.

Now the solution of the homogeneous equation can be expressed as $u^h(x) = \alpha u_1^h(x) + \beta u_2^h(x)$, where $u_1^h(x)$ and $u_2^h(x)$ results from the backward integration of

$$ a \frac{d^2 u_1^h}{dx^2} + b \frac{du_1^h}{dx} + cu_1^h = 0 \quad (38) $$

with $u_1^h(x=0) = 0$ and $\frac{du_1^h}{dx}|_{x=0} = 1$; and

$$ a \frac{d^2 u_2^h}{dx^2} + b \frac{du_2^h}{dx} + cu_2^h = 0 \quad (39) $$

with $u_2^h(x=0) = 1$ and $\frac{du_2^h}{dx}|_{x=0} = 0$.

Thus, finally, the general solution writes $u(x) = u^c(x) + \alpha u_1^h(x) + \beta u_2^h(x)$. The coefficient $\alpha$ and $\beta$ can be easily computed by enforcing the model boundary conditions:

$$ \beta = u_l \quad (40) $$

because $u^c(x=0) = u_1^h(x=0) = 0$ and $u_2^h(x=0) = 1$, and

$$ \alpha = \frac{u_r - u^c(x=L) - u_l \cdot u_2^h(x=L)}{u_1^h(x=L)} \quad (41) $$

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].

### 3.4 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 efficient way of circumventing this computational complexity.

For the sake of simplicity, a two dimensional thermal model define in $\Omega \subset \mathbb{R}^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 $\Omega_{rve}(X) = [0, L]^2$ can be extracted.

In the following, points in $\Omega_{rve}(X)$ are denoted by $\mathbf{x}$. We assume that the microstructure in each $\Omega_{rve}(X)$ consists of $M \times M$ cells, each one characterized by a value of the thermal conductivity.

If $h = \frac{L}{M}$, the cells boundaries are defined by the segments $(x_i = (i-1) \times h, y)$, $i = 1, \dots, M+1$ and $y \in [0, L]$; and $(x, y_i = (i-1) \times h)$, $i = 1, \dots, M+1$ and $x \in [0, L]$.

If we define the characteristic function $\chi_k^x(x)$, $k = 1, \dots, M$ by

$$ \chi_k^x(x) = \begin{cases} 1 & \text{if } x_k < x < x_{k+1} \\ 0 & \text{otherwise} \end{cases} \qquad (42) $$

and analogously $\chi_k^y(y)$, $k = 1, \dots, M$ by

$$ \chi_k^y(y) = \begin{cases} 1 & \text{if } y_k < y < y_{k+1} \\ 0 & \text{otherwise} \end{cases} \qquad (43) $$

any cell $C_{ij}$ with $i, j = 1, \dots, M$, can be defined by $C_{ij} = \chi_i^x(x) \times \chi_j^y(y)$.

Using this notation the thermal conductivity at the microscopic level can be expressed in the separated form:

$$ \mathbf{k}(\mathbf{x}) = \sum_{i=1}^{i=M} \sum_{j=1}^{j=M} \mathbf{k}_{ij} \chi_i^x(x) \times \chi_j^y(y) \qquad (44) $$

In what follows for the sake of simplicity we are assuming that the thermal conductivity is isotropic in each cell, assumption that allows writing

$$ k(\mathbf{x}) = \sum_{i=1}^{i=M} \sum_{j=1}^{j=M} k_{ij} \chi_i^x(x) \times \chi_j^y(y) \qquad (45) $$

Now, two possibilities can be considered for performing the homogenization:

1. The simplest one consists in applying the procedure described in Sect. 3. Thus, three steady state thermal models should be solved at each point X on the corresponding microscopic domain $\Omega_{rve}(X)$ where the homogenized thermal properties are searched.

2. The second possibility consists in introducing the cell thermal conductivities as extra-coordinates of the thermal model, which implies the following separated form of the temperature field

$$ T(x, y, k_{1,1}, k_{1,2}, \dots, k_{M,M}) \\ \approx \sum_{l=1}^{l=N} X_l(x) \cdot Y_l(y) \cdot \prod_{p=1}^{p=M} \prod_{q=1}^{q=M} \Xi_{pq}^l (k_{p,q}) \qquad (46) $$

Now, by solving the resulting thermal model only once, but define in a space involving $2+M \times 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 field 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 significant impact in the computational efforts needed to construct the separated representation. However, the number of terms involved in the separated representation, $N$, can increase significantly 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.

## 4 Advanced Non-concurrent Non-linear Homogenization

Until now only linear behaviors were addressed, however in practice the most serious difficulty appear as soon as non-linear behaviors are considered. In this section, we focus on this issue, and more precisely in two scenarios: The first 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.

### 4.1 History Independent Non-linear Behaviors

Using the notation introduced previously we assume that in a certain representative volume element $\Omega_{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 flux
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Concurrent homogenization of non-linear history independent thermal models</figcaption></figure>

$\mathbf{Q}(\mathbf{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. $\mathbf{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 flux that corresponds to the macroscopic heat flux

If one solves the thermal model in $\Omega_{rve}(\mathbf{X})$

$$ \left\{ \begin{array}{l} \nabla \cdot (\mathbf{k}(T) \cdot \nabla T(\mathbf{x})) = 0 \\ T(\mathbf{x} \in \partial\Omega_{rve}) = \mathbf{G} \cdot \mathbf{x} \end{array} \right. \qquad (47) $$

it is easy to prove that $\langle \nabla T(\mathbf{x}) \rangle = \mathbf{G}$.

The main difficult lies in the fact that as soon as the gradient of temperature changes at position $\mathbf{X}$ a new thermal model must be solved in the associated representative volume element $\Omega_{rve}(\mathbf{X})$. In the linear case, only a 3D thermal model for each space coordinate must be solved in each representative volume element $\Omega_{rve}(\mathbf{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 $\Omega_{rve}(\mathbf{X})$ could be solved for any macroscopic gradient $\mathbf{G}$. Thus, as soon as the macroscopic gradient is known, the macroscopic heat flux $\mathbf{Q}$ could be obtained without the necessity of solving a new thermal model. The main idea lies in introducing the components of $\mathbf{G}$ (denoted by $\mathbf{G}_1, \mathbf{G}_2, \mathbf{G}_3$) as extra-coordinates in the thermal model.

The main issue is how to transfer $\mathbf{G}$ from the boundary conditions to the partial differential equation itself. For this purpose, a new temperature field $\theta(\mathbf{x})$ is defined as:

$$ \theta(\mathbf{x}) = T(\mathbf{x}) - \mathbf{G} \cdot \mathbf{x} \qquad (48) $$

The thermal model (47) becomes:

$$ \left\{ \begin{array}{l} \nabla \cdot (\mathbf{k}(\theta, \mathbf{G}) \cdot \nabla (\theta(\mathbf{x}) + \mathbf{G} \cdot \mathbf{x})) = 0 \\ \theta(\mathbf{x} \in \partial\Omega_{rve}) = 0 \end{array} \right. \qquad (49) $$

whose solution is searched in a separated form

$$ \theta(x, y, z, \mathbf{G}_1, \mathbf{G}_2, \mathbf{G}_3) \\ \approx \sum_{i=1}^{N} X_i(x) \cdot Y_i(y) \cdot Z_i(z) \cdot \Gamma_i^1(\mathbf{G}_1) \cdot \Gamma_i^2(\mathbf{G}_2) \cdot \Gamma_i^3(\mathbf{G}_3) \qquad (50) $$

by applying the proper generalized decomposition as described in Sect. 2.

In mechanics, (for example non linear elasticity), this procedure suffices for determining the macroscopic stress $\Sigma$ from the macroscopic deformation tensor $\mathbf{E}$. This approach is sufficient because prescribing a displacement field $\mathbf{u} = \mathbf{E} \cdot \mathbf{x}$ on the boundary of $\Omega_{rve}$ ensures that the resulting microscopic displacement field verifies $\langle \epsilon \rangle = \mathbf{E}$. Obviously, a displacement given by $\mathbf{E} \cdot \mathbf{x} + \mathbf{U}_0$ (with $\mathbf{U}_0$ a uniform translation) could be prescribed on the boundary without modifying the solution because $\langle \epsilon \rangle = \mathbf{E}$ is still satisfactory 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 field is defined as:

$$ \tilde{\mathbf{u}}(\mathbf{x}) = \mathbf{u}(\mathbf{x}) - \mathbf{E} \cdot \mathbf{x} - \mathbf{U}_0 \qquad (51) $$

The introduction of the resulting expression of $\mathbf{u}(\mathbf{x})$ into the non-linear elastic problem, would make the translation field $\mathbf{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 through the material conductivity that usually depends on the tempera- ture itself. The same situation can be found in the Fick's dif- fusion law relating the concentration gradient and the chemi- cal species flu es through the so-called diffusion coefficien that is usually assumed depending on the concentration it- self. In all these cases the procedure just described for com- puting the macroscopic flu from the macroscopic gradi- ent (in our case **Q** from **G**) fails: Prescribing the tempera- ture **G** · **x** + T₀ on the boundary of Ω<sub>rve</sub> ensures condition ⟨∇T⟩ = **G**. However, the validity of the computed **Q** can- not be guaranteed because the computed temperature fiel depends on T₀, due to the dependence of the thermal con- ductivity on the temperature itself.

This problem of transferring the macroscopic tempera-
ture 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 equa-
tion 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:

$$
\theta(x, y, z, \mathbf{G}_1, \mathbf{G}_2, \mathbf{G}_3, T_0) \\
\approx \sum_{i=1}^{i=N} X_i(x) \cdot Y_i(y) \cdot Z_i(z) \cdot \Gamma_i^1(\mathbf{G}_1) \cdot \Gamma_i^2(\mathbf{G}_2) \cdot \Gamma_i^3(\mathbf{G}_3) \\
\cdot \Upsilon_i(T_0) \tag{52}
$$

4.2 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 mi-
croscopic one: (i) the gradient of the macroscopic tempera-
ture and (ii) the macroscopic temperature itself. The result-
ing microscopic problem is non-linear (the thermal proper-
ties depend on the microscopic temperature) and as just de-
scribed, one equation links the microscopic temperature to
the macroscopic one.

Even if the procedure is conceptually simple, its numeri-
cal implementation deserves certain difficultie since the so-
lution is nonlinear with respect to the gradient of the macro-

scopic temperature. Thus, many authors proposed a simpli-
fie modeling in which the thermal properties at the micro-
scopic level (that depend on the existing local temperature)
are frozen to values corresponding to the macroscopic tem-
perature. 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 mi-
croscopic temperature. This approach has been rigorously
justifie from the asymptotic expansion method for materi-
als with periodic microstructure in [42]. In a few words, in
this paper, the authors have shown that for a given 0-th or-
der temperature—i.e. macroscopic temperature—the effec-
tive conductivity of the non-linear temperature dependent
problem is equal to that of the linear problem with ther-
mal properties fi ed at the macroscopic temperature. There-
fore this paper provides a justificatio of an approach which
was widely applied in the context of computational homog-
enization of non-linear thermal models, see for example
[1, 28, 39].

Now, with the thermal properties frozen at the micro-
scopic scale, the non-linearity of the microscopic model dis-
appears. Thus, the microscopic solution only needs the pre-
scription of the macroscopic temperature gradient **G**(**X**) as
soon as the thermal conductivity is assumed known every-
where 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 mate-
rial thermal conductivity related to any macroscopic tem-
perature Θ(**X**). The solution of these microscopic problems
(T¹(**x**), T²(**x**) and T³(**x**) in Sect. 3.1.) for any macroscopic
temperature Θ suggests the separated representation:

$$
T^j(x, y, z, \Theta) \approx \sum_{i=1}^{i=N} X_i(x) \cdot Y_i(y) \cdot Z_i(z) \cdot \theta_i(\Theta) \\
j = 1, 2, 3 \tag{53}
$$

These solutions allow computing the homogenized con-
ductivity tensor for any value of the macroscopic tempera-
ture Θ. Through the dependence of the conductivity tensor
to the macroscopic temperature, this homogenization prob-
lem displays similarities with homogenization in nonlinear
elasticity, for which non-concurrent multiscale approaches
have been recently proposed [40, 41, 43]. In these refer-
ences, 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 pro-
vides automatically this interpolation and therefore enables
to compute in an elegant and efficien way the homogenized
behavior of non linear models.
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Concurrent homogenization of non-linear history dependent thermal models</figcaption></figure>

**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 $\nabla T^1$, $\nabla T^2$ and $\nabla T^3$ involved in the definition of K(X) with respect to the macroscopic temperature $\Theta$. 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.

## 4.3 History Dependent Non-linear Behaviors

In the previous section we assumed that the thermal model coefficient 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 extra-coordinates) 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 flux 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 \in \Omega_{rve}(X)$ depends on the thermal history at that position, i.e. $k = k(T(\tau))$, $0 \le \tau \le t$.

Now, as in the previous case, we should solve the thermal model define in $\Omega_{rve}(X)$, but in this case the thermal model must be integrated in the whole time interval $0 \le t \le t_{max}$:

$$ \frac{\partial T}{\partial t} - \nabla \cdot (k(T(t)) \nabla T) = 0 \quad (54) $$

with $T(x \in \partial \Omega_{rve}, t) = G(t) \cdot x + T_0(t)$.

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 finite element interpolation:

$$ G^i(t) = \sum_{j=1}^{j=M} N_j(t) \cdot G_j^i \quad (55) $$

where $G_j^i = G^i(t = t_j)$, $t_j = (j - 1) \times \Delta t$, and $N_j(t)$ the standard linear finite element shape functions.

We proceed in a similar way for approximating $T_0(t)$

$$ T_0(t) = \sum_{j=1}^{j=M} N_j(t) \cdot (T_0)_j \quad (56) $$

Now, the discrete values $G_j^i$, $i = 1, 2, 3$ and $j = 1, \dots, M$; and $(T_0)_j$, $j = 1, \dots, 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 \times 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:

$$ 
\begin{aligned}
& T(x, y, z, t, G_1^1, \dots, G_M^1, \dots, G_1^3, \dots, \\
& \qquad G_M^3, (T_0)_1, \dots, (T_0)_M) \\
& \approx \sum_{k=1}^{k=N} X_i(x) \cdot Y_i(y) \cdot Z_i(z) \cdot \Psi_i(t) \cdot \left( \prod_{j=1}^{j=3} \prod_{k=1}^{k=M} \Gamma_i^{j,k}(G_k^j) \right) \\
& \qquad \cdot \left( \prod_{k=1}^{k=M} \gamma_i^k((T_0)_k) \right)
\end{aligned} 
\quad (57) $$

*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 $\theta(\mathbf{x}, t) = T(\mathbf{x}, t) - G(t) \cdot \mathbf{x} - T_0(t)$.

## 4.4 History Dependent Non-linear Behaviors: A Simplified 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 macroscopic 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 \times 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.:

$$ T^j(x, y, z, t, \Theta_1, \dots, \Theta_M) $$

$$ \approx \sum_{i=1}^{i=N} X_i(x) \cdot Y_i(y) \cdot Z_i(z) \cdot \Psi_i(t) \cdot \left( \prod_{k=1}^{k=M} \theta_i^k(\Theta_k) \right) \quad (58) $$

that reduces the dimensionality of the problem with respect to problem (57) by removing $3 \times M$ coordinates.

# 5 Numerical Examples

## 5.1 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. 3.2), because it is used in Sect. 4 and in the numerical examples related to the procedures introduced in this section.

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 $\Omega_{vrte} = (0, 1) \times (0, 1)$ composed of $5 \times 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 $\Omega_{vrte}$ for two different boundary conditions: $T(\mathbf{x} \in \partial\Omega_{vrte}) = x$ in the first case and $T(\mathbf{x} \in \partial\Omega_{vrte}) = y$ in the second case. The resulting solutions are denoted $T^1(\mathbf{x})$ and $T^2(\mathbf{x})$ respectively. As described in Sect. 3.2 these solutions suffice for defining 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:

$$ T^j(x, y, k_{1,1}, k_{1,2}, \dots, k_{5,5}) \\
\approx \sum_{l=1}^{l=N} X_l(x) \cdot Y_l(y) \cdot \prod_{p=1}^{p=5} \prod_{q=1}^{q=5} \Xi_{pq}^l(k_{p,q}) \quad (59) $$

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 \times 5 = 27$) suffice to give access to the temperature field for any possible microstructure defined on this grid. Moreover, thanks to the PGD, the solution of these multidimensional models does not introduce major difficulties
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Microstructure composed of 5 × 5 cells. Each color denotes a value of the thermal conductivity assumed isotropic and constant in each cell: blue, green and red colors denote $k = 1, k = 1.5$ and $k = 2$ respectively</figcaption></figure>

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 defined by taking random values of the components of the conductivity tensor $k_{ij}$ that we denote by $k_{ik}^g$ (see Fig. 4). Then, the solution is obtained by particularizing the general PGD solution

$$T^1(x, y) \approx \sum_{l=1}^{l=N} X_l(x) \cdot Y_l(y) \cdot \prod_{p=1}^{p=5} \prod_{q=1}^{q=5} \Xi_{pq}^l (k_{p,q}^g) \quad (60)$$

$$T^1(x \in \partial \Omega_{rve}) = x$$

Temperatures $T^1(x, y)$ and $T^2(x, y)$ obtained using the PGD are compared with finite element solutions defined 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 means 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 finite 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 defined in a space of dimension 27 cannot be computed by using a standard mesh based discretization technique. It was also verified 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)$, $\Xi_{11}^i(k_{1,1})$ and $\Xi_{55}^i(k_{5,5})$ involved in the decomposition (59) for $N = 83$.

## 5.2 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 $\Omega_{rve} = (0, 1) \times (0, 1)$ is composed of two materials of different conductivities $k_1$ and $k_2$. The zone $\omega = (x_l, x_r) \times (y_d, y_t) \subset \Omega$ involves a material with conductivity $k_2$ whereas the material occupying the complementary region $\Omega - \omega$ involves a material of conductivity $k_1$ (see Fig. 7).

The region of conductivity $k_2$ can be expressed by the following separate form:

$$\omega = \chi^x(x) \times \chi^y(y) \quad (61)$$

<figure><img src="image_7.png" /><figcaption>**Fig. 5** $T(x) - x$ computed by applying both the PGD (left), and the finite element method (right) for the microstructure depicted in Fig. 1 and for the boundary condition $T(x \in \partial \Omega_{rve}) = x$</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Functions $X_i(x)$ (top-left), $Y_i(y)$ (top-right), $\Xi_{11}(k_{1,1})$ (down-left) and $\Xi_{55}(k_{5,5})$ (down-right)</figcaption></figure>

where both characteristic functions are defined by:

$$ \chi^x(x) = \begin{cases} 1 & \text{if } x_l < x < x_r \\ 0 & \text{otherwise} \end{cases} \qquad (62) $$

and

$$ \chi^y(y) = \begin{cases} 1 & \text{if } y_d < x < y_t \\ 0 & \text{otherwise} \end{cases} \qquad (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 justified in Sect. 4.1, their values are obtained considering that the microscopic temperature equals the macroscopic temperature $\Theta$ i.e. $k_1 = 10 + 9 \times \Theta$ and $k_2 = 90 + \Theta$, with $\Theta \in [0, 20]$. One can notice that this choice for the conductivities implies that when $\Theta = 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 \in \partial\Omega_{rve}) = x$ and $T(x \in \partial\Omega_{rve}) = y$, leading to the homogenized conductivity tensor. If one proceeds in this manner, as soon as the macroscopic temperature $\Theta$ evolves in time, the homogenized conductivity tensor should be recomputed for the thermal properties related to the new macroscopic temperature.

<figure><img src="image_9.png" /><figcaption>**Fig. 7** Two phase microstructure</figcaption></figure>

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, \Theta)$ and $T^2(x, y, \Theta)$ 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, \Theta)$ and $T^2(x, y, \Theta)$ are computed. Then,
the resulting temperature field for a given macroscopic tem-

perature $\Theta = \Theta^g$ (obtained by particularizing the general
solution, $T^1(x, y, \Theta = \Theta^g)$) is compared with the one com-
puted by applying the finite element method (with the ther-
mal properties defined at temperature $\Theta^g$). These results are
depicted in Fig. 8 (left PGD, right FE) proving the ability

<figure><img src="image_10.png" /><figcaption>**Fig. 8** $T(x) - x$ computed by applying the PGD (**left**) and the FEM (**right**) for different macroscopic temperatures $\Theta = 0$ (**top**), $\Theta = 10$ (**middle**) and $\Theta = 20$ (**down**), when the thermal model in $\Omega_{rve}$ was solved by assuming $T(x \in \partial\Omega_{rve}) = x$</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 9** Functions $X_i(x)$ (top-left), $Y_i(y)$ (top-right) and $\theta_i(\Theta)$ (down)</figcaption></figure>

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 significant functions $X_i(x)$, $Y_i(y)$ and $\theta_i(\Theta)$ ($i = 1, \dots, 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 $\Theta$.

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 extra-coordinate (the macroscopic temperature at the location of the RVE $\Omega_{rve}$), the necessity of solving a microscopic model for each value of the macroscopic temperature can be avoided.

## 5.3 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 fixed (equal to $\Theta(\tau)$, $0 \le \tau \le t$), following the simplified 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 $\Omega_{rve}$ could be computed and frozen.

Following the description given in Sect. 4.4, the approximation of $\Theta(t)$ in the whole time interval was interpolated from a polynomial of degree 4 that can be expressed from 5 nodal values of $\Theta$, denoted by $\Theta_i$, $i = 1, \dots, 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, \Theta_1, \dots, \Theta_M)$ and $T^2(x, y, \Theta_1, \dots, \Theta_M)$ according to (58). In this case, two problems define in a space involving x, y, t and the five temperatures $\Theta_i$, $i = 1, \dots, 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 fixed arbitrarily to one. Finally the five extra-coordinates related to the intermediate macroscopic temperatures $\Theta_i$, $i = 1, \dots, 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, \Theta_1, \dots, \Theta_M)$ and $T^2(x, y, \Theta_1, \dots, \Theta_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
<figure><img src="image_12.png" /><figcaption>**Fig. 10** $T(x) - x$ computed by applying the PGD for $\Theta_1 = 0, \Theta_2 = 5, \Theta_3 = 10, \Theta_4 = 5$ and $\Theta_5 = 0$ at times $t = 0.25$ (top-left), $t = 0.5$ (top-right), $t = 0.75$ (down-left) and $t = 1$ (down-right) when the thermal model in $\Omega_{rve}$ was solved by assuming $T(x \in \partial\Omega_{rve}) = x$</figcaption></figure>

decreases again linearly until reaching a null temperature at the final time. This problem, as soon as the thermal history is fixed, was also solved by means of the finit element method. To validate the proposed approach we depict 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, \Theta_1, \dots, \Theta_M)$ for the macroscopic thermal history just defined i.e. $T^1(x, y, \Theta_1 = 0, \Theta_2 = 5, \Theta_3 = 10, \Theta_4 = 5, \Theta_5 = 0)$. The corresponding solutions, computed by applying the finit element solution and enforcing $T(x \in \partial\Omega_{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 dependence. 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 sensitivity of the results to the regularity in that thermal dependence.

## 6 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 define 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 defined in a space involving $2 + 10 \times 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
<figure><img src="image_13.png" /><figcaption>**Fig. 11** $T(\mathbf{x}) - x$ computed by applying the FEM for $\Theta_1 = 0, \Theta_2 = 5, \Theta_3 = 10, \Theta_4 = 5$ and $\Theta_5 = 0$ at times $t = 0.25$ (top-left), $t = 0.5$ (top-right), $t = 0.75$ (down-left) and $t = 1$ (down-right) when the thermal model in $\Omega_{rve}$ was solved by assuming $T(\mathbf{x} \in \partial\Omega_{rve}) = x$</figcaption></figure>

temperature or macroscopic gradient), even in the transient case.

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(\mathbf{x})$ instead of $X_i(\mathbf{x}) \cdot Y_i(\mathbf{y}) \cdot Z_i(\mathbf{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 field of computational homogenization, rather than addressing models of practical interest.

## References

1. Alzina A, Toussaint E, Beakou A, Skoczen B (2006) Multiscale modelling of thermal conductivity in composite materials for cryogenic structures. Compos Struct 74(2):175-185

2. Ammar A, Mokdad B, Chinesta F, Keunings R (2006) A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. J Non-Newton Fluid Mech 139:153-176

3. Ammar A, Mokdad B, Chinesta F, Keunings R (2007) A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of com-
plex fluids Part II: transient simulation using space-time separated representations. J Non-Newton Fluid Mech 144:98-121

4. Ammar A, Normandin M, Daim F, Gonzalez D, Cueto E, Chinesta F (2010) Non-incremental strategies based on separated representations: Applications in computational rheology. Commun Math Sci 8(3):671-695

5. Arbenz P, van Lenthe GH, Mennel U, Müller R, Sala M (2008) A scalable multi-level preconditioner for matrix-free µ-finit element analysis of human bone structures. Int J Numer Methods Eng 73(7):927-947

6. Bóhm HJ (2009) A short introduction to basic aspects of continuum micromechanics. ISLB Report 208, TU Wien, Vienna, http://www.ilsbu.tuwien.ac.at/links/downloads/ilsbrep206.pdf

7. Chinesta F, Chaidron G (2001) On the steady solution of linear advection problems in steady recirculating flows. J Non-Newton Fluid Mech 98:65-80

8. Chaidron G, Chinesta F (2002) On the steady solution of nonlinear advection problems in steady recirculating flows. Comput Methods Appl Mech Eng 191:1159-1172

9. Chinesta F, Torres R, Ramón A, Rodrigo MC, Rodrigo M (2002) Homogenized thermal model in foods containing particulates. Int J Thermal Sci 41:1141-1150

10. Chinesta F, Chaidron G, Poitou A (2003) On the solution of the Fokker-Planck equations in steady recirculating flows involving short fiber suspensions. J Non-Newton Fluid Mech 113(2-3):97-125

11. Chinesta F, Ammar A, Lemarchand F, Beauchene P, Boust F (2008) Alleviating mesh constraints: model reduction, parallel time integration and high resolution homogenization. Comput Methods Appl Mech Eng 197(5):400-413

12. Chinesta F, Ammar A, Joyot P (2008) The nanometric and micro-metric scales of the structure and mechanics of materials revisited: an introduction to the challenges of fully deterministic numerical descriptions. Int J Multiscale Comput Eng 6(3):191-213

13. Chinesta F, Ammar A, Cueto E (2010) Recent advances in the use of the Proper Generalized Decomposition for solving multidimensional models. Arch Comput Methods Eng. doi:10.1007/S11831-010-9049-Y

14. Chinesta F, Ammar A, Cueto E (2010) Proper Generalized Decomposition of multiscale models. Int J Numer Methods Eng 83(8-9):1114-1132

15. de Borst R (2008) Challenges in computational materials science: multiple scales, multi-physics and evolving discontinuities. Comput Mater Sci 43(1):1-15

16. Feyel F (1999) Multiscale FE2 elastoviscoplastic analysis of composite structures. Comput Mater Sci 16(1-4):344-354

17. Fish J, Yuan Z (2005) Multiscale enrichment based on partition of unity. Int J Numer Methods Eng 62(10):1341-1359

18. Fish J (2006) Bridging the scales in nano engineering and science. J Nanopart Res 8(5):577-594

19. Geers MGD, Kouznetsova VG, Brekelmans WAM (2009) Multiscale computational homogenization: trends and challenges. J Comput Appl Math. doi:10.1016/j.cam.2009.08.077

20. Hughes TJR (1995) Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Comput Methods Appl Mech Eng 127(1-4):387-401

21. Jiang M, Jasiuk I, Ostoja-Starzewski M (2002) Apparent thermal conductivity of periodic two-dimensional composites. Comput Mater Sci 25(3):329-338

22. Kanit T, Forest S, Galliet I, Mounoury V, Jeulin D (2003) Determination of the size of the representative volume element for random composites: statistical and numerical approach. Int J Solids Struct 40(13-14):3647-3679

23. Kanit T, N'Guyen F, Forest S, Jeulin D, Reed M, Singleton S (2006) Apparent and effective physical properties of heterogeneous materials: representativity of samples of two materials

from food industry. Comput Methods Appl Mech Eng 195(33-36):3960-3982

24. Kanouté P, Boso DP, Chaboche JL, Schrefle BA (2009) Multiscale methods for composites: a review. Arch Comput Methods Eng 16(1):31-75

25. Kolda TG, Bader BW (2007) Tensor decompositions and applications. Technical Report SAND2007-6702, SANDIA National Laboratories

26. Ladeveze P (1999) Nonlinear computational structural mechanics. Springer, New York

27. Ladeveze P, Passieux JCh, Neron D (2009) The LATIN multiscale computational method and the proper orthogonal decomposition. Comput Methods Appl Mech Eng. doi:10.1016/j.cma.2009.06.023

28. Laschet G (2002) Homogenization of the thermal properties of transpiration cooled multi-layer plates. Comput Methods Appl Mech Eng 191(41-42):4535-4554

29. McVeigh C, Liu WK (2008) Linking microstructure and properties through a predictive multiresolution continuum. Comput Methods Appl Mech Eng 197(41-42):3268-3290

30. Michel JC, Suquet P (2004) Computational analysis of nonlinear composite structures using the nonuniform transformation fiel analysis. Comput Methods Appl Mech Eng 193(48-51):5477-5502

31. Mobasher-Amini A, Dureisseix D, Cartraud P (2009) Multi-scale domain decomposition method for large-scale structural analysis with a zooming technique: application to plate assembly. Int J Numer Methods Eng 79(4):417-443

32. Mokdad B, Pruliere E, Ammar A, Chinesta F (2007) On the simulation of kinetic theory models of complex fluids using the Fokker-Planck approach. Appl Rheol 17(2):1-14

33. Nouy A, Ladeveze P (2004) Multiscale computational strategy with time and space homogenization: a radial-type approximation technique for solving microproblems. Int J Multiscale Comput Eng 170(2):557-574

34. Nouy A (2009) Recent developments in spectral stochastic methods for the solution of stochastic partial differential equations. Arch Comput Methods Eng 16(3):251-285

35. Ostoja-Starzewski M (2006) Material spatial randomness: from statistical to representative volume element. Probab Eng Mech 21(2):112-132

36. Özdemir I, Brekelmans WAM, Geers MGD (2008) Computational homogenization for heat conduction in heterogeneous solids. Int J Numer Methods Eng 73(2):185-204

37. Sab K (1992) On the homogenization and the simulation of random materials. Eur J Mech A Solids 11(5):585-607

38. Strouboulis T, Zhang L, Babuška I (2004) p-version of the generalized FEM using mesh-based handbooks with applications to multiscale problems. Int J Numer Methods Eng 60(10):1639-1672

39. Telega JJ, Tokarzewski S, Galka A (2000) Effective conductivity of nonlinear two-phase media: homogenization and two-point Pade approximants. Acta Appl Math 61(1):295-315

40. Temizer I, Wriggers P (2007) An adaptive method for homogenization in orthotropic nonlinear elasticity. Comput Methods Appl Mech Eng 196(35-36):3409-3423

41. Temizer I, Zohdi T (2007) A numerical method for homogenization in non-linear elasticity. Comput Mech 40(2):281-298

42. Tokarzewski S, Andrianov I (2001) Effective coefficients for real non-linear and f'c titious linear temperature-dependent periodic composites. Int J Non-Linear Mech 36(1):187-195

43. Yvonnet J, Gonzalez D, He Q-C (2009) Numerically explicit potentials for the homogenization of nonlinear elastic heterogeneous materials. Comput Methods Appl Mech Eng 198(33-36):2723-2737

44. Zohdi T, Wriggers P (2005) An introduction to computational micromechanics. Springer, Berlin