M A Benaimeche

| J Yvonnet

| B Bary

-C He

A k-means clustering machine learning-based multiscale method for anelastic heterogeneous structures with internal variables

Keywords: Multiscale, FE2, k-means clustering, Nonlinear, Homogenization, Machine learning

A new machine-learning based multiscale method, called k-means FE 2 , is introduced to solve general nonlinear multiscale problems with internal variables and loading history-dependent behaviors, without use of surrogate models. The macro scale problem is reduced by constructing clusters of Gauss points in a structure which are estimated to be in the same mechanical state. A k-means clustering -machine learning technique is employed to select the Gauss points based on their strain state and sets of internal variables. Then, for all Gauss points in a cluster, only one micro nonlinear problem is solved, and its response is transferred to all integration points of the cluster in terms of mechanical properties. The solution converges with respect to the number of clusters, which is weakly depends on the number of macro mesh elements. Accelerations of FE 2 calculations up to a factor 50 are observed in typical applications. Arbitrary nonlinear behaviors including internal variables can be considered at the micro level. The method is applied to heterogeneous structures with local quasi-brittle and elastoplastic behaviors and, in particular, to a nuclear waste package structure subject to internal expansions.

INTRODUCTION

Machine learning, Artificial Intelligence and Data Sciences have recently attracted a lot of attention in the computational mechanics community. These techniques may offer new possibilities like reducing computational times, automating construction of models from data, or automating mappings between models and material characteristics or microstructural patterns.

One first example of use of machine learning in mechanics is the construction of models from data. For example, Artificial Neural Networks (ANN) have shown promising results to reproduce stress-strain curves of behaviors in complex configurations like high temperatures loadings [START_REF] Shakiba | Modeling the Effects of Cu content and deformation variables on the high-temperature flow behavior of dilute Al-Fe-Si alloys using an artificial neural network[END_REF] , design of steel structures 2 , vibrational behavior of structures [START_REF] Liu | Total solution for structural mechanics problems[END_REF][START_REF] Waszczyszyn | Neural networks in mechanics of structures and materials-new results and prospects of applications[END_REF] , or structural stability problems [START_REF] Rehmantahira | Artificial neural network prediction of buckling load of thin cylindrical shells under axial compression[END_REF] .

Another application importing techniques from machine learning is the recognition of patterns and the construction of a relationship between material topological features and some models. For examples in [START_REF] Rao | Three-dimensional convolutional neural network (3D-CNN) for heterogeneous material homogenization[END_REF][START_REF] Lißner | Data-driven microstructure property relations[END_REF] , machine learning was used to relate microstructures to their effective thermal or mechanical effective properties. In [START_REF] Ryckelynck | Mechanical dissimilarity of defects in welded joints via Grassmann manifold and machine learning[END_REF] Ryckelynck et al. used image recognition and non supervised machine learning techniques for classification issues and the construction of digital twins of material defects in nonlinear mechanical simulations.

More recently, Ortiz et al. [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF] developed a fully data-driven approach to replace constitutive law in finite elements simulations by direct use of data, by searching in a data base the closest data point to be used at the local Gauss points, then removing the use of empirical constitutive laws. In 10 , Carrara et al. constructed fracture models where modeling assumptions were removed and the governing equations stemming from variational principle were combined with discrete data points, leading to a model-free data-driven method for solving fracture problems.

Strategies using machine learning have also recently been developed to describe complex phenomena like dynamics trajectories or fluid sloshing without solving the partial differential equations [START_REF] Mozaffar | Deep learning predicts path-dependent plasticity[END_REF][START_REF] González | Learning non-Markovian physics from data[END_REF] .

Another exciting field of application of these techniques is multiscale modeling. One major difficulty in multiscale numerical simulation occurs when dealing with nonlinear materials. One direct approach to simulate the behavior of a nonlinear heterogeneous structure resorts to the so-called FE 2 technique [START_REF] Feyel | Multiscale FE 2 elastoviscoplastic analysis of composite structure[END_REF][START_REF] Feyel | FE 2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials[END_REF] , where two levels of Finite Element simulations are conducted, one at the scale of the whole structure, and the other at the microscale where nonlinear simulations are performed over Representative Volume Elements. The major drawback of such an approach is its computational cost, related to the unacceptable number of nonlinear calculations to be performed. Among the different strategies to reduce these calculations, one method, initiated in [START_REF] Yvonnet | Numerically explicit potentials for the homogenization of nonlinear elastic heterogeneous materials[END_REF][START_REF] Le | Computational homogenization of nonlinear elastic materials using neural networks[END_REF] , makes use of so-called data-driven approaches in which micro scale calculations are performed during an off-line stage, and which are then used as data in an on-line stage to reconstruct the macroscopic (effective) behavior. For this purpose, several techniques have been proposed, including interpolation methods [START_REF] Yvonnet | Numerically explicit potentials for the homogenization of nonlinear elastic heterogeneous materials[END_REF][START_REF] Fritzen | Two-stage data-driven homogenization for nonlinear solids using a reduced order model[END_REF] , Neural Networks [START_REF] Le | Computational homogenization of nonlinear elastic materials using neural networks[END_REF][START_REF] Lu | A data-driven computational homogenization method based on neural networks for the nonlinear anisotropic electrical response of graphene/polymer nanocomposites[END_REF][START_REF] Li | Multiscale computation on feedforward neural network and recurrent neural network[END_REF][START_REF] Logarzo | Smart constitutive laws: Inelastic homogenization through machine learning[END_REF][START_REF] Nguyen-Thanh | A surrogate model for computational homogenization of elastostatics at finite strain using high-dimensional model representation-based neural network[END_REF][START_REF] Wang | A multiscale multi-permeability poroplasticity model linked by recursive homogenizations and deep learning[END_REF][START_REF] Ghavamian | Accelerating multiscale finite element simulations of history-dependent materials using a recurrent neural network[END_REF] , Bayesian inference [START_REF] Wu | Bayesian inference of non-linear multiscale model parameters accelerated by a Deep Neural Network[END_REF] , Fourier series expansions [START_REF] Peigney | A Fourier-based machine learning technique with application in engineering[END_REF] , Gaussian process [START_REF] Rocha | On-the-fly construction of surrogate constitutive models for concurrent multiscale mechanical analysis through probabilistic machine learning[END_REF] , or self-clustering analysis [START_REF] Liu | Self-consistent clustering analysis: an efficient multi-scale scheme for inelastic heterogeneous materials[END_REF][START_REF] Liu | Microstructural material database for self-consistent clustering analysis of elastoplastic strain softening materials[END_REF] . In the related techniques, the off-line data collection is used in a regression process to construct an accurate surrogate model whose evaluation is of several orders of magnitude lower than performing one micro scale nonlinear calculation. A critical comparison of several regression techniques used in data-driven multiscale approaches can be found in [START_REF] Rocha | Micromechanics-based surrogate models for the response of composites: A critical comparison between a classical mesoscale constitutive model, hyper-reduction and neural networks[END_REF] . In [START_REF] Avery | A computationally tractable framework for nonlinear dynamic multiscale modeling of membrane woven fabrics[END_REF] , Avery et al. investigated and discussed several regression methods with ANN in homogenization problems of hyperelastic woven composites, and demonstrated its use in advanced dynamic/fluid structure applications. Recent advances of data driven techniques, including handling loading history-dependent behaviors like plasticity can be found in [START_REF] Wu | A recurrent neural network-accelerated multi-scale model for elasto-plastic heterogeneous materials subjected to random cyclic and non-proportional loading paths[END_REF][START_REF] Mozaffar | Deep learning predicts path-dependent plasticity[END_REF][START_REF] Ghavamian | Accelerating multiscale finite element simulations of history-dependent materials using a recurrent neural network[END_REF] . On-the fly construction of the surrogate model by probabilistic machine learning was proposed in [START_REF] Rocha | On-the-fly construction of surrogate constitutive models for concurrent multiscale mechanical analysis through probabilistic machine learning[END_REF] . Developments of neural networks techniques in FE [START_REF] Chojaczyk | Review and application of artificial neural networks models in reliability analysis of steel structures[END_REF] , including Feed-Forward and recurrent Neural Networks can be found in [START_REF] Li | Multiscale computation on feedforward neural network and recurrent neural network[END_REF][START_REF] Wu | A recurrent neural network-accelerated multi-scale model for elasto-plastic heterogeneous materials subjected to random cyclic and non-proportional loading paths[END_REF] . In [START_REF] Bhattacharjee | A nonlinear manifold-based reduced order model for multiscale analysis of heterogeneous hyperelastic materials[END_REF][START_REF] Bhattacharjee | A nonlinear data-driven reduced order model for computational homogenization with physics/pattern-guided sampling[END_REF] , both homogenization and localization of heterogeneous hyperelastic materials were performed using a digital database and the manifold-based nonlinear reduced order model (MNROM).

