Enhanced computational homogenization techniques for modelling size effects in polymer composites

Abstract

Several experimental investigations corroborate nanosized inclusions as being much more efficient reinforcements for strengthening polymers as compared to their microsized counterparts. The inadequacy of the standard first-order computational homogenization scheme, by virtue of lack of the requisite length scale to model such size effects, necessitates enhancements to the standard scheme. In this work, a thorough assessment of one such extension based on the idea of interface energetics is conducted. Systematic numerical experimentation and analysis demonstrate the limitation of the aforementioned approach in modeling mechanical behavior of composite materials where the filler material is much stiffer than the matrix. An alternative approach based on the idea of continuously graded interphases is introduced. Comprehensive evaluation of this technique by means of representative numerical examples reveals it to be the appropriate one for modeling nano-composite materials with different filler-matrix stiffness combinations.

1 Introduction

Incorporation of a secondary phase in the form of fibers or particles for the purpose of improving structural properties of polymers has been in practice for almost a century now. A daily life engineering application exemplifying this fact is the use of carbon black toughened rubbers in automobile tires. Glass-fiber-reinforced plastics, probably one of the most widely used polymer composites, present another classical instance. Furthermore, polymer composites, particulate composites in particular, find their usage in a variety of application fields ranging from aircraft structural components to biocompatible implants. See [31] and the references therein for a survey on applications of particulate composites.

A crucial development in the domain of polymer composites has been the inclusion of nanosized filler particles instead of the commonly adopted microsized reinforcements. This has lead to the development of the very exciting and challenging field of nano-composites. One of the principal factors for the recent prevalence of nanosized inclusions over the conventional microsized ones, as demonstrated by experimental studies [5, 6, 9, 11, 43], is their better structural reinforcement capability. Nanosized filler particles when well dispersed within the polymer matrix, are much more effective at improving the mechanical performance of the base polymer. In [11], specimens based on alumina-particle-reinforced vinyl ester resin, over an exhaustive range of particles diameters, ranging from $70\,\upmu \hbox {m}$ down to 15 nm, where subjected to uni-axial tensile testing. The effective elastic modulus of the composite material was found to increase significantly as the size of the filler particles decreased from the $\,\upmu \hbox {m}$ range to the nm range. Similar observations w.r.t the macroscopic elastic and viscoleastic properties were made in [5, 6], where nano-composites comprising of silica particles, with particle sizes ranging from 500 nm down to 15 nm, embedded in a poly methyl methacrylate matrix were tested. Some of the other morphological factors influencing the overall properties of a composite include the degree of interfacial adhesion between the particles and the matrix material, and the shape of the filler particles [10]. In this work, however, we focus our attention on modeling the influence of particle size in case of perfect spherical filler particles, while assuming perfect adhesion between the filler and the matrix phases.

An important aspect governing the properties of polymer composites is the formation of an interphase layer around the filler particles, as has been experimentally observed in the case of silica-elastomer nano-composites [4]. This phenomenon is much more dominant in nanosized filler particles owing to their extremely high specific surface area (SSA) [42]. The high SSA afforded by nano-particles provides many more opportunities for interaction with the matrix material for a prescribed volume fraction of filler material, thereby resulting in a much higher interphase volume fraction and consequently the aforementioned improvement. This phenomenon of greater strengthening through smaller (nanosized) particles as compared to larger (microsized) inclusions is generally referred to as the smaller is stronger size effect in literature, and will be abbreviated as size effect for the remainder of this paper.

In view of the size effect’s practical significance as described above, it needs to be considered during the computational modeling of nano-composites. Standard two-phase continuum mechanics based models fail to capture this size effect as they lack an intrinsic length scale which is essential to model this effect. Thus, appropriate enhancement of the standard model is required. Several attempts at modeling the size effect through analytical as well as computational approaches have been made in the past. These approaches could be broadly categorized into three classes as described below.

• I. These methodologies are based on explicit modeling of the variation in material properties within the interphase region. A naive approach in this regard is to consider the interphase region as a third phase with homogeneous material properties, and to vary the material properties over a range of admissible values in order to match the experimental results [8]. Another idea is to consider the interphase region to be formed of several layers, each of which has different material properties, as is done in developing the analytical model presented in [26]. The next level of sophistication in the modeling strategy would lead to an interphase region with continuously varying material properties. Such an approach was adopted in [27], where a small-strain model for computing the effective bulk modulus of a particulate composite with interphase was proposed.

• II. This class deals with approximation of the interphase region by means of a lower dimensional interface, primarily based on the theory of surface and interface elasticity by Gurtin and Murdoch [20, 29]. Herein, the interface, alike the bulk, is equipped with its own energetic structure and consequently, own tensorial stress and strain measures. Furthermore, the displacement across the interface is continuous while a jump in the stresses is allowed, culminating in a so-called elastic or coherent interface. An analytical study considering this effect during the computation of the overall behavior of solids containing nano-inhomogeneities is presented in [13]. In [49], a level set approach was employed to model the interface effects resulting in the size-dependent effective behavior of nano-composites. All of these studies were restricted to the small-strain setting. Applications involving interface/surface effects in the finite-deformation setting can be found in [23, 24].

• III. The last category comprises of techniques based on the so-called theory of generalized continua [14] as these possess an inherent length scale. See [15] and the references therein for details on the applications of generalized continuum models for modeling the size effect in heterogeneous materials. A major drawback to these techniques is the modeling effort resulting from the use of $C^1$ continuous finite elements and additional degrees of freedom since higher-order stress/strain measures are usually involved.

Asides, attempts have also been made to employ molecular dynamics (MD) simulations for understanding the variation in the elastic modulus within the interphase region [7, 33].

The myriad of possible applications of polymer nano-composites or heterogeneous materials in general has helped modeling their effective behavior garner a great deal of interest from the research community. Since the complete resolution of the heterogeneities for a specimen with complex micro-structure would be computationally prohibitive, the technique of computational homogenization [16, 17] is commonly employed in practice. Herein, a microscopic representative volume element (RVE) is attached to each point on the micro-scale and the macroscopic constitutive response is determined by means of averaging the response over the RVE, which in turn results from the solution of the micro-scale boundary value problem (BVP) (cf. Fig. 1). For applications of this technique to the modeling of nano-composites involving size effects, see [23, 40, 49].

The standard first-order computational homogenization approach, which is employed in this work, is not as it is, suitable for modeling the effective behavior of polymer nano-composites as it lacks a length scale necessary to model the size effect mentioned above. A straight-forward looking solution could be the second-order computational homogenization [16] scheme possessing the required length scale parameter. However, the involvement of higher order strains and stresses poses complications in terms of finite element implementation. Consequently, we consider enhancements to the first-order homogenization scheme for modeling size effects in polymer nano-composites, i.e. approaches belonging to classes I and II from the classification presented in Sect. 1. A systematic comparison of the interface energetics [24] approach and a novel continuously graded interphase approach is presented. Scenarios depicting varying stiffness of the matrix and the filler materials are presented. A finite deformation setting with a hyperelastic material model is considered.

This manuscript is organized as follows. We begin with a short review of computational homogenization in Sect. 2, followed by a concise summary of the fundamentals of the Interface-Energetics-Enhanced Computational Homogenization (IECH) approach (Sect. 3) including all the formulations essential for its finite element implementation. Next, we introduce in Sect. 4, the Gradient-Interphase-Enhanced Computational Homogenization (GICH) approach which is based on the concept of a continuously graded interphase region. A comprehensive evaluation of the two approaches, especially their applicability in modeling size effects in different composite materials, is presented by means of representative numerical examples in Sect. 5. Lastly, Sect. 6 summarizes the principal findings of this study and puts forward a few proposals for the future work.

Fig. 1

Schematic of the computational homogenization process depicting the BVPs at the two scales involved

2 Computational homogenization (CH)

A full-blown direct numerical simulation, at the continuum level, of an engineering component having dimensions typically of the order of centimeters and composed of a heterogeneous material such as a polymer nano-composite would ultimately turn out to be computationally prohibitive, even with the modern computing infrastructure. This necessitates some sort of averaging or homogenization, as it is commonly referred to in the literature, of the material properties. A simple analytical approach in this regard is the rule of mixtures, which essentially involves a weighted average of the constituent material properties in order to determine the composite material’s behavior. Some of the relatively more accurate analytical techniques include the Mori-Tanaka method and the Hashin-Strikman bounds, amongst others. The classical texts [34, 41] provide a generic, mathematically oriented introduction to the topic of analytical homogenization. See [30] for a detailed treatment on topic of micro-mechanics, the sub-field focusing on homogenization of mechanical properties for materials with heterogeneous micro-structure.

The numerics-based approach, i.e. computational homogenization [45, 46], on the other hand, allows for consideration of micro-structure details such as inclusion shape and dispersion of the filler particles. In this technique, the constitutive response at each point on the macro-scale is determined by subjecting a cut-out or the so-called RVE [3] of the material to a predefined macroscopic deformation gradient. The BVP with appropriate boundary conditions is solved at the micro-scale (cf. Fig. 1). The recent book [50] provides a detailed treatment on the topic of computational homogenization of mechanical properties. See also, [32, 39] for detailed, recent surveys on various homogenization techniques. The following sub-sections attempt to outline the underlying concepts of the standard first-order computational homogenization technique.

2.1 Kinematics and balance equations

Let $\varvec{\varphi }: \mathbf {X}\rightarrow \mathbf {x}$ be the deformation map describing the position in the deformed or spatial configuration of a point having position $\mathbf {X}$ in the undeformed or material configuration. The microscopic deformation gradient $\displaystyle \mathbf {F}(\mathbf {X}) :=\mathrm {Grad}\, \varvec{\varphi }$ is then given as the sum of the macroscopic deformation gradient $\mathbf {F}^{\mathrm {M}}$ and its fluctuation $\tilde{\mathbf {F}}(\mathbf {X}) :=\mathrm {Grad}\, \tilde{\varvec{\varphi }}$, i.e.

