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

HAL
open science

Machine learning-boosted nonlinear homogenization

Mikhael Tannous, Chady Ghnatios, Olivier Castelnau, Pedro Ponte

Castañeda, Francisco Chinesta

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

► To cite this version:

Mikhael Tannous, Chady Ghnatios, Olivier Castelnau, Pedro Ponte Castañeda, Francisco Chinesta. Machine learning-boosted nonlinear homogenization. Mechanics of Materials, 2025, 201, pp.105229. 10.1016/j.mechmat.2024.105229. hal-05161100

HAL Id: hal-05161100

https://hal.science/hal-05161100v1

Submitted on 14 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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.

Public Domain
# Machine learning-boosted nonlinear homogenization

Mikhael Tannous<sup>a, ∗</sup>, Chady Ghnatios<sup>b</sup>, Olivier Castelnau<sup>a</sup>, Pedro Ponte Castañeda<sup>c</sup>,
Francisco Chinesta<sup>a</sup>

Keywords:

Nonlinear material behavior

Homogenization

Machine learning

Digital twin

Two-phase composite

Polycrystal

Previous research has established nonlinear homogenization as an efficient technique for deriving macroscopic constitutive relations and field statistics in heterogeneous (i.e. composite) materials. This method involves optimal linearization of the nonlinear composite, resulting in a best linear comparison composite that shares identical microstructure and field statistics with the nonlinear material. However, the computational time associated with this method increases as the fidelity of the material representation improves, limiting its practical implementation in commercial finite element software for large-scale structural calculations in which a Representative Volume Element must be considered at each integration point. To overcome this limitation without sacrificing precision or efficiency, machine learning can be employed to develop a digital twin of the homogenization-based constitutive law. This approach enables real-time prediction of macroscopic material behavior while maintaining accuracy. The effectiveness of this approach has been demonstrated for two-phase composites with nonlinear power-law constitutive relations, and it has been successfully extended to model the complex three-dimensional behavior of viscoplastic polycrystals. In the latter case, a significant reduction in computational time has been achieved without compromising the precision of nonlinear homogenization method outputs.

# 1. Introduction

Introducing complex nonlinear heterogeneous material models into structural mechanics applications remains a difficult task from the computational viewpoint. In a large structure, the thermomechanical history experienced at different positions in the computation domain can differ substantially, with regions operating for example in the elastic linear domain and other regions experiencing irreversible deformation. Very often, the mechanical behavior depends on the whole loading history (e.g. strain hardening in metallic alloys) so that the irreversible deformation has to be computed with a constitutive relation that depends on the local thermomechanical state until the present time. This is for example the case for the computation of metal forming, car crashes, thermal shock, geophysical flows, etc. Thus, in a large structure the material behavior must be evaluated at each integration point. Even if equilibrium remains linear, it is global, whereas the behavior, nonlinear, remains local. However, one must compute and store the material's state at each integration point all along the simulation (loading history).

When the material exhibits a microstructure, the problem becomes
even harder. If the microstructure is random, such as porous materials,

polycrystalline aggregates or any random mixture of different mechan-
ical phases, a Representative Volume Element (RVE) of the material
must then be defined at each integration point. The RVE must be large
enough to describe accurately enough the microstructure in a statistical
way and the microstructure evolution along the thermomechanical
loading. For example, for polycrystalline aggregates, one typically has
to account for ~1000 different crystal lattice orientations. Full-field or
mean-field homogenization techniques operate on the RVE to estimate
or bound the actual effective behavior, stress/strain field statistics and
associated microstructure evolution from the actual stress/strain and
the experienced thermomechanical history.

In multiscale structure computations, the mechanical response at
the fine scale can be computed either using a full-field or a mean-field
techniques (Lebensohn et al., 2011), whereas the large scale (structure)
response is generally performed using the Finite Element (FE) method.
Published work include the coupling of two FE models (so-called
FE²) (e.g. Feyel and Chaboche, 2000), the coupling of the structure
FE model with the spectral (FFT) method of Moulinec and Suquet
(1994) (e.g. Spahn et al., 2014), and the coupling of the structure

<footnote><sup>a</sup> PIMM Lab, UMR CNRS 8006– Arts et Metiers Institute of Technology – Cnam, Paris, France</footnote>

<footnote><sup>b</sup> Mechanical Engineering Department, University of North Florida, Jacksonville, 32224, FL, USA</footnote>

<footnote><sup>c</sup> Department of Mechanical Engineering and Applied Mechanics, University of Pennsylvania, Philadelphia, USA</footnote>
FE model with mean-field homogenization models based on the self-consistent scheme (Blackman et al., 2017; Zecevic et al., 2017; Gu et al., 2017). In all these approaches, even if the use of mean-field homogenization techniques is much more efficient from the numerical point of view, the computation time is a real concern because complex calculations are needed at each integration point.

Machine-learning approaches have been applied to reproduce the effective behavior of polycrystalline material in a number of recently published papers. This include the work by Hasan et al. (2022), where the authors apply a ML-based approach to predict the elastic behavior of polycrystal made of cubic grains, assuming a uniform deformation throughout the aggregate (Taylor bound). More advanced estimation of the effective behavior, but still limited to linear behavior (conductivity and elastic stiffness) have been conducted by Dai et al. (2023) with the help of a full-field Finite Element method, limiting thus the size of the whole dataset to 5000 points due to significant computational efforts. A crystal plasticity based Finite-Element approach has also been used in McElfresh et al. (2022) to predict the hardening rate and yield stress in titanium polycrystalline aggregates, where four micro-structural parameters have been defined to describe the various micro-structures, leading to a sparse dataset (2000 points). Gu et al. (2024) make use of a neural network to predict the yield stress of polycrystalline aggregates as a function of micro-structural features, using a highly efficient uniform stress model (Reuss bound) for the building of a large dataset $2 \times 10^6$ points, but neglecting the potentially large impact of stress concentration. Liu et al. (2023) propose surrogate models for the elasto-viscoplastic behavior of polycrystalline aggregate, accounting for a crystal plasticity constitutive relation at the grain scale, for various 1D, 2D and 3D microstructures and dynamic loadings. In the case of 3D polycrystal, the effective behavior of the dataset is computed by a Taylor averaging approach, thus neglecting the effect of strain heterogeneities in the material. In Dai et al. (2021), a neural network has been trained on an earlier tangent extension of the self-consistent scheme for the prediction of the flow stress and texture evolution in polycrystalline copper.

Nonlinear homogenization is concerned with the extension of classical homogenization results for composites and polycrystals with linear properties, such as the bounds of Hashin and Shtrikman (1963) (HS) and the self-consistent (SC) estimates (Kröner, 1958) for composites and polycrystals with nonlinear (plastic, viscoplastic, etc.) properties. Early methods, such as the 'tangent' method of Hill (1965) and the 'secant' method (e.g., Tandon and Weng, 1988), provide estimates that, while reasonably accurate for weak nonlinearities, yield predictions for strongly nonlinear composites that are much stronger than numerical results obtained by means of full-field numerical simulations (see, for example, Ponte Castañeda and Suquet, 1998). In order to soften the response of the predictions for polycrystals an ad hoc method was proposed by Molinari et al. (1987), which gives improved results for certain types of polycrystals, but which can under-predict the effective flow stress of other types of polycrystals, especially for strong nonlinearities. Novel nonlinear homogenization methods yielding much more accurate predictions for the macroscopic response of these composites were introduced in the eighties and nineties. In particular, Ponte Castañeda (1991) developed new variational principles for the properties of suitably defined 'linear comparison composites' (LCCs), which allowed the direct conversion of estimates, such as the above-mentioned HS and SC estimates, for linear-elastic composites into corresponding estimates for viscoplastic composites. Corresponding variational linear comparison (VLC) were obtained by Nebozhyn et al. (2001) for viscoplastic polycrystalline materials. However, because of their bounding status, these estimates tended to yield estimates that were stronger than the corresponding numerical simulations. For this reason, improved 'second-order' variational estimates making use of more general linear comparison composites were introduced, including the so-called Fully Optimized Second Order (FOSO) method (Ponte Castañeda,

2015, 2016). Applications for model two- and three-dimensional viscoplastic composites were given by Furer and Ponte Castañeda (2018) and Kammer and Ponte Castañeda (2024), respectively, while results for viscoplastic polycrystalline materials have been given by Song and Ponte Castañeda (2018) with a first application in Castelnau et al. (2020).

