A Stochastic FE2 Data-Driven Method for Nonlinear Multiscale Modeling

A stochastic data-driven multilevel finite-element (FE2) method is introduced for random nonlinear multiscale calculations. A hybrid neural-network–interpolation (NN–I) scheme is proposed to construct a surrogate model of the macroscopic nonlinear constitutive law from representative-volume-element calculations, whose results are used as input data. Then, a FE2 method replacing the nonlinear multiscale calculations by the NN–I is developed. The NN–I scheme improved the accuracy of the neural-network surrogate model when insufficient data were available. Due to the achieved reduction in computational time, which was several orders of magnitude less than that to direct FE2, the use of such a machine-learning method is demonstrated for performing Monte Carlo simulations in nonlinear heterogeneous structures and propagating uncertainties in this context, and the identification of probabilistic models at the macroscale on some quantities of interest. Applications to nonlinear electric conduction in graphene–polymer composites are presented.


Introduction
Predicting the nonlinear behavior of materials from knowledge of their microstructure is a critical topic in engineering. For example, the development of 3D-printed micromaterials [1][2][3] or of nanomaterials [4,5] with nonlinear behaviors opens exciting opportunities for designing innovative functionalized and enhanced engineering systems. While linear effective properties of heterogeneous materials can be accurately estimated though either analytical [6,7] or numerical techniques [8], predicting the effective behavior of nonlinear materials requires more advanced techniques.
A direct but limited approach is the use of the representative volume element (RVE) to calibrate an empirical nonlinear model. A limitation of such techniques is the number of parameters to be calibrated for complex, nonlinear, or multiphysics problems. To more accurately describe the behavior of general nonlinear materials, the so-called multilevel finite-element (FE 2 ) method [9][10][11][12][13][14][15][16] or computational homogenization has been developed in recent years. In this approach, an RVE is associated to each Gaussian point of a finiteelement macrostructure, and a nonlinear problem must be solved in each integration point and for each iteration of the macrosolving procedure. The drawback of this method, however, is that it induces unaffordable computational times in practical applications.
Another idea, initiated in [27,28], is the use of so-called data-driven approaches in which microscale calculations are performed offline, and are then used as data in an online stage to reconstruct the macroscopic (effective) behavior. For this purpose, several techniques were proposed, including interpolation methods [27,29], neural networks [28,[30][31][32][33][34][35], Bayesian inference [36], Fourier series expansions [37], or Gaussian process regression [38]. In the related techniques, offline data collection is used in a regression process to construct an accurate surrogate model whose evaluation is several orders of magnitude lower than that performing one RVE nonlinear calculation. A critical comparison of several regression techniques used in data-driven multiscale approaches can be found in [39]. In [40], Avery et al. investigated and discussed several regression methods with ANN in homogenization problems of hyperelastic woven composites, and demonstrate its use in advanced dynamic or fluid structure applications. Recent advances of datadriven techniques, including handling history-dependent behaviors such as plasticity, can be found in [35,41,42]. On-the-fly construction of the surrogate model by probabilistic machine learning was proposed in [38]. Developments of neural-network techniques in FE 2 , including feed-forward and recurrent neural networks, can be found in [31,41]. In [43,44], a manifold-based nonlinear reduced-order model in tandem with a digital database was developed for the nonlinear multiscale analysis of hyperelastic structures involving neural networks, a kernel inverse/reconstruction map, and dimension reduction through an isomap.
Stochastic extensions of data-driven methods in multiscale applications are relatively new and unexplored. One of the first analyses in this context can be found in [45,46], where the NEXP method [27] was extended to stochastic problems. In these studies, stochastic parameters were introduced within the surrogate model using a separated representation-interpolation technique. Probability density functions related to the nonlinear macroscale problem were identified. In [47], a machine-learning strategy based on a three-dimensional convolutional neural network was introduced to evaluate the linear effective properties of random materials from geometrical descriptions of RVE. In [24], a framework for uncertainty quantification in a data-driven approach was proposed where self-consistent clustering analysis (SCA) [23,24] was used to reduce computational times in the learning step.
In this paper, the use of data-driven methods for heterogeneous nonlinear materials with uncertainties at both the micro-and the macroscale is addressed. Taking into account uncertainties in nonlinear multiscale methods implies (a) constructing a probabilistic surrogate macromodel from microcalculations, allowing for generating realizations of the macroresponse for a given macroloading; and (b) performing Monte Carlo simulations of the model at the macroscale to quantify uncertainties on the quantities of interest in the structure. In view of its immense computational requirements, direct use of FE 2 for stochastic nonlinear two-scale analysis is not possible. However, data-driven FE 2 approaches have comparable computational costs as compared to classical (one-scale) FEM calculations, and they open the route to developing stochastic two-scale nonlinear approaches. To the best of our knowledge, this problem remains relatively unexplored in the literature. A new stochastic data-driven approach based on RVE calculations was developed for taking into account random effects in nonlinear heterogeneous structures. First, preliminary RVE calculations were performed. These calculations include several microstructural features that varied, such as the distribution of heterogeneities and its volume fraction. Then, for each realization of the random microstructure, the space of macroscopic loading was sampled, and boundary conditions were prescribed on the RVE. Subsequently, the nonlinear problem was solved by FEM. This large database was used to construct a surrogate model whose inputs were the macroloading and the volume fraction, and its output was the macroscopic (homogenized) response. A new hybrid neural-network-interpolation (NN-I) surrogate model is proposed to provide an accurate response with a limited number of realizations. Once constructed, this model can be used within stochastic analysis of two-scale nonlinear structure calculations. At the macroscale, the volume fraction of heterogeneities is considered to be random here, and it was modeled as a stochastic field with given probabilistic characteristics. Then, during the macro-nonlinear resolution, solving the full nonlinear RVE was replaced by the proposed fast surrogate model, which allowed for performing hundreds of macro-non-linear calculations at the cost of classical FEM problems. As a result, statistical postprocessing can be performed on the macroquantities of interest, and probabilistic models could be identified.
The novelties of this paper are twofold. The first is the proposed neural-network-Interpolation FE 2 method, which is an extension of our previous neural-network FE2 method, developed in [28,30]. The NN-I scheme allows for modeling the stochastic spatial variability of the volume fraction in the frame of the FE 2 procedure, leading to the improved accuracy of the surrogate model when limited data are available. The second novelty is the application of this machine-learning method to nonlinear multiscale stochastic problems. Using the proposed approach, FE 2 calculations can be reduced by several orders of magnitude, allowing for Monte Carlo simulation on stochastic nonlinear multiscale structures. It is demonstrated for the first time that uncertainties can be propagated in this context, and probabilistic models can be identified.
The paper is organized as follows. Section 3 presents the equations of the nonlinear RVE problem, and the definitions of the input (macroelectric load) and output (homogenized electric flux) in the nonlinear composite. Section 4 introduces the hybrid neural network/interpolation scheme, and its construction using offline data on RVE is described. In Section 5, the present stochastic data-driven strategy is proposed. Lastly, numerical examples are presented in Section 6.

Brief Review of FE 2 Method for Nonlinear Conduction
The multilevel finite-element method [9,10], 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 [11][12][13][14][15][16], was introduced as a general multiscale method for solving nonlinear heterogeneous structural problems. The basic underlying idea is that two levels of finite elements must be 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 state (strain, electric field, etc.) are prescribed on the boundary of each RVE. After solving each nonlinear problem at each integration point, the appropriate macroscopic response (stress, electric flux), is averaged over the RVE and provided at the macrointegration 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 1).
For the sake of simplicity, a brief review of the method in a context of nonlinear conduction is presented. We consider a macroscopic structure associated with a domain Ω ⊂ R 3 , with a boundary ∂Ω. The assumption of scale separation is adopted (an extension of the method to second-order homogenization can be found in [14]). The microstructure was assumed to be characterized by an RVE associated with a domain Ω ⊂ R 3 , with boundary ∂Ω.
In the context of nonlinear electric conduction, electric field E(x) is related to the electric flux, or electric displacement j(x) by a nonlinear local constitutive relationship. Electric field E is defined by E(x) = −∇φ(x), where φ is the electric potential, ∇(.) is the gradient operator, and x is a material point within Ω. In the following, (.) notations denote macroscale quantities. For a given macroscopic electric field E, the RVE problem is to find φ(x), such that where ∇ · (.) is the divergence operator. The constitutive law is given by where F nl is a local nonlinear operator (specified in Section 3). The equilibrated electric field should satisfy where V is the volume of Ω. Equation (3) can be verified, e.g., by the following boundary condition: whereφ(x) is a periodic function over Ω. In the presence of imperfect interfaces and surface electric flux along interfaces (see [48]), the effective electric current J is defined according to where j s is a surface electric flux (see Section 3). In the so-called FE 2 method, the constitutive law J -E is unknown, but can be numerically obtained by solving a nonlinear problem over the RVE, detailed as follows (see Figure 1): Given E:

2.
Use a numerical method such as FEM with an iterative solver such as the Newton method to solve nonlinear Problems (1), (2), and (4) (see details in the following).

3.
Compute the spatial average of the electric flux over the RVE to obtain J.
In what follows, a detailed numerical implementation of a FE 2 problem in a context of nonlinear electric conduction is presented to better understand where Problems (1), (2), and (4) must be solved within finite-element calculation at the macroscopic scale. The macroscopic problem at the macroscale is given by with boundary conditions: where and Ω φ and ∂Ω J denote the Dirichlet and Neumann complementary boundaries, respectively.
In what follows, we assume J * n = 0. Then, the weak form corresponding to (6) is given by: Problem (8) is nonlinear. Then, a Newton method is employed to solve it. A first-order Taylor expansion of R(φ) gives where φ k is a solution provided at a previous iteration, and D ∆φ R(φ) is the Gateaux derivative, expressed by The corresponding linearized problem is given by with More details can be found in [48]. Classical FEM discretizing of (11) leads to linear system Then, the macroelectric potential is updated according to and (13) is solved until a convergence criterion is reached. In FE 2 , the main source of computational cost is the numerical evaluation of J and k, obtained by solving nonlinear RVE Problems (1), (2), and (4) at each Gaussian point. To address this issue, we introduce a data-driven approach where the estimation of J is provided by a neural-network-based surrogate model. Tangent matrix k can be computed by a perturbation approach as where e j is a unitary vector basis, and α << 1 a perturbation parameter.
Then, to compute macro-FEM nonlinear calculations, relationship J = F nl (E) is missing. In [30], we proposed a surrogate model to construct such a relationship using neural networks. In the present paper, this idea is extended to random microstructures, as detailed in the next section.