$$\begin{aligned} \mathbf {F}(\mathbf {X})&= \mathbf {F}^{\mathrm {M}} + \tilde{\mathbf {F}}(\mathbf {X}), \end{aligned}$$

(1)

$$\begin{aligned} \implies \varvec{\varphi }(\mathbf {X})&= \mathbf {F}^{\mathrm {M}} \cdot \mathbf {X}+ \tilde{\varvec{\varphi }}(\mathbf {X}), \end{aligned}$$

(2)

where $\mathbf {x}= \varvec{\varphi }(\mathbf {X})$ and $\tilde{\mathbf {x}} = \tilde{\varvec{\varphi }}(\mathbf {X})$ represent the position and its fluctuations in the spatial configuration while the gradient operator $\mathrm {Grad}\, \left\{ \bullet \right\}$ in the material configuration is defined as $\displaystyle \mathrm {Grad}\, \left\{ \bullet \right\} :=\frac{\partial \, \left\{ \bullet \right\} }{\partial \mathbf {X}}$. Further, the micro-scale BVP (without boundary conditions) reads

$$\begin{aligned} \mathrm {Div}\, \mathbf {P} = \mathbf {0}, \end{aligned}$$

(3)

with $\mathrm {Div}\, \mathbf {P}$ denoting the divergence of the microscopic Piola stress $\mathbf {P}$ in the material configuration, while the inertia and body forces acting at the micro-scale are neglected.

2.2 Constitutive model

The constitutive relation between the stresses and the strains is another important ingredient essential for solving the BVP of (3). In this work, we adopt a hyperelastic constitutive model of Neo-Hooke type for both the filler as well as the matrix phases. The bulk strain energy functional, decomposed into volumetric and isochoric parts, then reads

$$\begin{aligned} \Psi (\mathbf {C})&= \Psi _{\mathrm {vol}} + \Psi _{\mathrm {iso}} \end{aligned}$$

with

$$\begin{aligned} \Psi _{\mathrm {vol}}&= \frac{1}{4} \, \kappa \left[ J^2 - 1 - 2 \, \ln J \right] \end{aligned}$$

(4)

and

$$\begin{aligned} \Psi _{\mathrm {iso}}&= \frac{1}{2} \, \mu \left[ \displaystyle J^{ -\frac{2}{3}} \, \mathbf {C}: \mathbf {I}- 3 \right] , \end{aligned}$$

where $\kappa$ and $\mu$ denote the bulk and shear moduli, respectively. Further, $\mathbf {C}:={\mathbf {F}}^{t} \cdot \mathbf {F}$ and $\mathbf {I}$ are the right Cauchy-Green strain tensor and the second order unit tensor in the material configuration, respectively, while $J :=\det {\mathbf {F}}$ represents the determinant of the deformation gradient $\mathbf {F}$. Next, following the Coleman-Noll exploitation procedure (see for e.g. [21] for details), the second Piola-Kirchoff Stress is given as

$$\begin{aligned} \mathbf {S}= 2 \, \frac{\partial \Psi }{\partial \mathbf {C}} = \underbrace{2 \, \frac{\partial \Psi _{\mathrm {vol}}}{\partial \mathbf {C}}}_{\displaystyle =:\mathbf {S}_{\mathrm {vol}}} + \underbrace{2 \, \frac{\partial \Psi _{\mathrm {iso}}}{\partial \mathbf {C}}}_{\displaystyle =:\mathbf {S}_{\mathrm {iso}}} \end{aligned}$$

(5)

with the two parts being computed as

$$\begin{aligned} \mathbf {S}_{\mathrm {vol}} =&\ \frac{1}{2} \, \kappa \left[ J^2 - 1\right] \, {\mathbf {C}}^{-1}, \end{aligned}$$

(6)

and

$$\begin{aligned} \mathbf {S}_{\mathrm {iso}} =&\ \mu \, J^{ -\frac{2}{3}} \, \left[ \mathbf {I}- \frac{1}{3} \, \left[ \mathbf {C}: \mathbf {I}\right] \, {\mathbf {C}}^{-1} \right] \end{aligned}$$

(7)

respectively. Since an iterative scheme of Newton-Raphson type is usually employed to solve the system of nonlinear equations stemming from the finite element discretization, a consistent tangent modulus is required to achieve quadratic convergence, cf. Appendix 1.

2.3 Micro–macro transition

The macroscopic Piola stress $\mathbf {P}^{\mathrm {M}}$ is determined in accordance with the Hill-Mandel condition, which essentially dictates the equivalence of energies between the micro-scale and the macro-scale. Let $\displaystyle \delta \mathbf {F}^{\mathrm {M}}$ and $\delta \mathbf {F}$ be the macro- and microscopic virtual deformation gradients. Then, according to this energy equivalence, the macroscopic virtual work must be equal to the average over the RVE volume $\mathcal {V}$ of the microscopic virtual work, i.e.

$$\begin{aligned} \mathbf {P}^{\mathrm {M}} : \delta \mathbf {F}^{\mathrm {M}} = \frac{1}{|\mathcal {V} |} \int \limits _{\displaystyle \mathcal {V}} {\mathbf {P}: \delta \mathbf {F}} \ {\mathrm {dV}}. \end{aligned}$$

(8)

Inserting (1) in  (8) leads us to the so-called stress averaging theorem, which states that the macroscopic stress equals the average of its microscopic counter part, i.e.

$$\begin{aligned} \mathbf {P}^{\mathrm {M}} = \frac{1}{|\mathcal {V} |} \int \limits _{\displaystyle \mathcal {V}} {\mathbf {P}} \ {\mathrm {dV}} \end{aligned}$$

(9)

and the relation

$$\begin{aligned} \frac{1}{|\mathcal {V} |} \int \limits _{\displaystyle \mathcal {V}} {\mathbf {P}: \delta \tilde{\mathbf {F}}} \ {\mathrm {dV}} = 0, \end{aligned}$$

which on doing some basic calculus yields

$$\begin{aligned} ~&\frac{1}{|\mathcal {V} |} \int \limits _{\displaystyle \partial \mathcal {V}} {\delta \tilde{\mathbf {x}} \cdot \mathbf {t}} \ {\mathrm {dA}} = 0 \nonumber \\ \displaystyle \implies&\underbrace{\delta \tilde{\mathbf {x}} = 0 }_{\text {LBC}} \quad \text {or} \qquad \underbrace{\delta \tilde{\mathbf {x}}^{+} = \delta \tilde{\mathbf {x}}^{-} \ \text {and} \ \mathbf {t}^{+} = -\mathbf {t}^{-} }_{\text {PBC}}. \end{aligned}$$

(10)

In (10), the quantities $\delta \tilde{\mathbf {x}}$ and $\mathbf {t}$ represent the fluctuations in the virtual spatial coordinate and the tractions acting on the boundary of the RVE, respectively. The positive and negative superscripts denote the quantities living on the two opposite sides of the boundary. Thus, the Hill-Mandel condition also leads to different possible boundary conditions, namely linear boundary conditions (LBC) and periodic-antiperiodic boundary conditions (PBC), which could be applied to the micro-scopic RVE, cf (10). In this work, we employ PBCs as these result in a more realistic material response, cf. [28] for instance.

2.4 Variational formulation of the micro-scale BVP

The weak form or the so-called variational formulation is the start-point of the numerical solution process and is formulated in the following. Bringing together the balance equation from (3) and the PBC from (10), the strong form of the micro-scale BVP reads

$$\begin{aligned} \begin{aligned} \mathrm {Div}\, \mathbf {P}&= \mathbf {0}\quad \text {in} \quad \mathcal {B}_{\mathrm {0}}\\ \mathbf {x}^{+} - \mathbf {x}^{-}&= \mathbf {F}^{\mathrm {M}} \cdot \left[ \mathbf {X}^{+} - \mathbf {X}^{-}\right] \quad \text {on} \quad \partial \mathcal {B}_{\mathrm {0}} , \end{aligned} \end{aligned}$$

(11)

where Eq. (11)₂ represents an alternative formulation of the periodicity of fluctuations in (10)₂.

Next, the variational formulation is obtained by testing the strong form of the BVP, i.e. (11) with appropriate test function $\displaystyle \delta \varvec{\varphi }\in \varvec{\mathcal {H}}^{1}_{{per}}\left( \mathcal {B}_{\mathrm {0}}\right)$, and integrating the result over the corresponding domains. A straight forward application of the divergence theorem leads to the following variational formulation: Find $\displaystyle \varvec{\varphi }\in \varvec{\mathcal {H}}^{1}_{afper}\left( \mathcal {B}_{\mathrm {0}}\right)$ s.t.

$$\begin{aligned} \begin{aligned}&\int \limits _{\displaystyle \mathcal {B}_{\mathrm {0}}} {\mathbf {P}: \mathrm {Grad}\, \delta \varvec{\varphi }} \ {\mathrm {dV}} = 0 \qquad \forall \ \delta \varvec{\varphi }\in \varvec{\mathcal {H}}^{1}_{per}\left( \mathcal {B}_{\mathrm {0}}\right) , \end{aligned} \end{aligned}$$

(12)

where $\varvec{\mathcal {H}}^{1}_{{per}}$ represents the dim-dimensional (${dim=2,3}$) Sobolev space of vector valued functions possessing first weak derivative and the subscript per here denotes that the test functions are periodic w.r.t the boundary of the RVE, see [12, Chapter 3] for a detailed explanation. Similarly, the subscript afper in $\varvec{\mathcal {H}}^{1}_{{afper}}$ denotes that the trial functions are affine periodic w.r.t the boundary of the RVE, i.e. they satisfy relations of the form given by (11)₂.