Although very efficient from a numerical point of view, the coupling of advanced nonlinear homogenization scheme such as the FOSO method with a large scale FE structure model remains a challenge. However, at present, machine learning approaches offer a gateway for circumventing the just referred difficulties, by training a surrogate with rich synthetic data, which can then be used online under the stringent real-time constraints. Thus, thanks to these incipient data-driven technologies, one could consider the finest high-fidelity descriptions without impacting the computational complexity, even when addressing very large structures of industrial interest. This paper aims at showing the capabilities of machine learning approaches to construct a digital twin of accurate nonlinear methods.

The article is divided into four main sections including this introduction. As an initial step in Section 2, we aim to develop a digital twin capable of emulating a 2D nonlinear homogenization law. While the immediate computational gains may not be substantial, this endeavor serves as a testament to the potential of machine learning in this domain. Building upon this initial success, our objective in Section 3 is to extend this methodology to tackle a broader and more intricate homogenization challenge, that of a polycrystalline aggregate. The case of ice is used for illustrative purposes. In this latter scenario, the potential for significant computational time savings becomes apparent, underscoring the transformative impact of machine learning in accelerating complex computational tasks as highlighted in the conclusion (Section 4).

## 2. Digital twin of a 2D nonlinear homogenization law

### 2.1. Considered case

The FOSO nonlinear homogenization method by Ponte Castañeda (2016) offers a sophisticated means of accurately representing the macroscopic constitutive relations and field statistics within two-phase composite materials. This method provides a superior characterization of the physical properties of composites compared to the earlier variational linear comparison method (VLC) (Ponte Castañeda, 1991). However, this heightened accuracy is accompanied by an increased computational burden. Consequently, in this section, we explore the potential of machine learning to achieve comparable precision to that of Ponte Castañeda (2016), while benefiting from the computational advantages afforded by simplified approach of Ponte Castañeda (1991).

For simplicity, we first consider two-phase composites with isotropic viscoplastic phases and two-dimensional microstructures consisting of random (isotropic) distributions of aligned circular fibers, which are subjected to plane strain loading conditions (in the transverse plane). Both phases are power-law viscoplastic with nonlinearity exponent $n = \frac{1}{m}$ (where $0 \le m \le 1$ is the strain-rate sensitivity), and flow stresses $\sigma_{0}^{(1)}$ for the matrix phase and $\sigma_{0}^{(2)}$ for the fiber phase. The volume fraction of the fibers is denoted by $c$. The application of the VLC and FOSO nonlinear homogenization methods to this class of composites was considered in Furer and Ponte Castañeda (2018) making use of the linear homogenization estimates of Ponte Castañeda and Willis (1995) to estimate the effective properties of the LCC. The reader is referred to Furer and Ponte Castañeda (2018) for details, but the predictions of the VLC and FOSO methods can be expressed in terms of the effective or macroscopic flow stress $\tilde{\sigma}_{0}$ of the composite in the form

$$ \frac{\tilde{\sigma}_{0}}{\sigma_{0}^{(1)}} = f \left( n, c, \frac{\sigma_{0}^{(2)}}{\sigma_{0}^{(1)}} \right), \quad (1) $$
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Comparison between the VLC nonlinear homogenization method (Ponte Castañeda, 1991) and FOSO nonlinear homogenization method (Furer and Ponte Castañeda, 2018), as function of the nonlinearity for a given material state.</figcaption></figure>

where *f* represents a given function. It should be emphasized that *m* = 1 corresponds to linearly viscous response and that both nonlinear homogenization methods recover the linear result of Ponte Castañeda and Willis (1995) (which was used as input). At the other extreme, the value *m* = 0 corresponds to the strongly nonlinear limit of ideally plastic behavior. For values of 0 ≤ *m* < 1, the VLC estimate produces a rigorous upper bound for the corresponding FOSO estimate for the homogenized flow stress $\bar{\sigma}_0$ of the composite. It should also be noted that when the heterogeneity contrast parameter $\frac{\sigma_0^{(2)}}{\sigma_0^{(1)}}$ is equal to 1, the composite becomes a homogeneous material with effective property $\bar{\sigma}_0 = \sigma_0^{(1)} = \sigma_0^{(2)}$.

To illustrate the results, consider a straightforward example where $c = 0.20626$, $\sigma_0^{(1)} = 10$ and $\sigma_0^{(2)} = 2$. Fig. 1 juxtaposes the results obtained from the two methods, emphasizing the escalating disparity with decreasing *m*.

## 2.2. Elaborating a digital twin

In this section, a machine learning-based digital twin is developed to provide predictive solutions for Eq. (1). As a demonstrative first step, the only varying parameter from Eq. (1) is considered to be the nonlinearity parameter *n*. Appendix shows that accurate predictions can be achieved even with a simple neural network architecture, provided that the hyperparameters are carefully and appropriately selected.

However, for such a straightforward case where the only varying parameter is the nonlinearity parameter *n*, a polynomial regression could also replicate the results effectively. A more challenging situation arises when all three parameters of Eq. (1) are varied. To address this challenge, we generated data for each variable within a pertinent data range and leverage this dataset for training and testing a machine learning algorithm. In this context, the parameters were operated within the following ranges: *c* in [0.2, 0.7], $\frac{\sigma_0^{(2)}}{\sigma_0^{(1)}}$ in [0.2, 2], and *m* in [10⁻⁶, 1]. Note that *m* could not be set to zero, but we approached zero by setting a minimum value of 10⁻⁶. It is expected that, as we approach the lower bound of *m*, convergence of the FOSO nonlinear homogenization method becomes more challenging, leading to an increase in the computational cost.

Over this parameter range, generating as much data as needed is possible. However, considering only three parameters are involved, 10 values per parameter would sufficiently cover the data range, resulting in a dataset of 1000 points. Therefore, the training data set comprises 1000 data points from which 200 points were used as a validation data

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Neural Network architecture showing the number of layers, the number of neurons in each layer and the activation function used.</figcaption></figure>

set. However, for testing purposes, we generated 100,000 data points within the same range to evaluate the precision of the method over this significantly larger dataset.

For constructing the neural network, several architectures and activation functions were tested. Some activation functions deteriorate the convergence, such as, for instance, the ReLu activation function of Krizhevsky et al. (2012) and the linear activation function. On the other hand, the hyperbolic tangent (tanh) activation function of Lecun et al. (1998) leads to reduced convergence rate and lower precision level. The Softplus activation function (Dugas et al., 2000) proved to be the most suitable for this problem. One explanation could be that the FOSO nonlinear homogenization function as depicted in Fig. 1 resembles a logarithmic function which holds some similarities to the Softplus function which is defined as $\log(1+e^*)$. Fig. 2 shows the neural network architecture. We utilized 10 layers generating 853 parameters. The network architecture is made of:

* An input layer with 3 neurons corresponding to the three parameters of Eq. (1).

* Nine hidden layers with varying number of neurons as indicated in Fig. 2, all equipped with Softplus activation functions.

* An output layer corresponding to the flow stress ratio $\frac{\bar{\sigma}_0}{\sigma_0^{(1)}}$ with only one neuron and a linear activation function.

The learning rate was monitored to decrease over time, and convergence was achieved in 3500 epochs, resulting in a minimum least square loss function value as low as 10⁻⁷. Error on all training points was negligible. Subsequently, we examined the accuracy attained on the testing dataset.