Micro-Non-Linear Conduction Model for Graphene-Reinforced Composites
In this section, the nonlinear conduction model in graphene-reinforced polymer composites is defined. The nonlinear RVE problem is described as follows. The microstructure was assumed to be characterized by an RVE defined in a domain Ω ⊂ R 3 that contained N randomly distributed planar multilayer graphene sheets (see Figure 2). The graphene sheets were assumed to be aligned along the x − y plane. We chose this configuration for two reasons: (i) when samples made of graphene-reinforced polymer are obtained via injection molding, the graphene sheets can be aligned in the direction of the polymer flow [49]. Then, this configuration is representative of samples manufactured by the injection-molding process. Second, such an orientation induces strong anisotropy of the effective nonlinear conductive behavior of the material. The potential of the present approach to deal with such a challenging problem is then illustrated. To avoid meshing their thickness [48], the graphene sheets were modeled as highly conducting imperfect surfaces here [50]. The graphene surfaces with zero thickness are collectively denoted by Γ. The nonlinear behavior is related to the electric tunnelling effect here, which may be an explanation for the observed nonlinear behavior and low percolation thresholds in the nanocomposites (see [30,48]).
The energy of the system is defined by where density functions ω b and ω s are the bulk and surface density functions, respectively, expressed by In (18), E s (x) and j s (x) are the surface electric field and surface current density with respect to the graphene sheets, where E s = PE = −P∇φ, where P = I − n ⊗ n is a projector operator characterizing the projection of a vector along the tangent plane to Γ at a point x ∈ Γ, and n is the unit normal vector to Γ.
The nonlinear electric-conduction law including the tunneling effect in the polymer matrix is given by where d cut is a cut-off parameter, and k p 0 is the electric-conductivity tensor of the polymer matrix without tunneling effects. The distance function between graphene sheets d(x) is defined here as the sum of the two smallest distances between a point x in the polymer matrix and the two nearest-neighbor graphene sheets (see more details in [48]). Nonlinear tunneling current G was proposed by Simmons [51] as Above, Φ 0 and d denote barrier height and barrier width, respectively, e and m are the charge and the effective mass of electron, respectively, and h p is the Planck constant. Surface current density j s of graphene surface Γ is related to surface electric field E s [50] through where k s denotes the the surface conductivity of the graphene. This tensor can be related to the conductivity of the volume (multilayer) graphene as: where h is the thickness of the graphene sheet. Considering the constitutive equations above, and minimizing the electric power (17) with respect to the electric potential field, the weak form is obtained as where φ ∈ H 1 (Ω), φ satisfying the Dirichlet boundary conditions over ∂Ω and δφ ∈ H 1 (Ω), δφ = 0 over ∂Ω. The RVE is subjected to boundary conditions (4). The weak form (23) can be solved by the finite-element method.