3 Interface-energetics-enhanced computational homogenization (IECH)

The interface energetics [24, 25] approach, similar to other candidates of class-II (cf. Sect. 1), involves approximating the interphase region by a zero-thickness interface between the matrix and the filler materials. The underlying principle here is that of lower dimensional energetics, wherein, lower dimensional geometric entities, such as surfaces, interfaces between different materials, and curves formed by the intersection of material interfaces with surfaces are equipped with their own free energy densities, stress-strain measures and constitutive relations, and consequently contribute to the overall thermomechanics of the body. Focusing on energetic interfaces, these could be classified into different categories, depending on whether the jump in the deformation map $\varvec{\varphi }$ and/or in the traction, i.e. ${\llbracket \mathbf {P} \rrbracket \cdot \overline{\mathbf {N}}}$, acting on the interface is non-zero [25]. Here, the operator $\displaystyle \llbracket \left\{ \bullet \right\} \rrbracket$ represents the jump in the quantity $\left\{ \bullet \right\}$ across the interface, i.e. $\displaystyle \llbracket \left\{ \bullet \right\} \rrbracket :=\left\{ \bullet \right\} ^{+} - \left\{ \bullet \right\} ^{-}$, where the positive and negative superscripts denote the values on the two sides of the interface. Further, the normal to the interface is denoted by $\overline{\mathbf {N}}$. A perfect interface restricts the jumps in both these quantities and represents the situation evident in a standard two-phase continuum model. In this work, we employ the elastic interface model which is based on the theory of interface elasticity [29]. This model facilitates the modeling of the additional strengthening resulting from the smaller size of the filler particles by means of allowing a jump in the traction, i.e. ${\llbracket \varvec{\varphi } \rrbracket = \mathbf {0}, \llbracket \mathbf {P} \rrbracket \cdot \overline{\mathbf {N}}\ne \mathbf {0}}$ . Yet another possibility is the cohesive interface model which allows for opening of the interface between the two materials, since a jump in the displacement is allowed but not in the traction. This type of interface is primarily employed in modeling crack propagation such as in the cohesive zone model. A generalized interface model which allows for a jump in both the quantities, i.e. ${\llbracket \varvec{\varphi } \rrbracket \ne \mathbf {0}, \llbracket \mathbf {P} \rrbracket \cdot \overline{\mathbf {N}}\ne \mathbf {0}}$, is discussed in detail in [40].

Next, we summarize the overall working principle of the interface energetics approach to the size-dependent homogenization of heterogeneous materials. A detailed exposition of this topic can be found in [24, 25] and the references therein. In the following, interface quantities are distinguished from their bulk counterparts by means of an over-bar, i.e. for a bulk quantity $\left\{ \bullet \right\}$, the corresponding quantity for the interface shall be denoted as $\left\{ \overline{\bullet } \right\}$.

Fig. 2

: RVEs comprised of matrix material (shown in blue) and a circular inclusion (shown in red) and equipped with (a) an energetic interface or (b) a graded interphase of thickness t (shown in green). Here, the subscripts m, f and i refer to the matrix, filler and the interphase regions, respectively. (Color figure online)

3.1 Interface energetics: kinematics and balance equations

Consider a body $\mathcal {B}_{\mathrm {0}} \subset \mathbb {R}^{\dim }$, where$\dim =2,3$, composed of two different materials separated by an interface $\mathcal {I}_{\mathrm {0}}$ as depicted in Fig. 2a. Let $\overline{\mathbf {N}}= \overline{\mathbf {N}}\left( \overline{\mathbf {X}}\right)$ be the normal in the material configuration to the interface at a point P with coordinates $\overline{\mathbf {X}}$. The interface unit tensor in material configuration, denoted by $\overline{\mathbf {I}}$, is then defined as

$$\begin{aligned} \overline{\mathbf {I}}:=\mathbf {I}- \overline{\mathbf {N}}\otimes \overline{\mathbf {N}}, \end{aligned}$$

(13)

where $\mathbf {I}$ is the second order unit tensor introduced before. Furthermore, the interface gradient and divergence operators in the material configuration are defined as

$$\begin{aligned} \overline{\mathrm {Grad}}\, \left\{ \overline{\bullet } \right\}&:=\mathrm {Grad}\, \left\{ \overline{\bullet } \right\} \cdot \overline{\mathbf {I}}, \end{aligned}$$

(14)

and

$$\begin{aligned} \overline{\mathrm {Div}}\, \left\{ \overline{\bullet } \right\}&:=\overline{\mathrm {Grad}}\, \left\{ \overline{\bullet } \right\} : \overline{\mathbf {I}} \end{aligned}$$

(15)

respectively. It is noteworthy that $\overline{\mathbf {I}}$ is actually a projection tensor which projects the relevant physical quantities onto the interface. The projection process involves elimination of components normal to the interface. Consequently, for the elastic interface considered here, the interface deformation gradient $\overline{\mathbf {F}}$ is computed as

$$\begin{aligned} \nonumber \overline{\mathbf {F}}:=&\ \overline{\mathrm {Grad}}\, \overline{\varvec{\varphi }} \\ \nonumber =&\ \overline{\mathrm {Grad}}\, \varvec{\varphi } \qquad \left[ \overline{\varvec{\varphi }}= \varvec{\varphi }\big |_{\mathcal {I}_{\mathrm {0}}}\text { for elastic interface}\right] \\ \nonumber =&\ \mathrm {Grad}\, \varvec{\varphi } \cdot \overline{\mathbf {I}}\qquad [\text {using}~14] \\ \implies \overline{\mathbf {F}}=&\ \mathbf {F}\cdot \overline{\mathbf {I}}, \end{aligned}$$

(16)

where $\overline{\varvec{\varphi }}: \overline{\mathbf {X}}\rightarrow \overline{\mathbf {x}}$ is the interface deformation map describing the position in the spatial configuration of a point on the interface with position $\overline{\mathbf {X}}$ in the material configuration. The determinant $\overline{J}$ of the interface deformation gradient is defined as [44]

$$\begin{aligned} \overline{J}= \overline{\mathrm {Det}}\, \overline{\mathbf {F}} :=|\mathrm {COF}\, \mathbf {F} \cdot \mathbf {N}|, \end{aligned}$$

(17)

where ${\mathrm {COF}\, \mathbf {A} :=\mathrm {Det}\, \mathbf {A} \, {\mathbf {A}}^{-t}}$ denotes the cofactor of a second-order invertible tensor $\mathbf {A}$, which in the current scenario reads ${\mathrm {COF}\, \mathbf {F} = J \, {\mathbf {F}}^{-t}}$, since $J :=\mathrm {Det}\, \mathbf {F}$.

In a similar fashion, the interface counterpart of the right Cauchy-Green strain tensor $\mathbf {C}$, denoted by $\overline{\mathbf {C}}$, is computed as

$$\begin{aligned} \nonumber \overline{\mathbf {C}}:=&\ {\overline{\mathbf {F}}}^{t} \cdot \overline{\mathbf {F}}= \ {\left[ \mathbf {F}\cdot \overline{\mathbf {I}}\right] }^{t} \cdot \mathbf {F}\cdot \overline{\mathbf {I}}\qquad \left[ \text {using}~16\right] , \\[1ex] \nonumber =&\ {\overline{\mathbf {I}}}^{t} \cdot \left[ {\mathbf {F}}^{t} \cdot \mathbf {F}\right] \cdot \overline{\mathbf {I}}\\[1ex] \implies \overline{\mathbf {C}}=&\ \overline{\mathbf {I}}\cdot \mathbf {C}\cdot \overline{\mathbf {I}}, \end{aligned}$$

(18)

since ${\overline{\mathbf {I}}}^{t} = \overline{\mathbf {I}}$ by definition.

Similar to the balance of the linear momentum for the bulk (cf. (3)), it is also possible to derive the balance of linear momentum for the interface (see [24, 25] for details), which in the absence of distributed interface forces reads

$$\begin{aligned} \overline{\mathrm {Div}}\, \overline{\mathbf {P}} + \llbracket \mathbf {P} \rrbracket \cdot \overline{\mathbf {N}}= \mathbf {0}\quad \text {on} \quad \mathcal {I}_{\mathrm {0}}, \end{aligned}$$

(19)

where $\overline{\mathbf {P}}$ is the interface Piola stress tensor. Furthermore, the balance of angular momentum dictates the symmetry of the interface counter part of the second Piola-Kirchhoff stress tensor, i.e.

$$\begin{aligned} \begin{aligned}&{\overline{\mathbf {S}}}^{t} = \overline{\mathbf {S}}, \\[1ex]&\text {with} \quad \overline{\mathbf {S}}:=\overline{\mathbf {f}}\cdot \overline{\mathbf {P}}, \end{aligned} \end{aligned}$$

(20)

where $\overline{\mathbf {f}}$, as defined in (36), denotes the inverse deformation gradient for the interface.

3.2 Constitutive model

We consider for the interface a hyperelastic constitutive model, similar to that employed for the bulk in Sect. 2.2, with the interface free energy functional given as

$$\begin{aligned} \overline{\Psi }\left( \overline{\mathbf {C}}\right) =&\ \overline{\Psi }_1 + \overline{\Psi }_2 \end{aligned}$$

(21)

with

$$\begin{aligned} \overline{\Psi }_1 :=&\ \frac{1}{4} \, \overline{\kappa }\left[ \overline{J}^2 - 1 - 2 \, \ln \overline{J}\right] \end{aligned}$$

and

$$\begin{aligned} \overline{\Psi }_2 :=&\ \frac{1}{2} \, \overline{\mu }\left[ \overline{J}^{{-2} / {\overline{\dim }} } \, \overline{\mathbf {C}}: \overline{\mathbf {I}}- \overline{\dim }\right] \end{aligned}$$