However, the main challenge in the aforementioned techniques is that they are adapted only to some classes of constitutive behaviors and lack generality. More specifically, handling internal variables and loading-dependent behaviors still remains delicate. In addition, neural networks and related surrogate model techniques lack error estimation, especially when the parameters are out of the training parameter data set [START_REF] Rennie | New artificial neural networks for true triaxial stress state analysis and demonstration of intermediate principal stress effects on intact rock strength[END_REF] .

The objective of this paper is then to propose a new machine-learning based multiscale method, called k-means FE 2 , to solve general nonlinear multiscale problems with internal variables and loading-dependent behaviors, without use of surrogate models. For this purpose, we employ a new paradigm: instead of reducing the microscale problem calculations, we reduce the macro scale problem by constructing clusters of Gauss points in a structure which are assumed to be in the same mechanical state. A k-means clustering -machine learning technique is employed to select the Gauss points based on their strain state and sets of internal variable. Then, for all Gauss points in a cluster, only one micro nonlinear problem is solved, and its response is transferred to all points of the cluster. The operation is repeated until convergence is reached at both macro and micro scales. Drastic computational reductions can be achieved, and error can be controlled by gradually increasing the number of clusters at the macro scale. Note that the present method is totally different from the so-called Self-Clustering analysis proposed by Liu et al. [START_REF] Liu | Self-consistent clustering analysis: an efficient multi-scale scheme for inelastic heterogeneous materials[END_REF][START_REF] Liu | Microstructural material database for self-consistent clustering analysis of elastoplastic strain softening materials[END_REF] , which involves clustering algorithms at the micro scale, and which uses a Lippmann-Schwinger integral equation to accelerate micro equations. We summarize the main features of the proposed technique as follows:

• Convergence is obtained with respect to the number of clusters.

• Validation is made with respect to direct FE 2 solutions.

• Convergence is given with respect to the macro mesh size.

• Computational gains between 1 and 2 orders of magnitude can be achieved for similar accuracy when compared to full FE 2 solutions. Larger gains can be obtained if larger errors are accepted.

• A simple error analysis can be achieved by gradually increasing the number of clusters.

• Nonlinear behaviors can be arbitrary (plasticity, damage, etc.) with internal variables.

• The number of clusters weakly depends on the one of macro mesh elements.

The organization of the paper is as follows. In section 2, the basics ofclustering analysis are recalled. In section 3, the k-means FE 2 method is presented. In section 4, numerical examples are provided to illustrate the aforementioned features of the method and its potential use for industrial applications.

K-MEANS CLUSTERING

In this section, we briefly recall the basics of the k-means clustering method [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF] . This algorithm is used to classify data into groups, or clusters, with similar features. Clustering methods are widely used in data mining [START_REF] Fayyad | Advances in knowledge discovery and data mining[END_REF] , data compression [START_REF] Gersho | Vector quantization and signal compression[END_REF] , and pattern classification [START_REF] Duda | Pattern classification and scene analysis[END_REF] . Given a set of vectors

 = 1 , 2 , ..., , (1)
the k-means algorithms aims at defining clusters  , = 1, 2, ..., :

 = 1 , 2 , ..., , (2)
each containing vectors , where can be different for each cluster . We define the centroid of a cluster according to

= 1 ∑ ∈ (3)
and the Euclidean distance between one data point (vector) and the centroid of a cluster by

= ‖ -‖. (4)
where ‖ ‖ = √ ⋅ . Mathematically, the objective of the k-means algorithm is to find the sets  such that the following objective function

= ∑ ∑ ∈ ‖ -‖ 2 (5)
is minimized. The number is chosen by the user. The centroids are initialized by choosing arbitrary vectors ∈ , and assigning = at the first iteration. Then, the algorithm runs as follows [START_REF] Hartigan | AK-means clustering algorithm[END_REF] . At each iteration , for each data point ∈ , the Euclidean distance to the centroid of each cluster is computed according to (4). A data point (vector) is simply assigned to a cluster if its distance to the cluster is minimal among all clusters

 = | ‖ -‖ ≤ ‖ -‖, ∀ with ≠ . ( 6
)
The process goes on until a criterion is reached. Here, we choose the simple stopping criterion:

Δ = | ( ) -( -1) | = 0 (7)
where ( ) is given by ( 5) at iteration of the algorithm. The k-means algorithm is summarized in Algorithm 1.

An illustration in the context of a 3D linear finite element calculation is provided in Fig. 1. Here, the algorithm is applied to classify the Gauss points in the finite elements mesh according to their strain components values. A perforated plate under uniaxial tension is considered, whose dimensions, boundary conditions, and loading are described in section 4.2. Figs. 1a,1b, 1c show the distributions of the strain components 11 , 22 and 12 in the elements of the mesh. Here, as one Gauss point per element is used, the elements are associated to the Gauss points. After solving the FEM problem, the strain components are obtained in each element in a vector form as:

= 11 , 22 , 33 , 2 13 , 2 23 , 2 12 . ( 8
)
Here the number of clusters is chosen as = 5. Fig. 1d depicts the initial data points (elements) selected to initialize the clusters. Figs. 1e-1l show the positions of the clusters at different iterations. In this example, the convergence is achieved in 25 iterations. Note that the convergence rate of the algorithm depends on the number of clusters and on the number of vectors in  (here the number of elements in the FEM mesh).