Stochastic Nonlinear Machine-Learning Model
The objective of the present work was to construct a surrogate model relating macroscopic electric field E and volume fraction of graphene inclusions f to nonlinear macroscopic electric flux response J (see Figure 3). At the microscale, the microstructure was randomly distributed (see Figure 2). Here, due to the scale-separation assumption, it was expected that, despite the random nature of the microstructure, deterministic effective properties at the microscopic scale with respect to the microstructure would be obtained for a given volume fraction. At the macroscale, uncertainties where then only related to nonhomogeneous distributions of volume fractions. Then, it was assumed that, at the macroscale, the volume fraction was the only stochastic parameter.

Data Generation
We first define a set of K electric-field vectors, The values of E k were generated using Latin hypercube sampling (LHS) [52]. Then, we define a collection of microstructures as follows. A set of P volume fractions are defined, f α , α = 1, 2, . . . , P. For each volume fraction f α , N random microstructures satisfying volume fraction f α were generated and are denoted by Ω i α , i = 1, 2, . . . , N. Then, for each macroelectric field vector E k , each volume fraction f α and each realization of microstructure Ω i α , nonlinear problem (23) is solved by finite elements to obtain macroelectric displacement vector J k α,i . As discussed above, the scale-separation assumption allows for removing the stochastic nature of the microstructures at the macroscale. However, due to the RVE size and the random distribution of the inclusions, the outcome intensity of a given electric field is also stochastic. To this purpose, homogenization is performed using stochastic averaging, i.e., for each macroelectric field vector E k and each volume fraction f α , we compute the average over N microstructures realizations to obtain J

Construction of Neural-Network Surrogate Model
An issue in NN models is that a large set of data may be required to obtain good accuracy, especially for a large number of input parameters [28]. To overcome this limitation, we propose here a hybrid NN/interpolation surrogate model as follows.
First, for each volume fraction f α , α = 1, 2, . . . , P, used in the training dataset, we define a separate NN, denoted by N α , in order to construct the following relationship: Then, given E and for an arbitrary volume fraction f ∈ f 1 , f P a Lagrangian interpolation scheme is used to compute J(E, f ) as where N k ( f ) is the set of indices that includes only k out of P NNs, corresponding to the k volume fractions of f nearest to those in training dataset f 1 , f 2 , . . . , f P . In Equation (25), l j ( f ) are the Lagrangian basis polynomials. Here, k = 3 was employed where, as a result, polynomials l j ( f ) were second-order.
With this approach, a notion of locality is introduced in the interpolation scheme that leads to better overall predictions, especially in areas where relationship (E, f ) − J(E, f ) exhibits strong nonlinearity. A schematic of the surrogate-model construction is summarized in Figure 4.