respectively, with $\overline{\dim }$ representing the spatial dimension of the interface, i.e. $\overline{\dim }= 2$ if the bulk as modelled as 3D, and $\overline{\dim }=1$ for a 2D bulk (plane strain case). Further, $\overline{\kappa }$ and $\overline{\mu }$ denote the interface bulk and shear moduli respectively. It should be noted that the split of the strain energy density here is merely for the purpose of easing the process of stress and tangent computations and does not involve a volumetric-isochoric splitting as was the case for the bulk strain energy functional. Next, following the Coleman-Noll procedure as mentioned for the bulk, the interface second Piola-Kirchhoff stress tensor is given as

$$\begin{aligned} \overline{\mathbf {S}}= 2 \, \frac{\partial \overline{\Psi }}{\partial \overline{\mathbf {C}}} = \underbrace{2 \, \frac{\partial \overline{\Psi }_1}{\partial \overline{\mathbf {C}}}}_{\displaystyle =:\overline{\mathbf {S}}_1} + \underbrace{2 \, \frac{\partial \overline{\Psi }_2}{\partial \overline{\mathbf {C}}}}_{\displaystyle =:\overline{\mathbf {S}}_2}, \end{aligned}$$

(22)

with the two parts being computed as

$$\begin{aligned} \overline{\mathbf {S}}_1 =&\ \frac{1}{2} \, \overline{\kappa }\left[ \overline{J}^2 - 1\right] \, {\overline{\mathbf {C}}}^{-1}, \end{aligned}$$

(23)

and

$$\begin{aligned} \overline{\mathbf {S}}_2 =&\ \overline{\mu }\, \overline{J}^{\left[ {-2} / {\overline{\dim }}\right] } \, \left[ \overline{\mathbf {I}}- \frac{1}{\overline{\dim }} \, \left[ \overline{\mathbf {C}}: \overline{\mathbf {I}}\right] \, {\overline{\mathbf {C}}}^{-1} \right] \end{aligned}$$

(24)

respectively. For the formulation of the term ${\overline{\mathbf {C}}}^{-1}$ see Appendix 1. Furthermore, Appendix 1 summarizes the computation of the consistent interface tangent modulus as required for the finite element discretization based numerical solution process.

3.3 Micro-scale BVP and its variational formulation

For the RVE depicted in Fig. 2a, and subjected to a macroscopic deformation gradient $\mathbf {F}^{\mathrm {M}}$ under PBCs (cf. Sect. 2), the BVP to be solved, considering the contribution of the energetic interfaces reads

$$\begin{aligned} \begin{aligned} \mathrm {Div}\, \mathbf {P}&= \mathbf {0}\quad \text {in} \quad \mathcal {B}_{\mathrm {0}}\\ \overline{\mathrm {Div}}\, \overline{\mathbf {P}} + \llbracket \mathbf {P} \rrbracket \cdot \overline{\mathbf {N}}&= \mathbf {0}\quad \text {on} \quad \mathcal {I}_{\mathrm {0}},\\ \mathbf {x}^{+} - \mathbf {x}^{-}&= \mathbf {F}^{\mathrm {M}} \cdot \left[ \mathbf {X}^{+} - \mathbf {X}^{-}\right] \quad \text {on} \quad \partial \mathcal {B}_{\mathrm {0}} . \end{aligned} \end{aligned}$$

(25)

Note that the BVP (25) is essentially an extension of the standard case BVP of (11) by means of the balance Eq. (19) for the energetic interface.

Next, the variational formulation is obtained by testing the strong form of the BVP, i.e. (25)_1-2 with appropriate test functions $\displaystyle \delta \varvec{\varphi }\in \varvec{\mathcal {H}}^{1}_{{per}}\left( \mathcal {B}_{\mathrm {0}}\right)$ and $\displaystyle \delta \overline{\varvec{\varphi }}\in \varvec{\mathcal {H}}^{1}_{{per}}\left( \mathcal {I}_{\mathrm {0}}\right)$, respectively, and integrating the result over the corresponding domains. Some further manipulations involving divergence theorem, cf. [24], yield the weak form: Find $\displaystyle \varvec{\varphi }\in \varvec{\mathcal {H}}^{1}_{afper}\left( \mathcal {B}_{\mathrm {0}}\right)$ and $\displaystyle \overline{\varvec{\varphi }}\in \varvec{\mathcal {H}}^{1}_{afper}\left( \mathcal {I}_{\mathrm {0}}\right)$ s.t.

$$\begin{aligned} \begin{aligned}&\int \limits _{\displaystyle \mathcal {B}_{\mathrm {0}}} {\mathbf {P}: \mathrm {Grad}\, \delta \varvec{\varphi }} \ {\mathrm {dV}} + \int \limits _{\displaystyle \mathcal {I}_{\mathrm {0}}} {\overline{\mathbf {P}}: \overline{\mathrm {Grad}}\, \delta \overline{\varvec{\varphi }}} \ {\mathrm {dA}} = 0 \\ \qquad&\forall \ \delta \varvec{\varphi }\in \varvec{\mathcal {H}}^{1}_{per}\left( \mathcal {B}_{\mathrm {0}}\right) , \ \delta \overline{\varvec{\varphi }}\in \varvec{\mathcal {H}}^{1}_{per}\left( \mathcal {I}_{\mathrm {0}}\right) , \end{aligned} \end{aligned}$$

(26)

where the function spaces $\varvec{\mathcal {H}}^{1}_{{per}}$ and $\varvec{\mathcal {H}}^{1}_{{afper}}$ are as described in Sect. 2.4. Here, it has been assumed that the contributions that emerge from the energetic curves formed by energetic interfaces possibly intersecting with the boundary $\partial \mathcal {B}_{\mathrm {0}}$ cancel out. This assumption indeed always holds for the RVE homogenization problems dealt with in this work, since the filler particles and hence the energetic interfaces do not intersect the boundary of the RVE. Furthermore, even if a filler particle intersects the RVE, the contributions from the energetic curves on opposite sides of the RVE boundary will cancel out due to the fact that the RVE boundary is periodic as dictated by the application of PBCs.

3.4 Modified stress averaging

The presence of energetic interfaces influences the overall mechanical equilibrium of the RVE. Consequently, the inclusion of interface energetics in the computation homogenization scheme, cf. Sect. 2, culminates in an extended Hill-Mandel condition [23] wherein additional contributions due to the energetic interfaces are considered. The stress averaging theorem of (9) then takes the form

$$\begin{aligned} \nonumber \mathbf {P}^{\mathrm {M}}&= \frac{1}{|\mathcal {B}_{\mathrm {0}} |} \left[ \int \limits _{\displaystyle \mathcal {B}_{\mathrm {0}}} {\mathbf {P}} \ {\mathrm {dV}} + \int \limits _{\displaystyle \mathcal {I}_{\mathrm {0}}} {\overline{\mathbf {P}}} \ {\mathrm {dA}} \right] \\ \nonumber&= \frac{1}{|\mathcal {B}_{\mathrm {0}} |} \int \limits _{\displaystyle \mathcal {B}_{\mathrm {0}}} {\mathbf {P}} \ {\mathrm {dV}} + \underbrace{\frac{|\mathcal {I}_{\mathrm {0}} |}{|\mathcal {B}_{\mathrm {0}} |}}_{\displaystyle =:\frac{1}{S_l}} \, \frac{1}{|\mathcal {I}_{\mathrm {0}} |} \int \limits _{\displaystyle \mathcal {I}_{\mathrm {0}}} {\overline{\mathbf {P}}} \ {\mathrm {dA}}\\ \implies \mathbf {P}^{\mathrm {M}}= & {} lt;\mathbf {P}>_{_{\scriptstyle \mathcal {B}_{\mathrm {0}}}} + \frac{1}{S_l} \, <\overline{\mathbf {P}}>_{_{\scriptstyle \mathcal {I}_{\mathrm {0}}}}, \end{aligned}$$

(27)

where $|\mathcal {B}_{\mathrm {0}} |$ and $|\mathcal {I}_{\mathrm {0}} |$ denote the volume of the RVE and the surface area of the interface region, respectively, while the operator $<\left\{ \bullet \right\} >_{_{\scriptstyle }}$ represents the average of the quantity over the specified domain. Thus, the IECH approach includes the length scale $S_l$ with the unit of length, which is employed for capturing the desired size effect. This will be explained in detail in Sect. 5.1 by means of a numerical example.

Fig. 3

Variation in interphase properties with different values of the exponent n for the case of silica filler particles in epoxy matrix. Since the filler material, i.e. silica, is already much stiffer than the matrix, the value of the scaling factor $\alpha$ in (28) has been set to 1.0 here

4 Gradient-interphase-enhanced computational homogenization (GICH)

In contrast to the interface-energetics-based scheme described in the previous section, interphase-based approaches facilitate explicit modeling of the variation in material properties within the interphase region around the filler particles. In line with the experimental evidence presented in [4, 42], such a graded interphase would facilitate more accurate modeling of the near particle, and consequently overall behavior of polymer nano-composites. This, however, comes at the price of increased computational costs since the finite element mesh has to be sufficiently fine in order to appropriately resolve the interphase region geometry. Unlike the step-wise variation typically adopted, cf. [35], the Gradient-Interphase-Enhanced Computational Homogenization (GICH) approach introduced in this work allows for a continuous variation in the material properties within the interphase region, in the finite deformation setting. In the finite element sense, this means that the material properties at each quadrature point within the interphase region vary according to the point’s location relative to the filler particle.

A polynomial interpolation function, cf. [36, 48], is employed to determine the material properties within the interphase region. Consequently, the Young’s modulus $E_i\left( r\right)$ and the shear modulus $\mu _i\left( r\right)$ at a point within the interphase region, located at a distance $r_f \le r \le r_m$ from the filler particle, cf. Fig. 2b, are given by

