Effective properties of particulate nano-composites including Steigmann–Ogden model of material surface

Abstract

The objective of this work is inclusion of the Steigmann-Ogden interface in the Method of Conditional Moments to investigate the influence of surface effects on the effective properties of random particulate composites. The particular focus is centered on accounting for the surface bending stiffness. To this end, the notion of the energy-equivalent inhomogeneity developed for Gurtin–Murdoch interface is generalized to include the surface bending contribution. The crucial aspect of that generalization is identification of the formula defining energy associated with the surface bending. With the help of that formula, the real nano-particle and its surface are replaced by equivalent inhomogeneity with properties incorporating the surface effects. Closed-form expressions for the effective moduli of a composite with a matrix and randomly distributed spherical inhomogeneities are derived. The normalized shear moduli of nanoporous material as a function of void volume fraction is analyzed and evaluated in the context of other theoretical predictions.

1 Introduction

Interphases between the inhomogeneities and the matrix may have a very pronounced influence on the overall properties of composite materials. At the same time, their inclusion in mechanical (thermal, electric, etc.) analysis of those materials always entails additional complications whose level depends on the complexity of the interphase behavior and on the accuracy, with which that behavior is to be captured analytically. The interphases are typically three-dimensional continua but treating them as such is feasible only for simple geometry of the inhomogeneities and for simple loading conditions.