A K-MEANS FE 2

Brief review of the FE 2 method

The multilevel finite-element method [START_REF] Feyel | Multiscale FE 2 elastoviscoplastic analysis of composite structure[END_REF][START_REF] Feyel | FE 2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials[END_REF] , also called FE 2 in the literature, as it involves two levels of finite-element simulations, and independently proposed by several other authors and groups [START_REF] Ozdemir | Computational homogenization for heat conduction in heterogeneous solids[END_REF][START_REF] Geers | Multi-scale computational homogenization: Trends and challenges[END_REF][START_REF] Terada | A class of general algorithms for multi-scale analysis of heterogeneous media[END_REF][START_REF] Kouznetsova | Multi-scale constitutive modeling of heterogeneous materials with gradient enhanced computational homogenization scheme[END_REF][START_REF] Ghosh | A multilevel computational model for multi-scale damage analysis in composite and porous media[END_REF][START_REF] Geers | Multiscale Modeling of Microstructure-Property Relations[END_REF] , was introduced as a general multiscale method for solving nonlinear heterogeneous structural problems. The basic underlying idea is that two levels of finite element simulations are concurrently solved, one for each scale. At the macroscale, each integration point of the finite-element mesh is associated with a representative volume element (RVE). Boundary conditions depending on the macroscopic strain are prescribed on the boundary of each RVE. After solving each nonlinear RVE problem at each integration point, the macroscopic stress is averaged over the RVE and provided at the macro integration point. Then, the macroscopic constitutive law is available only through solving a nonlinear problem. These operations are repeated until convergence is reached at both scales (see Figure 2 (a)). We consider a macroscopic structure associated with a domain Ω ⊂ ℝ 3 , with a boundary Ω. The assumption of scale separation is adopted (an extension of FE 2 to second-order homogenization can be found in [START_REF] Kouznetsova | Multi-scale constitutive modeling of heterogeneous materials with gradient enhanced computational homogenization scheme[END_REF] ). The microstructure is assumed to be characterized by an RVE associated with a domain Ω ⊂ ℝ 3 , with boundary Ω. For a general nonlinear constitutive law, the strain ( ) is related to the stress ( ) by a nonlinear local constitutive relationship, where is a material point within Ω.

In the following, (.) notations will denote macroscale quantities. In the absence of body forces, for a given macroscopic strain and a set of microscopic internal variables defined by a vector ( ), the RVE problem is to find ( ), such that

∇ ⋅ ( ) = 0 ∀ ∈ Ω, (9)
where ∇ ⋅ (.) is the divergence operator. The constitutive law is given by

( ) =  ( ( ), ( )) . ( 10
)
where  is a local nonlinear operator. Under the assumption of sufficient displacement smoothness, the strain macroscopic strain is defined by

= 1 ∫ Ω ( ) Ω, (11)
where is the volume of Ω. Equation (11) holds, e.g. for the following boundary conditions (see e.g. [START_REF] Yvonnet | Computational Homogenization of Heterogeneous Materials with Finite Elements[END_REF] ):

( ) = + ̃ ( ) on Ω, (12)
where ̃ ( ) is a periodic function over Ω. Again upon a sufficient stress regularity, the macroscopic stress is defined by