$$\begin{aligned} \displaystyle E(r) = E_m + \left[ \frac{r_m - r}{r_m - r_f}\right] ^n \left[ \alpha E_f - E_m\right] \end{aligned}$$

(28)

and

$$\begin{aligned} \displaystyle \mu (r) = \mu _m + \left[ \frac{r_m - r}{r_m - r_f}\right] ^n \left[ \alpha \mu _f - \mu _m\right] , \end{aligned}$$

(29)

respectively. Here, the parameter $\alpha > 0$, also referred to as the interphase stiffness scaling factor later on, allows for stiffening of the interphase near the filler particle. This is required for the case of weaker inclusions in a stiffer matrix, as will be shown in Sect. 5.2. The equations [37, Table 4.1]

$$\begin{aligned} \displaystyle \nu _i(r) = \frac{E_i}{2 \, \mu _i} - 1 \end{aligned}$$

(30)

and

$$\begin{aligned} \displaystyle \kappa _i(r) = \frac{2 \, \mu _i\left[ 1 + \nu _i\right] }{3 \left[ 1 - 2\, \nu _i\right] } \end{aligned}$$

(31)

can then be employed to determine the corresponding value of the Poisson’s ratio $\nu _i\left( r\right)$ and the bulk modulus $\kappa _i \left( r\right)$ at a particular point. See Fig. 3 for an example of how the material properties vary within the interphase region. The influence of the interphase properties on the macroscopic material behavior is discussed in detail in Sect. 5.2 by means of concrete numerical examples.

Unlike the interface-energetics-based approach discussed in the previous section, the current technique involves no additional energetic response. Hence the BVP and consequently the weak form for the standard first-order computational homogenization scheme (cf. Sect. 2.4) are valid for the GICH approach.

Fig. 4

IECH size effect scaling: variation in the macroscopic Piola stress with filler particle size for the case rubber particles in epoxy matrix. Three different loading scenarios as described in Table 2 have been presented

5 Results and discussion

All numerical examples presented in the current study have been computed with an in-house simulation code written in C++ and based on the open source finite element library deal.II [1, 47]. Excessive, localized mesh refinement necessary to resolve the interphase region results in poor conditioning of the tangent stiffness matrix arising at each iteration of the Newton-Raphson scheme being used to solve the nonlinear problem. Hence, a direct linear solver shipped along with the deal.II library has been employed to solve the system of linear equations (SLE) arising during the iterative solution procedure. The unknown deformation map $\varvec{\varphi }$ has been approximated by means of bi-linear shape functions and a state of plane strain has been modeled for the 2D scenario. Further, for the purpose of finite element mesh generation, another in-house code RVEGen has been developed. Written in python and employing the open source framework GMSH [18] as a back end for mesh generation, this code currently allows for generation for finite element meshes for 2D RVEs consisting of randomly placed circular inclusions of arbitrary sizes. Furthermore, it also allows the user to control localized mesh refinement, especially in the interphase region, so as to yield the minimum number of degrees of freedom while ensuring sufficient resolution of the desired geometrical features. Special focus has been devoted to support for quadrilateral elements since our simulation code only supports quadrilateral (2D) / hexahederal (3D) elements.

In this section, we present a detailed comparison of the two approaches described in the previous sections by means of systematic numerical experimentation. In order to perform a holistic evaluation of the techniques at hand, two limit cases of material stiffness ratio between the matrix and the filler materials have been considered, i.e. the case of weak filler particles (rubber) in a stiffer matrix (epoxy) and the case of stiff filler particles (silica) in a weaker matrix (epoxy). Table 1 summarizes the requisite material properties. The symbols $E, \nu \ \text {and} \ \mu$ therein represent Young’s modulus, Poisson’s ratio and shear modulus, respectively. The values of the Young’s modulus and the Poisson’s ratio have been taken from the cited literature while those of shear modulus have been computed using Eq. (30). Furthermore, the two approaches, IECH and GICH have been examined for three different loading scenarios, namely uni-axial tension (T), shear (S) and volumetric loading (V), in order to validate their applicability to general loading conditions, cf. Table 2 for the corresponding macroscopic deformation gradients. The macroscopic Piola stress tensor $\mathbf {P}^{\mathrm {M}}$ or rather its relevant component, computed using the relation (27) for IECH and relation (9) for GICH, is considered as a measure of improvement in the overall elastic behavior and is thus compared against filler particle size. The term relevant here implies that the normal stress component $P^{\mathrm {M}}_{11}$ is considered for loading scenarios T and V, while the shear component $P^{\mathrm {M}}_{12}$ is studied for the loading scenario S. A smaller is stiffer size effect would then result in higher values of the relevant macroscopic stress component for smaller particle sizes.

Table 1 Material properties considered for numerical experimentation

Table 2 Values of the macroscopic deformation gradient $\mathbf {F}^{\mathrm {M}}$ for the different load cases considered in this study

Fig. 5

IECH size effect scaling: variation in the macroscopic Piola stress with filler particle size for the case silica particles in epoxy matrix. Three different loading scenarios as described in Table 2 have been presented

5.1 Size effect scaling for IECH

We consider the RVE of Fig. 2a having edge size $a_0$, containing a circular inclusion of radius $r_{f0}$, and thus having an initial length scale $\displaystyle {S_l^0 = {a_0^2}/{2\pi r_{f0}}}$, cf. (27). Upon scaling the RVE by a factor 1/K with $K>1$, the edge length and the filler particle radius become ${a = a_0/K}$ and ${r_f = r_{f0}/K}$ respectively. The length scale accordingly changes to ${S_l = {a^2}/{2\pi r_{f}}}$, thus yielding

$$\begin{aligned} \displaystyle \frac{S_l^0}{S_l}&= \frac{a_0^2}{2\pi r_{f0}} \, \frac{2\pi r_{f}}{a^2} = \frac{a_0^2}{2\pi r_{f0}} \, \frac{2\pi K r_{f0}}{K^2 \, a^2_0} = \frac{1}{K},\\ \implies S_l&= K \, S_l^0, \end{aligned}$$

while leaving the volume fraction of the filler material unaltered. Thus, with increasing values of the scaling factor K, we obtain a sequence of RVEs for which the length scale and hence the effectiveness of the energetic interface increases as the RVE and consequently the filler particle becomes smaller. This effectiveness is further influenced by the interface material parameters, cf. the interface strain energy functional of (21). It is noteworthy that for the 2D problem at hand, the interface is a curve, i.e. a one-dimensional entity, and thus only one material parameter is needed in order to characterize the deformation of the interface. Consequently, we set in (21) $\overline{\kappa }=0$ and vary $\overline{\mu }$ over a suitable range in order to study its impact on the overall deformation behavior of the RVE.

Figures 4 and 5 depict the size effect scaling, i.e. the variation in the relevant component of the macroscopic Piola stress tensor $\mathbf {P}^{\mathrm {M}}$ with respect to the filler particle radius $r_f$ for the two different stiffness extremes considered in this study. The interface shear modulus $\overline{\mu }$ is varied across a wide range of magnitudes starting from $|\overline{\mu } |=0$ until the point where ${|\overline{\mu } | > \max \left\{ |\mu _m |, |\mu _f |\right\} }$. Consequently, the subscripts E, R and S, in Figs. 4 and 5, correspond to epoxy, rubber and silica, respectively. This upper limit represents the scenario wherein the interface becomes almost rigid with respect to the two surrounding materials, such that a further increase in the interface stiffness would not help in improving the overall stiffness of the composite material. The scenario with $|\overline{\mu } |=0$ represents the case of no active energetic interface as is observed with the standard first-order homogenization technique. Here, no change in the macroscopic stress is exhibited upon reduction in particle size irrespective of the inclusions or the loadcase considered.

For the case of weak inclusions in a stiffer matrix, a smaller inclusion size results in a significantly stiffer macroscopic response as the strength of the interface, characterized here by $\overline{\mu }$, increases. Such behavior is observed for all the three loadcases considered, cf. Fig. 4a–c, respectively. This improvement in the overall stiffness can be understood by means of the stress contours presented in Fig. 6. The effectiveness of the energetic interface is influenced not only by its stiffness, but also the length scale $S_l$ or alternatively the particle size. Consequently, for very large particle sizes, no additional stiffening effect is observed by increasing the interface shear modulus as depicted by the left column of the stress contours. Contrarily, increasing the interface stiffness culminates in much higher stresses, especially for smaller particle sizes as can be observed from the right-most column. This behavior can be attributed to the fact that with an increase in interface’s stiffness and consequently its resistance to deformation, the neighboring matrix material needs to bear more stress and hence the microscopic stresses increases significantly. Moreover, for a highly active interface, such as that in the top row of Fig. 6, the stresses increase significantly with decreasing particle size.

Fig. 6

Contours of nodal values for the relevant Cauchy stress component $\sigma _{11}$ under uni-axial tensile loading, with varying particle sizes and interface strength for the case of rubber particles in epoxy matrix

Fig. 7

Deformed shape of the RVE subjected to shear loading for a chosen particle radius. The deformed shape has been superimposed on the undeformed one in order to visualize the relative displacement of the filler and the matrix

A similar pattern is also observed in the case of silica filler material as depicted through Fig. 5a–c. The change in magnitude of the relevant macroscopic stress component is, however quite minimal. Consider for instance the volumetric loading scenario wherein, for the case of silica particles, cf. Fig. 5c, the macroscopic stress only increases from 936 to 948 as the particle size reduces by four orders of magnitude. In stark contrast to this, for the case of rubber inclusions, the magnitude of the macroscopic stress more than doubles for the same reduction in filler particle radius, cf. Fig. 4c. This observation is corroborated by the fact that the deformation of the silica particle, which is already much stiffer than the epoxy matrix, is almost unaffected by the addition of a yet stiffer energetic interface. This can be further verified by the deformation contour under shear loading, cf. Fig. 7c, wherein the particle deformation owing to its high stiffness is almost negligible compared to that of the weaker matrix material. Accordingly, the applicability of IECH is confined to cases where the matrix material is much stiffer than the filler as is the case with rubber particles in epoxy matrix.

