A PGD-based homogenization technique for the resolution of nonlinear multiscale problems

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.

In engineering design, the numerical simulation of very large multiscale models is becoming increasingly important because of the need to describe realistic scenarios and derive tools to facilitate the virtual design of new structures. In the context of such huge calculations, multiscale approaches based on domain decomposition methods [1] or on model superposition [2,3] are widely used. Another popular family of strategies is based on computational homogenization, i.e. the use of PGD to achieve computational homogenization was investigated for thermal problems with materials whose microstructure can be described by a microscopic representative elementary volume (an assumption which is not required in our work). The issue of the multiscale analysis of structures using the PGD technique has also been addressed in [46] where the PGD method was used in a FE2 technique to solve the problems defined at the level of the RVE for any possible sets of geometry, material properties and boundary conditions. This work has been developed in the case of linear elasticity, which allowed to compute the homogenized behavior only once, in a preliminary phase that can be performed offline, whereas the aim of the present paper is to address nonlinear behaviors.
Once the algorithm has converged, the homogenized operator can be viewed as the homogenized tangent operator which, over the time interval and for each subdomain, links a small macrodisplacement perturbation to a small macroforce perturbation. In the case of viscoelastic materials, this operator characterizes the local macroscopic behavior of the structure.
The paper is organized as follows. Section 2 introduces the reference problem. The main features of the LATIN multiscale technique, including a brief description of the PGD method, are presented in Section 3. Section 4 focuses on the new algorithm and the determination of the PGD representation of the homogenized operator. Finally, Section 5 presents a numerical example which shows that the use of PGD in the macroproblem does not affect the convergence rate of the approach and leads to a reduction in the total computation time.