= 1 ∫ Ω ( ) Ω. (13)
In the so-called FE 2 method, the constitutive law -is unknown, but can be obtained numerically by solving a nonlinear problem over the RVE, as follows (see Figure 2 Given , ( ):

1. Prescribe boundary conditions (12) on Ω.

2.

Use a numerical method such as FEM with an iterative solver such as the Newton method to solve the nonlinear problems ( 9), (10), and (12) (see details in the following).

Compute the spatial average of the stress over the RVE to obtain .

In what follows, a detailed numerical implementation of the FE 2 method is presented to better understand where problems ( 9), (10), and (12) have to be solved within finite-element calculation at the macroscopic scale. In the absence of body forces, the problem at the macroscale is given by

∇ ⋅ = 0 in Ω, (14)
and completed with boundary conditions such as prescribed tractions * over Ω and prescribed displacemnts = * over Ω . The weak form associated with ( 14) is given by:

∫ Ω , ∶ ( ) Ω -∫ Ω * ⋅ Γ = ( ) = . ( 15
)
In a nonlinear context, Eq. ( 15) must be solved using an iterative method, e.g. a Newton method, as described in the following. The linearized problem related to (15) is given by:

∫ Ω ℂ , ∶ Δ ∶ Ω = -∫ Ω , ∶ ( ) Ω + ∫ Ω * ⋅ Γ, (16)
where and denote the displacement solution and the set of internal variables at the previous -th iteration, respectively, and ℂ = is the tangent elastic tensor. Solving (16), we obtain the correction to be carried out on the macroscopic displacements:

+1 = + Δ . ( 17
)
Problem ( 16) is solved until a convergence criterion is reached, i.e. when the norm of the residual is lower than a prescribed tolerance. Introducing an FEM discretization into (16), we obtain a linear system of equations in the form:

Δ = - ( 18
)
where Δ is the vector containing all nodal corrections of displacement components in the macroscopic mesh,

= ∫ Ω Ω (19)
and

= -∫ Ω Ω + ∫ Ω * Ω, (20)
where is the matrix form of ℂ , is the vector form of , and and are matrices of FEM shape functions and shape functions derivatives related to the macroscopic mesh. The components of ℂ can be evaluated numerically by perturbation using:

ℂ ( ) = ( ) ≃ + ( ) - (21)
with

( ) = 2 ⊗ + ⊗ ( 22
)
where is a small parameter such that << | | | | and , = 1, 2, 3 are unitary basis vectors. We can see that at each point of the macroscopic mesh, determining requires solving one nonlinear problem over the RVE while determining requires solving 6 additional calculations in 3D. We can provide an estimation of the total number of RVE calculations to be performed within a whole FE 2 simulation as follows:

≃ (1 + ) × × × × , (23)
where = 3 in 2D and = 6 in 3D, is the number of integration points per element in the macroscopic mesh, is the number of elements in the macroscopic mesh, is the average number of macroscopic Newton iterations before convergence at each time step, and is the number of time steps for the whole simulation (loading steps in a quasi-static simulation). As an illustration, considering a 3D mesh containing 100,000 elements (which is a rather coarse mesh for practical applications), = 1, = 4 and 10 time steps, = 10, we obtain = 28×10 6 nonlinear calculations on the RVE during the whole simulation. This simple example shows that for most applications, the classical FE 2 method is not applicable, unless very coarse meshes at both micro and macro levels are used. To drastically reduce this complexity, the k-means FE 2 method is proposed in the next section.

k-means FE 2

One bottleneck in the above FE 2 procedure is the total number of nonlinear problems associated with the number × of macroscopic integration points. In practical calculations, several RVEs are subjected to very close macroscopic strains . During a typical multiscale simulation, many redundant problems are then solved unnecessarily. Then, the idea is to identify clusters of integration points in the macro mesh according to the k-means algorithm described in section 2 and to associate a single RVE to each cluster  . The macroscopic strain is defined for the cluster  by:

= 1 ∑ ∈ ∑ ∈ (24)
where is the macroscopic strain at one integration point in the cluster  , is the volume of the element containing the integration point , and is its associated Gauss weight. Then, after solving the RVE problem, the corresponding macroscopic stress is affected to all integration points in  . The k-means algorithm has negligible cost as compared to the RVE calculations. Then, the new number of nonlinear RVE problems to be solved within the k-means FE 2 method can be evaluated as:

≃ (1 + ) × × × , ( 25
)
where is the average number of clusters along the whole simulation. If << × , a drastic time reduction can be achieved.

Clusters updating in presence of internal variables

At each loading step, the k-means algorithm is used to identify the clusters of macroscopic Gauss points in the macro mesh. It is worth noting that between two iterations of the macro Newton scheme, the number and distributions of clusters might evolve. In the case of nonlinear behavior with internal variables at the micro scale, precautions must be taken. The most general approach would be to include the internal variable state ( ) in the clustering analysis. However ( ) is a field defined in all point ∈ Ω of the RVE and its description would involve a high-dimensional vector. Due to this issue, we choose not to include the micro internal variable state ( ) in the cluster analysis, but to assign the RVEs that do not belong to the same cluster from one iteration to another to a new sub-cluster. Then, at one new macro Newton iteration, all RVEs belonging to one sub-cluster have close macroscopic strain and similar internal variables distributions ( ). In the following we define the clusters as the sets of RVEs grouped according to their macroscopic strain values as defined in section 2 and the sub-clusters as those grouped according to both strain and internal variables states ( ). Then, the actual number of RVE computations to be performed is associated with the number of sub-clusters.

To summarize:

• The clusters define regions of the macroscopic mesh with close strain states;

• Sub-clusters are regions of the macroscopic mesh with similar strain state and history of internal variables. Then, each sub-cluster is associated with one RVE calculation;

• The number of clusters is defined by the user and does not change through the iterations. At each iteration, the intersection between the previous and the actual sets of clusters defines the new sub-clusters, whose number may increase during the simulation.

Fig. 3 shows an illustrative example of the evolution of clusters and the number of RVEs (sub-clusters) in a schematic macroscopic structure, where each square represents a Gauss point. Fig. 3(a) displays the position of 3 initial clusters (red, blue and green) which are identical to the sub-clusters. Fig. 3(b)(c)(d) depicts an evolution of the 3 clusters at respective + 1, + 2, + 3-th iterations. The number of sub-clusters in the macroscopic structure has evolved, so that at the + 1-th iteration, cluster 1 contains 3 sub-clusters (3 RVEs) with different internal variable distributions, and cluster 3 comprises two sub-clusters (2 RVEs) with different internal variable distributions. At the end of the + 3-th iteration, the total number of sub-clusters is updated to 13. This procedure of updating the number of RVEs is integrated in the k-means FE 2 method. This approach has the drawback of splitting the clusters into sub-clusters at each change of cluster distributions. However, we will show in the numerical applications that this number still remains much lower than . It is important to note that the present approach does not depend on the number or type of internal variables within the RVE. As the generation of new sub-clusters is purely algorithmic (intersections of current and previous clusters), no a priori mechanical analysis is required. New generated sub-clusters then have the same strain history when they are created, but are subsequently submitted to different strain histories. As shown in the numerical examples, we can then apply this approach to arbitrary local history-dependent behaviors (damage, plasticity, etc.).

The k-means FE 2 procedure has been implemented in the finite element code Cast3m 47 . A summary of the k-means FE 2 is provided in Algorithm 2.

Additional remarks

• Note that for nonlinear behaviors without internal variables, the number of sub-clusters remains identical to the one of clusters and even larger gains of computational times can be achieved.

• The above approach can be combined with parallel computing to perform all RVE simulations associated to sub-clusters on different processors.

• If model reductions or surrogate model methods are available to replace the nonlinear RVE calculations, they can be combined with this framework.

NUMERICAL EXAMPLES

In this section, we present several numerical examples to demonstrate the following features of the method and its interests for nonlinear multiscale computations: (a) the convergence of the method with respect to the number of clusters; (b) the validation of the method with respect to direct FE 2 simulations; (c) the low dependence of the number of clusters on the macro mesh, (d) the computational gains compared with FE 2 method and (e) the potential of the method for arbitrary micro nonlinear behaviors with internal variables.

Material modeling at the micro scale

Small strains are assumed. We consider in the following microstructures composed of a nonlinear quasi-brittle matrix in which are distributed linear elastic or thermoelastic inclusions. As a result, the effective behavior of the material is nonlinear. It is computed by the multiscale method described in the previous sections. Linear and thermoelastic inclusions are described by the following constitutive law:

( ) = ℂ( ) ∶ ( -(Δℎ)) (26)
Algorithm 2 k-means FE 2

1: For = 1 ∶ (Loop over all time steps) 2:

= + 1 3:

Initialize:

+1 0 = 4:
While > (Macro Newton Loop).

5:

= + 1 (iterations) 6:

Input: ( +1 ) = 1, 2... .

7:

Define clusters  = 1, 2, ..., using Algorithm 1.

8:

Define sub-clusters  ; ( = 1 in the first time step, see section 3.3).

9:

Compute according to (24) .

10:

For = 1 ∶ (Loop over all clusters) .

11:

Input 12:
For = 1 ∶ (Loop over all sub-clusters in  )

13:

Solve the nonlinear problem ( 9)-( 12).

14:

Output: , .

15:

( +1 ) = , ∀ ∈  .
16: End.

( +1 ) = , ∀ ∈  .

21:

End. End.

27:

= +1 .

28: End.

where ℂ( ) is the elastic fourth-order tensor at point , and (Δℎ) is an eigenstrain whose evolution is assumed to be given as a function of a parameter Δℎ associated with the evolution of the system. The matrix is assumed to be quasi-brittle and described by the phase field model for fracture, which is briefly reviewed in Appendix A. Then, here, the field of internal variables within the RVE which induces the creation of new sub-clusters (see section 3.3) is the damage field, ( ) ≡ ( ).

Perforated plate in traction: local quasi-brittle behavior

The objective of this first example is to demonstrate the accuracy of the method and its convergence with respect to the number of clusters. We consider a plate of dimension 80 × 80 × 3 mm 3 with a hole of radius = 15 mm at its center. Due to symmetry, only a quarter of the structure is taken into account. The geometry and the boundary conditions are depicted in Fig. 4 (a). The mesh is composed of 401 linear tetrahedral elements (see Fig. 4 (b)). On the upper end of the sample, a displacement along the -direction is prescribed, while the degrees of freedom in the other directions are blocked.

In this example, the RVE is composed of aggregates embedded in a cement matrix. A view of aggregates (inclusions) is depicted in Fig. 5 (a). The mesh is generated by using the python script Combs [START_REF] Bourcier | Combs: open source python library for RVE generation[END_REF] operating in the CAD code Salome [START_REF]Salome software[END_REF] , see e.g. [START_REF] Benaimeche | Analytical and numerical prediction of mechanical properties of heterogeneous materials with expansive inclusions: Application to waste cementation[END_REF][START_REF] Bary | Analytical and 3D numerical analysis of the thermoviscoelastic behavior of concrete-like materials including interfaces[END_REF] and references therein. The mesh contains a total of 28042 tetrahedral elements. The volume fraction of inclusions is = 10%. The inclusions are assumed to be linear elastic with the following properties: = 80 , = 0.19, which denote the Young's modulus and the Poisson's ratio of inclusions, respectively. The matrix is assumed to be isotropic quasi-brittle and The k-means FE 2 method is applied for different numbers of clusters and its convergence is verified regarding the maximum von Mises stress in the macroscopic mesh with respect to the number of clusters in Fig. 6. We can note that 9 clusters are sufficient here to reach the convergence with regard to the chosen criterion, and that this converged solution is close to the reference FE 2 solution. In Fig. 7, we depict the force-displacement curve of the macro structure with respect to the number of clusters, again illustrating the convergence to the reference solution when increasing the number of clusters.

Figure 8 shows the evolution of the clusters during the simulation, respectively at the first loading step and at the final loading step. It is worth noting that as stated in section 3.3, the position of the clusters evolves during the simulation with the macro strain field.

We show the convergence of the local macro stress with respect to the number of clusters at the final loading step in Fig. 9, by depicting the von Mises stress fields in the structure for = 5, 6, 8 and comparing it with the reference FE 2 solution. Finally, we report in Table 1 the total number of nonlinear RVE problems to be solved using the k-means FE 2 for different numbers of clusters, as well as the corresponding computational times, and compare it to the full FE 2 method. We can note a drastic reduction of computational times, with a ratio of about 15 for 8 clusters, which corresponds to the converged solution in Fig. 7.

Total

3-point bending beam: local quasi-brittle behavior

In this next example, we show an additional feature of the method: the low dependence of the number of clusters with respect to the mesh size. For this purpose, we consider a homogeneous beam of dimensions 200 × 40 × 2 mm 3 in 3-point bending, as depicted in Fig. 10 (a). On ( = 0, = 0) mm, all degrees of freedom are blocked, whereas on ( = 200, = 0) mm the and displacements are blocked. A displacement in thedirection is prescribed at the middle of the beam. To test the influence of the macroscopic mesh refinement on the number of obtained clusters, 3 macroscopic meshes have been tested, including, respectively, 262, 431 and 786 tetrahedral elements (see Fig. 10). In this example, the RVE consists into a spherical elastic inclusion embedded into a quasi-brittle matrix. The dimensions of the RVE are 1 × 1 × 1 mm 2 . The inclusion is centered in the cubic RVE and has a radius = 0.14 mm. The mechanical properties of the inclusion and the matrix are the same than in the previous example. The k-means FE 2 method is applied for the different meshes and for different numbers of clusters. Convergence of the solutions with respect to the number of clusters for the different meshes is presented in Fig. 11, where the maximum von Mises stress is indicated at the final loading step of the simulations. For meshes 1 and 2, we have verified that the converged k-means solutions are in good agreement with the reference FE 2 solution. We did not perform this verification for the finest mesh, as the corresponding FE 2 computational costs was evaluated as too high (see table 2). As another observation, the total number of clusters required to converge to the reference solution does not increase significantly with the macro mesh size. This constitutes another important asset of the present method, even though this is obviously dependent on the regularity of the macro strain field. As another remark, it is worth pointing out that in the applications of the k-means FE 2 involving local damage, numerical regularization should be included at the macroscopic scale to guarantee convergence with respect to the macro mesh. For the sake of simplicity we did not implement such process here. Extensions with techniques like nonlocal damage model [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF] , gradient enhanced damage [START_REF] Peerlings | Gradient-enhanced damage for quasi-brittle materials[END_REF] , or gradient damage 54 could be considered in future studies. In addition, the existence of an RVE in the case of softening behavior is still a matter of debate in the literature and appropriate averaging techniques or analyses (see e.g. [START_REF] Vinh Phu | On the existence of representative volumes for softening quasi-brittle materials-a failure zone averaging scheme[END_REF] ) may be required in future extension of the present method in the case of local damage. As there is no macro regularization here, the 3 different meshes studied above are considered as 3 distinct cases and we do not intend to show a convergence with respect to the macro mesh in this example.

Figs. 12(a), (b) and (c) show the load-displacement curve of the beam subjected to 3-point bending for the 3 meshes and different numbers of clusters, and Fig. 12(d) shows the load-displacement curve of the 3 structures at convergence ( =8). According to these results the convergence of the model is reached for 8 clusters for the 3 meshes. The convergence with respect to the mesh refinement is also appreciated.

In Fig. 13, we illustrate the evolution of the clusters for the medium mesh when 5, 6 and 8 clusters are considered, for initial and final loading steps. As mentioned in section 3.3, the position of clusters changed in the macroscopic structure at the initial and final step, therefore the number of nonlinear RVE problem (sub-clusters) evolved.

In Fig. 14, we compare the macroscopic damage fields on the medium mesh for different numbers of clusters and the reference FE 2 solution. To better visualize the crack, a pseudo macro damage field is defined as = 1 -

√ ℂ 1111 ( ) ℂ 0 1111
, where denotes pseudo time (loading evolution). We can see a crack initiating in the middle of the beam and propagating from the bottom to the top. The shape of the damaged area depends on the number of clusters , and converges to the reference solution (FE 2 ) when = 9.

Finally, we compare the computational times on the same Linux machine as in the previous example in Table 2, as a function of RVE calculations and number of clusters. Results show that the number of RVEs and thus the computational times are directly impacted by the mesh density in FE 2 , while it affects in a minor way the computational times in k-means FE 2 . The results show that keeping 8 clusters (corresponding to a converged solution), a speed-up ratio of almost 30 in comparison with FE 2 can be achieved for a similar accuracy. For the finest mesh, we only estimate the computational times, as the duration of the full FE 2 is not affordable. Then, the present technique allows solving problems with mesh sizes which are too costly in direct FE 2 . Note that in the present work we have chosen damage as a local nonlinear behavior, but any other arbitrary local behaviors can be considered, like elastoplasticity or viscoplasticity.

As a final remark, if no FE 2 reference solution is available, the low computational costs of the method allow to use it for checking the convergence with respect to the number of clusters. In that case, the total computational costs must be calculated as the sum of the times related to several successive k-means analyses, with increasing numbers of clusters . For example, in the present case, considering the 3 successive cases = 5, = 6, = 8, the gains are still 5.4, 7.1 and 10.8 for meshes 1, 2 and 3, respectively, then yet an order of magnitude as compared to a direct FE 2 computation.

Number of clusters

Coarse mesh

Perforated plate in traction: local elastoplastic behavior

In this example, the k-means FE 2 method is employed to solve the same macroscopic structure than in section 4.2 with the same boundary conditions, but with an elastoplastic behavior within the RVE. At the microscopic scale, a porous RVE of dimensions (1 × 1 × 1 mm 3 ) is considered, with a pore in the center of radius r = 0.15 mm. The behavior of the matrix is elastoplastic with linear isotropic hardening under the 2 -flow theory. In that context, the constitutive law is given by: where is the plastic strain and where ℂ is the elastic tensor, assumed to be isotropic. The elastic domain is defined by

= ℂ ∶ ( -) (27)
0,0E+00 8,0E-04 1,6E-03 2,4E-03 3,2E-03 4,0E-03 0,0E+00 2,0E+06 4,0E+06 6,0E+06 8,0E+06 1,0E+07 Load (N
( , ) = 2 ( ) - - ≤ 0 (28)
where 2 ( ) = 3 2 ∶ , = -1 3 ( ), (.) being the trace operator, is the yield stress, is the linear hardening modulus and the cumulated equivalent plastic strain expressed by:

( ) = ∫ 0 ̇ ( ) (29)
with

̇ = 2 3 ̇ ∶ ̇ 1∕2
. The evolution of the plastic strain is given by the normality rule:

̇ = ̇ (30)
The mesh contains a total of 30452 tetrahedral elements. The mechanical properties of the matrix are as follows : = 15 , = 0.2, = 60 and = 1 . The k-means FE 2 method is applied for different numbers of clusters. Here, the field of internal variables within the RVE which induces the creation of new sub-clusters (see section 3.3) is the plastic strain field, ( ) ≡ ( ). The convergence is verified with respect to the maximum von Mises stress in the macroscopic structure for different number of clusters in Fig. 15(a). Convergence with respect to the number of clusters is reached for 9 clusters. In Fig. 15(b), the force-displacement curve of the macro-structure for different number of clusters is presented. Figure 16 shows the von Mises stress field in the macroscopic structure for = 7, 8, and 9 clusters and for the reference solution FE 2 at the final loading step, where again the convergence with respect to the number of clusters can be appreciated.

Table 3 summarizes the total number of nonlinear RVE problems solved by the k-means FE 2 method and the corresponding computation times for k = 6,7,8, and 9 clusters, as compared to a full FE 2 computation. A significant gain in terms of computational time is noticed, with a ratio of 17 for 9 clusters which corresponds to the converged solution.

Industrial example: waste package bahavior

The objective of this last example is to show the potential of the present method for industrial applications. We are concerned here with contaminated wastes generated by industrial processes, containing different types of materials including metallic particles. In order to immobilize these wastes, a solidification process may be used to aggregate them into packages that have to fulfill minimal mechanical characteristics. The choice of the matrix material mainly depends on the type and the behavior of the particles in contact with the matrix. We consider a waste package where the matrix is made of hydrated cement paste, used for its low porosity and low permeability to liquids and gases in the hardened state, and then for its resistance to aggressive environments. However, the presence of free-water in the pores induces a progressive corrosion of reactive metal inclusions [START_REF] Bart | Cement-based materials for nuclear waste storage[END_REF] . The corrosion reaction forms an expansive metal oxide layer around the metal inclusions, usually accompanied with a significant release of hydrogen gas [START_REF] Banos | Corrosion of uranium in liquid water under contained conditions with a headspace deuterium overpressure. Part 2: The ternary U+ H2O (l)+ D2 system[END_REF][START_REF] Haschke | Corrosion of uranium in air and water vapor: consequences for environmental dispersal[END_REF][START_REF] Stitt | In-situ, time resolved monitoring of uranium in BFS: OPC grout. Part 2: Corrosion in water[END_REF] . Moreover, the thickness of the developed corrosion product layer can be important with respect to the metal inclusion, and the transformation process may in theory be complete if reactants are available (the whole inclusion can be transformed into corrosion products). This volume expansion associated with the formation of the corrosion products (metal oxide) can generate spurious micro-cracks in the cement matrix [START_REF] Stitt | In-situ, time resolved monitoring of uranium in BFS: OPC grout. Part 2: Corrosion in water[END_REF][START_REF] Stitt | In-situ, time resolved monitoring of uranium in BFS: OPC grout. Part 1: Corrosion in water vapour[END_REF] .

The studied macro structure is depicted in Fig. 17 where only a quarter is considered due to symmetry. The constitutive material is composed of a matrix of cement paste in which are dispersed non-reactive inclusions (graphite) and metallic particles. This mixture is solidified within a cylindrical container (see Fig. 17). The thickness of the external metallic shell is 15 mm in the cylindrical area and 20 mm in flat bottom area.

The RVE corresponding to the mixture is depicted in Figs. 18. Three cases are considered: (a) a single metallic inclusion in a cement paste matrix with equivalent initial radius of 3.5 mm representing a volume fraction of 0.05 %; (b) 10% in volume of graphite and 0.05% of metallic inclusions and (c) 20% of graphite and 0.05% of metallic inclusions. The metallic inclusions have an equivalent radius = 2.5 mm in the last 2 cases. The graphite inclusions have polyhedral shapes. The size of the cubic RVE is 100 mm .

At the micro scale, the behavior of the graphite inclusions is assumed to be linear elastic, with properties provided in Table 4. The cement paste matrix is assumed to be quasi-brittle and described by the damage phase field model given in Appendix A. The behavior of uranium is linear thermoelastic, with linear properties described in Table 4 and with a constitutive law in the form:

where * = Δℎ is an eigenstrain associated with the swelling of the oxidized metallic inclusions, with = 0.182 m -1 (details on the calculation of such coefficient can be found e.g. in [START_REF] Benaimeche | Analytical and numerical prediction of mechanical properties of heterogeneous materials with expansive inclusions: Application to waste cementation[END_REF] ). In Eq. 31, Δℎ is the thickness of the corrosion product layer. In our simulations, Δℎ was taken as the loading parameter. When Δℎ increases, the local strains near the inclusions within the RVE create damage which modifies the macro behavior. In what follows, we study the evolution of the local stiffness within the structure for increasing values of the corrosion thickness using the proposed k-means FE 2 .

For illustrative purposes, Fig. 19 shows the macroscopic strain components with the effective expansion of * = 9.1 × 10 -7 .

Application 1

In this example the convergence of the k-means FE 2 with respect to the number of clusters in the macro structure is tested. The microstructure illustrated in figure 18 (a) is considered in the microscopic computations. The convergence of the k-means FE 2 method is verified by two different ways: (i) convergence with respect to the maximum von Mises stress within the structure at the final loading step and (ii) convergence with respect to the effective stiffness in chosen macro elements (when the value becomes quasi-constant as a function of Δℎ). Elements 1 , 2 , 3 , and 4 are considered (see Fig. 17 (b)). Fig. 20 shows the position of clusters in the structure for = 5, 6, 8 clusters at the final loading step. Fig. 21 shows the decrease of effective stiffnesses as a function of the thickness of the corrosion product layer Δℎ for different numbers of clusters. For elements 1 and 4 , convergence is verified for 6 clusters. A faster convergence in element 1 than in the element 4 is observed, which can be explained by the size of the cluster to which the element 1 belongs, which remains unchanged for the different values (see the red cluster in Fig. 20). As a comparison, an estimated computational time for this problem using the full FE 2 method would take 759 hours (45570 RVE calculations), while for = 8 clusters the k-means FE 2 only uses 8.9 hours (538 RVE calculations).

Application 2

In this second case, we consider both swelling metallic inclusions and graphite inclusions (Fig. 18). The RVEs meshes in Figs. 18 (b) and (c) contain 53706 and 73389 tetrahedral elements, respectively. Here = 8 clusters have been used. Fig. 22 shows the effective stiffness as a function of the thickness of the corrosion product layer Δℎ. When Δℎ = 0, the effective elastic moduli are the one of the undamaged material. A decrease of the effective stiffnesses is observed with the increase of the corrosion reaction in all elements of the package. In addition, the effective behavior of the material becomes progressively anisotropic because of the local distribution of graphite inclusions. Furthermore, a larger decrease is observed for all stiffnesses in the element 4 as compared to the elements 1 , 2 , and 3 in the two microstructures (b) and (c). The damage distribution of microstructures (b) and (c) in the elements 1 , 3 , and 4 is shown in Fig. 23, in which as expected, cracks appear around the metallic inclusions, and then they propagate within the microstructure. In addition, the damage is localized where the metallic inclusions are close to each other, and has a tendency to propagate in the direction of the nearest inclusion.

CONCLUSION

In this work, we have presented a new machine-learning based multiscale method, called k-means FE 2 , to solve general nonlinear multiscale problems with internal variables and loading history-dependent behaviors, without use of surrogate models.

In contrast with machine learning multiscale methods where a surrogate model is constructed to describe the micro nonlinear effective behavior, we reduce here the macro scale problem by constructing clusters of Gauss points in the structure which are assumed to be in the same mechanical state. A k-means clustering -machine learning technique is employed to select the Gauss points based on their strain state and sets of internal variable. Then, for all Gauss points in a cluster, only one micro nonlinear problem is solved, and its response in terms of effective mechanical properties is transferred to all points of the cluster. The operation is repeated until convergence is reached at both macro and micro scales. The method presents several advantageous features. First, its solution is convergent with respect to the number of clusters. This is expected since in the limit of one cluster for one single macro mesh element, we recover the direct FE 2 solution. We have shown in the different numerical examples that even for a low number of clusters, a very accurate solution can be achieved as compared to the full FE 2 solution taken as a reference. In addition, we have shown that the total number of clusters, which is directly related to the computational time, is only weakly dependent on the mesh refinement in the macro structure. In fact, it only depends on the regularity of the fields at the macro level. We have demonstrated accelerations of FE 2 calculations up to a factor 50 in typical applications involving local fracture within a Representative Volume Element. It is worth noting that arbitrary nonlinear behaviors can be considered at the micro level including internal variables and loading history-dependent behaviors. In contrast, in surrogate-based machine learning methods, handling loading history-dependent local behaviors is delicate. Finally, the convergence of the method with respect to the number of clusters allows determining the necessary number of clusters without any knowledge about the final solution by gradually increasing the number of clusters and checking the convergence. The method has been applied to heterogeneous structures with local quasi-brittle behavior described by damage, and to a nuclear waste storage structure of industrial interest. As a final remark, this framework can be complementary to parallel computing and to machine-learning surrogate models to reduce the RVE calculations for nonlinear multiscale calculations.

iteration 25 FIGURE 1 FIGURE 2

 2512 FIGURE 1 Application of the k-means algorithm for classifying the elements in a FEM linear calculation according to their strain components: (a), (b), (c): components of the strain tensor; (d): initialization of the clusters by choosing arbitrary elements in the structure; (e)-(l): positions of the clusters at iterations 1, 2,3,5,9,15,20 and 25 respectively.

  (a)):

FIGURE 3

 3 FIGURE 3 Schematic evolution of the clusters and number of RVE problems to be solved during a simulation. The square domains represent the elements of a macrcoscopic mesh. At first iteration of the Newton scheme (a), clusters are determined on the sole definition of macroscopic strains. At the second iteration (b), the clusters may be modified. The elements intersecting old and new clusters are splits into new "sub-clusters". Each sub-cluster is associated wih a nonlinear RVE calculation. For the next iterations (c) and (d), new sub-clusters can be formed.

FIGURE 4 FIGURE 5

 45 FIGURE 4 Perforated plate macro structure: (a): geometry and boundary conditions; (b) mesh.

FIGURE 6

 6 FIGURE 6 Convergence of the k-means FE 2 solution (von Mises stress) with respect to the number of clusters.

2 FIGURE 7

 27 FIGURE 7 Perforated plate: convergence of the load-displacement curve with respect to the number of clusters and reference FE 2 solution.

(a) 5 FIGURE 8

 58 FIGURE 8 Cluster positions for the medium mesh at the first step (a) (b) (c) , and the final step (d) (e) (f) in the macroscopic structure for = 5, 6, 8.

( a ) 2 FIGURE 9

 a29 FIGURE 9 von Mises stress in the macroscopic structure for = 5, 6, 8 clusters, k-means solutions (a), (b) and (c), and FE 2 reference solution (d).

FIGURE 10 3 -

 3 FIGURE 10 3-point bending beam; (a) geometry and boundary conditions; (b) coarse mesh; (c) medium mesh; (d) fine mesh.

FIGURE 11

 11 FIGURE 11 Convergence of the k-means FE 2 solutions with respect to the number of clusters for different meshes. The symbols indicate the maximum von Mises stress at the final loading step.

FIGURE 12

 12 FIGURE 12 (a)-(c): Load-displacement curves of the 3 meshes with different number of clusters and FE 2 reference solution (except for (c)) and (d): comparison of the curves for the three meshes at convergence with respect to the number of clusters ssumed to be obtained for 8 clusters.

(a) 5 FIGURE 13

 513 FIGURE 13 Cluster positions at the first step (a) (b) (c) , and at the final step (d) (e) (f) in the structure 2 for = 5, 6, 8.

2 FIGURE 14 = 5 ;

 2145 FIGURE 14 Macroscopic damage = 1 -√ 1111 ∕ 0 1111 for (a) = 5; (b) = 6, and (c) = 9 clusters; (d) reference (FE 2 ) solution.

2 FIGURE 15

 215 FIGURE 15 (a) Convergence of the k-means FE 2 solution (maximum von Mises stress) with respect to the number of clusters (b) convergence of the load-displacement curve with respect to the number of clusters and reference FE 2 solution.

2 FIGURE 16 von

 216 FIGURE 16 von Mises stress in the macroscopic structure for = 7, 8, 9 clusters, k-means solutions (a), (b) and (c), and FE 2 reference solution (d).

FIGURE 17

 17 FIGURE 17 Waste package structure: (a) geometry and (b) macro mesh.

FIGURE 18 3

 18 FIGURE18 3 RVEs in the waste package structure problem; a) a single metallic inclusion in a cement paste matrix; (b) graphite (10 % and metallic inclusion 0.05%) and (c) graphite (20 % and metallic inclusion 0.05 %).