Another noteworthy feature of the size effect scaling curves in Figs. 4 and 5 is the saturation of the size effect w.r.t the reducing filler particle radius, especially for higher interface stiffness and for stiffer filler material in particular, cf. 5b for instance. Moreover, the exact particle size regime where the size effect scaling is observed, i.e. in the micrometric or the nanometric range, is controlled by the units of the interface shear modulus $\overline{\mu }$. Thus, a comprehensive analysis of the interface effects necessitates deciphering the units of $\overline{\mu }$. We revisit the modified stress-averaging relation of (27) with a specific focus on the units of the involved quantities. Let us assume the unit of length to be nm. Further, the units of $\mu$ and hence those of $\mathbf {P}$, $\mathbf {P}^{\mathrm {M}}$ are taken to be $\displaystyle \left[ \hbox {N}/\hbox {mm}^{2}\right]$. The units of $\overline{\mu }$ and consequently those of $\overline{\mathbf {P}}$ are still unknown. From (27), we have

$$\begin{aligned} \displaystyle \left[ \frac{1}{|\mathcal {B}_{\mathrm {0}} |} \, \int \limits _{\displaystyle \mathcal {I}_{\mathrm {0}}} {\overline{\mathbf {P}}} \ {\mathrm {dA}}\right]&= \frac{1}{\left[ L^3\right] } \left[ \overline{\mathbf {P}}\right] \, \left[ L^2\right] = \frac{\hbox {N}}{\hbox {mm}^{2}}\nonumber \\ \implies \frac{1}{\left[ L\right] } \left[ \overline{\mathbf {P}}\right]&= \frac{\hbox {N}}{\hbox {mm}^2} \end{aligned}$$

(32)

where $\left[ L\right]$ is the dimension of length. Thus, from (32) is can be concluded that if ${\left[ L\right] = \hbox {nm}}$, then ${\displaystyle \left[ \overline{\mu }\right] = \hbox {N} / \hbox {km}}$. In a similar fashion, if ${\left[ L\right] = \hbox {mm}}$, then ${\displaystyle \left[ \overline{\mu }\right] = \hbox {N} / \hbox {mm}}$ in order to maintain consistency in the units of the macroscopic stress.

Fig. 8

GICH size effect scaling: variation in the macroscopic Piola stress with filler particle size for the case rubber particles in epoxy matrix. Three different loading scenarios as described in Table 2 and a fixed interphase thickness $t=5$ have been considered

Fig. 9

GICH size effect scaling: variation in the macroscopic Piola stress with filler particle size for the case silica particles in epoxy matrix. Three different loading scenarios as described in Table 2 and a fixed interphase thickness $t=5$ have been considered

The aptness of the IECH approach for modeling size effects in the case of weaker filler particles in stiff material is further restricted partially, since there occur geometrical instabilities, in the form of kinks under shear loading conditions as shown in Fig. 7, corresponding to the irregularities in the size effect scaling curves, cf. Fig. 4b. Such behavior is prominent for active interfaces, i.e. smaller particle sizes and high interface stiffness, cf. Fig. 7a, b, respectively. It is critical to note that the instabilities are not of numerical nature and are hence observed regardless of the employed finite element mesh. The appearance of kinks is indeed due to buckling of parts of the stiff interface, which are under excessive compression, for instance the top-right corner of the particle in Fig. 7b. This buckling behavior could be attributed to the absence of resistance from the surrounding filler material whose stiffness is negligible as compared to that of the active interface. Such a buckling behavior could sometimes also be useful, for instance in the study of wrinkling as required for modeling cortical folding [19]. A Silica filler particle, on the other hand, cf. Fig. 7c, owing to its much higher stiffness, does not deform much even in the presence of an active interface. Consequently no such instabilities are observed. Thus, it could concluded that an energetic elastic interface can only sustain limited compressive loads.

5.2 Size effect scaling for GICH

In order to decipher the size effect scaling behavior for the GICH approach, we consider a sequence of RVEs obtained by scaling the RVE of Fig. 2b similar to the previous sub-section. The key here is that the interphase thickness remains constant as the particle and hence RVE size increases. Consequently, the filler volume fraction remains constant while the so-called interphase volume fraction reduces significantly for larger particle sizes, which eventually leads to a much less effective interphase in terms of its impact on the overall behavior of the composite material. Further control over the effectiveness of the graded interphase is attained by means of the model parameters, namely the interpolation exponent n and the interphase stiffness scaling factor $\alpha$, as described previously in Sect. 4. It is worth mentioning that regardless of the spatial dimension of the problem, i.e. 2D or 3D, the GICH approach requires the same number of model parameters, as opposed to the IECH approach where the number of requisite constitutive parameters depends on the spatial dimension of the problem at hand. Furthermore, since all geometrical entities involved here are of the same spatial dimension, this approach is essentially free from the question of appropriate units for the physical quantities, which was not the case with the interface modeling approach. However, for the sake of consistency, all the lengths, i.e. the RVE edge length a, the filler radius $r_f$ and the interphase thickness t, are assumed to have the units of nm.

Fig. 10

Influence of the interphase thickness on size effect scaling for volumetric loading

Revisiting the relations (28) and (29), which describe the continuous variation in material properties within the interphase region, it is apparent that in case of weak filler particles in a stiff matrix, which corresponds to the case of rubber inclusions in epoxy matrix, a much higher value of the factor $\alpha$ (resulting in $\alpha \, E_f > E_m$) is necessary in order to achieve the desired stiffening effect, cf. Fig. 3. In this study, a value of $\alpha =20000.0$ has been employed for the case of rubber inclusions. Silica inclusions, on the other hand, being already much stiffer than the epoxy matrix, require no further interphase stiffness scaling and hence $\alpha =1.0$ suffices for such a composite material.

In the same vein as the previous subsection, we begin the discussion by introducing the size effect scaling curves as presented in Figs. 8 and 9, respectively for a fixed interphase thickness of $5\hbox { nm}$ and the two filler-matrix extremes described above. Furthermore, the quantities on the x- and y-axis also remain the same as for the previous approach, i.e. the relevant component of the macroscopic Piola stress $\mathbf {P}^{\mathrm {M}}$ is plotted against the decreasing filler particle radius $r_f$. The effectiveness of the continuously graded interphase is controlled by means of the interpolation exponent n, a very large value of which, i.e. $n \rightarrow \infty$, results in the interphase region having the same material properties as those of the matrix region, cf. (28) for instance. Such behavior is exhibited by the standard two-phase continuum model where no interphase effects and hence no additional stiffening are observed, as can be verified, for instance, by the curves corresponding to $n=10^5$ in Fig 8a–c for each of the different loading scenarios considered. A very small n, on the other hand, indicates the scenario where the interphase region has properties corresponding to those of filler material scaled with the interphase stiffness scaling factor $\alpha$. By virtue of the appropriate choice of $\alpha$ as explained above, the macroscopic Piola stress increases significantly for reduced particle sizes, especially for smaller values of the interpolation exponent. Quite significant size effect scaling is observed in case of rubber particles in epoxy matrix for all three loading scenarios, cf. Fig. 8a–c. Considering for instance the volumetric loading scenario depicted through Fig. 8c, the macroscopic stress component $P^{\mathrm {M}}_{11}$ more than doubles from its lowest value of $432\hbox {MPa}$, which corresponds to an inactive interphase, to $881\hbox {MPa}$ as the interphase becomes is fully activated.

A conspicuous increase in the macroscopic stresses is also observed for the case of silica particles in epoxy matrix as depicted in Fig. 9a–c. The macroscopic stresses scale, in the case of uni-axial tensile loading, from $563\hbox { MPa}$ to $621\hbox { MPa}$ as the interpolation exponent reduces from $1e^5$ to $1e^{-3}$, cf. Fig. 9a. This size effect scaling becomes even more pronounced as the thickness of the interphase region increases from $5\hbox { nm}$ to $10\hbox { nm}$, cf. Fig. 10. It is noteworthy that in the case of volumetric loading, the peak macroscopic stress increases from $1016\hbox { MPa}$ to $1131\hbox { MPa}$ for the case of silica fillers and even more so, for the case of rubber inclusions in epoxy matrix, i.e. from $882\hbox { MPa}$ to $1042\hbox { MPa}$. Since the rubber inclusions are much weaker than the epoxy matrix, increasing the thickness of the stiff interphase results in an even stiffer overall composite behavior. Hence, a thicker interphase culminates in greatly pronounced size effect scaling for rubber inclusions as can also be observed in a visual comparison of the graphs in Fig. 10a, b. A similar trend has also been observed for the other two loading scenarios not presented here.

Fig. 11

Influence of the interphase Poisson’s ratio variation on size effect scaling for volumetric loading

Another aspect of the graded interphase model is unravelled by studying the effect of keeping the interphase Poisson’s ratio constant with value $\nu _m$. Such an approach is usually adopted in literature (see [35] for instance) without a thorough critical reasoning, thus rendering this affect worthy of inspection. Figure 11 presents a comparison of the two aforementioned possibilities for variation in the interphase Poisson’s ratio. For the sake of brevity, only the results for volumetric loading are presented here since a similar pattern has also been observed for the other two loadcases. Interestingly, for the case of silica particles in epoxy matrix, a constant value of $\nu _i$ culminates in higher macroscopic stress as compared to the case where $\nu _i$ is varied according to (30). From Table 1, it is evident that $\nu _{\mathrm {E}} > \nu _{\mathrm {S}}$, thereby resulting in higher interphase Poisson’s ratio for the $\nu _i=\nu _m$ scenario and hence the higher macroscopic stresses. Contrarily, the constant Poisson’s ratio strategy leads to lower macroscopic stress for rubber particles in epoxy matrix, since $\nu _{\mathrm {E}} < \nu _{\mathrm {R}}$, cf. Fig. 11b. The difference in the macroscopic stress computed by these two strategies is particularly significant for intermediate values of the interpolation exponent such as $n=2.0$, even more so for the case of silica fillers, cf. Fig. 11a, b, respectively. The results seem to converge as $n \rightarrow 0$.