To illustrate the methodology, we initially constructed the solution for $c = 0.2026$ and $\frac{\sigma_0^{(2)}}{\sigma_0^{(1)}} = 0.2$ across 100 *m* values. These combinations had not been previously encountered by the network. The results are shown in Fig. 3, demonstrating that the predictions closely match FOSO nonlinear homogenization with high fidelity.

The same conclusion holds when examining Figs. 4 and 5. These figures explore variations in the volume fraction *c* and the heterogeneity contrast parameter $\frac{\sigma_0^{(2)}}{\sigma_0^{(1)}}$, respectively, when the nonlinearity coefficient *m* is set to 0.15.

Plotting the output of all 100,000 test cases in one figure yields Fig. 6. Here, the *x*-axis represents the FOSO nonlinear homogenization output $\frac{\bar{\sigma}_0}{\sigma_0^{(1)}}$, while the *y*-axis depicts the predicted values. Hence, the line *y* = *x* corresponds to the zero error data line. Upon inspection of Fig. 6, it is evident that all predictions lie on the first diagonal, with differences imperceptible to the naked eye unless the results are zoomed in.

To quantify the error, let us denote the output $\frac{\bar{\sigma}_0}{\sigma_0^{(1)}}$ by *y*. Thus, the error can be calculated as follows:
<figure><img src="image_5.png" /><figcaption>Fig. 3. Comparison of linear and nonlinear homogenization outputs with predicted values for a given material problem.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Comparison between the VLC nonlinear homogenization method (Ponte Castañeda, 1991), the FOSO nonlinear homogenization method (Furer and Ponte Castañeda, 2018) and the machine learning prediction as function of the volume fraction c for a nonlinearity coefficient m = 0.15.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Comparison between the VLC nonlinear homogenization method (Ponte Castañeda, 1991), the FOSO nonlinear homogenization method (Furer and Ponte Castañeda, 2018) and the machine learning prediction as function of the heterogeneity contrast parameter $\sigma_0^{(2)}/\sigma_0^{(1)}$ for a nonlinearity coefficient m = 0.15.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(a) Prediction versus analytical solution for all 100,000 test cases</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) Zoom on the indicated area of Figure 6a</figcaption></figure>

<caption>Fig. 6. Comparison of neural network predictions and nonlinear homogenization.</caption>

$$ \Delta y = \frac{|y_{\text{predicted}} - y_{\text{exact}}|}{|y_{\text{exact}}|} \times 100 $$

Here, $\Delta y$ represents the percentage error for each prediction value. After computing this value for the entire testing dataset, we find that the maximum error is 0.5%, as it would be expected when examining Fig. 6.

## 2.3. Additional use of the VLC bound

Despite machine learning's demonstrated capability in predicting nonlinear homogenized solutions, one may question whether this prediction could be further refined by incorporating not only the inputs but also the outcomes derived from VLC nonlinear homogenization. Essentially, we propose utilizing the VLC nonlinear homogenization method as a precursor to furnish a baseline from which the neural network can start its prediction of the FOSO nonlinear homogenized solution as done by Ghnatios et al. (2024). This approach aims not only to potentially enhance accuracy but also to expedite convergence. To investigate this route, we trained the neural network using identical architecture and dataset, augmenting the inputs with results from VLC nonlinear homogenization. Our analysis, limited to the first 100 epochs,
<figure><img src="image_10.png" /><figcaption>Fig. 7. Comparison of neural network performance with and without incorporating results from VLC nonlinear homogenization.</figcaption></figure>

yielded the results depicted in Fig. 7. Notably, by comparing Fig. 7 with Fig. A.23 we deduce that leveraging VLC nonlinear homogenization results facilitated convergence towards the values predicted by FOSO nonlinear homogenization. However, it is noteworthy that achieving comparable results without using the VLC results as input, and through further training necessitated over 1000 additional epochs. In addition, within just 200 epochs, the neural network incorporating VLC results was able to reach the same level of accuracy that required 3500 epochs on the network that is not equipped with VLC results. This dramatic reduction in training time underscores the efficiency of incorporating VLC outcomes, which significantly accelerates convergence. It is important to highlight that the 200 epochs required only around 3 min on a system powered by a 13th Gen Intel(R) Core(TM) i9-13950HX processor clocked at 2.20 GHz and equipped with 32 GB of memory, while the 3500 epochs took more than 53 min. Thus, the training time is substantially reduced, which presents a considerable advantage for offline training. However, it should be noted that no significant gains in online computational costs were observed—though importantly, no additional overhead was introduced during the online stage either.

## 2.4. Conclusions

It is evident that machine learning has effectively replicated the results obtained through FOSO nonlinear homogenization. This accomplishment brings forth two primary advantages. Firstly, FOSO nonlinear homogenization necessitates the use of an iterative solver like Newton-Raphson to compute the output. The convergence depends heavily on the initial guess provided and the numerical parameters. Even with a well-informed initial guess, the method struggles to converge as m approaches zero and for certain material combinations. In contrast, machine learning predictions are straightforward and offer high accuracy without relying on any iterative solver. The latter is only required for generating the data needed for training, a process that occurs once. Consequently, convergence issues are circumvented once the neural network has learned the material model. Secondly, while producing the 100,000 data points required approximately 7 s with the VLC nonlinear homogenization method and over two hours when employing the FOSO nonlinear homogenization technique, less than a second is sufficient to predict all the values using the trained network. Thus, a significant gain in computation time is achieved.

Given this success, our next endeavor in the following section is to extend the method to encompass a more complex material behavior law in 3D, specifically that of ice polycrystals. In this scenario, challenges emerge from the multitude of involved inputs and outputs, as well as

the nonlinearity of the underlying equations. However, the potential benefits are substantial, particularly in terms of the tremendous gains in computation time.

## 3. Digital twin of a polycrystalline aggregate: the case of ice

The behavior of viscoplastic polycrystalline materials is governed by the glide of lattice defects (dislocations) on specific crystallographic planes and along specific directions, in each grain, in response to the local stress state. Dislocation glide is known to be a highly nonlinear mechanism. Ice polycrystals are used here as an illustrative example of the capabilities of machine learning. Due to the crystallographic structure of the crystal (hexagonal lattice), only a limited number of dislocation glide planes and directions are possible, leading to a highly anisotropic constitutive relation at the grains scale (Schulson and Duval, 2009) and strong stress and strain concentration in specific grains. The microstructure of polycrystalline aggregates can be characterized by the morphologic texture (average grain shape) and the crystallographic texture (orientation distribution of the crystal lattice of the constituting grains). In polar ices, such as in Greenland and Antarctica, the ice crystals are randomly oriented (random texture) at the surface but preferred lattice orientations develop with increasing depth, due to the very large shear strain the ice is subjected to Castelnau et al. (1998), but also due to recrystallization processes (Chapelle et al., 1998). A consequence of both high grain anisotropy and strong crystallographic texture is that the effective constitutive relation of polar ices is strongly anisotropic. There have been attempts to account for this anisotropy in large scale flow models in order, for example, to better predict the mechanical response and in situ ice flow instabilities in Greenland and Antarctica that are expected with the climate warming. However, a state-of-the-art homogenization scheme such as the FOSO method is needed for such a complex material, while the coupling of the homogenization scheme with large scale flow model would require considerable computational power. In this context, machine learning could be leveraged to predict the complex ice behavior based on ice characteristics and relevant inputs. This could save the computational cost while preserving the accuracy of nonlinear homogenization techniques.

### 3.1. Material law parameters and outputs

In order to achieve this objective, we shall start by describing the ice structure and the different inputs that feed into the constitutive relation at the grain scale.