FIGURE 19

 19 FIGURE 19 Macroscopic strain components under effective expansion * = 9.1 × 10 -7 .

FIGURE 20

 20 FIGURE 20 Clusters position for ( = 5, 6, 8) at the final loading step.

FIGURE 21

 21 FIGURE 21 Effective stifnesses 1111 and 3333 as a function of Δℎ for the elements 1 and 4 (see Fig.17(b)) for different numbers of clusters K=5,6 and 8.

17

 17 FIGURE 21 Effective stifnesses 1111 and 3333 as a function of Δℎ for the elements 1 and 4 (see Fig.17(b)) for different numbers of clusters K=5,6 and 8.

  FIGURE 21 Effective stifnesses 1111 and 3333 as a function of Δℎ for the elements 1 and 4 (see Fig.17(b)) for different numbers of clusters K=5,6 and 8.

FIGURE 22

 22 FIGURE 22 Effective stifnesses 1111 , 2222 and 3333 as a function of Δℎ for the microstructures (b) and (c).

FIGURE 23

 23 FIGURE 23 Local damage in RVEs designated by 1 , 3 and 4 (see Fig.17(b) for the definition of these elements), for microstructures (b) and (c) (see Fig.18).

17

 17 FIGURE 23 Local damage in RVEs designated by 1 , 3 and 4 (see Fig.17(b) for the definition of these elements), for microstructures (b) and (c) (see Fig.18).

  FIGURE 23 Local damage in RVEs designated by 1 , 3 and 4 (see Fig.17(b) for the definition of these elements), for microstructures (b) and (c) (see Fig.18).