Furthermore, there occur instabilities for the case of rubber particles in epoxy matrix under shear loading, resulting in self-penetration, primarily for very large particle size and for very stiff interphases. This observation substantiates the missing (for $n=10^{-3}$) and incomplete curves (for $n \le 5.0$) in Fig. 8b. The cause of these instabilities is very much akin to those observed for the IECH approach with an identical filler-matrix combination, cf. Sect. 5.1. For very large particle sizes, i.e. with $t \ll a$, the interphase behaves similar to an interface and hence suffers from instabilities due to a lack of resistance from the surrounding filler material owing to its negligible stiffness as compared to that of the very stiff interphase. No such instabilities are observed, however, for cases where $t \not \ll a$ and thus the GICH technique is also suitable for modeling weak nanosized particles in a stiff matrix.

5.3 Comparison: IECH vs GICH

After having described the size effect modeling capabilities of the two approaches under different loading scenarios individually, and for different filler-matrix combinations, we now present a holistic comparison of IECH and GICH.

In contrast to the interface-energetics-based technique, the graded interphase model yields substantially better size effect scaling for the case of silica particles in epoxy matrix under all loading conditions considered, as can be observed in Fig. 12, There the size effect scaling curves resulting from IECH and GICH have been plotted on the same graph. Let us first examine the volumetric loading scenario in detail. Taking into account the case of a highly active interface/interphase, contrary to the meager scaling of macroscopic stresses from 936 to $948\hbox { MPa}$ in case of IECH, cf. Fig. 5c, a much more pronounced scaling of the macroscopic stresses from 936 to $1016\hbox { MPa}$ is observed for the case of GICH, cf. Fig. 9c. This observation clearly establishes the graded interphase modeling approach as the better suited candidate for capturing size effects in composite materials when the filler material is much stiffer than the matrix.

Fig. 12

Comparison of size effect scaling for the two approaches. Here, silica particles in epoxy matrix have been considered under volumetric loading

The other filler-matrix stiffness extreme studied in this paper, i.e. rubber in epoxy, on the other hand, leads to similar size effect scaling for the two approaches under uni-axial tension and volumetric loading as is evident by comparison of the relevant graphs in Figs. 4 and 8, respectively. Furthermore, as with IECH, the graded interphase approach also suffers from the problem of instabilities under shear loading, however for very large particle sizes as opposed to IECH, where the problem appears for small particle sizes which are of particular interest for nano-composite modeling. Consequently, the graded interphase approach is better suited for modeling nano-composites with weak inclusions in a stiffer matrix under generic loading conditions. The interface energetics technique, on the other hand, is only suited to loading scenarios where the compressive loads are limited.

Table 3 Number of DoFs in the 2D finite element system for the two approaches. Here, for the GICH approach, the number of DoFs corresponding to $t=5\hbox { nm}$ has been reported

Fig. 13

Structured hex meshes for the 3D RVEs with $a=100, r=25$ (and $t=5$), employed for IECH and GICH, respectively

In contrast to the interface-energetics-based approach, the size effect scaling curves for the graded interphase model do not converge to a particular value as the filler particle size decreases. This behavior is observed for all of the loading scenarios considered and irrespective of the relative stiffness between the matrix and the filler phases, cf. Figs. 5c and 9c for instance. A flattening of the size effect scaling curves for the GICH approach could be achieved by making the interphase thickness dependent upon the filler particle radius as the particle sizes reduced below a certain limit value. Another important aspect is that of the computational cost associated with each of these approaches. In the case of IECH, the finite element mesh, and hence the number of DoFs to be solved for, remains unchanged as the filler particle size increases. However, as depicted in Table 3, for GICH, the number of DoFs increases significantly with an increase in particle size, since the size of the interphase region to be resolved increases. The situation would become even worse for the case of 3D RVEs. Furthermore, the highly graded mesh refinement carried out near the interphase region in order to minimize the total number of DoFs culminates in a much higher condition number of the underlying tangent-stiffness matrix, thereby increasing further the computational effort required to solve the corresponding SLE. However, in case of small (nanosized) filler particles, a scenario of practical interest for modeling size effects in polymer nano-composites, the computational overheads of the graded interphase model are much less significant, cf. the first few rows of Table 3. Consequently, the GICH approach turns out to be a computationally feasible option for modeling size effects in polymer nano-composites.

5.4 Comparison in 3D

In order to authenticate the generality of the results presented in the above sub-sections, across spatial dimensions, we consider a representative example in 3D. A detailed study covering the entire length-scale spectrum, as conducted here for the 2D case, would be quite compute intensive, especially for GICH as described above. Consequently, we consider an RVE with a fixed geometric configuration, cf. Fig. 13, and focus on the influence of the strength scale parameters. We further restrict ourselves to the interesting scenario of silica filler particle in epoxy matrix, subjected to volumetric loading. As described previously in Sect. 5.1, two independent material parameters namely $\overline{\mu }$ and $\overline{\kappa }$ (or alternatively $\overline{\nu }$) are required to characterize an energetic interface in 3D. For the GICH approach, we set $\alpha = 1.0$ as before and study the influence of the polynomial exponent n on the macroscopic stress.

Figure 14 illustrates the variation of the relevant macroscopic Piola stress component, i.e. $P^{\mathrm {M}}_{11}$ with the respective strength scale parameters for the two approaches. Focusing on Fig. 14a, i.e. the results for IECH, the resulting stresses $P^{\mathrm {M}}_{11}$ for different admissible values of the interface material parameters $(\overline{\mu }, \overline{\nu })$ have been presented. Alike the 2D scenario, the IECH approach shows minimal scaling of the macroscopic Piola stress with respect to the increasing interface stiffness, cf. Figs. 14a and  5c, respectively. The results for the GICH approach, as expected from the 2D comparison of Sect. 5.3, are in sharp contrast to those for the interface energetics based technique. As depicted in Fig. 14b, a significant increase in the macroscopic stresses is observed as $n \rightarrow 0$. Furthermore, the influence of spatial variation of the interphase Possion’s ratio, is in full agreement with the observations made for the 2D scenario, i.e. a constant $\nu _i = \nu _m$ culminates in relatively higher macroscopic stress than those obtained when $\nu _i$ is computed using the relation (30), cf. Fig. 11a and  14b, respectively. Consequently, the graded interphase approach, despite its relatively higher computational cost, qualifies itself as the method of choice for modeling size effects in nano-composites even in 3D scenario.

Fig. 14

Influence of the respective strength scale parameters for the two approaches in 3D. A fixed micro-structure configuration, comprising of silica particle in epoxy matrix with $a=100,r=25$ (and $t=5$), and subjected to volumetric loading has been considered

6 Conclusions and outlook

Extensions to the standard first-order computational homogenization scheme, aiming at modeling size effects in polymer composites, have been presented, and their suitability for different filler-matrix stiffness combinations has been thoroughly evaluated. An assessment of the interface-energetics-based extension of the standard computational homogenization scheme, also referred to in this work as the IECH, approach has been presented. In view of the representative numerical examples studied in this paper, the realm of application for the IECH approach appears to be limited to cases where the filler material is much weaker than the matrix and further to a subset of the admissible loading conditions. This is primarily due to the fact that the energetic interface is not very effective at providing further stiffening effect in cases where the filler particle is already much stiffer than the matrix material. Furthermore, for the case of weaker filler particles, a highly active interface, particularly for smaller particle sizes, tends to buckle if the compressive loads are above a certain limit value.

An alternative approach considering a continuous variation in the material properties within the interphase region has been introduced in this paper and has been named the Gradient-Interphase-Enhanced Computational Homogenization (GICH). Extensive numerical experimentation and an in-depth analysis of the obtained results demonstrate the appropriateness of this technique for the mechanical modeling of nano-composites irrespective of the relative stiffness of filler material with respect to that of the matrix. This comes at a significantly higher computational cost however. Complete resolution of the interphase layer implies that the total number of DoFs in the finite element system is considerably larger than for the IECHs approach, even more so for 3D RVEs containing several inclusions. The high condition number of the stiffness matrix emanating from the excessive local refinement in the interphase region further raises the computational cost of this technique. For the practically relevant application, i.e. modeling of nano-composites, the requisite computational effort for GICH is not much higher than IECH. In a very recent work [38], an approach based on the idea of generalized weighted interfaces, and aiming to approximate the graded interphase’s effects through an generalized energetic interface with arbitrary interface location, has been presented.

As with all other parameter-based computational models, specification of the appropriate model parameters for GICH by means of validation with experimental observations would be the natural step forward. Moreover, it would be of potential merit to apply the technique to more complicated, random micro-structures, especially in 3D, as is already in progress. A further extension could be the use of such graded interphases for modeling fracture behavior of particulate composites.

Appendices

Appendix A: Formulation of ${\overline{\mathbf {C}}}^{-1}$

Here, the derivation of the term ${\overline{\mathbf {C}}}^{-1}$ is outlined, since it is required for computing the interface stresses, cf. (23) and (24), and the interface tangent modulus, cf (45) and (46).

Consider again the interface $\mathcal {I}_{\mathrm {0}}$ depicted in 2a, and let $\overline{\mathbf {n}}= \overline{\mathbf {n}}\left( \overline{\mathbf {x}}\right)$ be the spatial configuration normal at any point P on the interface. The interface unit tensor in the spatial configuration, denoted by $\overline{\mathbf {i}}$ is then defined in a similar fashion to its material counterpart, cf. (13), as

