6.2.1. Tunneling effect

Tunneling electric conductivity is a quantum phenomenon that allows electric conductivity across small isolating barriers, like thin polymer interphases between two highly conducting fillers. It has been shown that electric tunneling effect plays an important role in explaining very low percolation thresholds, as well as nonlinear electric conduction (Oskouyi et al. 2014, Sheng 1980) in such materials, especially when the characteristic distances between graphene sheets reduce to the order of nanometers. In the works of Otten and van der Schoot (2011) and Ambrosetti et al. (2010), the tunneling effect has been taken into account to analyze the percolation threshold of nanocomposites with polydisperse nanofillers, where the former focuses on carbon nanotubes and the latter extends to various particle shapes.

In a first approximation, the polymer–graphene nanocomposite can be seen as a juxtaposition of elementary graphene–polymer–graphene units (see Figure 6.1(a)). In order to understand the nonlinear electrical behavior of the nanocomposite, we first focus on the conduction mechanism of this unit. When the distance between two graphene sheets is nanoscaled, the polymer matrix can be seen as a dielectric film between two graphene electrodes. There is a strong analogy with a metal–insulator–metal device with tunnel junction. It is known that there are two types of conduction mechanisms in dielectric films (Hamann et al. 1988), which are interface-limited conduction mechanisms and bulk-limited conduction mechanism, respectively. The former depends on the electric properties at the polymer–graphene interface, while the latter depends on the electric properties of the polymer. In our work, we focus on solving the general problem and neglect the bulk-limited conduction mechanism, except Ohm’s law in varying polymers.

There are several interface-limited conduction mechanisms (Chiu 2014): (1) Schottky emission; (2) thermionic-field emission; (3) Fowler–Nordheim tunneling; (4) direct tunneling. Among them, the former two are only applied at high temperatures, when the electrons can get enough thermal energy to overcome the barrier so that the electronic current can flow in the conduction band. The latter two occur when the barrier is thin enough to permit its penetration by the quantum tunnel (Simmons 1971). To simplify the problem, we assume that the system is studied at low temperatures so that the thermal current can be neglected, and restrict the electron transportation between electrodes to the tunneling effect.

Figure 6.1 depicts the evolution of conduction band edge energy, E_c, which is at an equilibrium state across the graphene–polymer–graphene device in Figure 6.1(a).

For such a 1D model of two close graphene sheets separated by a polymer matrix, the electric constitutive law is expressed in the polymer in the form

[6.1]

where I denotes the electric current, E is the electric field and G (E, d) is a nonlinear function of the electric field E and the distance between the two graphene sheets d. An explicit formula for the electric tunneling effect through a potential square barrier was first derived by Simmons (1963) as:

[6.2]

Here, the tunneling current I is a nonlinear function of electric field E and depends on three parameters: the height Φ₀, width d of the barrier and the effective mass m of the electrons in graphene. This expression is obtained under the assumptions that the temperature is very low (i.e. T ≈ 0 K); the electric field E varies slowly over the width of the potential barrier (i.e. E ≈ V/d, and the rectangular barrier becomes trapezoidal under the effect of the electric field (see Figure 6.1(b))), and the trapezoidal barrier is approximated by the rectangular barrier with the same average height. This theory has been used and validated experimentally in the metal–insulator–semiconductor (Maity et al. 2016) and metal–insulator–metal capacitor field (Krishnan et al. 2008). Note that the current only comes from the electrons that cross the potential barrier. In Simmons’ model, the leakage current is not taken into account. This can be covered in the 3D extension of the conductivity model (see section 6.2.2). However, since we consider the bulk current in addition to the tunneling current, the leakage current will not affect the effective current in the composite.

6.2.2. Nonlinear electrical conduction model at the RVE scale

In this section, we describe a continuum RVE model for the graphene–polymer nanocomposite. At this scale, the graphene sheets are explicitly described at the scale of a quite large cluster of graphene nano platelets (see Figure 6.2) composed of a polymer matrix and randomly distributed and oriented graphene sheets, modeled as highly conducting imperfect interfaces (Yvonnet et al. 2008).

The electric tunneling effect between graphene sheets described below is taken into account in the continuum model.

In this model, the total electric power of the system, , is defined by

[6.3]

where Ω ⊂ ℝ³ is the domain defining the RVE, Γ denotes the surfaces associated with graphene collectively and the density functions ω^b and ω^s are expressed by

[6.4]

where E_i(x) = −∇_iø(x) = −ø_i(x) is the electric field, j(x) is the current density vector and ø(x) is the electric potential, where:

[6.5]

with as the electric conductivity tensor of the polymer matrix when neglecting the tunneling effect. In the above, the superscript s denotes surface quantities, for example, j^s is the surface current density. The surface electric field is defined with respect to its bulk counterpart as: with P_ij (x) = S_ij − n_i (x)n_j (x) as 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 local constitutive relationships relating j and j^s with E are defined in the graphene sheets by

[6.6]

where k^s denotes the surface electric conductivity of graphene and is dependent on the thickness t through:

[6.7]

In [6.7], k^g denotes the second-order electric conductivity tensor of the bulk graphite and n is the normal vector to the graphene sheet (for more details, see Lu et al. (2018), Yvonnet et al. (2008)).

[6.8]

where d_cut is a cut-off distance above which the tunneling effect can be neglected and G is a 3D extension of [6.2] defined by:

[6.9]

where Φ₀ is the energy barrier height that the electrons cross and h, e and m denote Plank’s constant, the charge of an electron and a material parameter, and d is a distance function defined in 3D (see more details in Lu et al. (2018)).

Minimizing [6.3] with respect to the displacement field, and using equations [6.5]–[6.8], we obtain the weak form, which can be solved by using the finite element method (see more details in Lu et al. (2018)):

[6.10]

where δø ∈ H¹ (Ω), δø = 0 over ∂Ω and ø ∈ H¹(Ω), ø satisfying the periodic boundary conditions over ∂Ω

[6.11]

and where is a periodic function over Ω, such as .

The effective electric conductivity tensor is defined as:

[6.12]

where is the effective current density whose components are expressed by:

[6.13]

and is the effective electric field given by:

[6.14]

6.3.1. Effect of barrier height on the percolation threshold

Using the above model, we study the effects of barrier height between graphene and various polymers on the nonlinear response of the nanocomposite. In our simulations, we have considered that the percolation threshold corresponds to a sharp variation of the effective conductivity above 10^–8 S/m. In Figure 6.3, the effective conductivity component is plotted as a function of the graphene volume fraction for the values Φ₀ = 0.17 eV, 0.3 eV, 1.0 eV and without tunneling effect. The aspect ratio is η = 50 and applied electric field . It should be noted that for small barrier height (0.17 eV and 0.3 eV), the computation results of electric conductivity, taking into account the tunneling effect, are much larger than the predictions without tunneling effect, while for Φ₀ = 1.0 eV, the effective conductivity characteristics exhibit no obvious difference, either with or without tunneling effect.

The percolation thresholds corresponding to Φ₀ =0.17, 0.3, 1.0 eV are 0.79 vol%, 0.92 vol% and 1.58 vol%, respectively. We can see from these simulations that the lower the barrier height is, the lower the percolation threshold. Experimental results reporting percolation thresholds for different polymer matrices and graphene types can be found in previous studies (Stankovich et al. 2006, Jiang et al. 2009, Ansari and Giannelis 2009, Liang et al. 2009).

6.3.2. Effect of graphene aspect ratio on the percolation threshold

Next, we use our numerical model to estimate the percolation threshold f* of the nanocomposite as a function of the graphene volume fraction and aspect ratio for graphene sheets. In Figure 6.4, the effective conductivity tensor component is computed for different aspect ratios η = 20, 50 and 100 as a function of volume fraction.

The side length of the graphene is fixed, and the various aspect ratios are obtained by changing the thickness of the graphene sheets. For each case, the values are averaged over 30 realizations of random distributions of graphene sheets within the RVE. As the behavior of the composite is nonlinear, the results of conductivity are presented for a fixed electric field . The barrier height is Φ₀ = 0.17 eV. According to our numerical simulations, the effective conductivity clearly depends on the aspect ratio, a larger η provides lower percolation threshold. The obtained percolation thresholds for η = 20,50, and 100 are 1.65 vol%, 0.79 vol%, and 0.33 vol%, respectively.

6.3.3. Effect of alignment of graphene sheets

Next, we evaluate the effect of alignment of graphene sheets on the effective conductivity of the composite. For this purpose, we consider, on the one hand, a microstructure with graphene sheets whose positions and orientations are randomly distributed, and, on the other hand, a microstructure where the positions of the graphene sheets are randomly distributed, but the orientation is fixed. Each point corresponds to the mean value over 30 realizations (see realization examples in Figure 6.5).

The results are presented in Figure 6.6. To clearly evidence the anisotropy, we have plotted the three components of the effective conductivity tensor. The parameters are E = 1.25 × 10^–3 V/nm, Φ₀ = 0.17 eV and η = 50.