Fig. 8 illustrates the dislocation glide planes and directions expected in ice single crystals. Slip on the basal plane family constitutes the main deformation mechanism and has been quite well characterized experimentally. Slip on prismatic and pyramidal planes should be much harder to activate and is often considered as possible secondary deformation mechanism in the literature. The shear-rate $\dot{\gamma}$ on those planes is a power-law function of the resolved shear stress $\tau$ and could be characterized by:

$$ \frac{\dot{\gamma}}{\dot{\gamma}_0} = \left(\frac{\tau}{\tau_0}\right)^n \qquad (2) $$

where $\dot{\gamma}_0$, $\tau_0$ and $n$ are the parameters of the constitutive relation, respectively a reference shear-rate, the critical resolved shear stress, and the stress sensitivity. Here, we choose $\dot{\gamma}_0 = 10^{-6} s^{-1}$ (typical strain-rate during laboratory experiment) for all slip systems. Data for the basal glide are quite well constrained by experimental results (Schulson and Duval, 2009): $n_{ba} = 2.0$, $\tau_{ba0} = 0.035$ MPa. Rheological data ($n_{pr}$, $n_{py}$, $\tau_{pr0}$ and $\tau_{py0}$) for prismatic and pyramidal slip are much less constrained, and therefore we will explore how machine learning can reproduce the effective behavior of polycrystalline aggregate whatever their value. The ratio $\tau_{pr0}/\tau_{ba0}$ or $\tau_{py0}/\tau_{ba0}$, which are both much larger than 1, determine the viscoplastic anisotropy of ice crystals.
<figure><img src="image_11.png" /><figcaption>Fig. 8. Possible glide planes and directions for ice single crystals.</figcaption></figure>

In this work, we also take into account the crystallographic texture of ice polycrystals in order to infer the capability of the neural network to reproduce the effective behavior of strongly anisotropic materials. We limit the study to transverse isotropic textures, as mostly encountered in ice sheets, and we use a parametrization similar to the one used in Maurel et al. (2015) with a single parameter θ describing the orientation distribution of crystal c-axes in space. The obtained textures are illustrated in Fig. 9. θ describes the degree of anisotropy of the polycrystal and ranges from -90° to 90°.

These five material parameters, together with the prescribed effective stress tensor, can be used to compute the effective strain-rate tensor according to Song and Ponte Castañeda (2018). The effective stress tensor $\bar{\sigma}$ is symmetric, and as the viscoplastic behavior guided by (2) does not depend on the first stress invariant, one can limit the stress inputs to five. The effective strain-rate tensor $\dot{\epsilon}$ has similar properties (isochoric deformations) leading, therefore, to five output values. Here, we use the reduced notation introduced by Lequeu et al. (1987) for second rank tensors:

$$ \sigma_1 = \frac{\sigma_{22} - \sigma_{11}}{\sqrt{2}}, \sigma_2 = \sqrt{\frac{3}{2}}\sigma_{33}, \sigma_3 = \sqrt{2}\sigma_{23}, \sigma_4 = \sqrt{2}\sigma_{13}, \sigma_5 = \sqrt{2}\sigma_{12} \quad (3) $$

(and similarly for the strain-rate tensor). Altogether, the anisotropic effective strain-rate of ice polycrystals can be written as

$$ \dot{\epsilon} = g(\bar{\sigma}, \theta, n_{pr}, n_{py}, \tau_{pr0}, \tau_{py0}), \quad (4) $$

where g designates a given function.

## 3.2. Data generation

The ice polycrystals nonlinear homogenization described in Song and Ponte Castañeda (2018) is solved using an iterative approach in order to get the effective strain-rate tensor $\dot{\epsilon}$ as described in Eq. (4). This computational step is facilitated by a FORTRAN executable, enhancing efficiency. However, individual computations typically require between 1.5 and 5 s on average. Notably, convergence challenges arise when approaching a cone angle θ of approximately 90°, resulting in a significant escalation in computational demands. Nevertheless, angles approaching 90° hold little physical significance, leading us to confine our focus to cone angles within the range of [-80°, 80°]. Specifically, emphasis is placed on the [-30°, 30°] interval due to its heightened relevance.

The non-linearity parameters $n_{pr}$ and $n_{py}$ typically fall within the range of [1, 5]. Similarly, values of $\tau_{pr0}$ and $\tau_{py0}$ are of interest within intervals of [0.2 MPa, 5 MPa] and [0.5 MPa, 50 MPa], respectively, as indicated by prior ice simulation studies (Castelnau et al., 2008). Analysis of the effective stress tensor $\bar{\sigma}$ reveals that, while individual stress components generally fluctuate within the [-2 MPa, 2 MPa] interval, they more frequently reside within a narrower range of [-0.2 MPa, 0.2 MPa]. Accordingly, we aim to augment data coverage within the [-2 MPa, 2 MPa] range, with a particular focus on the [-0.2 MPa, 0.2 MPa] subset.

Upon consideration of parameter ranges and associated computational costs, it becomes evident that generating a comprehensive dataset is impractical. With 10 input variables necessitating at least 10 values each to adequately cover their respective intervals, the resulting requirement for 10 billion simulations would entail a simulation time exceeding a century. Such a daunting prospect underscores the necessity for strategic data selection and increases the learning challenges for the machine learning techniques as discussed in Section 3.3.

The dataset, aimed at maintaining reasonable computational costs, should consist of approximately one hundred thousand data points. Our objective is to ensure comprehensive coverage across the entire parameter range to minimize the need for extrapolation when utilizing the machine learning model. To achieve optimal data selection, we adopt the Latin Hypercube selection method (LHS) of McKay et al. (1979), considering data points at interval limits and supplementing them with others within the intervals. We execute this process twice: once for values within the narrow interval and again for those within the broader range. On each occasion, we include boundary data points alongside fifty thousand additional points distributed using LHS. As a result, we attain a total of 102 048 data points, facilitating complete coverage of the parameter space with a focus on the narrower interval.

The data generation process required approximately 5 days using a 13th Gen Intel(R) Core(TM) i9-13950HX processor clocked at 2.20 GHz and equipped with 32 GB of memory.

## 3.3. Training a machine learning model

The collected data was stored in a CSV file for subsequent training purposes. Prior to initiating the training process, data normalization was imperative due to variations in input ranges and potential differences in output averages among strain rate components.

While various normalization techniques exist, the chosen approach ensures that both inputs and outputs are normalized within the range [0, 1]. This not only standardizes the data but also ensures its positivity, allowing for experimentation with different activation functions, some of which exclusively output positive values. For a given input or output value x, normalization is performed as follows:

$$ \|x\| = \frac{x - x_{min}}{x_{max} - x_{min}} \quad (5) $$

where $x_{min} = \min(x)$ and $x_{max} = \max(x)$ are computed over the entire dataset pertaining to x.

The dataset was divided into a training data set (including the validation data set) and a test data set, with 92% of the data dedicated to training, amounting to over 8000 data points for testing. This allocation provided a strong foundation for evaluating the network's ability to generalize to unseen data. The decision to use a 92-8 split, rather than the more common 80-20 split, was driven by two key factors. First, the large size of the data space means that even the 102 048 data points are not fully representative, and a larger training set is required to ensure effective learning and accurate predictions. Additionally, the relative scarcity of data helps mitigate overfitting. Nevertheless, 8% of the data was also set aside for validation, further reducing the amount available for training. Second, allocating 8% for testing ensures that this portion remains sufficiently representative of the vast parameter domain, allowing for a more reliable assessment of the model's predictive performance.