Nonlinear Stochastic Macroscale Calculations
Here, stochastic macroscopic structural problem is described. At the macroscale, it was assumed that there existed uncertainty in the local volume fraction f of the graphene sheets in the general case described by a 3D homogeneous Gaussian stochastic field in the x, y, and z axes. In particular, if x ≡ (x, y, z), then f (x) is considered to be of the form: In the above equation, µ and σ are the random field mean value and standard deviation, respectively, θ denotes the random outcome, and f 0 (x, θ) is a zero-mean unit variance Gaussian field with correlation structure R f 0 given by whereâ x ,â y andâ z are correlation length parameters along the x, y, and z axes, respectively. Next, an approximation of field f 0 can be obtained using the Karhunen-Loeve series expansion [53]. Specifically, let λ n , φ n denote the eigenvalues and eigenfunctions that satisfy the eigenvalue problem R f 0 (x 1 , x 2 )φ n (x 2 )dx 2 = λ n φ n (x 1 ), ∀n = 1, . . .. This is a Fredholm integral equation and it is typically solved using the finite-element method [53]. Then, f 0 can be written as with {z n } ∞ n=1 being a series of uncorrelated Gaussian random variables with zero mean and unit variance. In practice, the above series expansion is truncated after M KL terms, giving the following approximation of f 0 .
which yields, by virtue of Equation (26), Equation (30) allows for us to generate realizations of the field f (x, θ) by generating M KL -tuples (z 1 , . . . , z KL ) from their distribution. Subsequently, if we consider macrostructure Ω defined in Section 2 and the associated finite-element mesh, at each Gaussian point of element Ω e , e = 1, 2, . . . , N e with coordinate x, a random value of volume fraction f (x) can be assigned using (30). During the Newtonian procedure to solve the structural problem, for f and E given at each Gaussian point of the macromesh structure, the corresponding value of J is provided by the surrogate model (25) (see Figure 4). For one realization of the volume-fraction distribution generated by Equation (30), the cost of one two-scale nonlinear structural problem is drastically reduced with the present NN surrogate model, allowing for performing a large number of macrocalculations at a low cost to conduct statistics on quantities of interest in a structure.
Lastly, Monte Carlo simulations were performed on the macroscale problem by evaluating R realizations of macrostructures. For each realization r = 1, 2, . . . , R, the volume fraction was randomly generated in each Gaussian point by using Equation (30); in total, R nonlinear multiscale problems were solved using the above-described procedure. Lastly, statistics can be computed on quantities of interest using the R nonlinear FEM solutions. The overall procedure is summarized in Figure 5.

Data Collection
The data were obtained by performing preliminary calculations on the RVE described in Section 4.1. Eight different volume fractions were considered: f 1 = 0.53%, , f 2 = 0.66%, f 3 = 0.79%, f 4 = 0.92%, f 5 = 1.05%, f 6 = 1.19%, f 7 = 1.32%, and f 8 = 1.58% (see Figure 2). For each volume fraction, 15 realizations of random microstructures were generated except for the higher volume fraction, for which only 9 realizations were conducted. For each volume fraction and for each realization, 500 realizations of macroscopic electric field E were generated using Latin hypercube sampling. For each case, corresponding electric flux J was numerically computed by solving nonlinear Problem (23) on the RVE. The total number of solved nonlinear problems was then 57,000. All these calculations could be performed in parallel.