As expected, the numerical model clearly captures the anisotropic behavior of aligned graphene sheets (see Figure 6.6(b)). Another conclusion is that aligning the graphene sheets does significantly increase either the maximum effective conductivity in the direction, normal to the graphene sheets, or the percolation threshold, as compared to randomly oriented sheets. However, the increase in conductivity after the percolation threshold is sharper in the case of aligned graphene sheets.

Finally, in Figure 6.7 we depict the current density field in the polymer matrix when tunneling effect is taken into account for aligned and randomly oriented graphene sheets to show the percolation path of electric current in both of these configurations. The parameters are η = 50 and Φ₀ = 0.17 eV, f = 1.05 vol%. The electrical field is applied in the X-direction.

6.3.4. Comparison between numerical and experimental results

In the following, the obtained numerical results are compared with available experimental data by Stankovich et al. (2006) and Zhang et al. (2010), respectively. The details of the different material and numerical parameters can be found in Lu et al. (2018). The results presented in Figures 6.8 and 6.9 show encouraging agreement between experiments and numerical simulations, even though some discrepancies exist. These disagreements can be explained as follows.

First, the dispersion state of graphene is hard to control in the experimental preparation, and the aggregation of graphene in the microstructure may lead to the heterogeneity of the macroscopic property. However, in our present work, the numerical simulation is based on randomly distributed graphene sheets. Moreover, the numerical simulation ignores the defects in the nanocomposite, especially the ones on the interface between graphene and the polymer. It should also be noted that we focus on the simulation of nanocomposites with small graphene volume fraction due to the extremely small percolation threshold that we are interested in. The simulation of electric conductivity at high graphene volume fraction calls for a large computation cost and has quite a long computation time due to the refined meshes, so we did not take this further.

6.4.1. Construction of the surrogate model based on ANN: strategy

The key idea of this work is to construct a surrogate model based on neural networks in a multilevel finite element framework. In the present method, called neural network computational homogenization (NN-CH), we use a multi-layer perceptron (MLP) to construct the nonlinear relationship between the (input) effective electric field vector and the (output) effective electric current density vector using the model described in section 6.2.2. As this relationship is not known in advance, each set of values for requires a nonlinear finite element problem to be solved to obtain the corresponding , as defined in section 6.2.2. To avoid these costly calculations at the macro-scale, we propose here to define this relationship through the ANN surrogate model. For this purpose, a set of data consisting of N couples of vectors and their corresponding are provided to the neural network.

The procedure for constructing the ANN-based surrogate multiscale model is illustrated in Figures 6.10 and 6.11 and outlined in the following steps:

1. 1) generate N number of electric field vectors to be used for the MLP training;
2. 2) run N FEM nonlinear RVE calculations, for applied effective electric field vectors and compute the corresponding effective current densities , as described in section 6.2.2;
3. 3) train the MLP using as inputs and as the corresponding outputs;
4. 4) check for accuracy of the trained MLP;
5. 5) if the error (using appropriate metrics) is larger than a tolerance, add new input (through steps 1–2) to the set of data and perform steps 3–5 again.

6.4.2. Structural application

One advantage of the NN-CH approach is that once the surrogate model is constructed, many structural problems can be solved without new calculations on the RVE and taking advantage of the off-line calculation costs. In this next example, we consider the structure described in Figure 6.12. The side lengths of the domain are L_x = 16 μm, L_y =8 μm, L_z = 2 μm. The diameter of the hole in the middle of the cuboid is D = 4 μm.

A mesh of 1820 tetrahedral elements is constructed. The boundary condition is ø(x = 0 μm) = 0 V, ø(x= 16 μm) = 160 V on the corresponding end surfaces. The current density on the other boundary surfaces.

Results for an RVE of 0.53 vol% are shown in Figure 6.13. For NN-CH, the calculation took 267.45 s with 12 parallel workers, against 126 h for FE². A comparison of the different fields by FE² is provided in Figure 6.13, showing a very good agreement. Then, we show that accurate nonlinear models based on RVE can be used in structural engineering applications at reasonable costs.

6.5.1. Mechanical modeling

In this section, we use the theory of nonlinear continuum mechanics at finite strains, details of which can be found in previous studies (Gurtin et al. (2010), Tadmor et al. (2012), among others). We consider the RVE described in Figure 6.2. Here, the RVE in the reference configuration Ω₀ ∈ ℝ³, and the spatial configuration Ω_t ∈ ℝ³.