During the training phase, several machine learning techniques and architectures were explored, with varying degrees of success. Notably, while a single neural network exhibited the potential to predict all five outputs, it failed to achieve the desired level of precision. Consequently, five separate neural networks were trained, each tasked with predicting a specific component of the strain rate tensor. Despite their distinct targets, all networks shared identical architectures, activation functions, and characteristics. The neural network architecture that was used is the one illustrated in Fig. 10. It comprises 8 layers with a softplus activation function for each layer except the last one which has a linear
<figure><img src="image_12.png" /><figcaption>Fig. 9. c-axes pole figures of ice polycrystals. The opening angle $\theta$ is depicted by the grey arrows and indicated below the pole figures. The value $\theta = 0^{\circ}$ corresponds to randomly oriented c-axes (isotropic polycrystal), a positive $\theta$ value indicates that c-axes tend to cluster around the 3-axis, while a negative value indicates that c-axes tend to align parallel to the 1-2 plane. The RVE is constituted by 1000 grains.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. The neural network architecture used for each of the five output values of the strain rate tensor.</figcaption></figure>

activation function on its unique neuron. The loss function is a least square loss function. The learning rate starts at $10^{-3}$ and is reduced as a function of the number of epochs achieved. It reaches $5 \times 10^{-4}$ after 30,000 epochs. In fact, experiments showed that reducing the learning rate increased the convergence rate and contributed to a lower loss function.

Furthermore, the significance of network size became evident, with the optimal architecture illustrated in Fig. 10 generating a total of 913 parameters. Smaller networks struggled to achieve sufficient precision, while larger networks, surpassing 1800 parameters, initially converged quickly but plateaued at a suboptimal precision level.

The optimal network architecture of Fig. 10 consists of:

* An input layer with ten neurons corresponding to the parameters of Eqs. (3) and (4).

* Seven hidden layers with varying number of neurons as indicated on Fig. 10, all equipped with Softplus activation functions.

* An output layer corresponding to one of the five output values $\dot{\tilde{\epsilon}}_i$ ($i \in [1, 5]$) of the strain rate tensor with only one neuron and a linear activation function.

Moreover, as seen in Fig. 10, the neural network architecture features varying numbers of neurons across layers. The first layer contains 10 neurons, which aligns with the input layer, also of size 10. Since the network is designed to produce a single output, it is logical to progressively reduce the number of neurons in subsequent layers, leading to the final three layers having 10, 5, and 1 neuron, respectively. This gradual reduction facilitates a smooth transition from complex representations in the deeper layers to the single output required.

The intermediate layers, positioned between the input and output layers, contain more neurons because they play a crucial role in capturing and processing the underlying complexity of the data. Increasing the number of neurons in these middle layers allows the network to learn more intricate patterns and relationships in the dataset, enhancing its ability to model the non-linearities present in the problem.

By first expanding the network's capacity and then narrowing it towards the output, the architecture strikes a balance between representation power and efficient dimensionality reduction.

This design strategy follows a common neural network principle: deeper layers with more neurons can better extract features and process complex data, while gradually reducing the network size helps control overfitting and improve generalization without losing important learned representations.

Among the activation functions tested, Softplus (Dugas et al., 2000) emerged as the most effective. In details, the authors tested Sigmoid, ReLu, SeLu, tanh respectively of Glorot and Bengio (2010), Krizhevsky et al. (2012), Klambauer et al. (2017) and Lecun et al. (1998) as well as linear activation functions. While certain functions impeded convergence or prolonged the process, others, such as the hyperbolic tangent (tanh), were deemed acceptable but failed to maximize accuracy. Interestingly, softplus also demonstrated proficiency in predicting 2D nonlinear homogenization as demonstrated in Section 2.

It is noteworthy to mention that Random Forest surprisingly achieved similar performance while exhibiting significantly faster training times. Furthermore, it demonstrated the capability to predict all five outputs simultaneously. Random forests excel with large datasets as they reduce the need for interpolation by effectively capturing complex patterns. However, they are unsuitable for Finite Element Method (FEM) applications due to their non-smooth predictions and discontinuous second derivatives, which conflict with the smoothness requirements in FEM, particularly because FEM involves solving partial differential equations (PDEs) that demand smoothness and continuity of derivatives. In contrast, other machine learning techniques such as XG-Boost or Support Vector Machines struggled to learn from this dataset. This is why neural networks are mainly used in this article despite exhibiting larger training times compared to above mentioned techniques.

### 3.4. Analyzing the precision of machine learning predictions

In the subsequent analysis, we juxtapose the machine learning predictions against reference outputs derived from solving the nonlinear problem. Fig. 11 illustrates the outcomes pertaining to the first tensor value, $\dot{\tilde{\epsilon}}_1$. These values are obtained post-normalization reversal of Eq. (5), as represented by the formula:

$$ \dot{\tilde{\epsilon}}_1 = \|\dot{\tilde{\epsilon}}_1\| \times (\dot{\tilde{\epsilon}}_{1max} - \dot{\tilde{\epsilon}}_{1min}) + \dot{\tilde{\epsilon}}_{1min} $$

In Fig. 11, we plot the predicted values against the reference ones. Consequently, perfect alignment between prediction and reference would result in all data points lying along the $y = x$ axis, highlighted in black, indicating congruence with the reference values in our dataset. Additionally, the figure illustrates the predicted values corresponding to the training dataset in green and those associated with the test dataset in red. Notably, the precision achieved on the test dataset mirrors that of the training dataset, affirming the network's adeptness in learning from the data without succumbing to overfitting. To further elucidate this observation, Fig. 12 provides a magnified view of the central values from Fig. 11, where the majority of data points reside.
<figure><img src="image_14.png" /><figcaption>Fig. 11. Comparison between predicted and reference values for the first component of the strain rate tensor, $\dot{\epsilon}_1$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Zoomed-in depiction of the central values from Fig. 11.</figcaption></figure>

Furthermore, as discussed previously, emphasis was placed on generating additional data within the narrow interval corresponding to lower cone angles and stress tensor element values. Nonetheless, it is evident that the precision of predictions across the entire dataset remains consistent. The neural network demonstrates robustness in producing uniform precision across the entire range of data intervals, effectively capturing and generalizing the complexity inherent in the ice polycrystal law.

The precision achieved for $\dot{\epsilon}_1$ extends to the other components of the strain rate tensor. Figs. 13 and 14 illustrate the results for the fifth component, $\dot{\epsilon}_5$. Similar to $\dot{\epsilon}_1$, the quality of the results for $\dot{\epsilon}_5$ is apparent. However, a subtle distinction arises in this component: accuracy on certain interval limit values is marginally less pronounced. To quantify the error across all output values consistently, let us proceed with the analysis.

Both neural network training and random forest regression achieved low loss values and minimal percentage errors, whether on the training or test set. Plotting the reference versus predicted values for each case and assessing the relative error between the two yields Fig. 15 for $\dot{\epsilon}_1$ on the training dataset. Here, relative errors fall below 4% for the majority of cases, with only 29 instances exceeding an 8% error out of over 90 000 cases.

Examination of the test dataset reveals even higher precision than that achieved during training, as depicted in Fig. 16. Here, the relative error remains below 1.5% for most cases, with only 5 instances among more than 8000 exceeding an 8% error. This reaffirms our previous assertion regarding the neural network's ability to learn and generalize the polycrystalline material law.

This analysis extends to every output of the strain rate tensor, yielding similar results. For $\dot{\epsilon}_5$, although slightly lower precision was observed, Figs. 17 and 18 depict relative errors on the training and
<figure><img src="image_16.png" /><figcaption>Fig. 13. Comparison of predicted and reference values for the fifth component of the strain rate tensor, $\dot{\epsilon}_5$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. Zoomed-in view of the central values from Fig. 13.</figcaption></figure>

test datasets, respectively. While 36 cases exceeded an 8% error on the training dataset for $\dot{\epsilon}_5$, the majority exhibited errors of less than 4%. Notably, the few instances exceeding an 8% error display higher relative errors compared to other output variables. Despite this, these outliers remain exceptional considering the extensive dataset. A similar conclusion holds for the test dataset, where precision remains commendable.

Additionally, Fig. 20 compares reference values obtained for output $\dot{\epsilon}_2$, against predicted values, showcasing consistent results across both training and test datasets. Restricting this presentation to 1000 data points allows for a clearer visualization of the results.