Validation of Hybrid NN-Interpolation Surrogate Model
The accuracy of the proposed hybrid NN-interpolation surrogate model was first validated by comparing its response with full-field simulations on microstructures for different volume fractions. Regarding the characteristics of the trained neural networks, in all cases, one-hidden-layer architectures were considered with the optimal number of neurons varying for each case, as shown in Table 1. Moreover, the hyperbolic tangent function was employed as the activation function, and Levenberg-Marquardt as the optimizer in all NNs. The plotted curves were obtained as the average over the different realizations of the microstructure. Results are provided in Figure 6. For low volume fractions, the response was linear, while for larger volume fractions, the response was strongly nonlinear. In all cases, the surrogate model could accurately reproduce the effective nonlinear response of the material.
A validation of the interpolation procedure described in Section 4.2 is provided in Figure 7, where discrete data obtained by nonlinear FEM calculations on the RVEs are compared to the corresponding model predictions, computed using Equation (25) under various E x scenarios, with E y = E z = 0. The discrete data points obtained by FEM are denoted by marks, while the continuous interpolation with respect to the volume fraction is denoted by solid lines, which confirmed the good accuracy of this scheme.

Stochastic 2-Scale Nonlinear Structure Analysis
In this example, macroscopic stochastic nonlinear computations were performed using the procedure described in Section 5. In particular, 9 different Gaussian fields of volume fractions are investigated at the macroscale, where the studied macrostructure, described in Figure 8, was a plate with a central hole. The plate was subjected to potential boundary conditions such as Φ = Φ 1 on x = 0 and Φ = Φ 2 on x = L. A 3D mesh of 1934 elements is used to discretize the domain. Due to the low thickness of the structure, we assumed that the volume fraction did not vary in the z coordinate direction. Next, in order to define the aforementioned Gaussian fields, three different settings were first initialized: for Setting A, we set µ A = 0.9% and σ A = 0.11 µ A ; for Setting B, we set µ B = 1.05% and σ B = 0.19 µ B ; lastly, for Setting C, µ C = 1.05% and σ C = 0.38 µ C . Then, for each aforementioned setting, we considered a z =â y =â and assign three different values toâ, namely,â 1 = 6 µm,â 2 = 12 µm andâ 3 = 24 µm. A sample for each of these fields is illustrated in Figure 9. This figure indicates that an increase in the field standard deviation led to larger variations of volume fraction f along the spatial domain. Moreover, a small correlation-length parameter, such asâ = 6 µm, produced more "wavy" realizations, while for larger values (â = 12 and 24 µm) the realizations became smoother.
(a) Field A: µ A = 0.9%, σ A = 0.1%, For each of the 3 Gaussian distributions A, B, and C, we analyzed the 3 correlation lengthsâ 1 ,â 2 andâ 3 . For each case, we conducted 100 realizations. Then, in total, we conducted 900 FE 2 -NN simulations using the procedure described in Section 5. For each one, a stochastic distribution of volume fraction was generated in the elements using (30). The macroscopic quantity of interest is defined here as the average macroscopic flux in the domain Ω as where V is the volume of Ω. The convergence of the components of J * is depicted in Figure 10. In all cases, statistical convergence could be achieved. For the lowest average values f and standard deviation σ of the volume fractions (Cases 1-3 in Figure 10), correlation lengthâ did not have significant influence on the convergence rate. However, for larger values of f and σ, convergence could be much slower (e.g., Case 9 in Figure 10 ), where around 50 realizations are necessary to achieve convergence. This clearly illustrates the interest of the proposed surrogate-based multiscale method, where each realization is performed at the cost of a classical FEM simulation. In contrast, using standard FE 2 would not allow performing this kind of statistical analysis with available computer resources. (c)  Average distributions of local current densities over 100 realization are plotted in Figure 11 corresponding to distribution A and correlation lengthâ=6 µm. Clear anisotropy of the effective behavior induced by the aligned graphene sheets along the x − y plane can be appreciated. Comparing Figure 11a,b, we can see a clear difference in the magnitude of the J x and J z values, indicating that the effective conductivity in the z direction was much lower than that in the x − y plane. The present method could capture such anisotropic behavior in a nonlinear stochastic context. The evolution of the quantity of interest J * x was plotted with respect to the difference of the potential applied over macrostructure Φ 2 − Φ 1 in Figure 12. Various distributions of CNT volume fractions and different correlation lengths were taken into account for comparison. For each case, 100 realizations were computed, from which we obtained the average and deviation of J * x . For instance, in Figure 12a, correlation lengthâ = 6 µm is for all three different CNT volume-fraction distributions. The averaged value of J * x was independent on standard deviation σ of the Gaussian distribution, whereas the deviation of J * x increased slightly with increasing σ. The same phenomenon could also be observed in Figure 12b,c. Furthermore, by comparing Figure 12a-c, the increase in correlation length led to a tiny increase in the deviation of J * x , but had no effect on its averaged value. Lastly, in Figure 13, distributions of target values J * x , J * y and J * z are plotted for selected cases of the probabilistic models describing the distribution of the volume fraction in the macroscale. In Figure 14, the associated empirical cumulative distribution functions (ECDFs) are provided. These functions were identified from the histograms in Figure 13. These allow for a direct quantitative reading of key values of interest (minimum, maximum, mean, percentiles, etc.) regarding the macroscopic quantities. ECDFs also have the property of converging to the true CDF of the stochastic quantities of interest as the number of samples is increased [54]. Typically, an accurate estimate of a CDF would require a very large number of samples (>10 5 ); however, performing these many evaluations of nonlinear multiscale models would be computationally prohibitive. In this regard, the use of the proposed surrogate is the only viable approach to obtain reliable approximations of the CDFs under investigation. This demonstrates the potential of the present approach in constructing probabilistic models for macroquantities of interest in nonlinear multiscale models of random materials. deviation of J * x increases slightly with increasing σ. T observed in Figure 12b        (a) J * x , µ = 0.9%, σ = 0.11 µ,â = 24 µm; (b) J * x , µ = 1.05% , σ = 0.19 µ,â = 24 µm; (c) J * y , µ = 0.9%, σ = 0.11µ,â = 24 µm; (d) J * y , µ = 1.05%, σ = 0.19 µ,â = 24 µm; (e) J * z , µ = 0.9%, σ = 0.11 µ,â = 24; (f) J * z , µ = 1.05%, σ = 0.19 µ,â = 24 µm.