The two sides of the interface are denoted by and . The unit vector normal to the interface in the reference configuration is N(X). The displacement of the bulk, and the two sides of the interface are u, u^– and u⁺, respectively. The current positions x of the material particles at the initial position X are defined by x = X + u for the bulk. x^– and x⁺ are the current position for the two sides of interface. The graphene sheets are modeled as general imperfect interfaces (Chatzigeorgiou et al. 2017, Javili et al. 2017), satisfying [[u]] = u⁺ – u^– ≠ 0 and [[x]] ≠ 0; the traction through the graphene surfaces is also discontinuous, [[t]] = t⁺ – t^– ≠ 0, where t⁺ and t^– are, respectively, the traction on and associated with the normal N(X).

Here, we follow the theory of imperfect interface at finite strains, developed by Javili et al. (2017). The internal virtual work, δW_int, is given in reference configuration by the contributions of polymer bulk and graphene sheets, which are modeled by imperfect interfaces:

[6.15]

The polymer bulk contribution is

[6.16]

where P is the first Piola–Kirchoff stress, such as t = P · N and δF = ∇_x δu is the gradient of the virtual displacement Ju with respect to the Lagrangian coordinate X. It should be noted that in this chapter, we defined (A · b)_i = A_ijb_j , a.b = a_ib_i, (AB)_ij = A_ik B_kj.

The contribution of imperfect interfaces is given by

[6.17]

where P^(s) is surface first Piola–Kirchof stress and δF^(s) = ∇_X δu^s · denotes the surface gradient of the virtual displacement, where δu^s = 1/2(δu⁺ + δu^–) is a virtual displacement field on the interface and is the projector into the tangent plan of the interface in the reference configuration. {{t}} = 1/2 (t⁺ +1^–) is the average of traction. Note that the expression of internal virtual work with imperfect interfaces is true by constraining the interface motion to the midplane (Javili et al. 2017), that is, by imposing the definition of the interface displacement u^s = 1/2(u⁺ + u^–).