One can also check the particular case of isotropic ice polycrystals that are described by a cone angle $\theta = 0$ and for which $n_{pr} = n_{py} = \frac{1}{m}$. Fig. 19 shows a comparison of the FOSO nonlinear homogenization results with the machine learning prediction for the Von-Mises

equivalent strain rate value $\dot{\epsilon}_{VM} = \sqrt{\frac{2}{3} \times \sum_{i=1}^{5} (\dot{\epsilon}_i)^2}$ for $m$ varying between it is lower and upper limits. The different components of the stress tensor are set to 0.2 in this particular example. The material behavior is comparable to the power-law material behavior discussed in Section 2. The digital twin is able to predict this particular behavior with satisfactory precision.

These findings affirm the neural network's capability to capture the intricacies of the polycrystal effective behavior. Across the entire dataset, both observed and unseen by the network, high-error cases are conspicuously absent. This level of precision suffices for applications such as modeling the ice flow dynamic in natural ice masses.

Before concluding this section, it would be insightful to delve into the impact of data size on the quality of results. With a dataset comprising more than 100,000 data points at our disposal, we opted to conduct an investigation using a random forest regressor. The aim was
<figure><img src="image_18.png" /><figcaption>Fig. 15. Relative error on the training dataset for $\hat{\epsilon}_1$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16. Relative error on the test dataset for $\hat{\epsilon}_2$.</figcaption></figure>

to discern how the model's effectiveness varies with different sizes of training data. We initiated this exploration by randomly selecting 10,000 data points for training and evaluating its performance on the remaining data. Subsequently, we incremented the training dataset size by 10,000 points until it reached 90,000.

Fig. 21 illustrates the progression of mean squared error (MSE) on both the training and test datasets. It becomes evident that a dataset comprising solely 10,000 data points is insufficient for effective model training. However, as the dataset size increases, the error diminishes rapidly. Remarkably, beyond 80,000 data points, the MSE value on the training set closely aligns with that of the test set, exhibiting sluggish evolution thereafter. Hence, it appears that the dataset size utilized for training is optimal. While employing a larger dataset might marginally enhance accuracy, the associated computational effort to

generate additional data points might outweigh the potential gains, rendering it not particularly advantageous.

### 3.5. Conclusions

The work conducted in this section serves as a compelling demonstration of the power of machine learning in replicating intricate effective material behaviors. This achievement underscores similar findings to those elucidated in the 2D case discussed in Section 2.4. However, the complexity inherent in polycrystalline aggregates, encompassing five distinct outputs along with a multitude of parameters and inputs, presents a formidable challenge. Despite these intricacies, machine learning offers a distinct advantage due to its capacity to yield precise results with minimal computational overhead. In fact, machine learning
<figure><img src="image_20.png" /><figcaption>Fig. 17. Relative error on the training dataset for $\hat{\epsilon}_5$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 18. Relative error on the test dataset for $\hat{\epsilon}_5$.</figcaption></figure>

is able to accurately generate in less than a second what took 5 entire days to compute using FOSO method as mentioned in Section 3.2.

Moreover, as presented in Section 2.3, it is possible to train a machine learning algorithm using, in addition to the aforementioned inputs, outputs from a simpler estimation technique to learn the complex behavior of polycrystalline aggregates. This approach, which involves using machine learning to capture the difference between a simple estimation technique and a more sophisticated, computationally demanding one, may slightly enhance precision. However, this technique was not investigated in this section due to the fact that the computation of the VLC estimate is almost as demanding as that of the FOSO for polycrystals, and it is more expedient to use the FOSO estimate directly in this case. The potential gain in accuracy does not justify the additional computational expense.

The method presented in this section maximizes computational efficiency while preserving high accuracy. The reduction in computational

power requirements enables the handling of challenging physical phenomena, such as the creep flow observed in polar ice glaciers at physically representative scale, where traditional Finite Element Method (FEM) simulations would incur prohibitive computational costs. The next step along these lines would be to use machine learning for predicting microstructure evolution during polycrystal deformation. This point is left for future work.

## 4. Conclusions and perspectives

This article introduces a digital twin capable of accurately simulating the intricate behavior of heterogeneous materials. Initially, a neural network was employed to predict the behavior of nonlinear composite materials in two dimensions. This approach demonstrated remarkable efficiency, achieving accuracy with respect to FOSO nonlinear homogenization method while imposing minimal computational burden. Encouraged by this success, the methodology was extended to
<figure><img src="image_22.png" /><figcaption>Fig. 19. Comparison of FOSO nonlinear homogenization technique with the machine learning prediction for the particular case of isotropic polycrystals with $n_{pr} = n_{py} = n$. Results shown correspond to the Von-Mises equivalent strain rate value $\dot{\epsilon}_{VM}$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 20. Results obtained for $\dot{\epsilon}_2$ using 1000 data points.</figcaption></figure>

model the nonlinear homogenization of ice polycrystals in three dimensions with various microstructures and slip system resistances. Given the increased complexity and computational demands of 3D behavior, machine learning emerged as a valuable tool in crafting a faithful representation of the material's behavior. Notably, the predictive precision remained high while computational costs were substantially reduced. Furthermore, the neural network exhibited robust learning capabilities, enabling accurate predictions even with limited data—a critical advantage considering the large variable space.

However, the primary focus of this article has been on the constitutive behavior of materials, without any effort so far to integrate this behavior into a Finite Element Method (FEM) framework. Extending this work to an FEM context is one of the key motivations behind this

research. The main challenge arises from the fact that the material laws discussed here are formulated to compute strain rate based on known stress. However, conventional FEM solvers operate by deriving stress from strain or strain rate. Consequently, an essential step in a FEM application is the inversion of the constitutive model. Two important perspectives emerge from this challenge:

*   **Integration with FEM Solvers:** Testing the proposed methodology within an FEM environment is not straightforward and will require significant development effort. Specifically, it necessitates the creation of an adaptive FEM solver capable of handling these material laws. To address this, the authors are initiating a dedicated project aimed at implementing a complete ice polycrystal material law in an FEM context. This project will also offer the
<figure><img src="image_24.png" /><figcaption>Fig. 21. Evolution of Mean Squared Error (MSE) as a function of training dataset size.</figcaption></figure>

opportunity to refine the material law for better accuracy in
simulating ice flow in Antarctica. A more sophisticated model,
tailored to such complex physical scenarios, will be used.

*   **Inversion Stability and Convexity:** The inversion of the material laws described in this paper, such as the FOSO or VLC methodologies, is facilitated by the convexity of the mathematical expressions for the associated potentials. This convexity guarantees the uniqueness of the solution, thus ensuring the stability of the inversion process. Additionally, the machine learning (ML) approach presented in this work, within the studied parameter domain, is expected to perform well without causing divergence or computational instability. This is due to the close alignment of its predictions with those of the FOSO methodology. However, these assumptions will need to be validated once the methodology is fully implemented in an FEM solver. Should any issues arise, the dissipation potential $\overline{W}(\epsilon)$ could be employed, or the “potential” nature of the constitutive relations could be enforced to preserve solution uniqueness.

Nevertheless, care must be taken when extending this approach
beyond its training domain. Purely data-driven models, while ac-
curate within known data ranges, are unreliable for extrapolation.
Even small deviations from the training set can lead to significant
accuracy drops, potentially resulting in FEM solution divergence.
To mitigate this risk, it is advisable to cover the entire range of pa-
rameters during model training. In cases where out-of-distribution
data is encountered, fallback to the FOSO methodology is recom-
mended, despite the increased computational cost.

CRediT authorship contribution statement