TABLE 1

 1 Total number of nonlinear RVE problems and computational times on 32-cores LInux PC with 256 Go RAM with respect to the number of clusters for k-means and FE 2 solutions.

			number	Final computation
		of	calculated	time (hours)
		RVE problems
		(number	of
		sub-clusters)
	5	234		3.85 h
	6	291		4.87 h
	8	481		8.02 h
	9	540		9.01 h
	10	611		10 h
	FE 2 7426		123 h

TABLE 2

 2 Comparisons of computational costs for both k-means FE 2 and FE 2 methods for the 3 meshes and different numbers of clusters.

	Medium mesh	Fine mesh

TABLE 3

 3 Elastoplastic problem : Total number of nonlinear RVE problems and computational times on 32-cores Linux PC with 256 Go RAM with respect to the number of clusters for k-means and FE 2 solutions.

		Total	number	Final computation
		of	calculated	time (hours)
		RVE problems
		(number	of
		sub-clusters)
	7	873		13.7 h
	8	1062		17.87 h
	9	1155		19.02 h
	FE 2 19862		331 h

TABLE 4

 4 Elastic parameters of the different phases.

ACKNOWLEDGEMENTS

The financial support of Commissariat à l'Energie Atomique et aux Energies Alternatives (CEA) and of Orano, France, is gratefully aknowledged.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX A BRIEF REVIEW OF PHASE FIELD METHOD FOR QUASI-BRITTLE FRACTURE