Using the second Piola–Kirchhoff stress tensor S = F^–1P (respectively, the surface second Piola–Kirchhoff stress tensor, S^s = (F^s)^–1 P^s with , we finally obtain the following expression of internal virtual work (for a configuration in static equilibrium):

[6.18]

[6.19]

where δ∈ (respectively, δ∈^s) is the variation of the symmetric Green–Lagrange strain tensor, (respectively, the variation of the surface Green–Lagrange strain tensor, .

6.5.2. Constitutive laws

We choose the isotropic Saint Venant–Kirchhoff model for the bulk part, which is an extension of the linear elastic material model. The second Piola–Kirchoff stress is given by:

[6.20]

where the Lame’s coefficients are λ^(b) = 6890 MPa and μ^(b) = 680 MPa (Lu 2017). They are identified by deformation of MD simulation box of pure PE following classical procedures (Theodorou and Suter 1986).

As for the bulk, for sake of simplicity, we assume that the behavior of the imperfect interface is reversible. The surface elastic behavior of graphene is assumed to be isotropic along its tangent plane to the surface, such that the surface second Piola–Kirchoff stress is also given by the Saint Venant–Kirchhoff model:

[6.21]

where the surface Lame’s coefficient λ^(s) = 19.0 N.m^–1 and μ^(s) = 18.7 N.m^–1 are identified by the MD simulation (Lu 2017).

The expression of traction {{t}} is assumed to be aligned with the displacement jump [[u]] and is given by:

[6.22]

where g_cz (·) is the function identified by the MD simulations in equation [6.27].

6.5.3. Weak form of mechanical problem

In order to relate the micro stress and strain fields to the imposed macroscopic strain, the effective quantities are defined as in Javili et al. (2017):

[6.23]

and

[6.24]

where and are the effective deformation gradient and effective first Piola–Kirchoff stress, respectively. The effective Green-Lagrange strain tensor and the effective second Piola–Kirchoff stress are defined by:

[6.25]

To impose that the incremental internal virtual work at the micro-scale, δW_int, is equal to the macroscopic internal virtual work, , we chose periodic boundary conditions to satisfy the extended Hill–Mandel condition (Javili et al. 2017). The weak form associated with the mechanical problem at micro-scale is given by: find u ∈ H¹(Ω₀), satisfying the boundary conditions over δΩ₀ with ũ periodic, such that

[6.26]

for all .

We use FEM to discretize the solution space, linear tetrahedrons for the bulk part and linear triangles for interfaces. A Newton–Raphson procedure is used to solve this nonlinear problem step by step for small increments of effective strain gradient.

6.5.4. Identification of cohesive zone model

To study the separation in opening mode, graphene is moved in successive steps of 0.5Å along the Z direction by following a minimization procedure. The graphene atoms and bottom layer of the polymer are kept fixed (see Figure 6.14(a)). The separation process is depicted in Figure 6.14. The polymer chains undergo stretch at the beginning along the Z direction, then form highly oriented structures, called fibrils or nano-fibrils. Voids appear between the fibrils during the decohesion. This deformation mechanism observed during the simulation is similar to the nano-crazes of some semi-crystalline polymers, such as polybutene (Thomas et al. 2009, Detrez 2008). The size of the void grows along the separation direction and the extended chains slide along the graphene sheet to increase the fibrils, as described in the reference review paper by Kramer and Berger (1990). It should be noted that the separation is controlled by the chain desorption at the graphene surface by sliding, which is dominated by Van der Waals interactions.

The average force of polymer on graphene was monitored, from which we can get the normal traction force t_n of the cohesive zone as a function of the displacement of the graphene layer, [[u_n]], as shown in Figure 6.15. First, the force varies linearly with the displacement of the graphene sheet. Then, the curve bends to reach a maximum at 0.7 nm, called the yield threshold. This phase corresponds to the nano-fribils creation and to the cavity initiation. Once the yield threshold is crossed, the force decreases with the displacement of the graphene sheet. During this phase, the chains slip on the graphene sheet to feed the fibrils. It is interesting to note that the observed softening is related to the reduction of the contact area between the polymer chains and graphene. We propose fitting the MD results with the following empirical model:

[6.27]

The above model is incorporated in the mechanical model described in section 6.5.1 and is coupled with the nonlinear electrical model described in section 6.2.2.

6.5.5. Evolution of electrical properties under stretching of the composite

Using the above electro-mechanical model of nanocomposites, we study the impact of both stretching and decohesion at the graphene polymer interface on the electric conductivity. Indeed, the conductivity is controlled by the tunneling effect, which depends strongly on the distance between graphene. Therefore, we impose an macroscopic elongation, , on the RVEs, from . Due to the large computational times, only one RVE microstructure is randomly studied for each graphene volume fraction.

The deformed microstructures are stored for each 1% increment of deformation and the distance function, d(x), is updated. Introducing the new distance function, the electrical conductivities at different effective strain are shown in Figure 6.16. The boundary between the insulator and conductor is defined as 10^–8 S/m (Singh 2011), below which the material is supposed to be the insulator. In the other situation, it behaves like a conductor. Focusing on the electrical conductivity along the direction of deformation , in Figure 6.16(a) we can observe that the mechanical deformation has little effect on the electrical conductivity of the nanocomposites when the graphene volume fraction is below the percolation threshold, f < f_c = 0.52 vol%. When the graphene volume fraction is above the percolation threshold (f > f_c = 0.52 vol%), the electrical conductivity decreases with the applied elongation, but it should be noted that the nanocomposites remain as the conductor. However, if the graphene volume fraction is around the percolation threshold (f ≈ f_c = 0.52 vol%), a sharp decrease in the electrical conductivity can be seen when the nanocomposites are subjected to strain, which is regarded as a transition point from conductor to insulator. For instance, with 0.66 vol% graphene the transition point of the sample is , and with f = f_c = 0.52 vol% graphene it is . We note that this drop in conductivity does not occur for any curve above a volume fraction f ≥ 0.79 vol%. However, the elongation would not affect the electrical conductivities of the nanocomposites in the transverse directions and , which can be seen from Figures 6.16(b) and (c). This is probably a consequence of the fact that the simulations for each volume fraction have only been performed on one realization, and that some configurations might be more favorable to this effect. The results demonstrate that the graphene nanocomposite with certain graphene volume fraction (around percolation threshold) may lead to the existence of conductor to insulator behavior under finite strain. But it is also dependent on the microstructure of nanocomposite, and it only happens to those whose current path formed by graphene can be destroyed during the tension.

After observing this typical conductor-to-insulator transition for the composite with f = 0.66 vol% graphene by the proposed model, we compare the effective conductivity with the results that are estimated without considering the decohesion between graphene and the polymer matrix (i.e. we prescribe [[u]] = 0). It can be seen in Figure 6.17 that by neglecting the cohesive interface, the transition point increases from , which shows the important role of decohesion at the interface in predicting the piezoresistivity properties of polymer graphene nanocomposites. It is interesting to note that it is theoretically possible to design a composite that can go from conductor to insulator by varying the applied strain on the system. This transition can be induced mainly by the decohesion for weak interfaces, or only by strain for a stronger interface, but for a more important applied elongation.