**Mikhael Tannous:** Writing – original draft, Software, Methodology. **Chady Ghnatios:** Writing – review & editing, Validation. **Olivier Castelnau:** Writing – review & editing, Software. **Pedro Ponte Castañeda:** Writing – review & editing, Validation, Conceptualization. **Francisco Chinesta:** Writing – review & editing, Supervision, Methodology.

Declaration of competing interest

The authors declare that they have no known competing finan-
cial interests or personal relationships that could have appeared to
influence the work reported in this paper.

**Appendix. Machine learning setup and hyper parameters**

The primary objective of this appendix is to guide the reader
through the process of constructing and training a straightforward
neural network to predict the solution of Eq. (1). Additionally, it
highlights key technical choices made throughout this article, providing
insight into the methodology employed.

The simple neural network illustrated in Fig. A.22 is comprised of:

* One neuron in the input layer, corresponding to the nonlinearity parameter *m*

* A hidden layer with 32 neurons and employing a hyperbolic tangent activation function

* An output layer corresponding to the flow stress ratio $\frac{\sigma_0}{\sigma_0^{(1)}}$ with only one neuron and a hyperbolic tangent activation function

Moreover, the following setup and definitions are consistently ap-
plied not only in this simple example but throughout all applications
in this study. They have proven to be effective across various scenarios.

*   **Data Splitting:** The generated data was divided into three main subsets: training, validation, and testing. The specific percentages varied depending on the example. In cases with abundant data, at least 20% was allocated for testing, while the remaining data was used for training, with a minimum of 10% reserved for validation. The validation set was used to fine-tune the neural network and determine convergence. As this data has been indirectly seen by the network, it is treated as part of the training set when generating performance curves. In contrast, the test data, which was never seen by the network, was used solely to assess prediction accuracy.

*   **Loss Function:** Across all neural network architectures in this study, the minimal least squares error was used as the loss function.

*   **Convergence:** Convergence in this study was not defined by a fixed number of epochs. Instead, training continued until the loss dropped below a predefined threshold, typically $10^{-7}$ unless otherwise specified. This threshold had to be met on both the training and validation sets to prevent overfitting. Slow convergence was noted when the loss function decreased gradually, despite an increase in the number of epochs. Moreover, convergence was not announced and the training would also continue if the mean absolute error is above $10^{-5}$.

*   **Learning Rate:** For simpler cases, as illustrated in Fig. A.22, a constant learning rate was sufficient for achieving convergence. However, as complexity increased, it became necessary to reduce the learning rate as training progressed. In such instances, the initial learning rate was set to $10^{-3}$ and decreased over time to as low as $10^{-4}$ or lower. Higher learning rates impeded convergence.

*   **Optimizer:** The Adam optimizer was employed in this study due to its efficiency and ability to handle large datasets and noisy gradients. Adam, which stands for Adaptive Moment Estimation, combines the advantages of both the AdaGrad (Adaptive Gradient) and RMSProp (Root Mean Square Propagation) algorithms, making it particularly well-suited for non-stationary problems and sparse gradients (Kingma and Ba, 2014).

*   **Multiple Initializations:** For the applications of Section 3, the architecture was tested using at least three different initializations of the network weights to ensure the robustness and consistency of the results across different random starting points. This is particularly important when convergence is slow, as multiple initializations help assess whether slow convergence is due to suboptimal weight initialization or whether further tuning of the architecture or training parameters is needed.
<figure><img src="image_25.png" /><figcaption>**Fig. A.22.** A simple neural network architecture that generates the predictions of Fig. A.23.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. A.23.** Prediction of the nonlinear behavior of the scenario shown in Fig. 1 using a simple neural network.</figcaption></figure>

This latter network succeeded in predicting the nonlinear homogenized composite behavior of the given simple case as depicted in Fig. A.23. Although this setup appears straightforward, one might question why and how this particular architecture was chosen over alternative designs, and how the prediction accuracy can be both measured and validated. In this instance, there is only one parameter involved, with a single layer considered, the number of neurons in this layer must be selected. This selection can follow a grid search approach, where multiple layer sizes are tested, ranging from 8 to 96 neurons in increments of 8. The goal is to select the smallest layer size capable of providing adequate accuracy. Data size is also a crucial factor. It was evident that 10 data points were insufficient to train even this simple architecture, and a larger dataset was necessary as the layer size increased. Consequently, a dataset consisting of 100 data points, equally distributed across the data range, was chosen. Of this, 20% was used as a validation set to monitor overfitting and track the least square loss function as a function of the number of epochs. Training was considered complete when the loss function dropped below $10^{-7}$ for both the training and validation sets. The smallest layer size that achieved the desired accuracy was 32 neurons. Moreover, Fig. A.23 depicts predictions made over a thousand data points uniformly distributed across the *m* range. Importantly, 900 of these data points were never seen by the network during training, yet they were accurately predicted, as demonstrated by the curve.

## Data availability

No data was used for the research described in the article.

## References

Blackman, D., Boyce, D., Castelnau, O., Dawson, P.R., Laske, G., 2017. Effects of crystal preferred orientation on upper-mantle flow near plate boundaries: rheologic feedbacks and seismic anisotropy. Geophys. J. Int. 210 (3), 1481-1493.

Castelnau, O., Derrien, K., Ritterbex, S., Carrez, P., Cordier, P., Moulinec, H., 2020. Multiscale modeling of the effective viscoplastic behavior of Mg2SiO4 wadsleyite: bridging atomic and polycrystal scales. *Comp Rend Mec.* 10–11, 827–846.

Castelnau, O., Duval, P., Montagnat, M., Brenner, R., 2008. Elastoviscoplastic micromechanical modelling of the transient creep of ice. J. Geophys. Res. 113, http://dx.doi.org/10.1029/2008JB005751.

Castelnau, O., Shoji, H., Mangeney, A., Milsch, H., Duval, P., Miyamoto, A., Kuwabara, K., Watanabe, O., 1998. Anisotropic behavior of GRIP ices and flow in central greenland. Earth Planet. Sci. Lett. 160, 1–13.

Chapelle, S.D.L., Castelnau, O., Lipenkov, V., Duval, P., 1998. Dynamic recrystallization and texture development in ice as revealed by the study of deep ice cores in antarctica and greenland. J. Geophys. Res. 103, 5091–5105.

Dai, M., Demirel, M.F., Liu, X., Liang, Y., Hu, J.-M., 2023. Graph neural network for predicting the effective properties of polycrystalline materials: A comprehensive analysis. Comput. Mater. Sci. 230, http://dx.doi.org/10.1016/j.commatsci.2023.112461.

Dai, W., Wang, H., Guan, Q., Li, D., Peng, Y., Tomé, C.N., 2021. Studying the micromechanical behaviors of a polycrystalline metal by artificial neural networks. Acta Mater. 214, 117006. http://dx.doi.org/10.1016/j.actamat.2021.117006, URL: https://www.sciencedirect.com/science/article/pii/S1359645421003864.

Dugas, C., Bengio, Y., Bélisle, F., Nadeau, C., Garcia, R., 2000. Incorporating second-order functional knowledge for better option pricing. In: Leen, T., Dietterich, T., Tresp, V. (Eds.), Advances in Neural Information Processing Systems. Vol. 13, MIT Press.

Feyel, F., Chaboche, J.-L., 2000. FE2 multiscale approach for modelling the elasto-viscoplastic behaviour of long fibre SiC/Ti composite materials. Comput. Methods Appl. Mech. Engrg. 183 (3), 309–330. http://dx.doi.org/10.1016/S0045-7825(99)00224-8.

Furer, J., Ponte Castañeda, P., 2018. A symmetric fully optimized second-order method for nonlinear homogenization. ZAMM - Z. Angew. Math. Mech. 98 (2), 222–254. http://dx.doi.org/10.1002/zamm.201700065.

Ghnatios, C., Rodriguez, S., Tomezyk, J., Dupuis, Y., Mouterde, J., Silva, J.D., Chinesta, F., 2024. A hybrid twin based on machine learning enhanced reduced order model for real-time simulation of magnetic bearings. Adv. Model. Simul. Eng. Sci. 11 (3), http://dx.doi.org/10.1186/s40323-024-00258-2.