The phase field method [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | The Variational Approach to Fracture[END_REF][START_REF] Pham | The variational approach to damage: I. The foundations[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Buliga | Energy Minimizing Brittle Crack Propagation[END_REF][START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] is adopted here at the microscale to model the damage through microcracking within the microstructure. The formulations are briefly reviewed in this appendix.

For an elastic cracked body defined in a domain Ω ⊂ ℝ 3 containing sharp cracks denoted collectively as Γ, the total energy of the system is defined as:

where Ψ( , Γ) is the elastic strain density function and is the critical energy release rate in the sense of Griffith. In Eq. (A1), Ω denotes the portion of the boundary Ω where tractions * are prescribed. The above energy form can be replaced by a regularized one, given by:

where is a crack density function, whose model can be chosen among several possible forms, leading to a class of shapes for the regularized damage field near the crack (see e.g. [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] ). In this regularized framework, the cracks are no more described by surfaces but by a smooth field ( ). For an isotropic elastic solid, Ψ can be decomposed according to:

to only affect the damage to traction modes, and where ( ) is a degradation function such that (0) = 1, (1) = 0 and ′ (1) = 0 and Ψ + ( + ) and Ψ -( -) denote parts of the strain density related to tensile and compressive parts of the strain tensor, respectively, which are defined here following Miehe 67 by:

This form allows avoiding interpenetration when the cracks are closed without any special algorithm for auto-contact, which renders the implementation very simple. In (A4), and denote the elastic Lamé's constants in each phase. The operator ⟨.⟩ ± is defined as ⟨ ⟩ ± = ( ± | |)∕2, + is the tensile part while -is the compression part of the strain tensor obtained by the spectral decomposition:

in which v are the eigenvectors of the strain tensor . Other decompositions have been proposed as in [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] He | Closed-Form Coordinate-Free Decompositions of the Two-Dimensional Strain and Stress for Modeling Tension-Compression Dissymmetry[END_REF][START_REF] Nguyen | Implementation of a new strain split to model unilateral contact within the phase field method[END_REF] .

Variational approach to damage in tandem with regularization, called in the literature "phase field method", implies the minimization of the total energy with respect to the displacement field and the minimization of the energy with respect to the scalar field describing the crack surface in a smooth manner. This second minimization is subjected to an inequality constraint ̇ ≥ 0. To formulate this minimization problem in a simpler setting, a time-stepping  = 0 , 1 , ..., , +1 , ..., can be introduced. At each time step +1 , the problem is to find the displacement fields +1 and +1 such that

where  is a set of kinematically admissible fields. One possible algorithm to solve this problem is to use sequential solving of both minimization problems as

where ( ) is the Gateaux derivative, defined as:

The first equation (A7) defines the mechanical problem while the second one (A8) defines the phase field problem. These two problems are coupled as both involve the fields and .

Eq. (A7) can be developed according to:

where Ψ ( , ) = .

(A11)

For ( ) = (1 -) 2 + , where is a small scalar parameter introduced to avoid ill-conditioning when elements are fully damaged, and with Ψ defined as in (A3) we obtain:

It yields the classical weak form of the mechanical problem as follows: find ∈  = | = on Ω , ∈ 1 (Ω) such that:

for all ∈ 1 0 (Ω) = | = 0 on Ω , ∈ 1 (Ω) , where Ω denotes the portion of Ω where Dirichlet boundary conditions are prescibed. The Euler-Lagrange equation (strong form) associated with Eq. (A13) is given by:

where * are prescribed displacements. The first equation in (A8) can be developed as:

Choosing as

we obtain, after some calculations (see e.g. [START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF] ):

The associated Euler-Lagrange equations to (A17) are given by:

where Δ denotes the Laplacian operator. We follow the classical Miehe's approach [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] to enforce the irreversibility condition by introducing a history function  which substitutes Ψ + to handle loading and unloading and defined as:

It leads to:

and to the corresponding weak form: find ∈ 1 (Ω) such that

for all ∈ 1 (Ω). The problem can be solved by classical Finite Elements (see e.g. [START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF] ).