Conclusions
A stochastic data-driven multilevel finite-element (FE 2 ) method was proposed to solve nonlinear heterogeneous structures with uncertainties at both the micro-and the macrolevel. A hybrid neural-network-interpolation (NN-I) scheme was developed to improve the accuracy of NN surrogate models, allowing for the use of a lower number of representative volume element (RVE) nonlinear calculations, which serve as a database to train the neural networks. This NN-I surrogate model was used to develop a data-driven method for nonlinear heterogeneous conduction in a stochastic framework: uncertainties can be included on both the micro-and the macrolevel. More specifically, the drastic reduction in computational costs brought by the NN-I surrogate model allows Monte Carlo simulations of nonlinear heterogeneous structures. This framework was applied to propagate uncertainties in such nonlinear multiscale models, and demonstrated that it can be used to identify probabilistic models related to some quantities of interest at the macroscale in a fully nonlinear, anisotropic context. Author Contributions: Conceptualization, J.Y and V.P.; software, X.L. and L.P.; formal analysis, J.Y. and V.P.; investigation, X.L., J.Y., V.P., L.P. and I.K; writing-original draft preparation, J.Y., X.L. and V.P.; writing-review and editing, X.L., J.Y., V.P., L.P. and I.K.; supervision, J.Y.; project administration, J.Y.; funding acquisition, X.L. All authors have read and agreed to the published version of the manuscript.
Funding: Xiaoxin LU thanks the support from SIAT Innovation Program for Excellent Young Researchers 346 (E1G045).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

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