Glorot, X., Bengio, Y., 2010. Understanding the difficulty of training deep feedforward neural networks. In: Teh, Y.W., Titterington, M. (Eds.), Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics. In: Proceedings of Machine Learning Research, vol. 9, PMLR, Chia Laguna Resort, Sardinia, Italy, pp. 249–256.

Gu, T., Med, J.-R., Volpi, F., Castelnau, O., Forest, S., Herve-Luanco, E., LecouURIER, F., Proudhon, H., Renault, P.-O., Thilly, L., 2017. Multiscale modeling of the anisotropic electrical conductivity of architecture and nanostructured cu-nb composite wires and experimental comparison. Acta Mater. 141, 131–141.

Gu, Y., Stiles, C.D., El-Awady, J.A., 2024. A statistical perspective for predicting the strength of metals: Revisiting the hall-patch relationship using machine learning. Acta Mater. 266, 119631. http://dx.doi.org/10.1016/j.actamat.2023.119631.

Hasan, M., Mao, Y., Choudhary, K., Tavazza, F., Choudhary, A., Agrawal, A., Acar, P., 2022. Data-driven multi-scale modeling and optimization for elastic properties of cubic microstructures. Integrat. Mater. Manuf. Innov. 11 (2), 230–240. http://dx.doi.org/10.1007/s40192-022-00258-3.

Hashin, Z., Shtrikman, S., 1963. A variational approach to the theory of the elastic behaviour of multiphase materials. J. Mech. Phys. Solids 11 (2), 127–140.

Hill, R., 1965. Continuum micro-mechanics of elastoplastic polycrystals. J. Mech. Phys. Solids 13, 89–101.

Kammer, C., Ponte Castañeda, P., 2024. Fully optimized second-order estimates for the macroscopic behavior and field statistics of particle-reinforced viscoplastic composites. J. Mech. Phys. Solids 186, 105598. http://dx.doi.org/10.1016/j.jmps.2024.105598.

Kingma, D.P., Ba, J., 2014. Adam: A method for stochastic optimization. CoRR abs/1412.6980. URL: https://apisemanticscholar.org/CorpusID:6628106.

Klambauer, G., Unterthiner, T., Mayr, A., Hochreiter, S., 2017. Self-normalizing neural networks. In: Proceedings of the 31st International Conference on Neural Information Processing Systems. NIPS '17, Curran Associates Inc., Red Hook, NY, USA, pp. 972–981.

Krizhevsky, A., Sutskever, I., Hinton, G.E., 2012. ImageNet classification with deep convolutional neural networks. In: Pereira, F., Burges, C., Bottou, L., Weinberger, K. (Eds.), Advances in Neural Information Processing Systems. Vol. 25, Curran Associates, Inc..
Kröner, E., 1958. Berechnung der elastischen konstanten des vielkristalls aus den konstanten des einkristalls. Z. Phys. 151, 504–518.

Lebensohn, R.A., Ponte Castañeda, P., Brenner, R., Castelnau, O., 2011. Full-field vs. Homogenization methods to predict microstructure-property relations for polycrystalline materials. In: Ghosh, S., Dimiduk, D. (Eds.), Chapter 11 of Computational Methods for Microstructure-Property Relationships. Springer, pp. 393–441.

Lecun, Y., Bottou, L., Bengio, Y., Haffner, P., 1998. Gradient-based learning applied to document recognition. Proc. IEEE 86 (11), 2278–2324. http://dx.doi.org/10.1109/5.726791.

Lequeu, P., Gilormini, P., Montheillet, F., Bacroix, B., Jonas, J., 1987. Yield surfaces for textured polycrystals—I. Crystallographic approach. Acta Metall. 35 (2), 439–451. http://dx.doi.org/10.1016/0001-6160(87)90250-1.

Liu, B., Ocegueda, E., Trautner, M., Stuart, A.M., Bhattacharya, K., 2023. Learning macroscopic internal variables and history dependence from microscopic models. J. Mech. Phys. Solids 178, 105329. http://dx.doi.org/10.1016/j.jmps.2023.105329.

Maurel, A., Lund, F., Montagnat, M., 2015. Propagation of elastic waves through textured polycrystals: application to ice. Proc. R. Soc. A. 471, 20140988.

McElfresh, C., Roberts, C., He, S., Prikhodko, S., Marian, J., 2022. Using machine-learning to understand complex microstructural effects on the mechanical behavior of Ti-6Al-4V alloys. Comput. Mater. Sci. 208, 111267. http://dx.doi.org/10.1016/j.commatsci.2022.111267.

McKay, M.D., Beckman, R.J., Conover, W.J., 1979. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21 (2), 239–245, URL: http://www.jstor.org/stable/1268522.

Molinari, A., Canova, G.R., Ahzi, S., 1987. A self-consistent approach of the large deformation polycrystal viscoplasticity. Acta Metall. 35 (12), 2983–2994.

Moulinec, H., Suquet, P., 1994. A fast numerical method for computing the linear and nonlinear mechanical properties of composites. C. R. Acad. Sci. Paris 318 (IIb), 0417–1413.

Nebozhyn, M.V., Gilormini, P., Ponte Castañeda, P., 2001. Variational self-consistent estimates for cubic viscoplastic polycrystals: the effects of grain anisotropy and shape. J. Mech. Phys. Solids 49, 313–340.

Ponte Castañeda, P., 1991. The effective mechanical properties of nonlinear isotropic composites. J. Mech. Phys. Solids 39, 45–71. http://dx.doi.org/10.1016/0022-5096(91)90030-R.

Ponte Castañeda, P., 2015. Fully optimized second-order variational estimates for the macroscopic response and field statistics in viscoplastic crystalline composites. Proc. R. Soc. A 471, http://dx.doi.org/10.1098/rspa.2015.0665.

Ponte Castañeda, P., 2016. Stationary variational estimates for the effective response and field fluctuations in nonlinear composites. J. Mech. Phys. Solids 96, 660–682. http://dx.doi.org/10.1016/j.jmps.2016.06.010.

Ponte Castañeda, P., Suquet, P., 1998. Nonlinear composites. Adv. Appl. Mech. 34, 171–302. http://dx.doi.org/10.1016/S0065-2156(08)70321-1.

Ponte Castañeda, P., Willis, J.R., 1995. The effect of spatial distribution on the effective behavior of composite materials and cracked media. J. Mech. Phys. Solids 43, 1919–1951.

Schulson, E.M., Duval, P., 2009. Creep and Fracture of Ice. Cambridge University Press, http://dx.doi.org/10.1017/CBO9780511581397.

Song, D., Ponte Castañeda, P., 2018. Fully optimized second-order homogenization estimates for the macroscopic response and texture evolution of low-symmetry viscoplastic polycrystals. Int. J. Plast. 110, 272–293. http://dx.doi.org/10.1016/j.ijplas.2018.07.004.

Spahn, J., Andrä, H., Kabel, M., Müller, R., 2014. A multiscale approach for modeling progressive damage of composite materials using fast Fourier transforms. Comput. Methods Appl. Mech. Engrg. 268, 871–883. http://dx.doi.org/10.1016/j.cma.2013.10.017.

Tandon, G.P., Weng, G.J., 1988. A theory of particle-reinforced plasticity. J. Appl. Mech. 55 (1), 126–135. http://dx.doi.org/10.1115/1.3173618.

Zecevic, M., Beyerlein, I.J., Knezevic, M., 2017. Coupling elasto-plastic self-consistent crystal plasticity and implicit finite elements: Applications to compression, cyclic tension-compression, and bending to large strains. Int. J. Plast. 93, 187–211. http://dx.doi.org/10.1016/j.ijplas.2016.07.016.