*4.2. 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 <sup>N</sup> *<sup>α</sup>* , in order to construct the following relationship:

$$
\overline{\mathbf{J}}\_{\alpha}(\overline{\mathbf{E}}) = \mathcal{N}^{\alpha}(\overline{\mathbf{E}}).\tag{24}
$$

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

$$\overline{\mathbf{J}}(\overline{\mathbf{E}},f) = \sum\_{j \in N\_k(f)} l\_j(f)\overline{\mathbf{J}}\_j(\overline{\mathbf{E}}),\tag{25}$$

where *N<sup>k</sup>* (*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), *lj*(*f*) are the Lagrangian basis polynomials. Here, *k* = 3 was employed where, as a result, polynomials *lj*(*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.

**Figure 4.** Hybrid neural-network/interpolation surrogate model to describe macroscopic nonlinear behavior.

#### **5. 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:

$$f(\mathbf{x}, \theta) = \mu + \sigma f\_0(\mathbf{x}, \theta). \tag{26}$$

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 *Rf*<sup>0</sup> given by

$$\mathcal{R}\_{f\_0}(\mathbf{x}\_1, \mathbf{x}\_2) = \exp\left[-\left(\frac{|\mathbf{x}\_1 - \mathbf{x}\_2|}{\mathfrak{A}\_{\mathcal{X}}} + \frac{|y\_1 - y\_2|}{\mathfrak{A}\_{\mathcal{Y}}} + \frac{|z\_1 - z\_2|}{\mathfrak{A}\_{\mathcal{Z}}}\right)\right] \tag{27}$$

where *a*ˆ*x*, *a*ˆ*<sup>y</sup>* and *a*ˆ*<sup>z</sup>* are correlation length parameters along the *x*, *y*, and *z* axes, respectively.

Next, an approximation of field *f*<sup>0</sup> can be obtained using the Karhunen–Loeve series expansion [53]. Specifically, let *λn*, *φ<sup>n</sup>* denote the eigenvalues and eigenfunctions that satisfy the eigenvalue problem R *Rf*<sup>0</sup> (*x*1, *x*2)*φn*(*x*2)*dx*<sup>2</sup> = *λ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*<sup>0</sup> can be written as

$$f\_0(\mathbf{x}, \boldsymbol{\theta}) = \sum\_{n=1}^{\infty} \sqrt{\lambda\_n} z\_n(\boldsymbol{\theta}) \phi\_n(\mathbf{x}) \tag{28}$$

with {*zn*} ∞ *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 *MKL* terms, giving the following approximation of *f*0.

$$f\_0(\mathbf{x}, \boldsymbol{\theta}) \approx \sum\_{n=1}^{M\_{KL}} \sqrt{\lambda\_n} z\_n(\boldsymbol{\theta}) \phi\_n(\mathbf{x}) \,\tag{29}$$

which yields, by virtue of Equation (26),

$$f(\mathfrak{x}, \theta) \approx \mu + \sigma \sum\_{n=1}^{M\_{\rm KL}} \sqrt{\lambda\_n} z\_n(\theta) \phi\_n(\mathfrak{x}). \tag{30}$$

Equation (30) allows for us to generate realizations of the field *f*(*x*, *θ*) by generating *MKL*-tuples (*z*1, . . . , *zKL*) 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<sup>e</sup>* 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.

**Figure 5.** Stochastic nonlinear 2-scale procedure.

#### **6. Numerical Examples**

#### *6.1. 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* <sup>1</sup> = 0.53%, , *f* <sup>2</sup> = 0.66%, *f* <sup>3</sup> = 0.79%, *f* <sup>4</sup> = 0.92%, *f* <sup>5</sup> = 1.05%, *f* <sup>6</sup> = 1.19%, *f* <sup>7</sup> = 1.32%, and *f* <sup>8</sup> = 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.

### *6.2. 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.


**Table 1.** Characteristics of neural networks.

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<sup>x</sup>* scenarios, with *E<sup>y</sup>* = *E<sup>z</sup>* = 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.

**Figure 6.** Comparisons between direct simulations obtained by nonlinear FEM calculations on the RVE and the neuralnetwork–interpolation surrogate model: values of *J<sup>x</sup>* as a function of a unidirectional effective electric field *E<sup>x</sup>* , *E<sup>y</sup>* = *E<sup>z</sup>* = 0; (**a**) 0.53 vol%, (**b**) 0.79 vol%, (**c**) 1.19 vol%,(**d**) 1.58 vol%.

**Figure 7.** Comparisons between direct simulations obtained by nonlinear FEM calculations on the RVE and NN–I model under various *E<sup>x</sup>* ranging from 0.0050 to 0.0125 V/nm: values of *J<sup>x</sup>* as a function of the CNT volume fraction, *E<sup>y</sup>* = *E<sup>z</sup>* = 0.

#### *6.3. 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 Φ = Φ<sup>1</sup> on *x* = 0 and Φ = Φ<sup>2</sup> on *x* = *L*. A 3D mesh of 1934 elements is used to discretize the domain.

**Figure 8.** Structural problem: geometry, boundary conditions, and mesh.

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 *µ <sup>A</sup>* = 0.9% and *σ <sup>A</sup>* = 0.11 µ *<sup>A</sup>*; for Setting B, we set *µ <sup>B</sup>* = 1.05% and *σ <sup>B</sup>* = 0.19 µ *B* ; lastly, for Setting C, *µ <sup>C</sup>* = 1.05% and *σ <sup>C</sup>* = 0.38 µ *<sup>C</sup>*. Then, for each aforementioned setting, we considered *a*ˆ*<sup>z</sup>* = *a*ˆ*<sup>y</sup>* = *a*ˆ and assign three different values to *a*ˆ, namely, *a*ˆ<sup>1</sup> = 6 µm, *a*ˆ<sup>2</sup> = 12 µm and *a*ˆ<sup>3</sup> = 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 *a*ˆ = 6 µm, produced more "wavy" realizations, while for larger values (*a*ˆ = 12 and 24 µm) the realizations became smoother.

**Figure 9.** Sample realizations of three Gaussian fields A, B, and C for different correlation lengths *α*ˆ = 6, 12, 24 µm.

For each of the 3 Gaussian distributions A, B, and C, we analyzed the 3 correlation lengths *a*ˆ1, *a*ˆ<sup>2</sup> and *a*ˆ3. For each case, we conducted 100 realizations. Then, in total, we conducted 900 FE<sup>2</sup> -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

$$\mathbf{J}^\* = \frac{1}{\overline{V}} \int\_{\overline{\Omega}} \mathbf{\tilde{J}} d\Omega,\tag{31}$$

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 *a*ˆ 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<sup>2</sup> would not allow performing this kind of statistical analysis with available computer resources.

**Figure 10.** Averaged current-density components as a function of the number of realizations under various distributions of CNT volume fraction and different correlation lengths. (**a**) *J* ∗ *x* ; (**b**) *J* ∗ *y* ; (**c**) *J* ∗ *z* .

Average distributions of local current densities over 100 realization are plotted in Figure 11 corresponding to distribution A and correlation length *a*ˆ=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<sup>x</sup>* and *J<sup>z</sup>* 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.

**Figure 11.** Averaged current density ¯*J<sup>x</sup>* and ¯*J<sup>y</sup>* over 100 realizations calculated by the NN–I model for the composite structure for potential difference <sup>Φ</sup><sup>2</sup> − <sup>Φ</sup><sup>1</sup> = 144 V. The CNT volume fraction obeys distribution A with *<sup>µ</sup> <sup>A</sup>* = 0.9 vol%, *σ <sup>A</sup>* = 0.1%, and correlation length *a*ˆ=6 *µ*m: (**a**) *J<sup>x</sup>* -component: 3D view; (**b**) *J<sup>x</sup>* -component: plots along different planes; (**c**) *J<sup>y</sup>* -component: 3D view; (**d**) *J<sup>y</sup>* -component: plots along different planes.

The evolution of the quantity of interest *J* ∗ *<sup>x</sup>* was plotted with respect to the difference of the potential applied over macrostructure <sup>Φ</sup><sup>2</sup> − <sup>Φ</sup><sup>1</sup> 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 *a*ˆ = 6 µm is for all three different CNT volume-fraction distributions. The averaged value of *J* ∗ *<sup>x</sup>* 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* ∗ *<sup>y</sup>* and *J* ∗ *<sup>z</sup>* 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.

μ=0.9%,σ=0.1%,â=12 **Figure 12.** *Cont.*

μ=1.05%,σ=0.2%,â=12 μ=1.05%,σ=0.4%,â=12

0.46 0.47 0.48 0.49

0.062 0.064

0.06

0

0.06 0.062 0.064

0.5

1

J\*

(A/cm

)

2

x

1.5

2

2.5

96

μ=0.9%,σ=0.1%,â=24 μ=1.05%,σ=0.2%,â=24 μ=1.05%,σ=0.4%,â=24

96

2 - 1

0 24 48 72 96 120 144

2 - 1 (V)

0.46 0.47 0.48 0.49 (V)

120

2.2

2.3

144

120

144

0

0

0.5

0.5

0.06 0.061 0.062 0.063

0.06 0.061 0.062 0.063

J\*

J\*

2

2

x

x

1

1

1.5

1.5

2

2

2.5

2.5


2 - 1

0.465 0.47 0.475 0.48 0.485 0.49

0.465 0.47 0.475 0.48 0.485 0.49

μ=0.9%,σ=0.1%, â=6 μ=1.05%,σ=0.2%,â=6 μ=1.05%,σ=0.4%,â=6

μ=0.9%,σ=0.1%, â=6 μ=1.05%,σ=0.2%,â=6 μ=1.05%,σ=0.4%,â=6

96

96

0 24 48 72 96 120 144

0 24 48 72 96 120 144

(V)

(V)

120

120

144

144

**Figure 12.** Averaged current density *J* ∗ *<sup>x</sup>* and corresponding deviation as a function of potential difference <sup>Φ</sup><sup>2</sup> − <sup>Φ</sup><sup>1</sup> for various distributions of CNT volume fraction under different correlation lengths *a* ˆ. (**a**) *a* ˆ<sup>1</sup> = 6 µm; (**b**) *a* ˆ<sup>2</sup> = 12 µm; (**c**) *a* ˆ <sup>3</sup> = 24 µm. Color zones indicate ranges of values.

**Figure 13.** Histograms associated to probabilistic models describing distribution of **J** ∗ components at the macroscale. Values of **J** ∗ are reported for a fixed value of macroboundary condition <sup>Φ</sup><sup>2</sup> − <sup>Φ</sup><sup>1</sup> = <sup>144</sup>*V*; (**a**) *<sup>J</sup>* ∗ *x* , *µ* = 0.9%, *σ* = 0.11 *µ*, *a*ˆ = 24 µm; (**b**) *J* ∗ *x* , *µ* = 1.05% , *σ* = 0.19 *µ*, *a*ˆ = 24 µm; (**c**) *J* ∗ *y* , *µ* = 0.9%, *σ* = 0.11*µ*, *a*ˆ = 24 µm; (**d**) *J* ∗ *y* , *µ* = 1.05%, *σ* = 0.19 *µ*, *a*ˆ = 24 µm; (**e**) *J* ∗ *z* , *µ* = 0.9%, *σ* = 0.11 *µ*, *a*ˆ = 24; (**f**) *J* ∗ *z* , *µ* = 1.05%, *σ* = 0.19 *µ*, *a*ˆ = 24 µm.

**Figure 14.** Identified probabilistic models (empirical cumulated distribution functions) for generating distributions of **J** ∗ components at the macroscale. Values of **J** ∗ are reported for a fixed value of macroboundary condition <sup>Φ</sup><sup>2</sup> − <sup>Φ</sup><sup>1</sup> = 144 V; (**a**) *J* ∗ *x* , *µ* = 0.9%, *σ* = 0.11 *µ*, *a*ˆ = 24 µm; (**b**) *J* ∗ *x* , *µ* = 1.05% , *σ* = 0.19 *µ*, *a*ˆ = 24 µm; (**c**) *J* ∗ *y* , *µ* = 0.9%, *σ* = 0.11*µ*, *a*ˆ = 24 µm; (**d**) *J* ∗ *y* , *µ* = 1.05%, *σ* = 0.19 *µ*, *a*ˆ = 24 µm; (**e**) *J* ∗ *z* , *µ* = 0.9%, *σ* = 0.11 *µ*, *a*ˆ = 24; (**f**) *J* ∗ *z* , *µ* = 1.05%, *σ* = 0.19 *µ*, *a*ˆ = 24 µm.
