**2. Brief Review of FE<sup>2</sup> Method for Nonlinear Conduction**

The multilevel finite-element method [9,10], also called FE<sup>2</sup> in the literature, as it involves two levels of finite-element simulations, and independently proposed by several other authors and groups [11–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 <sup>Ω</sup> <sup>⊂</sup> <sup>R</sup><sup>3</sup> , 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 <sup>Ω</sup> <sup>⊂</sup> <sup>R</sup><sup>3</sup> , 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

$$\nabla \cdot \mathbf{j}(\mathbf{x}) = \mathbf{0} \,\forall \mathbf{x} \in \Omega,\tag{1}$$

where ∇ · (.) is the divergence operator. The constitutive law is given by

$$\mathbf{j(x)} = \mathcal{F}^{\mathrm{nl}}(\mathbf{E(x)}).\tag{2}$$

where <sup>F</sup>*nl* is a local nonlinear operator (specified in Section 3). The equilibrated electric field should satisfy

$$\overline{\mathbf{E}} = \frac{1}{V} \int\_{\Omega} \mathbf{E}(\mathbf{x}) d\Omega,\tag{3}$$

where *V* is the volume of Ω. Equation (3) can be verified, e.g., by the following boundary condition:

$$\boldsymbol{\phi}(\mathbf{x}) = -\overline{\mathbf{E}} \cdot \mathbf{x} + \tilde{\boldsymbol{\phi}}(\mathbf{x}) \text{ on } \partial \Omega,\tag{4}$$

where *φ*˜(**x**) is a periodic function over Ω.

**Figure 1.** Schematic of classical FE<sup>2</sup> method for nonlinear heterogeneous conduction problem (adapted from [8]).

In the presence of imperfect interfaces and surface electric flux along interfaces (see [48]), the effective electric current **J** is defined according to

$$\mathbf{\tilde{J}} = \frac{1}{V} \left( \int\_{\Omega} \mathbf{j}(\mathbf{x}) d\Omega + \int\_{\Gamma} \mathbf{j}^{s}(\mathbf{x}) d\Gamma \right), \tag{5}$$

where **j** *s* is a surface electric flux (see Section 3). In the so-called FE<sup>2</sup> 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**:


In what follows, a detailed numerical implementation of a FE<sup>2</sup> 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

$$\nabla \cdot \overline{\mathbf{J}} = 0 \text{ in } \overline{\Omega}, \tag{6}$$

with boundary conditions:

$$
\overline{\boldsymbol{\phi}} = \overline{\boldsymbol{\phi}}^\* \text{ on } \partial \overline{\Omega}\_{\boldsymbol{\phi}\prime} \ \overline{\mathbf{J}} \cdot \mathbf{n} = \overline{\boldsymbol{\zeta}}\_n^\* \text{ on } \partial \overline{\Omega}\_{\boldsymbol{l}\prime} \tag{7}
$$

where and Ω*<sup>φ</sup>* and *∂*Ω*<sup>J</sup>* denote the Dirichlet and Neumann complementary boundaries, respectively.

In what follows, we assume *J* ∗ *<sup>n</sup>* = 0. Then, the weak form corresponding to (6) is given by:

$$\int\_{\overline{\Omega}} \overline{\mathbf{J}}(\overline{\Phi}) \cdot \nabla(\delta \overline{\Phi}) d\Omega = \overline{\mathcal{R}}(\overline{\Phi}) = 0. \tag{8}$$

Problem (8) is nonlinear. Then, a Newton method is employed to solve it. A first-order Taylor expansion of *R*(*φ*) gives

$$
\overline{\mathcal{R}}(\overline{\phi}^k + \Delta \overline{\phi}) \simeq \overline{\mathcal{R}}(\overline{\phi}^k) + D\_{\Delta \overline{\phi}} \overline{\mathcal{R}}(\overline{\phi}^k), \tag{9}
$$

where *φ k* is a solution provided at a previous iteration, and *D*∆*φR*(*φ*) is the Gateaux derivative, expressed by

$$D\_{\Delta \overline{\phi}} \overline{\mathcal{R}}(\overline{\phi}^k) = \left[ \frac{d}{d\alpha} (\overline{\mathcal{R}}(\overline{\phi} + \alpha \Delta \overline{\phi})) \right]\_{\alpha = 0}. \tag{10}$$

The corresponding linearized problem is given by

$$D\_{\Delta \overline{\phi}} \overline{R}(\overline{\phi}^k) = -\overline{R}(\overline{\phi}^k), \tag{11}$$

with

$$D\_{\Delta\overline{\phi}}\overline{\mathcal{R}}(\overline{\phi}^k) = -\int\_{\overline{\Omega}} \overline{\mathbf{k}}(\phi^k) \nabla(\Delta\overline{\phi}) \cdot \nabla(\delta\phi) d\Omega. \tag{12}$$

More details can be found in [48]. Classical FEM discretizing of (11) leads to linear system

$$
\overline{\mathbf{K}}\_T(\overline{\phi}^k)\Delta\overline{\phi} = -\mathbf{R}(\overline{\phi}^k). \tag{13}
$$

Then, the macroelectric potential is updated according to

$$
\overline{\phi}^{k+1} = \overline{\phi}^k + \Delta \overline{\phi} \tag{14}
$$

and (13) is solved until a convergence criterion is reached. In FE<sup>2</sup> , 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

$$\left(\overline{k}\_T\right)\_{i\overline{j}}(\overline{\mathbf{E}}) = \frac{\partial \overline{f}}{\partial \overline{E}} \simeq \frac{\overline{J}\_i(\overline{\mathbf{E}} + \delta \overline{\mathbf{E}}^{(j)}) - \overline{J}\_i(\overline{\mathbf{E}})}{\mathfrak{a}} \tag{15}$$

with

$$
\delta \overline{\mathbf{E}}^{(j)} = \mathfrak{a} \mathbf{e}\_{j\prime} \tag{16}
$$

where **e***<sup>j</sup>* 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.

#### **3. 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 <sup>Ω</sup> <sup>⊂</sup> <sup>R</sup><sup>3</sup> 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.

**Figure 2.** Realizations of microscopic RVE with various graphene volume fractions: (**a**) 0.53 vol%, (**b**) 0.66 vol%, (**c**) 0.79 vol%, (**d**) 0.92 vol%, (**e**) 1.05 vol%, (**f**) 1.19 vol%, (**g**) 1.32 vol%, (**h**) 1.58 vol% [30].

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

$$\mathcal{W} = \int\_{\Omega} \omega^b(\mathbf{x}) d\Omega + \int\_{\Gamma} \omega^s(\mathbf{x}) d\Gamma\_\prime \tag{17}$$

where density functions *ω<sup>b</sup>* and *ω<sup>s</sup>* are the bulk and surface density functions, respectively, expressed by

$$
\omega^b(\mathbf{x}) = \frac{1}{2}\mathbf{j}(\mathbf{x}) \cdot \mathbf{E}(\mathbf{x}) \,, \quad \omega^s(\mathbf{x}) = \frac{1}{2}\mathbf{j}^s(\mathbf{x}) \cdot \mathbf{E}^s(\mathbf{x}) . \tag{18}
$$

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** *<sup>s</sup>* <sup>=</sup> **PE** <sup>=</sup> <sup>−</sup>**P**∇*φ*, where **<sup>P</sup>** <sup>=</sup> **<sup>I</sup>** <sup>−</sup> **<sup>n</sup>** <sup>⊗</sup> **<sup>n</sup>** is a projector operator characterizing the projection of a vector along the tangent plane to Γ at a point **<sup>x</sup>** ∈ <sup>Γ</sup>, and **<sup>n</sup>** is the unit normal vector to <sup>Γ</sup>.

The nonlinear electric-conduction law including the tunneling effect in the polymer matrix is given by

$$\mathbf{j} = \begin{cases} \mathbf{k}\_0^p \mathbf{E} & \text{if } d(\mathbf{x}) \le d\_{cut\prime} \\ \mathcal{G}(\mathbf{E}, d) \frac{\mathbf{E}}{|\mathbf{E}|} & \text{if } d(\mathbf{x}) > d\_{cut\prime} \end{cases} \tag{19}$$

where *dcut* 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

$$\begin{split} \mathcal{G}(E,d) &= \frac{2.2e^3 E^2}{8\pi h\_p \Phi\_0} \exp(-\frac{8\pi}{2.96 h\_p e E} (2m)^{\frac{1}{2}} \Phi\_0^{\frac{3}{2}}) \\ &+ [3 \cdot \frac{(2m\Phi\_0)^{\frac{1}{2}}}{2d}] (e/h\_p)^2 Ed \exp[-(\frac{4\pi d}{h\_p})(2m\Phi\_0)^{\frac{1}{2}}]. \end{split} \tag{20}$$

Above, Φ<sup>0</sup> and *d* denote barrier height and barrier width, respectively, *e* and *m* are the charge and the effective mass of electron, respectively, and *h<sup>p</sup>* is the Planck constant. Surface current density **j** *<sup>s</sup>* of graphene surface Γ is related to surface electric field **E** *s* [50] through

$$\mathbf{j}^s(\mathbf{x}) = \mathbf{k}^s \mathbf{E}^s,\tag{21}$$

where **k** *<sup>s</sup>* denotes the the surface conductivity of the graphene. This tensor can be related to the conductivity of the volume (multilayer) graphene as:

$$\mathbf{k}^s = h\mathbf{S}^\*, \quad \mathbf{S}^\* = \mathbf{k}^\mathcal{S} - \frac{(\mathbf{k}^\mathcal{S}\mathbf{n}) \otimes (\mathbf{k}^\mathcal{S}\mathbf{n})}{\mathbf{k}^\mathcal{S} : (\mathbf{n} \otimes \mathbf{n})}. \tag{22}$$

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

$$\int\_{\Omega} \mathbf{j}(\phi) \cdot \nabla (\delta \phi) d\Omega - \int\_{\Gamma} \mathbf{P} \nabla \phi \cdot \mathbf{k}^s \mathbf{P} \nabla (\delta \phi) d\Gamma = 0,\tag{23}$$

where *<sup>φ</sup>* <sup>∈</sup> *<sup>H</sup>*<sup>1</sup> (Ω), *<sup>φ</sup>* satisfying the Dirichlet boundary conditions over *<sup>∂</sup>*<sup>Ω</sup> and *δφ* <sup>∈</sup> *<sup>H</sup>*<sup>1</sup> (Ω), *δφ* = 0 over *∂*Ω. The RVE is subjected to boundary conditions (4). The weak form (23) can be solved by the finite-element method.

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

**Figure 3.** Local model: effective flux **J** depends nonlinearly on the macroscopic electric field **E**, volume fraction, and local random distribution of phases.

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.

#### *4.1. Data Generation*

We first define a set of *K* electric-field vectors, **E** *<sup>k</sup>* = *E*1, *E*2, *E*<sup>3</sup> *k* , (*k* = 1, 2, . . . , *K*). The values of **E** *<sup>k</sup>* 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 <sup>α</sup>* were generated and are denoted by Ω*<sup>i</sup> α* , *i* = 1, 2, . . . , *N*.

Then, for each macroelectric field vector **E** *k* , each volume fraction *f <sup>α</sup>* and each realization of microstructure Ω*<sup>i</sup> α* , 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** *k <sup>α</sup>* = <sup>1</sup> *<sup>N</sup>* ∑ *N i*=1 **J** *k α*,*i* .