To cover more complex situations, most notably composites involving many interacting inhomogeneities, some effort has been invested to develop various simplified models of interphases ([1,2,3,4,5,6,7,8,9,10] among other). The most practical and popular of them are (arguably) the Gurtin–Murdoch [11, 12] model and the spring layer model (e.g., [1, 4, 5, 13,14,15,16]. The Gurtin–Murdoch and related models of surface elasticity have been used to study beams, plates, and shells, e.g. Miller and Shenoy [17], Altenbach and Eremeyev [18] among others.

The generalization of Gurtin–Murdoch model was proposed by Steigmann and Ogden [19, 20] who introduced the resistance of the surface to both stretch and bending. It means that the surface energy in the Steigmann-Ogden model includes both the surface strain tensor and the surface curvature tensor. The Steigmann–Ogden model has been used in Chhapadia et al. [21], Mohammadi and Sharma [22] to study bending of nano-sized cantilever beams. In these works the Steigmann–Ogden constants are determined by using combination of atomistic simulations and a simple continuum model.

In Javili et al. [23, 24], dell’Isola and Seppecher [25], dell’Isola et al. [26] it is demonstrated that the higher-gradient theories could entail surface tensors of stresses and couple stresses, as well as other stress resultants.

Within the Toupin–Mindlin formulation [27,28,29] of the strain gradient elasticity the mathematical study of static and dynamic boundary value problems with surface stresses described by Steigmann-Ogden model was presented in Eremeyev and Lebedev [30], Eremeev [31]. Similar variational formulation was used in the case of statics of classical elasticity with the Steigmann-Ogden model in Zemlyanova and Mogilevskaya [32]. The effect of curvature-dependent interfacial energy was also studied by Gao et al. [33] for finite deformation. The effective moduli of nanocomposites were analyzed in Gao et al. [34]. The effective properties of the isotropic particulate composites with Steigmann-Ogden interface are derived in Zemlyanova and Mogilevskaya [35]. Cylinder of finite length with Steigmann-Ogden interface was studied in Nazarenko et al. [36].

The main goal of this work is to show that the new concept of energy-equivalent inhomogeneity (EEI), recently presented in [37,38,39], permits direct evaluation of effective properties of nanomaterials, which includes Steigmann-Ogden interface model. In addition, we aim at evaluating the effectiveness of the proposed approach by means of numerical examples and comparisons with results obtained using other techniques, also those based on the Lurie solution for sphere [40] in which all governing equations in the inhomogeneity/interphase system are satisfied exactly. Lurie’s solution was also used for analysis of the effective behavior of the composites with spherical layered particles in [41, 42].

To illustrate that the presented notion of EEI can be used in combination with any method of evaluating the effective properties of composites, our choice is the method of conditional moments (MCM) rather than the self-consistent approach that is used more commonly.

This article is organized as follows. In the next section the governing equations of the problem of effective properties of random composites with interface is presented. Section 3 deals with the notion of energy-equivalent inhomogeneities. The main contribution of the proposed approach, which is determination of properties of the equivalent spherical inhomogeneity with Steigmann-Ogden interface, is presented in this Section in general terms, with the development quantifying those contributions relegated to Appendices 1 and 2. Even though the choice of the method used to develop the effective properties of an entire composite is in this work a secondary issue, to make it self-contained in the Sect. 4 a brief description of the MCM and its application to the problems with interphases are provided. The effective properties of particulate composite with Steigmann-Ogden interface are evaluated. Numerical results and comparisons are presented in Sect. 5, while discussion of the approach and conclusions are contained in the last Sect. 6.

2 Governing equations

Consider a representative macro-volume V consisting of a matrix with randomly distributed nano-inhomogeneities. Under conditions of uniform loading the macroscopic stress ${\overline{\mathbf{\sigma }}}$ and strain ${\overline{\mathbf{\varepsilon }}}$ are connected by following relations:

$$ {\overline{\mathbf{\sigma }}} = {\mathbf{C}}^{*} {\mathbf{:}}{\overline{\mathbf{\varepsilon }}}, $$

(1)

where ${\mathbf{C}}^{*}$ is the effective stiffness tensor, and the overbar denotes the operation of the statistical averaging.

For linear elastic materials the problem of finding the effective stiffness tensor requires the solution of the following set of equations.

• Equations of equilibrium:

  $$ \nabla \cdot \varvec{\upsigma} \left( {\mathbf{x}} \right) = 0, $$

  (2)

• Hooke’s law:

  $$ {{\varvec{\upsigma}}}{(}{\mathbf{x}}{)} = {\mathbf{C}}{\kern 1pt} {(}{\mathbf{x}}{)}{\mathbf{:}}{{\varvec{\upvarepsilon}}}{(}{\mathbf{x}}{),} $$

  (3)

• Linear kinematic relation:

  $$ {{\varvec{\upvarepsilon}}}\left( {\mathbf{x}} \right) = {\text{sym}}\left( {\nabla {\mathbf{u}}\left( {\mathbf{x}} \right)} \right), $$

  (4)

where ${\mathbf{x}}$ is the position vector of a micro-point.

In Eqs. (2)–(4) $\,{{\varvec{\upsigma}}}{(}{\mathbf{x}}{)}$ and ${{\varvec{\upvarepsilon}}}{(}{\mathbf{x}}{)}$ are the stress and strain tensors in the bulk material (matrix or inhomogeneity), ${\mathbf{u}}{(}{\mathbf{x}}{)}$ is the displacement vector; $\, \nabla$ is the three dimensional nabla operator; “·” identifies the single dot product of two tensors; the fourth-order tensor of elastic parameters ${\mathbf{C}}{(}{\mathbf{x}}{)}$ is a random, statistically homogeneous function of coordinates with a finite scale of correlation and linked to the inclusion and to the matrix properties via

$$ {\mathbf{C}}{(}{\mathbf{x}}{)} = {\mathbf{C}}_{1} H\left( {z{(}{\mathbf{x}}{)}} \right) + {\mathbf{C}}_{2} H\left( { - z{(}{\mathbf{x}}{)}} \right), $$

(5)

where $H$ is the Heaviside function and ${\mathbf{C}}_{1}$ and ${\mathbf{C}}_{2}$ denote the values of the tensors of elastic moduli in the inhomogeneities and in the matrix, respectively. The function $z{(}{\mathbf{x}}{)}$ is any function satisfying the following requirements:

$$ \begin{gathered} z{(}{\mathbf{x}}{)} > 0,\quad {\text{if}}\quad {\mathbf{x}} \in V_{1} \hfill \\ z{(}{\mathbf{x}}{)} = 0,\quad {\text{if}}\quad {\mathbf{x}} \in S_{I} \hfill \\ z{(}{\mathbf{x}}{)} < 0,\quad {\text{if}}\quad {\mathbf{x}} \in V_{2} , \hfill \\ \end{gathered} $$

(6)

where $V_{1}$ and $V_{2}$ are the domains of the inhomogeneities and the matrix, respectively and $S_{I}$ is the surface of inhomogeneities.

When the surface (interface) effects are described by the Steigmann-Ogden model [19, 20], Eqs. (2)–(4) need to be supplemented by the equilibrium equations an boundary conditions at the interface $S_{I}$ between the matrix and the nano-inhomogeneities. These relations are derived within the Toupin–Mindlin formulation of the strain gradient elasticity in Eremeev [31] and have following form:

displacements continuity on $S_{I}$

$$ \left[\kern-0.15em\left[ {{\mathbf{u}}{(}{\mathbf{x}}{)}} \right]\kern-0.15em\right]{\kern 1pt}_{{S_{I} \, }} = {\mathbf{0}}{;} $$

(7)

Surface equilibrium conditions on $S_{I}$

$$ \left[\kern-0.15em\left[ {{\varvec{\sigma}}{(}{\mathbf{x}}{)}} \right]\kern-0.15em\right]{\kern 1pt}_{{S_{I} \, }} \cdot {\mathbf{n}}\left( {\mathbf{x}} \right) = \, \nabla_{{S_{I} }} \cdot \left[ {{{\varvec{\upsigma}}}_{S} \left( {\mathbf{x}} \right) - \left( {\nabla_{{S_{I} }} \cdot {\mathbf{M}}_{S} \left( {\mathbf{x}} \right)} \right){\mathbf{n}}\left( {\mathbf{x}} \right)} \right] - 2H{\mathbf{n}}\left( {\mathbf{x}} \right) \cdot \left( {\nabla_{{S_{I} }} \cdot {\mathbf{M}}_{S} \left( {\mathbf{x}} \right)} \right){\mathbf{n}}\left( {\mathbf{x}} \right). $$

(8)

The unit vector ${\mathbf{n}}$ is normal to $S_{I}$. It is assumed, that at each interface ${\mathbf{n}}$ points away from the inhomogeneity. The square brackets indicate the jump of the field quantities across the interface, defined as their value on the side towards, which vector ${\mathbf{n}}$ is pointing minus their value on the side from which it is pointing;$\, \nabla_{{S_{I} }}$ is the surface nabla operator; $2H = {\text{tr}}{\mathbf{B}}{(}{\mathbf{x}}{)}$ is main curvature; ${\mathbf{B}}{(}{\mathbf{x}}{)} = - \nabla_{{S_{I} }} {\mathbf{n}}{(}{\mathbf{x}}{)}$ is the curvature tensor; the surface stress tensor ${{\varvec{\upsigma}}}_{S}^{{}}$ [12] is defined as

$$ {{\varvec{\upsigma}}}_{S} \left( {\mathbf{x}} \right) = \tau_{0} \mathop {{\mathbf{I}}_{S} }\limits^{2} + 2\left[ {\mu_{S} - \tau_{0} } \right]{{\varvec{\upvarepsilon}}}_{S} \left( {\mathbf{x}} \right) + \left[ {\lambda_{S} + \tau_{0} } \right]{\text{tr}}\left( {{{\varvec{\upvarepsilon}}}_{S} \left( {\mathbf{x}} \right)} \right)\mathop {{\mathbf{I}}_{S} }\limits^{2} + \tau_{0} \nabla {\mathbf{u}}\left( {\mathbf{x}} \right), $$

(9)

where ${{\varvec{\upvarepsilon}}}_{S}$ is the interface/surface strain tensor, $\mathop {\mathbf{I}_{S}}\limits^{2}$ represents the second-rank identity tensor in the plane tangent to the surface, $\tau_{0}^{{}}$ is the magnitude of the deformation-independent (residual) surface/interfacial tension (assumed “hydrostatic” and constant in Gurtin–Murdoch model), $\lambda_{S}^{{}}$, $\mu_{S}^{{}}$ are surface Lamé constants, while $\nabla_{S} {\mathbf{u}}{(}{\mathbf{x}}{)}$ denotes the surface gradient of the interface displacement field.

The surface couple stress tensor ${\mathbf{M}}_{S}$, which describes surface bending [19, 20] has following form:

$$ {\mathbf{M}}_{S} \left( {\mathbf{X}} \right) = 2\mu_{B} {{\varvec{\upkappa}}}_{S} \left( {\mathbf{X}} \right) + \lambda_{B} {\text{tr}}\left( {{{\varvec{\upkappa}}}_{S} \left( {\mathbf{X}} \right)} \right)\mathop {{\mathbf{I}}_{S} }\limits^{2}. $$

(10)

Here $\lambda_{B}$, $\mu_{B}$ are additional material parameters describing the bending stiffness of the material surface. The surface strain tensor ${{\varvec{\upvarepsilon}}}_{S}$ and the bending strain measure (tensor of changes of curvature) ${{\varvec{\upkappa}}}_{S}$

$$ {{\varvec{\upvarepsilon}}}_{S} \left( {\mathbf{x}} \right) = {\text{sym}}\left( {\mathop {{\mathbf{I}}_{S} }\limits^{2} \left( {\mathbf{X}} \right) \cdot \nabla_{S} {\mathbf{u}}\left( {\mathbf{x}} \right)} \right); $$

(11)

$$ {{\varvec{\upkappa}}}_{S} = {\text{Sym}}\left( {\mathop {{\mathbf{I}}_{S} }\limits^{2} \left( {\mathbf{x}} \right)\nabla_{S} \vartheta \left( {\mathbf{x}} \right)} \right), $$

(12)

in which $\vartheta {(}{\mathbf{x}}{)}$ is displacement of the end of vector ${\mathbf{n}}{(}{\mathbf{x}}{)}$ due to rotation of the surface

$$ \vartheta \left( {\mathbf{x}} \right) = \nabla_{{S_{I} }} \left( {{\mathbf{n}}\left( {\mathbf{x}} \right) \cdot {\mathbf{u}}\left( {\mathbf{x}} \right)} \right) + {\mathbf{B}}\left( {\mathbf{x}} \right) \cdot {\mathbf{u}}\left( {\mathbf{x}} \right). $$

(13)

The system of differential Eqs. (2)–(6), (7) can be transformed to the system of integral equations with the help of Green’s function ${\mathbf{G}}{(}{\mathbf{x}}{)}$, which is defined by the following boundary-value problem

$$ div\;\left( {{\mathbf{C}}_{c} {\mathbf{:}}\nabla {\mathbf{G}}{(}{\mathbf{x}}{)}} \right) +\updelta {(}{\mathbf{x}}{)}\,\mathop {\mathbf{I}}\limits^{2} = {\mathbf{0}} \ldots {\mathbf{G}}{(}{\mathbf{x}}{)}\left| {_{\infty } } \right. = {\mathbf{0}}, $$

(14)

where ${\mathbf{C}}_{c}$ is the constant tensor describing elasticity of the selected reference medium, $\updelta \left( {\mathbf{x}} \right)$ denotes the Dirac delta function and $\mathop {\mathbf{I}}\limits^{2}$ is the identity tensor of rank two. Then the fluctuations in the displacement field within the entire region V are described by the following formula (see [37, 43]:

$$ {\bf{u}}^{0} \left( {\mathbf{x}} \right) = \int\limits_{{V_{y} }} {{\bf{G}\left( {x - y} \right)} \cdot div\left( {{\bf{C}}^{0} \left( {\bf{y}} \right):\varepsilon \left( {\bf{y}} \right) - \varvec{\upbeta}} \right)dV_{y} } - \oint\limits_{{S_{y} }} {{\bf{G}}\left( {\bf{{x - y}}} \right) \cdot \nabla_{{S_{I} }} \cdot \left[ {\sigma_{S} \left( {\mathbf{x}} \right) + \left( {\nabla_{{S_{I} }} \cdot {\bf{M}}_{S} \left( {\mathbf{x}} \right)} \right){\bf{n}}\left( {\mathbf{x}} \right)} \right]ds_{y} } , $$

(15)

where ${{\varvec{\upbeta}}}$ is an arbitrary constant.

We consider the macro-volumes and macro-surfaces as infinite (as in physical experiments, they need to be considerably lager than the dimensions of inhomogeneities) and the boundary conditions have the form:

$$ {\mathbf{u}}^{0} \left( {\mathbf{x}} \right)\left| {_{\infty } } \right. = \bf{0}. $$

(16)

The linear kinematic relations of Eq. (4) combined with Eq. (15) and with the Gauss theorem leads to the following stochastically non-linear integral equations for the random strain field (see Nazarenko et al. [43] for more details)

$$ {{\varvec{\upvarepsilon}}}\left( x \right) = {\overline{\mathbf{\varepsilon }}} + {\mathbf{K}}\left( {{\mathbf{x}} - {\mathbf{y}}} \right) * \left[ {{\mathbf{C}}^{0} \left( {\mathbf{y}} \right):{{\varvec{\upvarepsilon}}}\left( {\mathbf{y}} \right)} \right] - {\text{sum}}\left\{ {\nabla_{x} \oint\limits_{{S_{I} }} {{\mathbf{G}}\left( {{\mathbf{x}} - {\mathbf{y}}} \right) \cdot \nabla_{{S_{I} }} \cdot \left[ {\sigma_{S} \left( {\mathbf{x}} \right) + \left( {\nabla_{{S_{I} }} \cdot {\mathbf{M}}_{S} \left( {\mathbf{x}} \right)} \right){\mathbf{n}}\left( {\mathbf{x}} \right)} \right]dS_{Y} } } \right\}. $$

(17)

The operator ${\mathbf{K}}{(}{\mathbf{x}} - {\mathbf{y}}{)}$ acts according to

$$ {\mathbf{K}}\left( {{\mathbf{x}} - {\mathbf{y}}} \right) * {{\varvec{\uppsi}}}\left( {\mathbf{y}} \right) = \int\limits_{V} {{\text{sym}}\left( {\nabla_{x} \left( {\nabla_{x} {\mathbf{G}}\left( {{\mathbf{x}} - {\mathbf{y}}} \right)} \right)} \right):\left( {{{\varvec{\uppsi}}}\left( {\mathbf{y}} \right) - {\overline{\mathbf{\psi }}}} \right)dV_{y} } , $$

(18)

while ${\overline{\mathbf{\varepsilon }}}$ and ${\overline{\mathbf{\psi }}}$ are mean (expectation) values of ${{\varvec{\upvarepsilon}}}{(}{\mathbf{x}}{)}$ and ${{\varvec{\uppsi}}}{(}{\mathbf{y}}{)}$. ${\mathbf{C}}^{0} \left( {\mathbf{y}} \right)$ in Eq. (17) is the fluctuations in elastic constants and is defined as

$$ {\mathbf{C}}^{0} \left( {\mathbf{y}} \right) = {\mathbf{C}}\left( {\mathbf{y}} \right) - {\mathbf{C}}_{c} . $$

(19)

The issues related to selection of the tensor ${\mathbf{C}}_{c}$ were discussed in [37, 43,44,45] and the Reader is referred to these articles.

The way to account for the surface effects described by Eq. (17) is rigorous. However, evaluation of the surface integral containing $\left[ {{{\varvec{\upsigma}}}_{S} {(}{\mathbf{x}}{)} + (\nabla_{{S_{I} }} \cdot {\mathbf{M}}_{S} {(}{\mathbf{x}}{)}){\mathbf{n}}{(}{\mathbf{x}}{)}} \right]$ in the right-hand side of Eq. (17) is a rather complex problem for the case of arbitrary loading. In Nazarenko et al. [43], this was done only under some approximating assumptions, valid only for overall volumetric deformations. To bypass that difficulty an approach based on MCM in combination with a notion of the EEI was proposed in [37,38,39]. In this approach the contribution of the surface integral in Eq. (17) is accounted for by a proper adjustment of the bulk properties of the inhomogeneities, and the material is then analyzed as standard composite, i.e. involving no surfaces possessing their own mechanical characteristics.

According with the above remarks, in the EEI approach Eq. (17) is replaced by the equation for the energy-equivalent system (see [37]

$$ {{\varvec{\upvarepsilon}}}\left( {\mathbf{x}} \right) = {\overline{\mathbf{\varepsilon }}} + {\mathbf{\rm K}}\left( {{\mathbf{x}} - {\mathbf{y}}} \right)*\left[ {{\tilde{\mathbf{C}}}^{0} \left( {\mathbf{y}} \right):{{\varvec{\upvarepsilon}}}\left( {\mathbf{y}} \right)} \right], $$

(20)

in which the surface integral of Eq. (17) is incorporated into ${\tilde{\mathbf{C}}}^{0} \left( {\mathbf{x}} \right)$ as a result of changing properties of the inhomogeneities from ${\mathbf{C}}_{1}$ to ${\mathbf{C}}_{{{\text{eq}}}}$. Detailed developments leading to ${\tilde{\mathbf{C}}}^{0} \left( {\mathbf{x}} \right)$ are presented in Sect. 3. Tensor ${\tilde{\mathbf{C}}}^{0} \left( {\mathbf{x}} \right)$ is defined as:

$$ {\tilde{\mathbf{C}}}^{0} \left( {\mathbf{x}} \right) = {\tilde{\mathbf{C}}}\left( {\mathbf{x}} \right) - {\mathbf{C}}_{C} \;{\text{and}}\;{\tilde{\mathbf{C}}}\left( {\mathbf{x}} \right) = {\mathbf{C}}_{{{\text{eq}}}} H\left( {z\left( {\mathbf{x}} \right)} \right) + {\mathbf{C}}_{2} H\left( { - z\left( {\mathbf{x}} \right)} \right). $$

(21)

Equation (20) is solved in an averaged sense by the MCM to extract the average strains. As described subsequently, this leads to the determination of average stresses and the effective elastic constants.

3 Energy-equivalent inhomogeneity

3.1 The concept of energy-equivalent inhomogeneity

The idea of the EEI is simple: given that the interfaces in Steigmann-Ogden [19, 20] model are coherent (no jump in displacements across the interface), the surface strains of the interface can be related to the strains in the bulk material of the nano-inhomogeneity (or of the matrix). Furthermore, since the two-point approximation of the MCM—typically used in analysis based on that approach—effectively implies that the strains within the nano-inhomogeneities are assumed constant (for details see [43, 38]), it is convenient to relate the interface strains to those of the inhomogeneities (not those of the matrix, where the strain field is not assumed constant). Under those conditions it is reasonable to expect that the overall stiffness characteristics of the analyzed nano-materials should remain essentially the same whether the stiffness provided by the interface is treated independently (as done in [43] or lumped together with the original stiffness of the nano-ihomogeneity. With such modified stiffnesses of all nano-inhomogeneities, combining their original properties and those of their interfaces, the composite may be considered as standard, in which there is no need for an independent inclusion of the surface effects.

The advantage of the approach based on the notion of EEI is that it bypasses all the technical difficulties caused by the presence of independently treated interfaces, encountered for example in Nazarenko et al. [43]. All the formulas for the effective properties of standard random heterogeneous materials (i.e. those not involving interfaces), that can be relatively easily developed using MCM [44, 45], become then directly applicable to nano-materials (i.e. involving interfaces), as long as the properties of inhomogeneities are properly modified to include the surface effects.

To find the modified properties of the EEI it is postulated that, for arbitrary—but homogeneous—deformation, those modified properties render the elastic energy of the EEI that is equal to the sum of the energies of the unmodified inhomogeneity and the energy of its interface. The postulate made here yields:

$$ \frac{1}{2}V_{1} \left( {{\overline{\mathbf{\varepsilon }}}:{\mathbf{C}}_{{{\text{eq}}}} :{\overline{\mathbf{\varepsilon }}}_{1} } \right) = \frac{1}{2}V_{1} \left( {{\overline{\mathbf{\varepsilon }}}:{\mathbf{C}}_{1} :{\overline{\mathbf{\varepsilon }}}_{1} } \right) + E_{S} . $$

(22)

Here ${\mathbf{C}}_{{{\text{eq}}}}$ is the stiffness tensor of the modified inhomogeneity to be found, ${\mathbf{C}}_{1}$ is the stiffness tensor of the original inhomogeneity; $V_{I}$ is the volume of the inhomogeneity and $S_{I}$ is its surface; ${\overline{\mathbf{\varepsilon }}}_{1}$ is the constant strain tensor within the volume of the inhomogeneity; $E_{S}$ is the surface energy.

In the case of Steigmann-Ogden interface Eqs. (7)–(13) the surface energy can be represented as

$$ E_{S} = U_{T} + U_{B} , $$

(23)

where $U_{T}$ and $U_{B}$ are the energies related to the surface tension and the surface bending

$$ U_{T} = \frac{1}{2}\int\limits_{{S_{I} }} {\int\limits_{{S_{I} }} {\left[ {2\overline{\mu }_{S} \varepsilon_{S} :\varepsilon_{S} + \overline{\lambda }_{S} {\text{tr}}\left( {\varepsilon_{S} } \right)^{2} + \tau_{0} \nabla_{S} {\mathbf{u}}:\nabla_{S} {\mathbf{u}}} \right]dS} } ; $$

(24)

$$ U_{B} = \frac{1}{2}\int\limits_{S} {\left[ {2\mu_{b} {{\varvec{\upkappa}}}:{{\varvec{\upkappa}}} + \lambda_{B} {\text{tr}}\left( {{\varvec{\upkappa}}} \right)^{2} } \right]dS} . $$

(25)

In Nazarenko et al. [39], it is shown, that for spherical and isotropic inhomogeneities, and for isotropic surface properties described by Gurtin–Murdoch model of material surface expressed in Eqs. (7)–(9) if ${\mathbf{M}}_{S} {(}{\mathbf{x}}{)} = 0$, the stiffness tensor ${\mathbf{C}}_{{{\text{eq}}}}$ of equivalent inhomogeneity is also isotropic and its bulk and shear moduli are defined as

$$ K_{{{\text{eq}}}} = K_{1} + \hat{K}_{T} ;\;{\upmu }_{{{\text{eq}}}} = {\upmu }_{1} + {\hat{\mu }}_{T} , $$

(26)

where $K_{1}$ and $\mu_{1}$ are the constants (bulk and shear moduli) of the original inhomogeneity, while $\hat{K}_{T}$ and $\hat{\mu }_{T}$ are (see details in [39]

$$ \hat{K}_{T} = 2\frac{{\left[ {2{\overline{\mu }}_{S} + 2{\overline{\lambda }}_{S} + \tau_{0} } \right]^{{}} }}{3a},\;\hat{\mu }_{T} = \frac{{7\overline{\mu }_{S} + \overline{\lambda }_{S} + 5\tau_{0} }}{5a}, $$

(27)

with $\overline{\lambda }_{S} = \lambda_{S} + \tau_{0}$, $\overline{\mu }_{S} = \mu_{S} - \tau_{0}$ appearing as a result of the surface tension contribution in Eq. (22) and Eq. (8); $a$ is radius of spherical inhomogeneity. If ${\mathbf{M}}_{S} \ne 0$ in Eqs. (7)–(9), the tensor ${\mathbf{C}}_{{{\text{eq}}}}$ should be also isotropic and its bulk and shear moduli can be defined as

$$ K_{{{\text{eq}}}} = K_{1} + \hat{K}_{T} + \hat{K}_{B} ;\;\mu_{{{\text{eq}}}} = \mu_{1} + \hat{\mu }_{T} + \hat{\mu }_{B} , $$

(28)

where $\hat{\mu }_{B}$ and $\hat{K}_{B}$ are the additional contribution of the surface bending.

The development, when the term ${\mathbf{M}}_{S}$ of Eq. (8) is neglected, was presented in [39]. Here the focus of evaluation of the properties of EEI is on the contribution of the surface bending. Use of the complete Eq. (8) in analysis may turn out to be important in some practical applications, where bending of surface should be accounted (e.g. [46]. Inclusion of the complete Eq. (8) within the framework of the EEI is outlined in the next subsection, with some supporting derivations presented in the related "Appendix 1".

3.2 Contribution of the surface bending to the energy of equivalent inhomogeneity

3.2.1 Evaluation of the surface energy related to the bending

The overarching idea pursued here is the same as that described in the preceding subsection. The only outstanding issue that needs to be addressed is how to determine the surface contribution in Eq. (25) in order to account for the presence of the surface bending in Eq. (8). To this end, the tensor of curvature changes will be evaluated first.

It is assumed, that the strains ${\overline{\mathbf{\varepsilon }}}_{1}$ that an inhomogeneity is subjected to are constant. Under those conditions the displacements in of the surface of that inhomogeneity can be expressed as

$$ {\mathbf{u}}\left( {\xi^{\Lambda } } \right) = {\overline{\mathbf{\varepsilon }}}_{1} \cdot {\mathbf{r}}\left( {\xi^{\Lambda } } \right), $$

(29)

where ${\mathbf{r}}\left( {\xi^{\Lambda } } \right)$ is the position vector of a point on that surface which is locally parameterized by $\xi^{\Lambda }$, $\Lambda \in \left\{ {\,1\,,\;2} \right\}$. Consequently, (cf. [47]

$$ \nabla_{S} {\mathbf{u}} = \left( {{{\varvec{\upvarepsilon}}}_{1} \cdot {\mathbf{r}}} \right)_{,\Delta } \otimes {\mathbf{G}}^{\Delta } = {\overline{\mathbf{\varepsilon }}}_{1} \cdot \left( {{\mathbf{r}}_{,\Delta } \otimes {\mathbf{G}}^{\Delta } } \right) = {\overline{\mathbf{\varepsilon }}}_{1} \cdot \left( {{\mathbf{G}}_{\Delta } \otimes {\mathbf{G}}^{\Delta } } \right) = {\overline{\mathbf{\varepsilon }}}_{1} \cdot \mathop {\mathbf{I}}\limits^{2}_{S} . $$

(30)

Here ${\mathbf{G}}_{\Delta } = {\mathbf{r}},_{\Delta }$ are the vectors of the natural basis associated with the parametrization $\xi^{\Delta }$ (tangent to the surface) and ${\mathbf{G}}^{\Delta }$ are the vectors of the dual, or reciprocal, basis also tangent to the surface) satisfying the condition ${\mathbf{G}}_{\Delta } \cdot {\mathbf{G}}^{\Lambda } = \delta_{\Delta }^{\Lambda }$ with $\delta_{\Delta }^{\Lambda }$ being the “Kronecker delta”.

The tensor of curvature changes is determined as

$$ {{\varvec{\upkappa}}} = {\text{sym}}\left( {\mathop {\mathbf{I}_{S}}\limits^{2} \cdot \nabla_{S} \vartheta } \right), $$

(31)

with

$$ \vartheta = {{\varvec{\upomega}}}_{N} \cdot {\mathbf{n}}, $$

(32)

where

$$ {{\varvec{\upomega}}}_{N} = - \left( {{\mathbf{n}} \cdot \nabla_{S} {\mathbf{u}}} \right) \otimes {\mathbf{n}}. $$

(33)

Considering Eq. (30) $\vartheta$ can be defined as

$$ \vartheta = - {\mathbf{n}} \cdot {\overline{\mathbf{\varepsilon }}}_{1} \cdot \mathop {\mathbf{I}_{S}}\limits^{2} = - {\overline{\mathbf{\varepsilon }}}_{1} :{\mathbf{n}} \otimes \mathop {\mathbf{I}_{S}}\limits^{2} , $$

(34)

which gives

$$ \nabla_{S} \vartheta = - {\overline{\mathbf{\varepsilon }}}_{1} :\nabla_{S} \left( {{\mathbf{n}} \otimes \mathop {\mathbf{I}_{S}}\limits^{2} } \right); $$

(35)

$$ {{\varvec{\upkappa}}} = - {\text{sym}}\left[ {\mathop {\mathbf{I}_{S}}\limits^{2} \cdot \left( {{\overline{\mathbf{\varepsilon }}}_{1} :\nabla_{S} \left( {{\mathbf{n}} \otimes \mathop {\mathbf{I}_{S}}\limits^{2} } \right)} \right)} \right] $$

(36)

The above formula indicate that ${\overline{\mathbf{\varepsilon }}}_{1}$ contracted with the first 2 vectors of $\nabla_{S} \left( {{\mathbf{n}} \otimes \mathop {\mathbf{I}_{S}}\limits^{2} } \right)$ and $\mathop {\mathbf{I}_{S}}\limits^{2}$ operate on the third vector of dyadic product in $\nabla_{S} \left( {{\mathbf{n}} \otimes \mathop {\mathbf{I}_{S}}\limits^{2} } \right)$. This means that multiplication by $\mathop {\mathbf{I}_{S}}\limits^{2}$ eliminate the ${\mathbf{n}} \otimes {\mathbf{G}}_{\Lambda } \otimes {\mathbf{n}} \otimes {\mathbf{G}}^{\Delta }$ (${\mathbf{G}}_{\Lambda } \bot {\mathbf{n}}$) and the remaining two parts are unchanged. So

$$ {{\varvec{\upkappa}}} = - {\text{sym}}\;\left[ {{\mathbf{I}}_{S} \cdot \left( {\nabla_{S} \left( {{\mathbf{n}} \otimes {\mathbf{I}}_{S} } \right)} \right)^{T} :{\overline{\mathbf{\varepsilon }}}_{1} } \right] $$

$$ \begin{aligned} {{\varvec{\upkappa}}} = & - {\text{sym}}\left[ {{\mathbf{I}}_{S} \cdot \left( {\nabla_{S} \left( {{\mathbf{n}} \otimes {\mathbf{I}}_{S} } \right)} \right)^{T} :{\overline{\mathbf{\varepsilon }}}_{1} } \right] \\ & = - {\text{sym}}\left[ {\left( { - B_{\Delta }^{\Pi } {\mathbf{G}}_{\Pi } \otimes {\mathbf{I}}_{S} \otimes {\mathbf{G}}^{\Delta } + B_{\Lambda \Delta } {\mathbf{n}} \otimes {\mathbf{n}} \otimes {\mathbf{G}}^{\Lambda } \otimes {\mathbf{G}}^{\Delta } } \right):{\overline{\mathbf{\varepsilon }}}_{1} } \right]. \\ \end{aligned} $$

(37)

Evaluation of the components of the tensors $B_{\Lambda \Delta }$ and $B_{\Delta }^{\Pi }$ is illustrated in "Appendix 1", where curvature tensors for spherical inhomogeneity of radius $a$ are given in Table 1.

Table 1 Curvature tensors for sphere of radius $a$

Considering values $B_{\Lambda \Delta }$ and $B_{\Lambda }^{\Delta }$ from Table

$$ {{\varvec{\upkappa}}} = - {\text{sym}}\left[ {\frac{1}{a}\left( {{\mathbf{G}}^{\Lambda } \otimes {\mathbf{G}}^{1} \otimes {\mathbf{G}}_{1} \otimes {\mathbf{G}}_{\Lambda } + {\mathbf{G}}^{\Lambda } \otimes {\mathbf{G}}^{2} \otimes {\mathbf{G}}_{2} \otimes {\mathbf{G}}_{\Lambda } - } \right.} \right.\left. {\left. {{\mathbf{G}}^{1} \otimes {\mathbf{G}}_{1} \otimes {\mathbf{n}} \otimes {\mathbf{n}} - {\mathbf{G}}^{2} \otimes {\mathbf{G}}_{2} \otimes {\mathbf{n}} \otimes {\mathbf{n}}} \right):{\overline{\mathbf{\varepsilon }}}_{1} } \right]. $$

(38)

Given that for ${\mathbf{G}}_{\Delta } \cdot {\mathbf{G}}_{\Lambda } = 0$ if $\Delta \ne \Lambda$, the curvilinear coordinates $\xi^{\Lambda }$ introduced to parameterize the surface of the inhomogeneities are orthogonal, the unit vectors ${\overline{\mathbf{G}}}_{\Delta } = \frac{{{\mathbf{G}}_{\Delta } }}{{\left| {\left. {{\mathbf{G}}_{\Delta } } \right|} \right.}} = {\overline{\mathbf{G}}}^{\Delta } = \frac{{{\mathbf{G}}^{\Delta } }}{{\left| {\left. {{\mathbf{G}}^{\Delta } } \right|} \right.}}$ can be substituted for ${\mathbf{G}}_{\Delta }$ and ${\mathbf{G}}^{\Delta }$ in the above equation and it can be rewritten in the form:

$$ \begin{aligned} {{\varvec{\upkappa}}} = & - {\text{sym}}\left[ {\frac{1}{a}\left( {\,\left[ {\,\left( {{\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} + } \right.} \right.{\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{2} + {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{1} + } \right.} \right. \\ & + {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} - {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\mathbf{n}} \otimes {\mathbf{n}} - \left. {\left. {{\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\mathbf{n}} \otimes {\mathbf{n}}} \right):{\overline{\mathbf{\varepsilon }}}_{1} } \right]. \\ \end{aligned} $$

(39)

Accounting for symmetry of ${\overline{\mathbf{\varepsilon }}}_{1}$ tensor of curvature changes can be defined as

$$ \begin{aligned} {{\varvec{\upkappa}}} = & - \frac{1}{a}\left[ {\left( {{\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} + } \right.} \right.\left( {{\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{1} + {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{2} } \right) \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{2} \\ & + {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} - {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\mathbf{n}} \otimes {\mathbf{n}} - \left. {\left. {{\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\mathbf{n}} \otimes {\mathbf{n}}} \right):{\overline{\mathbf{\varepsilon }}}_{1} } \right]. \\ \end{aligned} $$

(40)

Then the surface strain energy in the case of Steigmann-Ogden model of interface is defined as (see for details [31, 32])

$$ E_{S} = \frac{1}{2}\oint\limits_{{S_{I} }} {\left[ {2\overline{\mu }_{S} {{\varvec{\upvarepsilon}}}_{S} :{{\varvec{\upvarepsilon}}}_{S} + \overline{\lambda }_{S} tr\left( {{{\varvec{\upvarepsilon}}}_{S} } \right)^{2} + \tau_{0} \nabla_{S} {\mathbf{u}}:\nabla_{S} {\mathbf{u}} + 2\mu_{B} {{\varvec{\upkappa}}}:{{\varvec{\upkappa}}} + \lambda_{B} \left( {{\text{tr}}{{\varvec{\upkappa}}}} \right)^{2} } \right]dS} . $$

(41)

The first three terms of the above integrand are identical with the surface strain energy given by Gurtin–Murdoch model and properties of equivalent inhomogeneities related to these terms are determined in Nazarenko et al. [37,38,39]. The last two terms of Eq. (41) represent the surface strain energy related to the surface bending. In the next section, the working formula for the properties of the equivalent inhomogeneity is presented.

3.2.2 Constitutive tensor of the energy-equivalent inhomogeneity

Considering Eq. (40) last two terms of Eq. (41) are given as

$$ \begin{gathered} {{\varvec{\upkappa}}}:{{\varvec{\upkappa}}} = \frac{1}{{a^{2} }}{\overline{\mathbf{\varepsilon }}}_{1} :\left[ {{\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} + {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} + 2{\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{2} - {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\mathbf{n}} \otimes {\mathbf{n}}} \right. \hfill \\ \left. { - {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\mathbf{n}} \otimes {\mathbf{n}} - {\mathbf{n}} \otimes {\mathbf{n}} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} - {\mathbf{n}} \otimes {\mathbf{n}} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} + 2{\mathbf{n}} \otimes {\mathbf{n}} \otimes {\mathbf{n}} \otimes {\mathbf{n}}} \right]:{\overline{\mathbf{\varepsilon }}}_{1} ; \hfill \\ \end{gathered} $$

(42)

$$ \begin{gathered} {\text{tr}}\left( {{\varvec{\upkappa}}} \right)^{2} = \frac{1}{{a^{2} }}{\overline{\mathbf{\varepsilon }}}_{1} :\left[ {{\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} {\mathbf{ + \overline{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} {\mathbf{ + \overline{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} } \right. \hfill \\ \left. { - 2{\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\mathbf{n}} \otimes {\mathbf{n}} - 2{\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} {\mathbf{n}} \otimes {\mathbf{n}} - 2{\mathbf{n}} \otimes {\mathbf{n}} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} - 2{\mathbf{n}} \otimes {\mathbf{n}} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} + 4{\mathbf{n}} \otimes {\mathbf{n}} \otimes {\mathbf{n}} \otimes {\mathbf{n}}} \right]:{\overline{\mathbf{\varepsilon }}}_{1} \hfill \\ \end{gathered} $$

(43)

The surface energy of Eq. (41) is a sum of the surface tension and the surface bending Eqs. (23)–(25) and we focus only on the latter

$$ U_{B} = \frac{1}{2}\oint\limits_{S} {\left[ {2\mu_{B} {{\varvec{\upkappa}}}:{{\varvec{\upkappa}}} + \lambda_{B} {\text{tr}}\left( {{\varvec{\upkappa}}} \right)^{2} } \right]dS} = U_{{\mu_{B} }} + U_{{\lambda_{B} }} , $$

(44)

where $U_{{\mu_{B} }} = \frac{1}{2}\oint\limits_{S} {\left[ {2\mu_{B} {{\varvec{\upkappa}}}:{{\varvec{\upkappa}}}} \right]} dS$ is defined as

$$ \begin{aligned} U_{{\mu_{B} }} = & \frac{1}{{2a^{2} }}\oint\limits_{{S_{I} }} {\left[ {{\overline{\mathbf{\varepsilon }}}_{1} :\left( {{\overline{\mathbf{G}}}_{1} } \right. \otimes } \right.} \,{\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} + {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \\ & + 2{\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{2} - {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\mathbf{n}} \otimes {\mathbf{n}} \\ & - {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\mathbf{n}} \otimes {\mathbf{n}} - {\mathbf{n}} \otimes {\mathbf{n}} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} - {\mathbf{n}} \otimes {\mathbf{n}} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} + 2{\mathbf{n}} \otimes {\mathbf{n}} \otimes {\mathbf{n}} \otimes \left. {\mathbf{n}} \right){:}\left. {{\overline{\mathbf{\varepsilon }}}_{1} } \right]d\,S \\ \end{aligned} $$

(45)

and $U_{{\lambda_{B} }} = \oint\limits_{S} {\left[ {\lambda_{B} {\text{tr}}\left( {{\varvec{\upkappa}}} \right)^{2} } \right]dS}$ is

$$ \begin{aligned} U_{{\lambda_{B} }} = & \frac{1}{{2a^{2} }}\oint\limits_{{S_{I} }} {\left[ {{\overline{\mathbf{\varepsilon }}}_{1} :\left( {{\overline{\mathbf{G}}}_{1} } \right. \otimes } \right.} {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} + {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} + {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \\ & + {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} - 2{\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\mathbf{n}} \otimes {\mathbf{n}} \\ & - 2{\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\mathbf{n}} \otimes {\mathbf{n}} - 2{\mathbf{n}} \otimes {\mathbf{n}} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} - 2{\mathbf{n}} \otimes {\mathbf{n}} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} + 4{\mathbf{n}} \otimes {\mathbf{n}} \otimes {\mathbf{n}} \otimes \left. {\mathbf{n}} \right){:}\left. {{\overline{\mathbf{\varepsilon }}}_{1} } \right]d\,S. \\ \end{aligned} $$

(46)

These last formulas Eqs. (45), (46) can be put in Eq. (44) and the following form of the surface energy yields

$$ {\text{E}}_{S} = U_{T} + U_{{\mu_{B} }} + U_{{\lambda_{B} }} \, = {\overline{\mathbf{\varepsilon }}}_{1} :\left( {{\mathbf{K}}_{T} + {\mathbf{K}}_{{\mu_{B} }} + {\mathbf{K}}_{{\lambda_{B} }} } \right):{\overline{\mathbf{\varepsilon }}}_{1} , $$

(47)

where

$$ U_{T} = {\overline{\mathbf{\varepsilon }}}_{1} :{\mathbf{K}}_{T} :{\overline{\mathbf{\varepsilon }}}_{1} ;\;U_{{\mu_{B} }} = {\overline{\mathbf{\varepsilon }}}_{1} :{\mathbf{K}}_{{\mu_{B} }} :{\overline{\mathbf{\varepsilon }}}_{1} ;\;U_{{\lambda_{B} }} = {\overline{\mathbf{\varepsilon }}}_{1} :{\mathbf{K}}_{{\lambda_{B} }} :{\overline{\mathbf{\varepsilon }}}_{1} . $$

(48)

The last result for $E_{S}$ and the energy equivalence expressed by Eq. (22) leads to the following formula for the effective moduli of equivalent inhomogeneity:

$$ {\mathbf{C}}_{{{\text{eq}}}} = {\mathbf{C}}_{1} + \frac{1}{{V_{I} }}\left( {{\mathbf{K}}_{T} + {\mathbf{K}}_{{\mu_{B} }} + {\mathbf{K}}_{{\lambda_{B} }} } \right) $$

(49)

Then the problem of properties of equivalent inhomogeneities is reduced to evaluation of the components of the above tensors ${\mathbf{K}}_{{\mu_{B} }}$ and ${\mathbf{K}}_{{\lambda_{B} }}$, which is illustrated in "Appendix 2" Eqs. (98), (99). Contribution of the surface bending in this case is

$$ \hat{K}_{B} = 0,\;\hat{\mu }_{B} = \frac{{15\mu_{B} + 9\lambda_{B} }}{{5a^{3} }}, $$

(50)

where $\lambda_{B}$, $\mu_{B}$ are additional material parameters describing the bending stiffness of the material surface Eq. (9). In the presence of ${\mathbf{M}}_{S}$ in Eqs. (7)–(13), the tensor ${\mathbf{C}}_{{{\text{eq}}}}$ is also isotropic and its constants can be defined as

$$ K_{{{\text{eq}}}} = K_{1} + \hat{K}_{T} ;\;\mu_{{{\text{eq}}}} = \mu_{1} + \hat{\mu }_{T} + \hat{\mu }_{B} , $$

(51)

where $\hat{\mu }_{T}$ and $\hat{K}_{T}$ are defined in Eq. (26).

4 Effective properties

4.1 Basic features of MCM

Details of the evaluation of the effective stiffness tensor for composites with randomly distributed spherical particles are presented in Nazarenko et al. [38, 44]. Applying the methodology described in those papers, the following expression for effective moduli of composite with perfect interphase is obtained:

$$ {\mathbf{C}}^{ * } = {\overline{\mathbf{C}}} + {\mathbf{C}}_{1} {\mathbf{C}}_{3} {\mathbf{:}}{\kern 1pt} \left[ {\mathop {\mathbf{I}}\limits^{4} - {\mathbf{L}}{\mathbf{:C}}^{\prime } } \right]^{ - 1} {\mathbf{:}}\left[ {{\mathbf{C}}_{2} {\mathbf{L}}{\mathbf{:C}}_{3} } \right], $$

(52)

where ${\overline{\mathbf{C}}}$, ${\mathbf{C}}_{3}$ and ${\mathbf{C}}^{\prime }$ are

$$ {\overline{\mathbf{C}}} = \sum\limits_{k = 1}^{2} {c_{k} } {\mathbf{C}}_{k} ,{\mathbf{C}}_{3} = {\mathbf{C}}_{1} - {\mathbf{C}}_{2} ,{\mathbf{C}}^{\prime } = c_{1} {\mathbf{C}}_{2} + c_{2} {\mathbf{C}}_{1} - {\mathbf{C}}_{c} , $$

(53)

while ${\mathbf{C}}_{c}$ is the constitutive tensor for the “reference medium” whose selection is specified subsequently. It is shown Nazarenko et al. [45] that the tensor ${\mathbf{L}}$ in Eq. (52) coincides with the negative classical Hill tensor ${\mathbf{P}}$ (${\mathbf{P}} = {\mathbf{S}}:{\mathbf{C}}_{{2}}^{ - 1}$, with ${\mathbf{S}}$ being the Eshelby tensor; c.f. [48]. The difference is that the Hill tensor ${\mathbf{P}}$ is defined using properties of the matrix material whereas tensor ${\mathbf{L}}$ is related to the stiffness of the reference medium ${\mathbf{C}}_{c}$.

$$ {\mathbf{L}} = 2b\mathop {\mathbf{I}}\limits^{4} + a\mathop {\mathbf{I}}\limits^{2} \otimes \mathop {\mathbf{I}}\limits^{2} , $$

(54)

with

$$ a = \frac{{{\uplambda }_{c} + {\upmu }_{c} }}{{15{\upmu }_{c} [{\uplambda }_{c} + 2{\upmu }_{c} ]}},\;b = - \frac{{3{\uplambda }_{c} + 8{\upmu }_{c} }}{{30{\upmu }_{c} [{\uplambda }_{c} + 2{\upmu }_{c} ]}}, $$

(55)

and with ${\uplambda }_{c}$, ${\upmu }_{c}$ being the Lamé constants of the reference medium defined according to the rule Nazarenko et al. [37]

$$ {\uplambda }_{c} = \left\{ {\begin{array}{*{20}l} {c_{{1}} {\uplambda }_{{1}} + c_{{2}} {\uplambda }_{{2}} ,\quad \quad \quad \quad \quad \;\;{\text{if}}\;\;{\uplambda }_{{1}} \le {\uplambda }_{{2}} } \hfill \\ {\left[ {c_{{1}} \left( {{\uplambda }_{{1}} } \right)^{ - 1} + c_{{2}} \left( {{\uplambda }_{{2}} } \right)^{ - 1} } \right]^{ - 1} ,\;\;{\text{if}}\;\;{\uplambda }_{{1}} \ge {\uplambda }_{{2}} } \hfill \\ \end{array} } \right., $$

(56)

$$ {\upmu }_{c} = \left\{ {\begin{array}{*{20}l} {c_{{1}} {\upmu }_{{1}} + c_{{2}} {\upmu }_{{2}} ,\quad \quad \quad \quad \;\,\quad {\text{if}}\;\;{\upmu }_{{1}} \le {\upmu }_{{2}} } \hfill \\ {\left[ {c_{{1}} \left( {{\upmu }_{{1}} } \right)^{ - 1} + c_{{2}} \left( {{\upmu }_{{2}} } \right)^{ - 1} } \right]^{ - 1} ,\;\;{\text{if}}\;\;{\upmu }_{{1}} \ge {\upmu }_{{2}} } \hfill \\ \end{array} } \right. $$

(57)

4.2 Application of MCM to problems with interfaces

The effective properties of the composite with Steigmann–Ogden model of interfaces can be obtained by MCM from Eqs. (52)–(57), valid for perfect interfaces, if the properties of the spherical particles (${\mathbf{C}}_{1}^{\,}$) is replaced with those of equivalent inhomogeneity (${\mathbf{C}}_{{{\text{eq}}}}$) and, evaluated in Sect. 3.2.3. This leads to the following expression

$$ {\mathbf{C}}^{*} = {\mathbf{\tilde{\overline{C}}}} + c_{1} {\tilde{\mathbf{C}}}_{3} {\kern 1pt} {\mathbf{:}}{\kern 1pt} \left[ {\mathop {\mathbf{I}}\limits^{4} - {\tilde{\mathbf{L}}}{\mathbf{:\tilde{C}}}{^{\prime}}} \right]^{ - 1} {\mathbf{:}}\left[ {c_{2} {\tilde{\mathbf{L}}}{\mathbf{:\tilde{C}}}_{3} } \right], $$

(58)

where ${\mathbf{\tilde{\overline{C}}}}$, ${\tilde{\mathbf{C}}}_{3}$, ${\tilde{\mathbf{L}}}$, and ${\tilde{\mathbf{C}}}{^{\prime}}$ are determined in accordance with Eqs. (53)-(55) if ${\mathbf{C}}_{1}^{\,}$ is replaced by ${\mathbf{C}}_{{{\text{eq}}}}$.

The scalar expression for the effective bulk and shear moduli of the composite can be obtained from a tensorial formula (52), cf. [37, 38]:

$$ K^{*} = c_{1} {\text{K}}_{{{\text{eq}}}} + c_{2} K_{2} - \frac{{c_{1} c_{2} \left[ {{\text{K}}_{{{\text{eq}}}} - K_{2} } \right]^{\,2} }}{{c_{1} K_{2} + c_{2} {\text{K}}_{{{\text{eq}}}} + {\raise0.7ex\hbox{$4$} \!\mathord{\left/ {\vphantom {4 3}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$3$}}{\tilde{\mu }}_{c} }}, $$

(59)

$$ {\upmu }^{*} = c_{1} {\upmu }_{{{\text{eq}}}} + c_{{2}} {\upmu }_{{2}} + \frac{{{4}c_{1} c_{{2}} \tilde{b}\,\left[ {{\upmu }_{{{\text{eq}}}} - {\upmu }_{{2}} } \right]^{\,2} }}{{{1} - {4}\tilde{b}\,\left[ {c_{{2}} {\upmu }_{{{\text{eq}}}} + c_{1} {\upmu }_{{2}} - {\tilde{\mu }}_{c} } \right]^{{}} }}, $$

(60)

where

$$ \tilde{b} = - \frac{{3\left[ {\tilde{K}_{c} + 2{\tilde{\mu }}_{c} } \right]}}{{10{\tilde{\mu }}_{c} [3\tilde{K}_{c} + 4{\tilde{\mu }}_{c} ]}} $$

(61)

and $\tilde{K}_{c}$, $\tilde{\mu }_{c}$ are follows:

$$ \tilde{K}_{c} = \left\{ {\begin{array}{*{20}l} {c_{1} {\text{K}}_{{{\text{eq}}}} + c_{{2}} K_{{2}} ,\quad \quad \quad \quad \quad \;\;{\text{if}}\;\;{\text{K}}_{{{\text{eq}}}} \le K_{{2}} } \hfill \\ {\left[ {c_{1} \left( {{\text{K}}_{{{\text{eq}}}} } \right)^{ - 1} + c_{{2}} \left( {K_{{2}} } \right)^{ - 1} } \right]^{ - 1} ,\;\;{\text{if}}\;\;{\text{K}}_{{{\text{eq}}}} \ge K_{{2}} } \hfill \\ \end{array} } \right., $$

(62)

$$ {\tilde{\mu }}_{c} = \left\{ {\begin{array}{*{20}l} {c_{1} {\upmu }_{{{\text{eq}}}} + c_{{2}} {\upmu }_{{2}} ,\quad \quad \quad \quad \quad \;\;{\text{if}}\;\;{\upmu }_{{{\text{eq}}}} \le {\upmu }_{{2}} } \hfill \\ {\left[ {c_{1} \left( {{\upmu }_{{{\text{eq}}}} } \right)^{ - 1} + c_{{2}} \left( {{\upmu }_{{2}} } \right)^{ - 1} } \right]^{ - 1} ,\;\;{\text{if}}\;\;{\upmu }_{{{\text{eq}}}} \ge {\upmu }_{{2}} } \hfill \\ \end{array} } \right. $$

(63)

In the above formulas $K_{{2}}$ and ${\upmu }_{{2}}$ are the bulk and shear moduli of the matrix, whereas ${\text{K}}_{{{\text{eq}}}}$ and ${\upmu }_{{{\text{eq}}}}$ are determined in Eqs. (50), (51), (26).

5 Numerical comparisons and discussion

As a numerical example, we consider a composite material consisting of silver matrix with the properties $E_{2} = 50\;{\text{GPa}}$ and $\nu_{2} = 0.37$ containing spherical cavities ($\mu_{1} = K_{1} = 0$). The free surface properties are those presented by Mohammadi and Sharma [22]. In their article surface properties are determined for surface [100]: $\tau_{0} = 0.3701\;{{\text{N}} \mathord{\left/ {\vphantom {{\text{N}} {\text{m}}}} \right. \kern-\nulldelimiterspace} {\text{m}}}$; $\mu_{S} = - 2.6948\;{{\text{N}} \mathord{\left/ {\vphantom {{\text{N}} {\text{m}}}} \right. \kern-\nulldelimiterspace} {\text{m}}}$; $\mu_{B} = 12.3 \cdot 10^{ - 19} {\text{Nm}}$. In the article, results are presented only for shear moduli. Assuming that relationship between Lamé parameters $\lambda$ and $\mu$ for bulk material and Lamé parameters for surface tension and for surface bending are the same we added missing parameters using relationship between Lamé parameters for bulk silver ${\raise0.7ex\hbox{$\lambda $} \!\mathord{\left/ {\vphantom {\lambda \mu }}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{$\mu $}} = 2.85$. So, in our calculation we take hypothetical values of $\lambda_{S} = - 7.69\;{{\text{N}} \mathord{\left/ {\vphantom {{\text{N}} {\text{m}}}} \right. \kern-\nulldelimiterspace} {\text{m}}}$; $\lambda_{B} = 35 \cdot 10^{ - 19} {\text{N}} \cdot {\text{m}}$.

The variation of the normalized shear modulus $\mu^{*} /\mu_{2}$ with the void volume fractions calculated by the MCM in combination with EEI approach for the spherical cavities of radius $a = 5$ nm with Steigmann-Ogden model of interface is shown in Fig. 1 (black solid line). For comparison, the normalized shear modulus for the same material obtained accounting only for parameters of Gurtin–Murdoch interface model ($\mu_{B} = 0$, $\lambda_{B} = 0$) (dash dot line) and for classical case without surface energy (short dash line) are also shown in Fig. 1. It is seen that difference between the normalized shear moduli obtained with accounting for Steigmann-Ogden and for Gurtin–Murdoch models of interface (when bending surface parameters are neglected) differ insignificantly what is illustrated in Fig. 2, which is a zoomed section of Fig. 1. It means that for case of interface whose thickness is vanishingly small, the Gurtin–Murdoch model of interface is dominant and influence of the surface bending on effective properties is negligibly small. At the same time, the normalized shear moduli obtained with accounting only for Gurtin–Murdoch and for classical case without surface energy differ essentially, especially for high pore volume fraction. As expected, this numerical illustration indicates that the influence of the additional surface bending parameters of the Steigmann-Ogden model is insignificant for the shear moduli. It can be explained by the fact, that as shown by Benveniste and Miloh [1] the problem of an inhomogeneity with a thin interphase layer can be reduced to that of an imperfectly bonded inhomogeneity with the interface conditions described by one of seven distinct regimes, which depend on the elastic properties of the interphase layer and the components of composite. The Gurtin–Murdoch model corresponds to the membrane type interface regime, while the Steigmann-Ogden model corresponds to the inextensible shell regime. The surface bending stiffness of the membrane type interface is negligible, and the results presented in the manuscript support this. However, there may be a special case where the bending stiffness of the surface should be taken into account (e.g. [46].

Fig. 1

Dependence of shear modulus $\mu^{*} /\mu_{2}$ for nanoporous silver on void volume fraction $c_{1}$; radius of a spherical cavity $a = 5$(nm)

Fig. 2

Dependence of shear modulus $\mu^{*} /\mu_{2}$ for nanoporous silver on void volume fraction $c_{1}$; radius of a spherical cavity $a = 5$ (nm) (zoomed section of Fig. 1)

As a second numerical result, we consider porous material with spherical cavities of radius $a = 5$ nm and following parameters adopted after Zemlyanova and Mogilevskaya [35]: $v_{2} = 0.3$, $\tau_{0} = 0\;{{\text{N}} \mathord{\left/ {\vphantom {{\text{N}} {\text{m}}}} \right. \kern-\nulldelimiterspace} {\text{m}}}$; $\frac{{\mu_{S} }}{{\mu_{2} a}} = {0}{\text{.030156}}$; $\frac{{\lambda_{S} }}{{\mu_{2} a}} = {0}{\text{.060312}}$; $\frac{{5\mu_{B} + 3\lambda_{B} }}{{\mu_{2} a^{3} }} = {0}{\text{.00028382}}$.

The variation of the normalized shear modulus $\mu^{*} /\mu_{2}$ with void volume fractions calculated by the MCM in combination with the EEI approach (MCM EEI black solid line) for Steigmann-Ogden model of interface and for classical case without surface energy (MCM dash dot line) is shown in Fig. 3. For comparison analogical results for the normalized effective shear modulus with Steigmann-Ogden model of interface and for the classical case obtained in Zemlyanova and Mogilevskaya [32, 35] (short dash dot line and short dash line) are presented in Fig. 3 as well. The results of Zemlyanova and Mogilevskaya [35] were obtained: for equivalent inhomogeneity on the base of exact Lurie’s solution [40] for spherical particle and for effective properties of entire composite on the base of Maxwell homogenization scheme [49]. As it is seen, tendencies and numerical values for the results calculated by different methods demonstrate that the effective shear modulus of considered porous material presented in Zemlyanova and Mogilevskaya [35] are higher than those determined by EEI approach in combination with MCM. It can be explained by the conceptual differences in the MCM and Maxwell scheme, which in the case of material with spherical cavities is identical to Mori–Tanaka method (MTM) [48]. It is known, that for porous material MTM shows upper bound for effective moduli. For confirmation of this explanation the normalized effective shear modulus with Steigmann-Ogden model of interface and for the classical case were calculated by EEI approach in combination with MTM (MTM EEI short dash dot line and MTM dot line). These results and those of Zemlyanova and Mogilevskaya [35] are shown in Fig. 4. It is seen that the shear moduli calculated by MTM in combination with EEI approach are identical with the both obtained on the base of Maxwell scheme in combination with exact Lurie’s solution. It means, that the contribution of Steigmann-Ogden interface to effective properties of entire composite is identical for the both cases. It can be considered as additional verification of accuracy of EEI approach where properties of equivalent inhomogeneity are identical with those obtained on the base of Lurie’s solution. This is a very positive sign for the present approach given that in the definition of equivalent inhomogeneity based on Lurie’s solution all governing equations in the inhomogeneity/interphase system are satisfied exactly. At the same time, the effective properties of entire composite calculated by MCM are lower than upper bound calculated by MTM.

Fig. 3

Dependence of shear modulus $\mu^{*} /\mu_{2}$ for nanoporous material on void volume fraction $c_{1}$ (MCM EI and Zemlyanova et al. [35]); radius of a spherical cavity $a = 5$ (nm)

Fig. 4

Dependence of shear modulus $\mu^{*} /\mu_{2}$ for nanoporous material on void volume fraction $c_{1}$$c_{1}$ (MTM EI and Zemlyanova et al. [35]); radius of a spherical cavity $a = 5$ (nm)

Remark

It should be noted that properties of equivalent spherical inhomogeneity with spring layer interphase obtained by EEI approach [38] have been compared with those calculated using the equivalent inhomogeneity defined on the base of Lurie’s solution [38]. It was observed, that both solutions are virtually identical for all volume fraction of inhomogeneities.

6 Conclusions

A mathematical model, employing the concept of the EEI in combination with the MCM [37, 39, 50], has been generalized to introduce the surface effects described by the Steigmann-Ogden model [19, 20] derived within the strain gradient elasticity [31]. A particular focus was centered on accounting for the surface bending contribution in the definition of the EEI. The model was also used in combination with MCM to determine the effective properties of materials with randomly distributed nano-particles with Steigmann-Ogden model of interface. Effective shear modulus of spherical EEI with Steigmann-Ogden model of interface is presented in a closed-form and compared with those obtained on the base of Lurie’s solution for sphere.

The properties of the EEI are determined based on the derived definition of surface energy, which includes the surface tension and the surface bending, assuming the uniform state of strains within the inhomogeneity. This assumption is particularly suitable in the context of the MCM used here. The restriction of the MCM to the two-point approximation adopted here is tantamount to the assumption that deformation of the inhomogeneities is uniform. This is in perfect agreement with the way the EEI is defined in Sect. 3.

As a numerical illustration of the presented approach, nanoporous silver is studied. It is shown that the equivalent bulk modulus of spherical equivalent inhomogeneity does not depend on the surface bending parameters.

Shear moduli of nanoporous silver has been analyzed for varying volume content of the nano-cavities for Steigmann-Ogden model of interface and for classical case without accounting for the surface effects. The size effect introduced due to contribution of the residual stresses, elasticity and bending of the matrix/nanoparticles interface in nanoporous silver is accounted for in the expressions for effective bulk and shear moduli of the composite. It has been shown that the contribution of the surface bending to the shear moduli is insignificant.

As a second numerical example, the effective shear moduli of the porous material, evaluated on the basis of MCM in combination with EEI approach was considered in context of comparison with numerical results available in the literature and, is in a good agreement with one calculated in Zemlyanova and Mogilevskaya [32, 35]. It is interesting to note that, in spite of the difference between expressions for the shear modulus of equivalent inhomogeneity based on EEI and that obtained on the base of Lurie’s solution for sphere, the effective properties of the composite calculated using those two definitions of the equivalent inhomogeneity are virtually identical for all volume fraction of voids. This indicates that the basic assumptions underlying the proposed approach are quite sound, as the approach based on Lurie’s solution exactly fulfills all governing equations within the original inhomogeneity and the interphase.

To conclude it is worth mentioning that the definition of the EEI is general and can be used in the case of inhomogeneities of other shapes than spherical, e.g., ellipsoidal or cylindrical (e.g. [36]. It can be very naturally combined with the MCM and appears to be potentially amenable for inclusion of other than Gurtin–Murdoch or Steigmann-Ogden interface models. Therefore, the MCM with combination of the EEI can be applied for analysis of materials containing inhomogeneities with more complex geometric and mechanical characteristics. The important characteristic of the proposed approach is its ability to provide closed-form expressions for the effective properties of nano-composites. Closed-form results are important, especially if the influence of different problem parameters needs to be analyzed.

Appendices

Appendix 1: Components of curvature tensors for sphere

Let’s assume that the surface of interest is locally parameterized by $\xi^{\Lambda }$, $\Lambda \in \left\{ {\,1\,,\;2} \right\}$, that is the position vector of a point on that surface is expressed as ${\mathbf{r}}\left( {\xi^{\Lambda } } \right)$. Then, one can define a couple of vectors ${\mathbf{G}}_{\Lambda }$

$$ {\mathbf{G}}_{\Lambda } \; = \;\frac{{\partial {\mathbf{r}}}}{{\partial \xi^{\Lambda } }}\; \equiv \;{\mathbf{r}},_{\Lambda } , $$

(64)

which forms the vector basis in the linear space tangent to the surface S, called the natural basis. Another basis in the same tangent space, denoted by ${\mathbf{G}}^{{\Delta }}$ and called dual or reciprocal, is defined via the following orthogonality condition

$$ {\mathbf{G}}_{{\Lambda }} \cdot {\mathbf{G}}^{{\Delta }} \; = \;\delta^{\Delta }_{\Lambda } , $$

(65)

where the “·” represents the “dot” (or “inner”) product of vectors and $\delta^{\Delta }_{\Lambda }$ is the Kronecker “delta”. The bases ${\mathbf{G}}_{\Lambda }$ and ${\mathbf{G}}^{{\Delta }}$ are functions of $\xi^{\Lambda }$ and their derivatives can be expressed by the well-known Gauss-Weingarten formulas (see Itskov [47], for example). For the natural basis these formulas are (cf. Equation (64) for notations)

$$ {\mathbf{G}}_{\Lambda } ,_{\Sigma } \; = \;\Gamma_{\Lambda \Sigma }^{\Delta } {\mathbf{G}}_{\Delta } + B_{\Lambda \Sigma } {\mathbf{n}}\; \equiv \;\;\Gamma_{\Lambda \Sigma }^{1} {\mathbf{G}}_{1} \; + \;\Gamma_{\Lambda \Sigma }^{2} {\mathbf{G}}_{2} + B_{\Lambda \Sigma } {\mathbf{n}}, $$

(66)

with a unit vector n normal to the surface. Here (as shown in the above equation) an index repeated in the subscript and superscript position implies summation, $\Gamma_{\Lambda \Sigma }^{\Delta } = {\mathbf{G}}_{\Lambda } ,_{\Sigma } \cdot {\mathbf{G}}^{\Delta }$ are the so-called Christoffel symbols (of the second kind) and the components of the local curvature tensor are

$$ B_{\Lambda \Sigma } = {\mathbf{G}}_{\Lambda } ,_{\Sigma } \cdot {\mathbf{n}}{\mathbf{.}} $$

(67)

Equation (66) together with Eq. (64) imply that $B_{\Lambda \Sigma } = B_{\Sigma \Lambda }$ whereas definition of the Christoffel symbols and Eq. (64) imply the following symmetry property $\Gamma_{\Lambda \Sigma }^{\Omega } = \Gamma_{\Sigma \Lambda }^{\Omega }$. The analogical formulas for the derivatives of vectors of the dual basis are

$$ {\mathbf{G}}^{\Delta } ,_{\Sigma } \; = \; - \Gamma_{\Lambda \Sigma }^{\Delta } {\mathbf{G}}^{\Lambda } - B^{\Delta }_{\Sigma } {\mathbf{n}}\; \equiv \;\; - \Gamma_{1\Sigma }^{\Delta } {\mathbf{G}}^{1} \; - \;\Gamma_{2\Sigma }^{\Delta } {\mathbf{G}}^{2} - B^{\Delta }_{\Sigma } {\kern 1pt} {\mathbf{n}}, $$

(68)

where $B^{\Delta }_{\Sigma }$ are the so-called mixed components of the local curvature tensor.

The curvature tensor B (components of which appear in Eqs. (66) and (67)) can be represented in several ways shown below (as well as many other ways)

$$ {\mathbf{B}} = B_{\Delta \Lambda } {\mathbf{G}}^{\Delta } \otimes {\mathbf{G}}^{\Lambda } = B^{\Delta \Lambda } {\mathbf{G}}_{\Delta } \otimes {\mathbf{G}}_{\Lambda } = B_{\Delta }^{\Lambda } {\mathbf{G}}^{\Delta } \otimes {\mathbf{G}}_{\Lambda } = B^{\Delta }_{\Lambda } {\mathbf{G}}_{\Delta } \otimes {\mathbf{G}}^{\Lambda } . $$

(69)

In the above equation double summation is implied and the (indexed) coefficients multiplying the dyadic are various components of the tensor ${\mathbf{B}}$. They all can be different, but they are related to each other by transformation formulas involving the so-called Gram matrices related to the natural or dual bases. Those matrices are defined as follows

$$ G_{\Delta \Lambda } = {\mathbf{G}}_{\Delta } \cdot {\mathbf{G}}_{\Lambda } ,\;G^{\Delta \Lambda } = {\mathbf{G}}^{\Delta } \cdot {\mathbf{G}}^{\Lambda } . $$

(70)

The relationship between various components of the curvature tensor B one can present the following

$$ B_{\Delta }^{\Lambda } = B_{\Delta \Sigma } G^{\Sigma \Lambda } . $$

(71)

In the case of spherical inhomogeneity of radius $a$, the position vector ${\mathbf{r}}$ of a point on the surface of the inhomogeneity, and the corresponding unit vector n, normal to that surface, may be expressed as follows

$$ {\mathbf{r}} = a\left[ {\begin{array}{*{20}c} {\cos \phi \sin \theta } \\ {\sin \phi \sin \theta } \\ { - \cos \theta } \\ \end{array} } \right],\;{\mathbf{n}} = \left[ {\begin{array}{*{20}c} {{\text{cos}}\,\phi {\kern 1pt} \,\sin \,\theta } \\ {\sin \,\phi \;\sin \,\theta } \\ { - {\text{cos}}\,\theta } \\ \end{array} } \right],\;0 \le \varphi \le 2\pi ,\;0 \le \theta \le \pi . $$

(72)

The natural basis is

$$ {\mathbf{G}}_{1} = r_{,\phi } = a\,\left[ {\begin{array}{*{20}c} { - {\text{sin}}\,\theta \,{\text{sin}}\,\phi {\kern 1pt} \,} \\ {{\text{sin}}\,\theta \,{\text{cos}}\,\phi \;} \\ {0} \\ \end{array} } \right],\;{\mathbf{G}}_{2} = r_{,\theta } = a\,\left[ {\begin{array}{*{20}c} {\cos \,\phi {\kern 1pt} \,\cos \,\theta } \\ {\sin \,\phi \;\cos \,\theta } \\ {\sin \,\theta } \\ \end{array} } \right],\;{\mathbf{G}}_{3} = {\mathbf{G}}^{{3}} = {\mathbf{n}} = \left[ {\begin{array}{*{20}c} {{\text{cos}}\,\phi {\kern 1pt} \,\sin \,\theta } \\ {\sin \,\phi \;\sin \,\theta } \\ { - {\text{cos}}\,\theta } \\ \end{array} } \right], $$

(73)

$$ {\mathbf{G}}_{1,1} = r_{,\phi \phi } = a\,\left[ {\begin{array}{*{20}c} { - {\text{sin}}\,\theta \,\cos \,\phi {\kern 1pt} \,} \\ {{\text{ - sin}}\,\theta \,\sin \;} \\ {0} \\ \end{array} } \right],\;{\mathbf{G}}_{2,1} = r_{,\theta \phi } = a\,\left[ {\begin{array}{*{20}c} { - \sin \,\phi {\kern 1pt} \,\cos \,\theta } \\ {\cos \,\phi \;\cos \,\theta } \\ {\sin \,\theta } \\ \end{array} } \right],\;{\mathbf{G}}_{2,2} = r_{,\theta \theta } = a\,\left[ {\begin{array}{*{20}c} { - \cos \,\phi {\kern 1pt} \,\sin \,\theta } \\ { - \sin \,\phi \;\sin \,\theta } \\ {cos\theta } \\ \end{array} } \right] $$

(74)

The dual basis is defined as

$$ {\mathbf{G}}^{{1}} = \frac{1}{a}\,\left[ {\begin{array}{*{20}c} { - {{{\text{sin}}\,\phi } \mathord{\left/ {\vphantom {{{\text{sin}}\,\phi } {{\text{sin}}\,\theta }}} \right. \kern-\nulldelimiterspace} {{\text{sin}}\,\theta }}{\kern 1pt} \,} \\ {{{{\text{cos}}\,\phi } \mathord{\left/ {\vphantom {{{\text{cos}}\,\phi } {{\text{sin}}\,\theta }}} \right. \kern-\nulldelimiterspace} {{\text{sin}}\,\theta }}\;} \\ {0} \\ \end{array} } \right],\;{\mathbf{G}}^{{2}} = \frac{1}{a}\,\left[ {\begin{array}{*{20}c} {\cos \,\phi {\kern 1pt} \,\cos \,\theta } \\ {\sin \,\phi \;\cos \,\theta } \\ {\sin \,\theta } \\ \end{array} } \right],\;{\mathbf{G}}_{3} = {\mathbf{G}}^{{3}} = {\mathbf{n}} = \left[ {\begin{array}{*{20}c} {{\text{cos}}\,\phi {\kern 1pt} \,\sin \,\theta } \\ {\sin \,\phi \;\sin \,\theta } \\ { - {\text{cos}}\,\theta } \\ \end{array} } \right]. $$

(75)

Then the curvature tensors for spherical inhomogeneity of radius $a$ are given:

See Table 1

Appendix 2: Properties of the energy-equivalent inhomogeneity accounting for surface bending

For illustration of some technical details, ${\mathbf{K}}_{{\mu_{B} }}$ and ${\mathbf{K}}_{{\lambda_{B} }}$ of Eqs. (45)–(49) are evaluated in this Appendix. In addition to ${\mathbf{K}}_{{\mu_{B} }}$ and ${\mathbf{K}}_{{\lambda_{B} }}$, the contribution of surface tension to properties of equivalent inhomogeneity includes other term ${\mathbf{K}}_{T}$ present in Eq. (49), however evaluation of this term is presented in Nazarenko et al. [39].

Assuming that inhomogeneities are spheres of radius $a$ and using spherical coordinate system Eq. (72), ${\mathbf{K}}_{{\mu_{B} }}$ of Eqs. (45), (47), (48) is described by

$$ \begin{aligned} K_{{\mu_{B} }} = & \frac{{3\mu_{B} }}{{2\pi a^{3} }}\int\limits_{0}^{\pi } {\int\limits_{0}^{2\pi } {\left[ {{\overline{\mathbf{G}}}_{1} \otimes } \right.} } \,{\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} + {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} + 2{\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{2} - {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\mathbf{n}} \otimes {\mathbf{n}} \\ & - {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\mathbf{n}} \otimes {\mathbf{n}} - {\mathbf{n}} \otimes {\mathbf{n}} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} - {\mathbf{n}} \otimes {\mathbf{n}} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} + 2{\mathbf{n}} \otimes {\mathbf{n}} \otimes {\mathbf{n}} \otimes \left. {\mathbf{n}} \right]\sin \theta \,d\phi \,d\theta ; \\ \end{aligned} $$

(76)

$$ \begin{aligned} K_{{\mu_{b} \left[ {1111} \right]}} = & \frac{{3\mu_{B} }}{{2\pi r^{3} }}\int\limits_{0}^{\pi } {\int\limits_{0}^{2\pi } {\left[ {{\text{sin}}\,\phi^{4} } \right. + } } \cos \,\phi^{4} \cos \,\theta^{4} + 2{\text{sin}}\,\phi^{2} \cos \,\phi^{2} \cos \,\theta^{2} - 2{\text{sin}}\,\phi^{2} \cos \,\phi^{2} \sin \,\theta^{2} \\ & - 2\cos \,\phi^{4} \cos \,\theta^{2} \sin \,\theta^{2} + 2\cos \,\phi^{4} \sin \,\left. {\theta^{4} } \right]\sin \theta \,d\phi \,d\theta = 4\frac{{\mu_{B} }}{{a^{3} }}; \\ \end{aligned} $$

(77)

$$ K_{{\mu_{B} \left[ {2222} \right]}} = K_{{\mu_{B} \left[ {3333} \right]}} = 4\frac{{\mu_{B} }}{{\pi a^{3} }}; $$

(78)

$$ \begin{aligned} K_{{\mu_{B} \left[ {1122} \right]}} & = \frac{{3\mu_{B} }}{{2\pi a^{3} }}\int\limits_{0}^{\pi } {\int\limits_{0}^{2\pi } {\left[ {{\text{sin}}\,\phi^{2} \cos \,\phi^{2} } \right. + } } {\text{sin}}\,\phi^{2} \cos \,\phi^{2} \cos \,\theta^{4} - 2{\text{sin}}\,\phi^{2} \cos \,\phi^{2} \cos \,\theta^{2} - {\text{sin}}\,\phi^{4} \sin \,\theta^{2} \\ & - \cos \,\phi^{4} \sin \,\theta^{2} - 2\sin \,\phi^{2} \cos \,\phi^{2} \cos \,\theta^{2} \sin \,\theta^{2} + 2\sin \,\phi^{2} \cos \,\phi^{2} \sin \,\left. {\theta^{4} } \right]\sin \theta \,d\phi \,d\theta = - 2\frac{{\mu_{B} }}{{a^{3} }} \\ \end{aligned} $$

(79)

$$ K_{{\mu_{B} \left[ {1133} \right]}} = K_{{\mu_{B} \left[ {2233} \right]}} = - 2\frac{{\mu_{B} }}{{a^{3} }}; $$

(80)

$$ \begin{aligned} K_{{\mu_{B} \left[ {1212} \right]}} = & \frac{{3\mu_{B} }}{{2\pi a^{3} }}\int\limits_{0}^{\pi } {\int\limits_{0}^{2\pi } {\left[ {{\text{sin}}\,\phi^{2} \cos \,\phi^{2} } \right. + } } {\text{sin}}\,\phi^{2} \cos \,\phi^{2} \cos \,\theta^{4} + 2{\text{sin}}\,\phi^{4} \cos \,\theta^{2} + 2{\text{sin}}\,\phi^{2} \cos \,\phi^{2} \sin \,\theta^{2} \\ & - 2\sin \,\phi^{2} \cos \,\phi^{2} \cos \,\theta^{2} \sin \,\theta^{2} + 2\sin \,\phi^{2} \cos \,\phi^{2} \sin \,\left. {\theta^{4} } \right]\sin \theta \,d\phi \,d\theta = 4\frac{{\mu_{B} }}{{a^{3} }} \\ \end{aligned} $$

(81)

$$ K_{{\mu_{B} \left[ {1313} \right]}} = K_{{\mu_{B} \left[ {2323} \right]}} = 4\frac{{\mu_{B} }}{{a^{3} }}; $$

(82)

$$ \begin{aligned} K_{{\mu_{B} \left[ {2112} \right]}} = & \frac{{3\mu_{B} }}{{2\pi a^{3} }}\int\limits_{0}^{\pi } {\int\limits_{0}^{2\pi } {\left[ {{\text{sin}}\,\phi^{2} \cos \,\phi^{2} } \right. + } } {\text{sin}}\,\phi^{2} \cos \,\phi^{2} \cos \,\theta^{4} - 2{\text{sin}}\,\phi^{2} \cos \,\phi^{2} \cos \,\theta^{2} + 2{\text{sin}}\,\phi^{2} \cos \,\phi^{2} \sin \,\theta^{2} \\ & - 2\sin \,\phi^{2} \cos \,\phi^{2} \cos \,\theta^{2} \sin \,\theta^{2} + 2\sin \,\phi^{2} \cos \,\phi^{2} \sin \,\left. {\theta^{4} } \right]\sin \theta \,d\phi \,d\theta = 2\frac{{\mu_{B} }}{{a^{3} }} \\ \end{aligned} $$

(83)

$$ K_{{\mu_{B} \left[ {3113} \right]}} = K_{{\mu_{B} \left[ {3223} \right]}} = 2\frac{{\mu_{B} }}{{a^{3} }}. $$

(84)

Consequently, the contribution of ${\mathbf{K}}_{{\mu_{B} }}$ to bulk and shear moduli of the equivalent inhomogeneity are

$$ \hat{K}_{{\mu_{B} }} = K_{{\mu_{B} \left[ {1111} \right]}} + K_{{\mu_{B} \left[ {1122} \right]}} + K_{{\mu_{B} \left[ {1133} \right]}} = 4\frac{{\mu_{B} }}{{a^{3} }} - 2\frac{{\mu_{B} }}{{a^{3} }} - 2\frac{{\mu_{B} }}{{a^{3} }} = 0; $$

(85)

$$ \hat{\mu }_{B} = \frac{1}{2}\left[ {K_{{\mu_{B} \left[ {1212} \right]}} + K_{{\mu_{B} \left[ {2112} \right]}} } \right] = 3\frac{{\mu_{B} }}{{a^{3} }}. $$

(86)

${\mathbf{K}}_{{\lambda_{B} }}$ of Eqs. (46), (47), (48) is described by

$$ \begin{aligned} K_{{\lambda_{B} }} = & \frac{{3\lambda_{B} }}{{4\pi a^{3} }}\int\limits_{0}^{\pi } {\int\limits_{0}^{2\pi } {\left[ {{\overline{\mathbf{G}}}_{1} \otimes } \right.} } \,{\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} + {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} + {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \\ & + {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} - 2{\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\mathbf{n}} \otimes {\mathbf{n}} - 2{\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\mathbf{n}} \otimes {\mathbf{n}} \\ & - 2{\mathbf{n}} \otimes {\mathbf{n}} \otimes {\overline{\mathbf{G}}}_{1} \otimes {\overline{\mathbf{G}}}_{1} - 2{\mathbf{n}} \otimes {\mathbf{n}} \otimes {\overline{\mathbf{G}}}_{2} \otimes {\overline{\mathbf{G}}}_{2} + 4{\mathbf{n}} \otimes {\mathbf{n}} \otimes {\mathbf{n}} \otimes \left. {\mathbf{n}} \right]\sin \theta \,d\phi \,d\theta \\ \end{aligned} $$

(87)

$$ \begin{aligned} K_{{\lambda_{B} \left[ {1111} \right]}} = & \frac{{3\lambda_{B} }}{{4\pi a^{3} }}\int\limits_{0}^{\pi } {\int\limits_{0}^{2\pi } {\left[ {{\text{sin}}\,\phi^{4} } \right. + } } \cos \,\phi^{4} \cos \,\theta^{4} + 2{\text{sin}}\,\phi^{2} \cos \,\phi^{2} \cos \,\theta^{2} - 4{\text{sin}}\,\phi^{2} \cos \,\phi^{2} \sin \,\theta^{2} \\ & - 4\cos \,\phi^{4} \cos \,\theta^{2} \sin \,\theta^{2} + 4\cos \,\phi^{4} \sin \,\left. {\theta^{4} } \right]\sin \theta \,d\phi \,d\theta = \frac{12}{5}\frac{{\lambda_{B} }}{{a^{3} }}; \\ \end{aligned} $$

(88)

$$ K_{{\lambda_{B} \left[ {2222} \right]}} = K_{{\lambda_{B} \left[ {3333} \right]}} = \frac{12}{5}\frac{{\lambda_{B} }}{{a^{3} }}; $$

(89)

$$ \begin{aligned} K_{{\lambda_{B} \left[ {1122} \right]}} & = \frac{{3\lambda_{B} }}{{4\pi a^{3} }}\int\limits_{0}^{\pi } {\int\limits_{0}^{2\pi } {\left[ {{\text{sin}}\,\phi^{2} \cos \,\phi^{2} } \right. + } } {\text{sin}}\,\phi^{2} \cos \,\phi^{2} \cos \,\theta^{4} + {\text{sin}}\,\phi^{4} \cos \,\theta^{2} \\ & + \cos \,\phi^{4} \cos \,\theta^{2} - 2{\text{sin}}\,\phi^{4} \sin \,\theta^{2} - 2\cos \,\phi^{4} \sin \,\theta^{2} \\ & - 4\cos \,\phi^{2} \sin \,\phi^{2} \cos \,\theta^{2} \sin \,\theta^{2} + 4\cos \,\phi^{2} \sin \,\phi^{2} \sin \,\left. {\theta^{4} } \right]\sin \theta \,d\phi \,d\theta = - \frac{6}{5}\frac{{\lambda_{B} }}{{a^{3} }}; \\ \end{aligned} $$

(90)

$$ K_{{\lambda_{B} \left[ {1133} \right]}} = K_{{\lambda_{B} \left[ {2233} \right]}} = - \frac{6}{5}\frac{{\lambda_{B} }}{{a^{3} }}; $$

(91)

$$ \begin{aligned} K_{{\lambda_{b} \left[ {1212} \right]}} & = \frac{{3\lambda_{b} }}{{4\pi r^{3} }}\int\limits_{0}^{\pi } {\int\limits_{0}^{2\pi } {\left[ {{\text{sin}}\,\phi^{2} \cos \,\phi^{2} } \right. + } } {\text{sin}}\,\phi^{2} \cos \,\phi^{2} \cos \,\theta^{4} \\ & - 2{\text{sin}}\,\phi^{2} \cos \,\phi^{2} \cos \,\theta^{2} + 4{\text{sin}}\,\phi^{2} \cos \,\phi^{2} \sin \,\theta^{2} \\ & - 4\cos \,\phi^{2} \sin \,\phi^{2} \cos \,\theta^{2} \sin \,\theta^{2} + 4\cos \,\phi^{2} \sin \,\phi^{2} \sin \,\left. {\theta^{4} } \right]\sin \theta \,d\phi \,d\theta = \frac{9}{5}\frac{{\lambda_{B} }}{{a^{3} }}; \\ \end{aligned} $$

(92)

$$ K_{{\lambda_{B} \left[ {1313} \right]}} = K_{{\lambda_{B} \left[ {2323} \right]}} = \frac{9}{5}\frac{{\lambda_{B} }}{{a^{3} }}; $$

(93)

$$ \begin{aligned} K_{{\lambda_{B} \left[ {2112} \right]}} & = \frac{{3\lambda_{B} }}{{4\pi a^{3} }}\int\limits_{0}^{\pi } {\int\limits_{0}^{2\pi } {\left[ {{\text{sin}}\,\phi^{2} \cos \,\phi^{2} } \right. + } } {\text{sin}}\,\phi^{2} \cos \,\phi^{2} \cos \,\theta^{4} \\ & - 2{\text{sin}}\,\phi^{2} \cos \,\phi^{2} \cos \,\theta^{2} + 4{\text{sin}}\,\phi^{2} \cos \,\phi^{2} \sin \,\theta^{2} \\ & - 4\cos \,\phi^{2} \sin \,\phi^{2} \cos \,\theta^{2} \sin \,\theta^{2} + 4\cos \,\phi^{2} \sin \,\phi^{2} \sin \,\left. {\theta^{4} } \right]\sin \theta \,d\phi \,d\theta = \frac{9}{5}\frac{{\lambda_{B} }}{{a^{3} }}; \\ \end{aligned} $$

(94)

$$ K_{{\lambda_{B} \left[ {3113} \right]}} = K_{{\lambda_{B} \left[ {3223} \right]}} = \frac{9}{5}\frac{{\lambda_{B} }}{{a^{3} }}. $$

(95)

The contribution of ${\mathbf{K}}_{{\lambda_{B} }}$ to bulk and shear moduli of the equivalent inhomogeneity are

$$ \hat{K}_{{\lambda_{B} }} = K_{{\lambda_{B} \left[ {1111} \right]}} + K_{{\lambda_{B} \left[ {1122} \right]}} + K_{{\lambda_{B} \left[ {1133} \right]}} = \frac{12}{5}\frac{{\lambda_{B} }}{{a^{3} }} - \frac{6}{5}\frac{{\lambda_{B} }}{{a^{3} }} - \frac{6}{5}\frac{{\lambda_{B} }}{{a^{3} }} = 0; $$

(96)

$$ \hat{\mu }_{B} = \frac{1}{2}\left[ {K_{{\lambda_{B} \left[ {1212} \right]}} + K_{{\lambda_{B} \left[ {2112} \right]}} } \right] = \frac{9}{5}\frac{{\lambda_{B} }}{{a^{3} }}. $$

(97)

The contribution of surface bending to bulk and shear moduli of the equivalent inhomogeneity Eq. (28) is

$$ \hat{K}_{B} = \hat{K}_{B} + \hat{K}_{B} = 0; $$

(98)

$$ \hat{\mu }_{B} = \hat{\mu }_{B} + \hat{\mu }_{B} = 3\frac{{\mu_{B} }}{{a^{3} }} + \frac{{9\lambda_{B} }}{{5a^{3} }} = \frac{{15\mu_{B} + 9\lambda_{B} }}{{5a^{3} }}. $$

(99)