Multiscale Analysis of Elastic Properties of Nano-Reinforced Materials Exhibiting Surface Effects. Application for Determination of Effective Shear Modulus

: This work concerns a multiscale analysis of nano-reinforced heterogeneous materials. Such materials exhibit surface effects that must be taken into account in the homogenization procedure. In this study, a coherent imperfect interface model was employed to characterize the jumps of mechanical properties through the interface region between the matrix and the nanoﬁllers. As the hypothesis of scale separation was adopted, a generalized self-consistent micromechanical scheme was employed for the determination of the homogenized elastic moduli. An explicit calculation for the determination of effective shear modulus is presented, together with a numerical application illustrating the surface effect. It is shown that the coherent imperfect interface model is capable of exploring the surface effect in nano-reinforced materials, as demonstrated experimentally in the literature.


Introduction
Composite or composite material consists of at least two elementary phases [1][2][3][4]. This is an artificial assembly, achieved by exploiting the capacity of adhesion between the constituents-namely the matrix and the reinforcements. The composite (or the reinforced matrix) exhibits many mechanical, thermal, and electrical properties, which are beneficial in comparison to those of the pure matrix. Thanks to these various reinforcements, composites are increasingly being used in industry. When one, two or three of the dimensions of the reinforcement are less than 100 nanometers, the resulting material is called a nanocomposite [5][6][7]. Nanocomposites are a novel field of industrial research. Because of the nanometric dimensions of the reinforcement, nanocomposites are highly complex to manufacture, and therefore also to model mechanically. Polymer-based nanocomposites are very frequently used in industry due to the relative ease with which they can be produced using several manufacturing methods, and also the low cost of this type of material. For example, nanocomposites based on thermoplastic polyolefin are commonly used for manufacturing in the automotive industry [8,9]. In this study, we focus on nanocomposites with a polymer matrix reinforced by nanoparticles.
Perhaps the most widely discussed aspect of nanocomposites is the existence of a disturbed region of the matrix surrounding the nanoinclusions, called the interface region, due to the nature of adhesion between the elementary phases in the nanocomposite. This observation has been made in the literature in a range of experimental [10,11] and numerical [12,13] works. At very small scales (in this instance, nanometric), atomic bonds become the characteristics used to describe interactions between the nanofillers and the matrix, or between different nanofiller particles. Thus, multiscale modeling of such a material must take these atomic characteristics into account. A difficulty often encountered in modeling such a material is that the macroscopic properties are highly dependent on the coupling between several parameters, such as the size effect or the aggregation of reinforcements [14].
Modeling of the interface region is crucially important to predict the effective properties of the heterogeneous materials [15,16]. From a general standpoint, the interface region is of non-zero volume fraction in the material, and contains different mechanical properties with the matrix and the inclusion (see [17,18]). Therefore, in general, its thickness may be constant or fluctuating, and its mechanical properties may be homogeneous or non-homogeneous [19][20][21][22][23][24]. In the literature, for approximation of the effective properties of materials with an interface region, there are two main categories of approaches to simulate this region. The first is to approach the effective properties using interface models, where we assume that the interface region is zero volume fraction in the material. The basic idea in constructing interface models is to replace the interface region with an interface of zero thickness (i.e., zero volume fraction) by stipulating that the jumps of physical and mechanical properties through this interface be almost equal to the jumps of those same properties through the interface region, with a fixed error (see [17,18,[25][26][27][28]). In this work, we are interested in the jumps of mechanical properties (i.e., displacement, strain and stress fields, respectively), also known as discontinuities (see [29]).
This study is devoted to a multiscale analysis of the elastic properties of particulate nanocomposites presenting surface effects, together with an application for the determination of effective shear modulus. The paper is organized as follows. Section 2 is devoted to the modeling of the mechanical properties of the interface region using an interface model. The micromechanical homogenization scheme of nano-reinforced materials is introduced in Section 3. Section 4 is given over to explicit calculation for the determination of the effective shear modulus. Finally, Section 5 presents the numerical application, showing surface effects in nano-reinforced materials through the interface model.

Geometric Description of Interface
We use the term interface to denote a surface (or a transition zone) separating two bodies or two parts of a body. In order to describe the topology of the interface denoted by Γ, we consider a solid Ω ⊂ R 3 made up of two materials (obeying any behavior laws), respectively occupying the subdomains Ω (1) and Ω (2) , assumed to be separated by the interface Γ such that (see Figure 1), n(x) Subdomains Ω (1) and Ω (2) separated by interface Γ.
Mathematically, the interface Γ can be defined by the following level-set function f : where the function f is classically assumed to be twice differentiable. The unit vector normal to Γ at point x, denoted by n(x), is then defined by with ∇ f (x) = 0, ∀x ∈ Γ.

Perfect Elastic Interface
In the case of so-called perfect interfaces, the displacement field and stress vector are assumed to be continuous at the interface (see [26,30,31]). We can specify the displacement field u(x) and stress field σ(x) of the domain Ω by introducing and We then write the continuity hypotheses for the point x ∈ Γ • Continuity of the displacement field •

Continuity of the stress vector
[[σ(x)]] n(x) = σ (2) By applying Hadamard's relation for a vector function (see [26]), we obtain the jump of the displacement gradient at the interface and using the mechanical relation ε = 1 2 (∇u + ∇u T ), we then deduce the jump of the strain field at the interface Γ We introduce the surface strain tensor ε s by and by introducing the result of Equation (10), we obtain where the tensorial relation , with A being a second-order tensor and a and b being two vectors, is introduced in Gu [30], we then find that [[ε s ]] = 0, using the property of two second-order projection operators, t T = t, and t n = 0 (see [30]). This result shows that the surface strains are continuous across the interface Γ. We introduce the following tensorial relation where P is the operator defined by Equation (A4). We can then rewrite the result of Equation (13) [[ε s ]] = t ε (2) This is the first relation of continuity at the interface Γ for the strain field. Then, the continuity of the stress vector implies that where P ⊥ is the operator defined by Equation (A5). Finally, the relationships of continuity across the perfect interface Γ are More information on projection operators and elastic moduli can be found in Appendixes A and B.

Imperfect Elastic Interface
In several situations, the perfect interface hypothesis is no longer adopted (see [26]). The most important reason is that the energy of the interface is not negligible compared to the energy of the domain (see [31]). Because of this phenomenon, there is a major increase in effective properties. In this section, in the wake of the previous section's presentation of the problem of the surface elasticity so that the jumps of mechanical values through the interface are almost equal to the jumps of the same values through the interface region, we present a general and systematic approach based on asymptotic analysis to establish imperfect interface models (see [29,30]). Figure 2a shows a simple depiction of the three-phase configuration, the interface region of thickness t(x), occupying the domain Ω * (0) , separating the two sub-domains Ω * (1) and Ω * (2) (see [30]). The basic idea to construct an interface model is to replace the area of Ω * (0) in the three-phase configuration with an interface Γ of zero thickness in the two-phase configuration by stipulating that the jumps of mechanical properties through the interface Γ are almost equal to the jumps of the same properties through the surface Γ 1 and Γ 2 , with a fixed error (see [26]). Γ is called an imperfect interface. The two-phase configuration is represented by Figure 2b, where at the interface Γ, there are discontinuities of displacement, strain and stress fields. These discontinuities are called interface conditions (see [17]); they are added to the classical homogenization equations to solve the problem of heterogeneous material with imperfect interfaces.
In the literature, there are several types of imperfect interface models, such as the imperfect linear spring interface (see [28]), and coherent imperfect interface (see [17,18,26,31]). These two types are classified essentially by discontinuities of mechanical fields across the imperfect interface. We would like to focus here on the coherent imperfect interface model. To describe the mechanical properties through the coherent imperfect interface, we consider a domain of heterogeneous material having a coherent elastic interface between the matrix and the inclusion. We then introduce the displacement fields u s , strain ε s , and stress σ s of the coherent interface. Mechanical properties across the coherent imperfect interface Γ are represented below, using the terminology introduced in [17,18,31]. • Continuity of the displacement field Similarly, when the displacement field is continuous at the interface, we obtain the continuity relationship of surface deformation at the interface using Equation (17) • Discontinuity of the stress vector (see [31]) with div s being the surface divergence operator, and σ s the stress tensor of the interface. The stress vector is no longer continuous at the interface. The following relation represents this discontinuity (see Equation (18)) These relations are used in the same way as the interface conditions in the section discussing the homogenization scheme.

Homogenization of Heterogeneous Materials Presenting an Interface
It is assumed that there is a separation of scale, such that the homogenization operator becomes deterministic. Therefore, a classical (pseudo-) analytical (see [31]) or numerical (see [32]) model could be used to homogenize the microstructure with an imperfect interface. In this work, a pseudo-analytical estimation system was applied to solve the homogenization problem.

Description of Microstructure
The microstructure in question is presented in Figure 3. Each phase-the inclusion (denoted by i), the matrix (denoted by m), and the interface (denoted by s)-is composed of linear, homogeneous and isotropic material. The elastic constants of the inclusion, the matrix and the interface are denoted by (κ i , µ i ), (κ m , µ m ), and (κ s , µ s ), respectively. In this study, the inclusion occupies the volume of radius r 1 and the matrix occupies the volume, the radius of which ranges from r 1 to r 2 .
where in the isotropic case, the elasticity tensor c m or c i could be specified: with p = m(matrix) or p = i(inclusion).

Interface Condition
The elasticity model of the interface is characterized by the following equation where in the isotropic case, the elastic constants of the interface are such that (λ s , µ s ), the constitutive law is written as follows (see [18]) where 1 is the identity tensor in two-dimensional space. We represent the constitutive law in a matrix form as follows (see [31]) with the relationship κ s = 2(λ s + µ s ) (see [17]). We subsequently detail the discontinuity of the stress vector at the interface represented by Equation (21). We introduce the result of Povstenko's work (see [33]). By studying the balance of the interface Γ, we obtain the explicit form of the discontinuity (see [17,31] for more details): where the surface divergence operator div s is calculated in spherical coordinates by the following formula (see [26]) This interface condition is added to the classical homogenization equations in the following sections.

Generalized Self-Consistent Scheme
We propose, in this work, to use the generalized self-consistent scheme to solve the homogenization problem. This scheme is detailed explicitly in the work of Christensen and Lo [34] to solve the thermomechanical problem. They described the model as a heterogeneous sphere integrated into an equivalent infinite medium V, this medium being the effective domain. The heterogeneous sphere was made up of the matrix and inclusion. The interface between the matrix and the equivalent medium was considered perfect. The interface between the matrix and the inclusion was considered to be coherent imperfect. The diagram of the generalized self-consistent system is shown in Figure 4, where the letter "i" presents the inclusion, "m" the matrix, and "e" the equivalent effective medium.
The equivalent effective medium of elastic moduli κ eff and µ eff is characterized by the isotropic constitutive law σ e = c e : ε e . Below, we represent the self-consistent condition to find the effective moduli of the material. For the domain V of free energy U 0 , we impose the homogeneous conditions on the boundary in strain and in stress, respectively. We set u 0 as the displacement field imposed on the boundary corresponding to the strain condition, and t 0 the imposed stress vector corresponding to the stress condition. First, if we imply the condition in homogeneous strain with the boundary of the form and after having integrated the heterogeneous sphere in the equivalent medium, we obtain the energy of the system On the other hand, if we imply the homogeneous condition in stress to the boundary and in the same way, we find the energy of the system We then obtain the self-consistent condition from two Equations (33) and (35) ∂V t u 0 − t 0 u dS = 0.
This is the condition for determining the effective modulus of the heterogeneous material, κ eff and µ eff , using the generalized self-consistent system. We will apply it in the next section to determine the effective shear modulus.

Boundary Conditions for Determination of Effective Shear Modulus
In this part, we represent the homogeneous strain condition at the boundary to determine the effective shear modulus of the material. In this case, as pure shear is imposed, the homogeneous strain boundary condition ∂V is represented by where γ 0 is a constant. In spherical coordinates, this condition is written By imposing the self-consistent condition, we obtain the effective shear modulus of the material, which is a function of κ s , µ s , κ i , µ i , κ m , µ m , and r 1 , r 2 : The explicit calculation to determine the effective shear modulus is presented in the next section.

Explicit Calculation to Determine Effective Shear Modulus
In this section, we present the steps to solve the homogenization problem seeking the effective shear modulus of the material. Figure 5 shows the flowchart of the proposed approach to determine the effective modulus of the heterogeneous material. Details of the calculation are presented in the following. The form of the solution for the displacement field of each phase is conjectured from the form of the boundary condition (see [31]). By applying the homogeneous strain condition to the boundary presented by Equation (38), we propose the form of the displacement solution

Displacement field of inclusion, matrix and equivalent phases
where the index (p) means that we are dealing the inclusion phase (p = 1 or i), matrix (p = 2 or m), or equivalent medium (p = 3 or e), respectively, and f (p) and g (p) are unknown functions to be solved. Then, we derive the displacement field to obtain the strain field Then, by Hooke's law, we obtain the stress field of each phase We then obtain the equilibrium equations for each phase, in spherical coordinates These equilibrium equations give a system of second-order differential equations. We solve this system by indicating the precise form of the two functions f and g for each phase, and search for the constants appearing by the boundary conditions and the interface conditions. The system of differential equations obtained is of the form where a (p) 2 are the constants as a function of the moduli λ p and µ p of each phase To solve this system, we consider its characteristic equation of the form: with We then find that m (p) and n (p) are constant in the isotropic case. They are independent of the elastic moduli. Then the solutions of Equation (51) are −4, −2, 1 and 3. We therefore deduce the form of the two functions f and g for each phase: Below, we seek to simplify the solutions. For the equivalent medium, when the radius r tends towards infinity, it is necessary to satisfy the homogeneous strain condition imposed on the boundary; then, we obtain A = γ 0 . The form of the solutions for the equivalent medium is therefore: Similarly, when the radius r tends toward zero, we consider that f and g are finite. Hence, we obtain the form of solutions for inclusion We still have to solve a system of eight equations with eight unknown constants. The conditions at the interface are, at r = r 1 : where the first two equations are the conditions of continuity of the displacement field, described by Equation (19). The following strain conditions are deduced from Equation (20), the surface part of the strain tensor is continuous across the interface. Finally, the relations of stress discontinuities are represented in Equation (21), with the constitutive law of the interface being given in Equation (28).
At the interface at r = r 2 , as this interface between the matrix and the equivalent medium is perfect, we have We obtain a system of linear equations, which allows us to find the constants sought. This system is of the form where the vector of the unknowns b then the vector c c = [0 0 0 0 − r 2 10 and the matrix A is given by with the components of matrix A are functions of the elastic moduli of each phase, and of the geometry of the field. Finally, according to the terminology of the self-consistent condition presented in Section 3.3, we obtain an equation that gives the effective shear modulus µ eff as a function of κ s , µ s , κ i , µ i , κ m , µ m , and r 1 , r 2 .
Equation (68) is solved by using the consistency condition for the equivalent medium presented previously in Section 3.3. The effective shear modulus of the heterogeneous material is calculated as the positive root of a quadratic equation derived from the above-mentioned consistency condition.

Input Parameters for Numerical Application
Considering the nanoinclusion of elastic properties E i = 40 GPa and ν i = 0.22, and the matrix with E m = 4 GPa and ν m = 0.37 (typical values for polymeric materials reinforced with silica nanoparticles [12,35]), we calculate the elastic moduli (κ i , µ i ), and (κ m , µ m ) by the relation We obtain κ i = 23.8095 GPa, µ i = 16.3934 GPa and κ m = 5.1282 GPa, µ m = 1.4599 GPa. In the following, we offer the choice of surface elastic properties κ s and µ s for numerical application. The elastic surface moduli have already been calculated by several researchers, by the method of stress simulation at the interface (see [36]), or by atomistic simulation (see [37]), or by simulating a fine or stiff interphase using the asymptotic approach (see [17,26,29,38]). In this case, we use the approximation proposed in [17,26,29,38] for a coherent imperfect interface. We propose a Young's modulus ratio E s E m = 1000, the surface Poisson's ratio ν s = 0.25, and the constant thickness interface h = 0.001 nm = 10 −12 m. Such a very thin and very stiff layer between the matrix and the nanoinclusion is considered in order to apply the linearly elastic interface stress model earlier proposed by Benveniste [29]. We present the formula that gives the values of elastic surface moduli as a function of E s , ν s and h: We obtain κ s = 4.2667 N/m and µ s = 2.6667 N/m. In the following section, we offer a numerical simulation of the surface effect on the effective shear modulus using the micromechanical system.

Surface Effect on Effective Shear Modulus
Using the input parameters introduced above, we propose here to illustrate the surface effects on the effective shear modulus (see Equation (68)). We use the notation µ eff imperfect.interface and µ eff perfect.interface for the effective shear modulus using imperfect and perfect interfaces, respectively. Figure 6 shows the surface effects of the heterogeneous material containing nanoinclusions. We can see that there is an increase in effective properties, and the smaller the radius of the fillers, the greater the difference in the property. If the radius of inclusion r 1 ≤ 5 nm, this increase is very marked. If the radius of inclusion r 1 > 50 nm, the difference in effective properties is very small, and we can therefore consider that there is no longer any surface effect. In other words, a classical solution can be applied in this case.
On the other hand, by using the present approach, if the radius of inclusion tends to an infinitesimal value, the specific area of the reinforcement phase is then increased to an infinite value. Such an augmentation leads to an increase in surface interactions between the matrix and the reinforcement, and therefore an increase in the effective properties. However, this point should be revisited from both experimental and physical points of view at the smallest scales. For instance, in Wang and Ye [39], the radial displacement field due to surface effect has been derived and showed that the effect of an imperfect interface is only profound within the region near the interface. From an experimental point of view, such surface effect on the effective elastic properties has been clearly demonstrated, especially in the case of polymer-based particulate nanocomposites. A typical example can be found in [40], where a polysiloxane matrix was reinforced by silica nanoparticles of three different sizes −15, 35 and 60 nm, respectively. The surface effect on the Young's modulus of nanocomposite was clearly observed by comparing the results for the nanofillers with diameters of 15 and 35 nm, respectively. For instance, the effective Young's modulus was about 4.3 and 3.95 GPa, using nanofillers of 15 and 35 nm, respectively (at the same volume fraction of reinforcement of 10%), where the Young's modulus of the pure polymer was 3 GPa. However, such surface effect was not observed with 60 nm particles. A possible explanation is that for this size of reinforcement (i.e., larger than 50 nm in diameter), the specific area of the nanoparticles is decreased, leading to a decrease in surface interactions and therefore a decrease of the macroscopic properties. Finally, other experimental results on surface effect are available such as [41][42][43][44].

Conclusions
In this work, we are interested in the multiscale analysis of heterogeneous materials presenting surface effects. A model of the coherent imperfect elastic interface was used, together with a generalized self-consistent micromechanical scheme for determining the effective elastic properties of the nano-reinforced materials. In this study, the calculation for effective shear modulus was detailed. Numerical application showed that the imperfect interface model was capable of taking the surface effect as a function of nanofiller size into account.
In further studies, the surface energy should be discussed, to describe whether or not the surface elasticity tensor is positive-definite. In addition, both experiments and numerical simulations (for instance, in Molecular Dynamics) demonstrated the existence of a perturbed area at the boundary between the inclusion and the matrix phase. From a probabilistic standpoint, the elastic properties of that perturbed area randomly fluctuate from point to point. Consequently, a probabilistic investigation should be considered in future research for describing the surface effect in nano-reinforced materials.
Funding: This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Conflicts of Interest:
The authors declare that there is no conflict of interest.

Appendix A. Projection Operators
Consider the two projection operators defined by n(x) = n(x) ⊗ n(x), (A1) where i denotes the second-order symmetric identity tensor. These operators n(x) and t(x) represent the projections in the normal direction n(x) and in the tangential direction to Γ at point x, respectively. Note that n(x) and t(x) are two orthogonal and complementary projectors, such that n 2 = n, t 2 = t, n T = n, t T = t, n t = t n = 0, n + t = i, where T represents the transpose operator. From these two projection operators (of the second order), we can define the following fourth-order projection operators with index notation: where these two projection operators satisfy the following properties P 2 = P , P ⊥ 2 = P ⊥ , P T = P , P ⊥ T = P ⊥ , P : P ⊥ = P ⊥ : P = 0, P + P ⊥ = I, with I being the fourth-order identity tensor. These properties show that P and P ⊥ are orthogonal and complementary.