$$\begin{aligned} \overline{\mathbf {i}}:=\mathbf {i}- \overline{\mathbf {n}}\otimes \overline{\mathbf {n}}, \end{aligned}$$

(33)

where $\mathbf {i}$ denotes the second-order unit tensor in the spatial configuration. Further, in line with (14), the interface gradient operator in the spatial configuration is defined as

$$\begin{aligned} \overline{\mathrm {grad}}\, \left\{ \overline{\bullet } \right\}&:=\mathrm {grad}\, \left\{ \overline{\bullet } \right\} \cdot \overline{\mathbf {i}}, \end{aligned}$$

(34)

with $\displaystyle \mathrm {grad}\, \left\{ \bullet \right\} :=\frac{\partial \, \left\{ \bullet \right\} }{\partial \mathbf {x}}$.

Next, let ${\displaystyle \varvec{\phi }: \mathbf {x}\rightarrow \mathbf {X}}$ and ${\displaystyle \overline{\varvec{\phi }}: \overline{\mathbf {x}}\rightarrow \overline{\mathbf {X}}}$ be the inverse deformation maps for the bulk and the interface, respectively. The inverse deformation gradients for the bulk and the interface are then defined as

$$\begin{aligned} \mathbf {f}:=&\ \mathrm {grad}\, \varvec{\phi } = {\mathbf {F}}^{-1} \end{aligned}$$

(35)

and

$$\begin{aligned} \overline{\mathbf {f}}:=&\ \overline{\mathrm {grad}}\, \overline{\varvec{\phi }}\nonumber \\ =&\ \overline{\mathrm {grad}}\, \varvec{\phi } \qquad \left[ \overline{\varvec{\phi }}= \varvec{\phi }\big |_{\mathcal {I}_{\mathrm {0}}}\text { for elastic interface}\right] \nonumber \\ =&\ \mathrm {grad}\, \varvec{\phi } \cdot \overline{\mathbf {i}}\qquad [\text {using}~(34)]\nonumber \\ \implies \overline{\mathbf {f}}=&\ \mathbf {f}\cdot \overline{\mathbf {i}}\qquad [\text {using}~35] \end{aligned}$$

(36)

respectively. The spatial unit normal $\mathbf {n}$ needed to compute $\overline{\mathbf {i}}$ using (33) is related to the material one, denoted as $\mathbf {N}$, by the relation [44]

$$\begin{aligned} \mathbf {n}= \mathbf {N}\cdot \mathbf {f}\, {|\mathbf {N}\cdot \mathbf {f}|}^{-1}. \end{aligned}$$

(37)

Finally, the inverse ${\overline{\mathbf {C}}}^{-1}$ of the interface right Cauchy-Green strain tensor is defined as

$$\begin{aligned} {\overline{\mathbf {C}}}^{-1} :=&\ \overline{\mathbf {f}}\cdot {\overline{\mathbf {f}}}^{t} = \ \left[ \mathbf {f}\cdot \overline{\mathbf {i}}\right] \cdot {\left[ \mathbf {f}\cdot \overline{\mathbf {i}}\right] }^{t} \qquad \left[ \text {using}~(36)\right] ,\nonumber \\ =&\ \mathbf {f}\cdot \overline{\mathbf {i}}\cdot {\overline{\mathbf {i}}}^{t} \cdot {\mathbf {f}}^{t}, \end{aligned}$$

(38)

and since by defintion of $\overline{\mathbf {i}}$ (cf. (33)) it holds that ${\overline{\mathbf {i}}}^{t} = \overline{\mathbf {i}}$ and $\overline{\mathbf {i}}\cdot \overline{\mathbf {i}}= \overline{\mathbf {i}}$, thereby yielding

$$\begin{aligned} \implies {\overline{\mathbf {C}}}^{-1} =&\ \mathbf {f}\cdot \overline{\mathbf {i}}\cdot {\mathbf {f}}^{t} \end{aligned}$$

(39)

which is the required expression.

B Consistent tangent modulus

1.1 B.1 Tangent modulus for bulk

A consistent tangent modulus as required for quadratic convergence of the Newton-Raphson scheme employed to solve the nonlinear system of equations is presented here. The tangent modulus is given as

$$\begin{aligned} \mathbb {C}= 2 \, \frac{\partial \mathbf {S}}{\partial \mathbf {C}} = \underbrace{2 \, \frac{\partial \mathbf {S}_{\mathrm {vol}}}{\partial \mathbf {C}}}_{\displaystyle =:\mathbb {C}_{\mathrm {vol}}} + \underbrace{2 \, \frac{\partial \mathbf {S}_{\mathrm {iso}}}{\partial \mathbf {C}}}_{\displaystyle =:\mathbb {C}_{\mathrm {iso}}}, \end{aligned}$$

(40)

with the two parts computed as

$$\begin{aligned} \mathbb {C}_{\mathrm {vol}} =&\ \kappa \, \left[ J^2 \, {\mathbf {C}}^{-1} \otimes {\mathbf {C}}^{-1} + \left[ J^2 - 1\right] \, \frac{\partial {\mathbf {C}}^{-1} }{\partial \mathbf {C}} \right] \end{aligned}$$

(41)

and

$$\begin{aligned} \mathbb {C}_{\mathrm {iso}} =&\ -\frac{2}{3} \,\mu \, J^{ -\frac{2}{3}} \, \left[ \mathbf {I}\otimes {\mathbf {C}}^{-1} + {\mathbf {C}}^{-1} \otimes \mathbf {I}\right] \nonumber \\&-\frac{2}{3} \,\mu \, J^{-\frac{2}{3}} \left[ \mathbf {C}: \mathbf {I}\right] \left[ \frac{\partial {\mathbf {C}}^{-1}}{\partial \mathbf {C}} - \frac{1}{3} \, {\mathbf {C}}^{-1} \otimes {\mathbf {C}}^{-1} \right] \\ \nonumber \end{aligned}$$

(42)

respectively. Here, the fourth order tensor $\displaystyle \frac{\partial {\mathbf {C}}^{-1} }{\partial \mathbf {C}}$ can be computed using the relation

$$\begin{aligned} \frac{\partial {\mathbf {A}}^{-1}}{\partial \mathbf {A}} = -\frac{1}{2} \left[ {\mathbf {A}}^{-1} \, \overline{\otimes }\, {\mathbf {A}}^{-1} + {\mathbf {A}}^{-1} \, \underline{\otimes }\, {\mathbf {A}}^{-1}\right] , \end{aligned}$$

(43)

where the non-standard tensor products are defined as

$$\begin{aligned} \mathbf {A}\, \overline{\otimes }\, \mathbf {A}= A_{IK} \, A_{LJ} \, \mathbf {E}_I \otimes \mathbf {E}_J \otimes \mathbf {E}_K \otimes \mathbf {E}_L \end{aligned}$$

and

$$\begin{aligned} \mathbf {A}\, \underline{\otimes }\, \mathbf {A}= A_{IL} \, A_{KJ} \, \mathbf {E}_I \otimes \mathbf {E}_J \otimes \mathbf {E}_K \otimes \mathbf {E}_L, \end{aligned}$$

with $\mathbf {A}$ denoting an invertible, second-order tensor, while $\mathbf {E}_I$ represents the coordinate system unit vector (see, for e.g. [21] for details).

1.2 B.2 Tangent modulus for interface

The consistent tangent modulus for the interface, denoted by $\overline{\mathbb {C}}$, is given as

$$\begin{aligned} \overline{\mathbb {C}}= 2 \, \frac{\partial \overline{\mathbf {S}}}{\partial \overline{\mathbf {C}}} = \underbrace{2 \, \frac{\partial \mathbf {S}_1}{\partial \overline{\mathbf {C}}}}_{\displaystyle =:\overline{\mathbb {C}}_1} + \underbrace{2 \, \frac{\partial \mathbf {S}_2}{\partial \overline{\mathbf {C}}}}_{\displaystyle =:\overline{\mathbb {C}}_2}, \end{aligned}$$

(44)

where the two parts can be determined using

$$\begin{aligned} \overline{\mathbb {C}}_1 =&\ \overline{\kappa }\, \left[ \overline{J}^2 \, {\overline{\mathbf {C}}}^{-1} \otimes {\overline{\mathbf {C}}}^{-1} + \left[ \overline{J}^2 - 1\right] \, \frac{\partial {\overline{\mathbf {C}}}^{-1} }{\partial \overline{\mathbf {C}}} \right] \end{aligned}$$

(45)

and

$$\begin{aligned} \overline{\mathbb {C}}_2= & {} \ -\frac{2}{\overline{\dim }} \,\overline{\mu }\, \overline{J}^{\left[ {-2} / {\overline{\dim }}\right] } \, \left[ \overline{\mathbf {I}}\otimes {\overline{\mathbf {C}}}^{-1} + {\overline{\mathbf {C}}}^{-1} \otimes \overline{\mathbf {I}}\right] \nonumber \\&-\frac{2}{\overline{\dim }} \,\overline{\mu }\, \overline{J}^{\left[ {-2} / {\overline{\dim }}\right] } \left[ \overline{\mathbf {C}}: \overline{\mathbf {I}}\right] \left[ \frac{\partial {\overline{\mathbf {C}}}^{-1}}{\partial \overline{\mathbf {C}}} - \frac{1}{\overline{\dim }} \, {\overline{\mathbf {C}}}^{-1} \otimes {\overline{\mathbf {C}}}^{-1} \right] ,\nonumber \\ \end{aligned}$$

(46)

respectively. Here, the fourth order tensor $\displaystyle \frac{\partial {\overline{\mathbf {C}}}^{-1} }{\partial \overline{\mathbf {C}}}$ can be computed using (43).