Description of the reference problem
Let us consider, under the assumption of small perturbations, the quasi-static and isothermal evolution of a structure defined in time-space domain I Â X, where I ¼ ½0; T and X 2 R D (D ¼ 2; 3 being the spatial dimension). Boundary @X is assumed to be divided into two non-overlapping parts @ 1 X and @ 2 X such that @ 1 X [ @ 2 X ¼ @X and @ 1 X \ @ 2 X ¼ ;. Surface forces F d are applied over boundary @ 2 X and a displacement u d is prescribed over @ 1 X. The structure is also subjected to a body force f d (see Fig. 1).
Let e be the strain associated with displacement u and let r be the Cauchy stress. The strain can be divided into an elastic strain e e and an inelastic strain e p so that e ¼ e e þ e p . The elastic strain e e is related to the stress through the state law: where K is Hooke's tensor. The inelastic strain is related to the stress through the evolution law: where B is a differential operator which can be time-dependent and nonlinear. In the numerical examples of Section 5, two behaviors will be investigated: a linear viscoelastic behavior _ e p ¼ K À1 =g r and a nonlinear time-dependent behavior _ e p ¼ K À1 =ðEgðtÞÞTrðrÞr. All these quantities are defined at each time-space point ðt; xÞ 2 I Â X and assigned initial values at time t ¼ 0. Thus, the state of the structure is defined by the set s ¼ ð_ e p ; rÞ.
Such a representation of a material's behavior is typical of a viscoelastic constitutive law expressed using what is called a ''normal formulation'' [28]. Other types of constitutive laws, such as elastic-viscoplastic behavior, which requires the introduction of internal variables, can be found in [47,48]. Dealing with such laws does not require any change in the strategy which will be developed in this paper.

The LATIN multiscale method
In this section, only the main characteristics of the LATIN multiscale method are reviewed. Additional details can be found in A. The interested reader can also refer to [37].

Multiscale description of the interface unknowns
Two different descriptions of the interface unknowns (displacements and forces) are introduced: a macroscale description (denoted with superscript M ) and a microcomplement (denoted with superscript m ). Over an interface U EE 0 between subdomains X E and X E 0 , one has: The choice of the spaces of the macrodisplacements and macroforces, denoted respectively W M EE 0 and F M EE 0 , is arbitrary, but has a strong influence on the scalability of the method. A major point of the strategy, which grants it its multiscale character, is that the set of the macroforces is required to verify the transmission conditions at the interfaces a priori at each iteration. This assumption is called the ''admissibility of the macroforces'' and is detailed in Appendix A.1.

The LATIN solver
The governing equations of the problem are divided into two groups: the linear equations and the local equations. First, let A d be the space of the solutions of the linear equations (defined globally in the spatial domain) whose elements s ¼ ðs E Þ X E &X verify: the admissibility condition (6) for every X E & X : s E 2 A dE the admissibility of the macroforces (see (A.4) and (A.5) for details) Then, let C be the space of the solutionsŝ ¼ ðs E Þ X E &X of the local (possibly nonlinear) equations whose elements satisfy: the evolution law (2) for every X E & X : 8ðt; xÞ 2 I Â X E , _ e pE ¼ Bðr E Þ the interface behavior (7) over every interface U EE 0 The solution of the problem s ref ¼ A d \ C is built by using the LATIN method. This method is an iterative scheme between two stages called the ''local stage'' and the ''linear stage'', repeated alternatively until convergence. At Iteration n þ 1, in the local stage detailed in Appendix A.2, one seeks a solutionŝ nþ1=2 of the local nonlinear equations in space C; in the linear stage detailed in Appendix A.3, one seeks a solution s nþ1 of the global linear equations in space A d . The introduction of search directions enables one to switch back and forth between space C and space A d . The convergence of the algorithm is checked using the criterion described in Appendix A.4. This two-search-direction algorithm is quite similar to that described in [49] for non-overlapping domain decomposition methods in linear elasticity.

The linear stage: a series of microproblems and a macroproblem
The linear stage detailed in Appendix A.3 leads to a series of microproblems defined in each time-space subdomain I Â X E along with a macroproblem defined over the entire time interval I and for the whole set of interfaces. These problems are defined in the following.

The series of microproblems
Introducing the weak form of search direction (A.8) and (A.10) into Eq. (A.12), one gets the following microproblem for each subdomain X E : This macroproblem boils down to a linear homogenized problem in space and time defined over the whole set of interfaces and the entire time interval.

Resolution of the linear stage
The linear stage consists in solving a series of microproblems (9) defined in each time-space subdomain I Â X E along with macroproblem (12) defined over the entire time interval I and for the whole set of interfaces, leading to Lagrange multiplier _ W M E . One can note that parallel resolution of these microproblems is possible, which reduces the computation cost.

Proper Generalized Decomposition
As described in Appendix A.3, many microproblems defined on the subdomain level must be solved during the linear stage. This stage can be very time-consuming and requires a significant amount of storage. This difficulty is addressed thanks to the introduction of the Proper Generalized Decomposition (PGD) technique. PGD was originally introduced into the LATIN method as part of the solver under the name ''radial approximation'' [28]. It has been shown in previous works that PGD leads to a significant reduction in computation time [36,37,44].
The main idea consists in representing a generic function of time and space q ¼ qðt; xÞ using separated variables: where k k ðtÞ are functions of time alone and K k ðxÞ are functions of space alone. The PGD representation is introduced into the LATIN scheme in order to describe the unknowns of microproblems (9). At the beginning of an iteration, if a reduced-order basis (ROB) is known, it is used to build a first reduced-order model to solve the new microproblem. If the quality of the solution is insufficient, the ROB is enriched by adding new PGD functions using a greedy power algorithm based on the minimization of a functional built from the verification of the search direction (A.8) and (A.10). Further details on the introduction of PGD into the microproblems and examples showing the ability of the method to reduce the computation cost of the microproblems can be found in [37,44,48]. The novelty of the present work lies in the reduction of the computation and storage cost of the homogenized operator by extending the use of the PGD technique to its construction. whose solution is now sought using a PGD representation of stress tensor r E in which not only time and space, but also the variable which defines the Lagrange multiplier are separated: The objective of this choice is to build the homogenized operator at each time t, at each point x and for each value of the Lagrange multiplier defined by c E . The triplets fA k ; B k ; C k g of PGD representation (23) are obtained by means of Power Algorithm 1, which is based on the rewriting of Problem (22) as the minimization of the error: with the introduction of the norms: M and m are symmetric, positive definite operators (which, in the following numerical example, will be chosen as M ¼ H À1 and m ¼ h À1 ) and C is the space of parameters c E . The calculation of norms (25) requires the evaluation of integrals over the space C of definition of parameters c E . The bounds of this space are a priori unknown and, at this stage in our development, are chosen to be large enough for r E to be known for a range of c E which is sufficient to satisfy Macroproblem (12). The optimization of this choice will be the subject of future developments, but it has little impact on the global computation cost.
Finally, the procedure enables one to build the homogenized operator L F E which maps The corresponding operation can be written formally as: where triplets ða k ; b k ; c k Þ are deduced easily from triplets ðA k ; B k ; C k Þ by using the fact that F M E is the macropart of r E n E at the boundary of X E (with n E being the outer unit normal to boundary @X E ). Let us note that, due to the linearity of (22) with respect to c E , functions c k are also linear in c E .

Discussion of the computation and storage costs
At enrichment step M of Algorithm 1, the construction of a new triplet fA M ; B M ; C M g given fA k ; B k ; C k g k¼1;...;MÀ1 requires the minimization of a problem which couples time, space and macroquantities. This problem is solved using a fixed-point algorithm with ' max iterations. Thus, step M involves ' max resolutions of D-dimensional spatial problems (D ¼ 2; 3) for function B M ; ' max resolutions of 1-dimensional ODEs for function A M , and ' max resolutions of 1-dimensional problems for function C M . The complexity of this procedure lies almost entirely in the ' max D-dimensional steady-state problems (compared to which the cost of the 1-dimensional problems is negligible). In our numerical tests, the initialization of A 0 ðtÞ and C 0 ðc E Þ did not affect the convergence of the algorithm significantly and, in practice, only a few iterations of the algorithm were necessary to converge. Therefore, the calculation of a new triplet was linked to the resolution of rd ' max spatial microproblems and the global computation cost of the homogenized operator was in the order of that of Mrd ' max spatial microproblems.
In the examples presented in the next section, following [37], the number of iterations was arbitrarily set to a small value (' max ¼ 2). Even if the convergence of the fixed-point algorithm is not reached, the quality of the decomposition is enhanced by the next triplets. It is known that the convergence can be slower than in the case of two variables decompositions but, in the numerical examples investigated in the paper, this issue does affect the global convergence rate of the strategy. Further investigations should be necessary for other material behaviors. The number M of triplets required to describe the homogenized operator was small. In the case of linear behavior, such as in viscoelasticity, the tangent search direction H ¼ @B=@rjr is constant. Therefore, fixed search directions H and h were used throughout the iterations and the homogenized operator needed to be calculated only once at the beginning of the iterative process. Usually, in that case, M ¼ 2 triplets are enough.
Conversely, in the case of nonlinear behavior, parameters H and h are recalculated throughout the iterations and the homogenized operator must be updated [48]. In practice, the search directions are recalculated only for the first iterations, then remain fixed when their variations become negligible. However, we did not use this technique in our work in order to simplify the discussion about the computation cost in the general case. Therefore, M was set to be equal to the number of iterations of the LATIN method which were necessary to reach a given error. This number can also be decreased significantly if one reuses the already known triplets as a reduced-order basis when performing a given iteration: then, the values of B k and C k are fixed and the recalculation of time function A k alone suffices to find a better representation of the operator without solving any D-dimensional problem. Then, one introduces a criterion in order to add a new triplet if, and only if, the given reduced-order basis turns out to be insufficient.
The global cost of the Mrd ' max spatial problems required to build the homogenized operator must be compared to the cost of the classical procedure described in Section 4.1. Let us recall that this procedure consists in solving rdðN þ 1Þ Eq. (15) using a time-space PGD-separated form of the unknowns. Such a resolution is very close to that described in Algorithm 1 and the same computation cost evaluation applies. This cost is in the order of that of M' max spatial problems for each of the rdðN þ 1Þ equations, which corresponds to the global cost of MrdðN þ 1Þ' max D-dimensional spatial problems. Therefore, the advantage of our procedure is obvious because its cost is independent of the number of time steps. In terms of storage cost, our procedure requires the storage of only the M triplets (whose cost is associated mainly with the spatial functions) compared to the storage of N þ 1 spatial operators (one for each of the times considered) with the classical procedure.
Relation (29) means that once the algorithm has converged L W E can be viewed as the local macrobehavior near the calculated solution, defined over I ¼ ½0; T and in each subdomain X E . This behavior establishes a relation between a small perturbation in terms of macrodisplacements and a small perturbation in terms of macroforces.
In the particular case of a linear problems, such as in viscoelasticity, the tangent search direction is independent of the solution being sought. Therefore, it is calculated only once at the beginning of the algorithm and then remains constant throughout the iterations. Thus, the previous interpretation is valid throughout the iterations.
If search direction H is not tangent to subspace C, the algorithm bears some similarities to a quasi-Newton scheme and no obvious physical interpretation can be given to the homogenized operator.

Numerical examples
Now, the strategy presented in the previous sections will be validated by introducing some numerical tests in order to show that the PGD representation of the homogenized operator guarantees the same order of convergence as the standard technique while leading to a significant reduction in the computation time when the operator needs to be updated.
Under the assumption of plane strain, we considered the evolution over I ¼ ½0; T (T ¼ 2 s) of a two-dimensional L-shaped structure made of 12 heterogenous cells (see Fig. 4). Each cell was a 0.5 m Â 0.5 m square with an inclusion of radius 0.15 m. The matrix and the inclusion were made of two different materials denoted respectively Material 1 (shown in light gray) and Material 2 (shown in dark gray). The structure was clamped at the bottom and subjected along its right side to a prescribed pressure F d ðtÞ ¼ F 0 sin 2pt=T x with F 0 ¼ 100 MPa. The domain was divided into 12 subdomains, each corresponding to a cell. In two dimensions, the spatial macrobasis consists of d ¼ 4 vectors. The number of interfaces to be taken into account for the construction of the homogenized operator was r ¼ 16. Each subdomain was discretized into linear triangular finite elements, leading to a total of 3675 DOFs (see Fig. 5(a)). The time interval I was discretized into N time steps. A Euler implicit integration scheme was used. We studied the influence of N on the cost of the resolution assuming that the convergence of the algorithm is achieved when the error indicator g defined in Appendix A.4 becomes less than 10 À5 . Since our example is relatively simple, it was possible to build a reference solution using a direct incremental approach and to plot the actual error g ref defined in Appendix A.4.

Linear viscoelastic behavior
First, we assumed for Material 1 and Material 2 linear viscoelastic behavior, such that: where K i denotes the Hooke's tensor of Material i. The Poisson's ratio was the same for both materials (m 1 ¼ m 2 ¼ 0:3) but the Young's moduli and the viscosity parameters were different: E 1 ¼ 50 GPa, E 2 ¼ 200 GPa and g 1 ¼ 10 s, g 2 ¼ 1000 s. Fig. 5(b)-(d) shows the stress components at time t ¼ 3T=4. In order to validate these results, a first comparison with a classical direct incremental scheme is given by Fig. 6, which shows the horizontal displacement of interface point P (see Fig. 4). One can observe that the solution obtained with the direct incremental method (dotted line) and the solution obtained using our new strategy (continuous line) match perfectly.
Figs. 7 and 8 enable one to study the convergence of the procedure. The continuous line represents the convergence curve of our approach (LATIN method + PGD representation of the homogenized operator). The convergence of the standard approach (LATIN method + standard homogenized operator [37]) is represented by circles. The two convergence curves    match perfectly, which shows that the introduction of PGD into the representation of the homogenized operator did not affect the convergence of the approach. If not enough triplets had been used in the PGD description, a bad representation of the homogenized operator would have been obtained, leading to a poor convergence rate. Now, in order to compare the two approaches in terms of their performance, we will consider the computation cost of generating the homogenized operator. Let us observe that, due to the linearity of the problem in this example, the search direction operators ðH; hÞ are independent of time and remain constant throughout the iterations, which enabled the homogenized operator to be calculated only once at the beginning of the analysis. The CPU costs were compared in terms of the number of linear spatial microproblems (defined on the subdomain scale) which had to be solved to build the homogenized operator. In our example, the 12 subdomains were identical and had the same number of DOFs. Thus, all the spatial microproblems had the same computation cost, denoted t:u. (time unit).
The following study focuses on the influence of the number N of time steps in the discretization. For the sake of generality, we assumed the temporal discretization to be non-uniform, thus requiring the resolution of a linear system at each time step. (In the particular case of a uniform discretization, the linear system would have been factorized only once at the beginning.) The results are summarized in Table 1, which shows that the number of linear systems to be solved is independent of the number of time steps when using our technique, but increases proportionately to N þ 1 when using the standard method. The gain achieved with the new method is significant, especially when the number N of time steps increases.

Nonlinear time-dependent behavior
Now, in order to study the behavior of the new procedure in the case of a nonlinear problem, let us consider that the behavior of Material 1 and Material 2 is described as: where K i and E i are the same as before, but the viscosity parameters are time-dependent: g 1 ¼ g 2 ¼ g 0 sin t=2p with g 0 ¼ 100 s.
The search direction operators ðH; hÞ could have been considered to be independent of time and constant throughout the iterations (using, for example, time-averaged pseudo-tangent operators), leading to a unique homogenized operator, but this would have potentially resulted in a poor convergence rate due to the nonlinearity of the problem. In order to guarantee rapid convergence, the search direction operators needed to be updated to recalculate the tangent search direction. As mentioned earlier, practical criteria could be set up in order to avoid updating the operators when their variations are small, but this feature was not introduced in this study so the discussion on the computation cost would be simpler. Therefore, we assumed that the time-dependent homogenized operator had to be recalculated at each iteration.
Two variants of the strategy were considered: the first variant, hereafter labeled '(constant)', consisted in calculating the homogenized operator only once at the beginning; the second variant, hereafter labeled '(tangent)', consisted in updating the operator at each iteration.
The results are shown in Figs. 9 and 10. Fig. 9 shows the comparison of the convergence rates of these two variants without and with the new PGD-based construction of the homogenized operator. As expected, the use of an updated tangent search direction improved the convergence rate of the method and, like in the linear behavior case, the PGD representation did not affect the convergence of the strategy. However, the more interesting point concerns the performance of the approach in terms of the computation cost of the homogenized operator. Using the same notations as in Table 1, Tables 2  and 3 demonstrate the gain obtained using our strategy with either a constant or a tangent (updated at each iteration) search direction.
Concerning Table 2, one can observe that the computation cost of the standard technique is the same as in the linear case. Indeed, since the homogenized operator is calculated only once at the beginning of the algorithm, its computation time depends only on the number of time steps. When the PGD representation is introduced, the computation cost is higher than in the linear case because a few more triplets are necessary to describe the operator due to this complex behavior. However, the gain is still significant and goes up to about 7 in the case of 40 time steps.
Conversely, in Table 3, one can see that in the case of a tangent search direction the operator is recalculated at each iteration and, therefore, the number of iterations needed to reach convergence affects the computation cost. The gain can go up to about 18 in the case of 40 time steps. Fig. 11 summarizes these results by showing the evolution of the gain afforded by the PGD-based homogenization technique compared to the standard approach as a function of the number of time steps. The continuous line corresponds to a fixed homogenized operator and the dashed line corresponds to a variable tangent search direction.     Table 3 Comparison of the computation costs in terms of the number of spatial microproblems (t.u.) to be solved in order to build the homogenized operator (using a tangent operator).  Let us recall that the objective of this numerical example was to evaluate the PGD-based construction of the homogenized operator. That is the reason why the performance comparison has been presented only in terms of the number of global resolutions (spatial microproblems) which are necessary to build the operator. In our example, the number of subdomains was small and the construction of the homogenized operator represented only 30% of the total computation time. In more complex situations, such as in the problems addressed in [37] in which the number of subdomains is large, the computation cost of the homogenized operator can account for most of the total computation time. In that case, a gain of 18 in the calculation of the homogenized operator would be all the more significant in terms of the overall cost of the analysis.

Conclusions
This paper is based on the LATIN multiscale computational strategy, which enables one to use the PGD model order reduction technique to reduce the computation cost. The multiscale nature of the strategy is based on the separation of the interface unknowns according to a macroscale and a microscale. The macroscale is used to represent the homogenized behavior of the structure and to accelerate the convergence of the iterative algorithm. In previous works based on this strategy, the PGD technique was used to reduce the computation cost of the problems defined on the microscale and to build a reducedorder model of the microscopic behavior.
The novelty of this paper is that it extends the use of the PGD representation to the construction of the homogenized operator. The numerical example shows that the convergence rate of the approach is unaffected by this new PGD representation, but a significant reduction in computation cost can be achieved, especially in the case of nonlinear behavior. The proposed strategy enables a reduced-order model of the homogenized behavior of structures made of nonlinear materials to be calculated offline and at reasonable cost. This homogenized behavior, which characterizes the local macroscopic behavior of the structure and is defined over the whole time interval and in each subdomain, establishes a relation between a small perturbation in terms of macrodisplacements and a small perturbation in terms of macroforces. Once this homogenized behavior has been calculated and stored, it can be used online to calculate the response of the structure rapidly for new loading cases.
The choice of the spaces of the macrodisplacements and macroforces, denoted respectively W M EE 0 and F M EE 0 , is arbitrary, but has a strong influence on the scalability of the method. In practice, in the spatial domain, the macropart is chosen to be the linear part of the forces and displacements, which guarantees the scalability of the domain decomposition method in space, but other possible choices are discussed in [51]. Scalability in time is more difficult to achieve and requires an adaptive generation of the temporal macrobasis. A detailed analysis of this question was presented in [50]. Here, the separation of the macroscale from the microscale is carried out only in space, which means that the temporal microscale and macroscale are identical. The extension to the multiscale temporal case is straightforward.
Once spaces W M EE 0 and F M EE 0 have been defined, the macropart of the displacement field can be defined as: and the macropart of the force field can be defined as: The spaces of the interface quantities can be extended to the other interfaces N E of subdomain X E as: The extension to the whole set of the subdomains of X leads to W M and F M . The micropart of the displacements and forces can be derived from (8).
Over the interfaces, the macroforces are assumed to verify the transmission conditions:

ðA:5Þ
This assumption is called the ''admissibility of the macroquantities''. The corresponding subspace is defined as F M ad . Now let us introduce the space W M ad of the displacements which are continuous at the interfaces and equal to the macropart of prescribed displacement u d over @ 1 X. The spaces whose elements have their macroparts in W M ad and F M ad are denoted W ad and F ad .

A.2. The local stage
At iteration n + 1, given s n 2 A d , the local stage leads to an approximation of the solutionŝ nþ1=2 in space C such that ðŝ nþ1=2 À s n Þ belongs to search direction E þ . Skipping the subscripts for the sake of clarity, in each subdomain X E & X; E þ is defined by: Hðr E À r E Þ þ ð_ e pE À _ e pE Þ ¼ 0 where H and h are symmetric, positive definite operators which have an influence on the convergence rate of the algorithm, but do not affect the solution. The solutionŝ of the local stage is the ''intersection'' of spaces E þ and C. This leads to a set of equations which are defined locally in space.

A.3. The linear stage
The linear stage consists, given the solutionŝ nþ1=2 2 C of the local stage, in finding a solution s n 2 A d such that ðs nþ1 Àŝ nþ1=2 Þ belongs to search direction E À . Skipping the subscripts, E À is defined by: which is written in a weak sense in order to take into account the admissibility of the macroforces (A.4) and (A.5).
In [28], it was shown that an optimal choice of H in terms of the convergence rate consists in defining the manifold E À as the vector space which is tangent to C inŝ. This requires the calculation of H ¼ @B=@rjr at each iteration (or, at least, every few iterations) as in a classical Newton algorithm. h can be viewed as a ''viscosity'' effect at the interface. The choice of this parameter is discussed in [37].
The admissibility of the macroforces F Ã 2 F ad in (A.9), which expresses that macropart F MÃ belongs to F M ad , is enforced by introducing a Lagrange multiplier _ W