Nonlinear elastic moduli of composite materials with nonlinear spherical inclusions dispersed in a nonlinear matrix

TL;DR

This paper develops a theoretical method to evaluate the nonlinear elastic properties of composite materials with nonlinear spherical inclusions in a nonlinear matrix, providing analytical formulas and numerical verification.

Contribution

It introduces an analytical approach to determine both linear and nonlinear elastic moduli of composites with nonlinear inclusions, including explicit formulas for effective Murnaghan moduli.

Findings

Effective Murnaghan moduli depend linearly on constituent moduli.

Analytical formulas for composite moduli are derived.

Numerical modeling confirms theoretical results.

Abstract

A theory is developed for evaluation of nonlinear elastic moduli of composite materials with nonlinear inclusions dispersed in another nonlinear material (matrix). We elaborate a method aimed for determination of elastic parameters of a composite: its linear elastic moduli (second-order elastic constants) and nonlinear elastic moduli, which are known as the Murnaghan moduli (third-order elastic constants). We find an analytical form for the effective Murnaghan moduli of a composite with spherical inclusions. The effective moduli depend linearly on Murnaghan moduli of constituents. The results obtained have been verified through numerical modeling using the finite element method.

Tables

Table 1. Table 1: Elastic moduli of polycarbonate and polystyrene (GPa) ( Hughes and Kelly ( 1953 ); Krüger et al. ( 1991 ); Dreiden et al. ( 2011 ) ).

	$K$	$\mu$	$l$	$m$	$n$
Polycarbonate	3.93	0.84	$-50.0$	$-12.2$	$-32.0$
Polystyrene	4.20	1.50	$-18.9$	$-13.3$	$-10.0$

Table 2. Table 2: Elastic moduli of a composite with the volume fraction c = 0.4 % 𝑐 percent 0.4 c=0.4\% for different number of finite elements N k subscript 𝑁 𝑘 N_{k} . Relative moduli of the form Δ ​ K / c = ( K eff − K 0 ) / c Δ 𝐾 𝑐 subscript 𝐾 eff subscript 𝐾 0 𝑐 \Delta K/c=(K_{\rm eff}-K_{0})/c are given in GPa units. Extrapolated values are denoted by the infinity symbol. The last line of the table was calculated using the final analytical expression ( 55 ).

$k$	$N_k$	$\Delta K/c$	$\Delta \mu /c$	$\Delta l/c$	$\Delta m/c$	$\Delta n/c$
0	4736	0.25682	0.50149	28.066	$-3.180$	3.342
1	37888	0.25647	0.49276	28.050	$-3.152$	2.452
2	303104	0.25638	0.49048	28.046	$-3.140$	2.223
$\infty$	extrapolation	0.25635	0.48972	28.045	$-3.137$	2.146
–	theory	0.25635	0.48971	28.045	$-3.135$	2.145

Equations

$$\Pi =\frac{\lambda +2\mu }{2}{\mathrm{I}}_1^2-2\mu {\mathrm{I}}_2+\frac{l+2m}{3}{\mathrm{I}}_1^3-2m{\mathrm{I}}_1{\mathrm{I}}_2+n{\mathrm{I}}_3,$$

$${\mathrm{I}}_1=\mathrm{tr}\hat{\epsilon };{\mathrm{I}}_2=\frac{1}{2}\left((\mathrm{tr}\hat{\epsilon }{)}^2-\mathrm{tr}{\hat{\epsilon }}^2\right);{\mathrm{I}}_3=\det \hat{\epsilon },$$

$$\varepsilon_{i\!j}(\vec{x})=\frac{1}{2}\Big{(}u_{i\!j}(\vec{x})+u_{ji}(\vec{x})+u_{ki}(\vec{x})u_{kj}(\vec{x})\Big{)},$$

$$u_{ij}(\vec{x})=\frac{\partial u_i(\vec{x})}{\partial x_j}$$

$$\Pi (\vec{x})=\frac{1}{2}{C}_{ijkl}(\vec{x})u_{ij}(\vec{x})u_{kl}(\vec{x})+\frac{1}{3}{N}_{ijklmn}(\vec{x})u_{ij}(\vec{x})u_{kl}(\vec{x})u_{mn}(\vec{x})+O(u^4).$$

$$C_{ijkl}=\lambda {\delta }_{ij}{\delta }_{kl}+\mu \left({\delta }_{il}{\delta }_{jk}+{\delta }_{ik}{\delta }_{jl}\right),$$

$$\displaystyle\begin{multlined}N_{i\!jkhts}=\frac{\lambda}{2}\left(\delta_{hk}\delta_{it}\delta_{js}+\delta_{hs}\delta_{i\!j}\delta_{kt}+\delta_{hj}\delta_{ik}\delta_{st}\right)+\frac{\mu}{2}\Big{(}\delta_{hs}\delta_{it}\delta_{jk}\\ +\delta_{ht}\delta_{ik}\delta_{js}+\delta_{hs}\delta_{ik}\delta_{jt}+\delta_{hj}\delta_{it}\delta_{ks}+\delta_{hj}\delta_{is}\delta_{kt}+\delta_{hi}\delta_{js}\delta_{kt}\Big{)}\\ +l\delta_{hk}\delta_{i\!j}\delta_{st}+\frac{m}{2}\Big{(}\delta_{hj}\delta_{ik}\delta_{st}+\delta_{hi}\delta_{jk}\delta_{st}+\delta_{hk}\delta_{it}\delta_{js}+\delta_{hk}\delta_{is}\delta_{jt}\\ +\delta_{ht}\delta_{i\!j}\delta_{ks}+\delta_{hs}\delta_{i\!j}\delta_{kt}-2\delta_{hk}\delta_{i\!j}\delta_{st}\Big{)}\\ +\frac{n}{8}\Big{(}\delta_{ht}\delta_{is}\delta_{jk}+\delta_{hs}\delta_{it}\delta_{jk}-2\delta_{hi}\delta_{jk}\delta_{st}+\delta_{ht}\delta_{ik}\delta_{js}-2\delta_{hk}\delta_{it}\delta_{js}\\ +\delta_{hs}\delta_{ik}\delta_{jt}-2\delta_{hk}\delta_{is}\delta_{jt}-2\delta_{ht}\delta_{i\!j}\delta_{ks}+\delta_{hj}\delta_{it}\delta_{ks}-2\delta_{hj}\delta_{ik}\delta_{st}\\ +\delta_{hj}\delta_{is}\delta_{kt}-2\delta_{hs}\delta_{i\!j}\delta_{kt}+\delta_{hi}\delta_{js}\delta_{kt}+4\delta_{hk}\delta_{i\!j}\delta_{st}+\delta_{hi}\delta_{jt}\delta_{ks}\Big{)},\end{multlined}N_{i\!jkhts}=\frac{\lambda}{2}\left(\delta_{hk}\delta_{it}\delta_{js}+\delta_{hs}\delta_{i\!j}\delta_{kt}+\delta_{hj}\delta_{ik}\delta_{st}\right)+\frac{\mu}{2}\Big{(}\delta_{hs}\delta_{it}\delta_{jk}\\ +\delta_{ht}\delta_{ik}\delta_{js}+\delta_{hs}\delta_{ik}\delta_{jt}+\delta_{hj}\delta_{it}\delta_{ks}+\delta_{hj}\delta_{is}\delta_{kt}+\delta_{hi}\delta_{js}\delta_{kt}\Big{)}\\ +l\delta_{hk}\delta_{i\!j}\delta_{st}+\frac{m}{2}\Big{(}\delta_{hj}\delta_{ik}\delta_{st}+\delta_{hi}\delta_{jk}\delta_{st}+\delta_{hk}\delta_{it}\delta_{js}+\delta_{hk}\delta_{is}\delta_{jt}\\ +\delta_{ht}\delta_{i\!j}\delta_{ks}+\delta_{hs}\delta_{i\!j}\delta_{kt}-2\delta_{hk}\delta_{i\!j}\delta_{st}\Big{)}\\ +\frac{n}{8}\Big{(}\delta_{ht}\delta_{is}\delta_{jk}+\delta_{hs}\delta_{it}\delta_{jk}-2\delta_{hi}\delta_{jk}\delta_{st}+\delta_{ht}\delta_{ik}\delta_{js}-2\delta_{hk}\delta_{it}\delta_{js}\\ +\delta_{hs}\delta_{ik}\delta_{jt}-2\delta_{hk}\delta_{is}\delta_{jt}-2\delta_{ht}\delta_{i\!j}\delta_{ks}+\delta_{hj}\delta_{it}\delta_{ks}-2\delta_{hj}\delta_{ik}\delta_{st}\\ +\delta_{hj}\delta_{is}\delta_{kt}-2\delta_{hs}\delta_{i\!j}\delta_{kt}+\delta_{hi}\delta_{js}\delta_{kt}+4\delta_{hk}\delta_{i\!j}\delta_{st}+\delta_{hi}\delta_{jt}\delta_{ks}\Big{)},$$

$$K=\lambda +\frac{2\mu }{3}.$$

$$u_{i}(\vec{x})\big{|}_{\partial\Omega}=u^{0}_{i\!j}x_{j},$$

$${\overline{u}}_{ij}=\frac{1}{V}{\int }_{\Omega }{u}_{ij}(\vec{x})d\vec{x}.$$

$$\begin{gathered} {\overline{u}}_{ij}=\frac{1}{V}{\int }_{\Omega }{u}_{ij}(\vec{x})d\vec{x}=\frac{1}{V}{\int }_{\partial \Omega }{u}_i(\vec{x})n_{j}dS= \\ \frac{1}{V}{\int }_{\partial \Omega }{u}_{ik}^0{x}_k{n}_{j}dS=\frac{1}{V}{\int }_{\Omega }{u}_{ik}^0{x}_{kj}d\vec{x}=u_{ij}^0. \end{gathered}$$

$$\overline{\Pi }=\frac{1}{V}\int \Pi (\vec{x})d\vec{x}.$$

$$\overline{\Pi }=\frac{1}{2}{C}_{ijkl}^{\mathrm{eff}}{\overline{u}}_{ij}{\overline{u}}_{kl}+\frac{1}{3}{N}_{ijklmn}^{\mathrm{eff}}{\overline{u}}_{ij}{\overline{u}}_{kl}{\overline{u}}_{mn}+O\left({\overline{u}}^4\right),$$

$$u_{ij}(\vec{x})=u_{ij}^L(\vec{x})+u_{ij}^N(\vec{x}),$$

$$\displaystyle u_{i}^{L}\big{|}_{\partial\Omega}=u_{i\!j}^{0}x_{j},$$

$$\displaystyle u_{i}^{N}\big{|}_{\partial\Omega}=0.$$

$$\overline{\Pi}=\frac{1}{V}\int\Bigl{(}\frac{1}{2}C_{i\!jkl}(\vec{x})u_{i\!j}^{L}(\vec{x})u_{kl}^{L}(\vec{x})+C_{i\!jkl}(\vec{x})u_{i\!j}^{L}(\vec{x})u_{kl}^{N}(\vec{x})\\ +\frac{1}{3}N_{i\!jklmn}(\vec{x})u_{i\!j}^{L}(\vec{x})u_{kl}^{L}(\vec{x})u_{mn}^{L}(\vec{x})\Bigl{)}d\vec{x}.$$

$$\overline{\Pi}=\frac{1}{V}\int\Bigl{(}\frac{1}{2}C_{i\!jkl}(\vec{x})u_{i\!j}^{L}(\vec{x})u_{kl}^{L}(\vec{x})+\frac{1}{3}N_{i\!jklmn}(\vec{x})u_{i\!j}^{L}(\vec{x})u_{kl}^{L}(\vec{x})u_{mn}^{L}(\vec{x})\Bigl{)}d\vec{x}.$$

$${\Pi }^{\infty }=\frac{1}{2}{C}_{ijkl}^0{u}_{ij}^{\infty }{u}_{kl}^{\infty }+\frac{1}{3}{N}_{ijklmn}^0{u}_{ij}^{\infty }{u}_{kl}^{\infty }{u}_{mn}^{\infty },$$

$$\overline{\Pi}=\frac{1}{V}\int_{\Omega}\Pi(\vec{x})d\vec{x}=\Pi^{\infty}+\frac{1}{V}\int_{\Omega}\Bigl{(}\Pi(\vec{x})-\Pi^{\infty}\Bigr{)}d\vec{x}.$$

$$\overline{\Pi}=\Pi^{\infty}+\Bigl{<}\Pi^{1}(\vec{x})-\Pi^{\infty}\Bigr{>}_{\vec{x}},$$

$${\Pi }^1(\vec{x})=\frac{1}{2}{C}_{ijkl}(\vec{x})u_{ij}^1(\vec{x})u_{kl}^1(\vec{x})+\frac{1}{3}{N}_{ijklmn}(\vec{x})u_{ij}^1(\vec{x})u_{kl}^1(\vec{x})u_{mn}^1(\vec{x}),$$

$$\Bigl{<}f(\vec{x})\Bigr{>}_{\vec{x}}=\frac{N}{V}\Biggl{\{}\int_{\vec{x}\in\omega_{1}}f(\vec{x})d\vec{x}+(1-c)\int_{\vec{x}\notin\omega_{1}}f(\vec{x})d\vec{x}\Biggr{\}},$$

$$\overline{u}_{i\!j}=\frac{1}{V}\int_{\Omega}u_{i\!j}(\vec{x})d\vec{x}=u_{i\!j}^{\infty}+\Bigl{<}u_{i\!j}^{1}(\vec{x})-u_{i\!j}^{\infty}\Bigr{>}_{\vec{x}}.$$

$$u_{ij}^1(\vec{x})=L_{ijkl}(\vec{x})u_{kl}^{\ast },$$

$$L_{ijkl}(\vec{x})=S_{ijkl}(\vec{x})+L_{ijkl}^{\infty }.$$

$$\hat{L}^{\infty}=\bigl{(}\hat{I}-\bigl{(}\hat{C}^{0}\bigr{)}^{-1}\hat{C}^{1}\bigr{)}^{-1}-\hat{S}^{0}.$$

$$\begin{gathered} S_{ijkl}(\vec{x})=\frac{1}{8\pi (1-{\nu }_0)}\frac{{\partial }^4\Psi }{\partial x_i\partial x_j\partial x_k\partial x_l}-\frac{{\nu }_0}{1-{\nu }_0}\frac{{\delta }_{kl}}{4\pi }\frac{{\partial }^2\Phi }{\partial x_i\partial x_j} \\ -\frac{1}{8\pi }\left({\delta }_{il}\frac{{\partial }^2\Phi }{\partial x_j\partial x_k}+{\delta }_{ik}\frac{{\partial }^2\Phi }{\partial x_j\partial x_l}+{\delta }_{jl}\frac{{\partial }^2\Phi }{\partial x_i\partial x_k}+{\delta }_{jk}\frac{{\partial }^2\Phi }{\partial x_i\partial x_l}\right), \end{gathered}$$

$$\Phi (\vec{x})={\int }_{{\omega }_1}\frac{1}{|\vec{x}-{\vec{x}}_1|}d{\vec{x}}_1,$$

$$\Psi (\vec{x})={\int }_{{\omega }_1}|\vec{x}-{\vec{x}}_1|d{\vec{x}}_1.$$

$$\begin{gathered} {\Pi }^1(\vec{x})=\frac{1}{2}{C}_{ijkl}(\vec{x})L_{ij{i}^{\prime }{j}^{\prime }}(\vec{x})L_{kl{k}^{\prime }{l}^{\prime }}(\vec{x})u_{i^{\prime }{j}^{\prime }}^{\ast }{u}_{k^{\prime }{l}^{\prime }}^{\ast } \\ +\frac{1}{3}{N}_{ijklmn}(\vec{x})L_{ij{i}^{\prime }{j}^{\prime }}(\vec{x})L_{kl{k}^{\prime }{l}^{\prime }}(\vec{x})L_{mn{m}^{\prime }{n}^{\prime }}(\vec{x})u_{i^{\prime }{j}^{\prime }}^{\ast }{u}_{k^{\prime }{l}^{\prime }}^{\ast }{u}_{m^{\prime }{n}^{\prime }}^{\ast }. \end{gathered}$$

$$u_{ij}^{\infty }=L_{ijkl}^{\infty }{u}_{kl}^{\ast }.$$

Full text

††journal: International Journal of Solids and Structures

\geometry

margin=0.55in

1 Introduction

Composite materials are increasingly used in everyday life (Phil and Soutis (2014); Ishikawa et al. (2018); Toozandehjani et al. (2018)). In this paper, we consider composites from the point of view of their macroscopic mechanical properties. A number of techniques have been developed for strengthening of composite materials and for improvement of their elastic properties (Fawaz and Mittal (2015)). These techniques provide a variety of different types of composite materials which can be described by different theoretical models. Laminated composites can be considered as materials with piecewise homogeneous properties (Moulana and Nariboli (1976); Gula and Samsonov (2017)). Composites containing various kinds of particles, usually small-size and evenly distributed, can be regarded in most of models as homogeneous materials (Colombo and Giordano (2011); Giordano (2017)). Fiber-reinforced composites and those with large-size particles (Mortazavian and Fatemi (2015); Giordano (2013)) exhibit intermediate properties (Zaferani (2018)). In macro-mechanical models, all the above-mentioned composites except laminated ones can be considered as homogeneous materials with modified properties.

The use of composite materials allows to achieve unique physical and mechanical properties. Under certain conditions, composites can provide better properties than those of the matrix material and inclusions separately (Koratkar et al. (2005); Uddin and Sun (2008)). Therefore, the determination and prediction of elastic properties of such composite materials are highly desirable for various applications.

In the first approximation, elastic properties of a material are governed by the linear elastic moduli, which characterize elastic stresses of a solid at small strains. Since linear moduli give a quadratic contribution to the elastic energy, they are also called the second-order moduli. Linear elastic properties of composite materials have been extensively studied both in theory and experiments and expressions for second-order moduli of elasticity have been derived (Eshelby (1957); Halpin and Kardos (1976)).

However, as strain increases, nonlinear elastic properties start to contribute to material behaviour. Nonlinear elasticity of a material is governed by the third-order elastic moduli. These nonlinear moduli are responsible, among other things, for the formation and propagation of nonlinear strain solitary waves (solitons for brevity) (Samsonov (2001)).

The first mathematical description of such waves was made by Korteweg and de Vries (Korteweg and De Vries (1895)). At the end of the last century, the theory developed by A. M. Samsonov described the propagation of bulk solitons in solids (see Samsonov (2001) and references therein) and became a basis for first successful experiments on generation and observation of bulk strain solitons in nonlinearly elastic solid waveguides (Dreiden et al. (1988)).

The basic parameter governing soliton formation and propagation in a solid is a nonlinearity coefficient comprising a combination of linear and nonlinear elastic moduli of the material (Samsonov (2001)). Thus, for correct prediction and description of soliton formation and behaviour in a composite waveguide, one should know the effective nonlinear elastic moduli of the composite. The assessment of these values is thus in high demand.

There are several models of materials that take into account nonlinear corrections. In these models both elastic (Murnaghan (1937)) and inelastic (Mooney (1940); Cohen and Dai (1993)) materials are analyzed. We will consider only elastic deformations, which means that after an applied deformation, the body returns to its original non-deformed state. Therefore, we will use the Murnaghan’s theory where the potential energy density is expressed in the form:

[TABLE]

where $I_{1},I_{2},I_{3}$ are invariants of the strain tensor $\hat{\varepsilon}$:

[TABLE]

$l,m,n$ are Murnaghan (third-order) moduli and $\lambda$, $\mu$ are Lame (second-order) parameters. Elastic moduli are coefficients in the power expansion of the potential energy density, and the set of these values defines elastic properties of the material. Third-order elastic moduli show the deviation from the linear elasticity of solids.

To date, several experimental studies have been carried out on measurements of the third-order elastic moduli. However, values measured for the same material can vary considerably in different experiments (Renaud et al. (2016); Winkler and McGowan (2004)). Moreover, measurement error is often huge and can reach 100% (Lang and Gupta (2011); Hughes and Kelly (1953); Korobov et al. (2013)). Therefore, a theoretical approach is needed for an accurate evaluation of nonlinear elastic moduli of composites.

There have been several theoretical studies performed on the subject. The Eshelby’s theory (Eshelby (1957)) for a single ellipsoidal inclusion was generalized to the case of nonlinear elastic properties of inclusion (Sevostianov and Vakulenko (2001); Tsvelodub (2000, 2004)). Macroscopic elastic properties of composites with linear matrix and nonlinear inclusions was obtained by S. Giordano and L. Colombo (Giordano et al. (2009); Colombo and Giordano (2011); Giordano (2017)). In practice, stiffer inclusions are often introduced into a softer matrix, for example, to increase the stiffness of the resulting nanocomposite. Therefore, during deformation of such a composite material, the matrix undergoes greater deformation than inclusions. Accordingly, in this case, the nonlinear properties of the matrix play a more essential role than the nonlinear properties of the inclusions.

In this paper, we consider the general case when both the inclusions and the matrix possess nonlinear elastic moduli with an arbitrary ratio between the stiffness of the inclusions and the matrix. We also assume the applicability of the nonlinear theory of elasticity. In other words, we assume that the inclusions are of a sufficiently large size, at which the molecular structure of the substance does not play a significant role.

The paper is organized as follows. In Section 2 the nonlinear elasticity problem for the composite is formulated. For a given linear and nonlinear elastic properties of composite constituents, we provide a procedure to find effective values of macroscopic moduli. In Section 3 linear iterations are used to solve the nonlinear problem. In Section 4 the Eshelby’s theory is applied to the linearized problem. In Section 5 we find the effective linear and nonlinear moduli of a composite with spherical inclusions. In Section 6 the obtained formulas are verified using the finite element method. Finally, in Section 7 we discuss the results obtained.

2 Formulation of the problem

Let us consider a certain volume $\Omega$, where a composite consisting of a matrix and inclusions is located. To determine effective elastic moduli of the composite, let us apply some external deformation to it. The following strain tensor can be written for a point $\vec{x}$ in equilibrium:

[TABLE]

where

[TABLE]

is the distortion tensor (strain gradient).

Since we are interested in nonlinear properties of the composite, the nonlinear term in the definition of the strain tensor (3) cannot be considered small. We assume that deformations are finite, and cannot be considered small enough to disregard the nonlinear term in (3). We assume also that deformations are elastic, which means that after unloading, the body returns to its original state.

At each point $\vec{x}$, the potential energy density $\Pi(\vec{x})$ depends on the strain tensor $\varepsilon_{i\!j}(\vec{x})$ and on the values of elastic moduli of the material at that point. In general, this dependence can be represented as an expansion in powers of $u_{i\!j}(\vec{x})$:

[TABLE]

For homogeneous and isotropic body, linear and nonlinear stiffness tensors have the following form:

[TABLE]

where $l,m,n$ are Murnaghan moduli and $\lambda$, $\mu$ are Lame parameters. Further, we will mainly use a pair of linear elastic moduli $K$, $\mu$, since the resulting equation has a simpler form in these notations. Here $K$ is the bulk elastic modulus, which is expressed through $\lambda$, $\mu$ in the form:

[TABLE]

In a general case, the elastic properties of the matrix and inclusions are different. For both materials we have different sets of linear and non-linear elastic moduli: $K_{0},\mu_{0},l_{0},m_{0},n_{0}$ for the matrix and $K_{1},\mu_{1},l_{1},m_{1},n_{1}$ for the inclusions. Thus, volumes $\omega_{i}$ containing inclusions are characterized by elastic moduli tensors $\hat{C}^{1}$ and $\hat{N}^{1}$. The remaining volume $\Omega\setminus(\omega_{1}\cup\omega_{2}\cup...)$ containing matrix material is characterized by elastic moduli tensors $\hat{C}^{0}$ and $\hat{N}^{0}$. To obtain effective macroscopic elastic moduli of the composite, we need to average the elastic moduli over the composite volume $\Omega$ in a special way.

To apply a deformation to the composite, we assume that the boundary conditions have a form

[TABLE]

where the boundary deformation $u^{0}_{i\!j}$ is given and does not depend on $\vec{x}$. In the case of a homogeneous medium in equilibrium, we would have a homogeneous deformation, which is described by a displacement $u_{i}(\vec{x})=u^{0}_{i\!j}x_{j}$. However, for a composite or any other inhomogeneous medium, the displacement $u_{i}(\vec{x})$ is not a linear function of $\vec{x}$ coordinate. Consider a volume-average distortion tensor

[TABLE]

Note that the equality $\overline{u}_{i\!j}=u^{0}_{i\!j}$ is satisfied for a nonlinear elastic medium as well. Indeed,

[TABLE]

For each applied deformation (defined by $\overline{u}_{i\!j}=u^{0}_{i\!j}$), the average elastic energy of the composite is

[TABLE]

It can be represented as an expansion in powers of $\overline{u}_{i\!j}$:

[TABLE]

which defines the effective elastic moduli tensors $\hat{C}^{\rm eff}$ and $\hat{N}^{\rm eff}$. Since the composite is homogeneous and isotropic at the macroscopic scale, the effective elastic moduli tensors $\hat{C}^{\rm eff}$ and $\hat{N}^{\rm eff}$ have the same general form as $\hat{C}$ and $\hat{N}$ (Eqs. (6) and (14)) with effective elastic moduli $\lambda_{\rm eff}$, $\mu_{\rm eff}$, $l_{\rm eff}$, $m_{\rm eff}$, $n_{\rm eff}$.

Equation (20) can be solved both analytically and numerically. In Sections 3–5 we perform an analytical averaging of the elastic energy and find an analytical expression for the effective elastic moduli. In Section 6 we find numerical values of effective moduli using the finite element method.

3 Reduction to a linear problem and assumptions

Let us solve the problem iteratively by the method of successive approximations. The distortion tensor $u_{i\!j}(\vec{x})$ can be written as

[TABLE]

where $u_{i\!j}^{L}(\vec{x})$ is the solution to the linear problem and $u_{i\!j}^{N}(\vec{x})\ll u_{i\!j}^{L}(\vec{x})$ is a small nonlinear contribution. Then the boundary conditions take the form:

[TABLE]

Using Eqs. (5) and (21) and taking into account that $u_{i\!j}^{N}(\vec{x})\ll u_{i\!j}^{L}(\vec{x})$, the average density of potential energy can be written in the form:

[TABLE]

Integrating the second term in parts and applying the Gauss theorem (Giordano (2013)), we note that

[TABLE]

Thus, in order to find the average density of potential energy with the required accuracy, it is necessary to calculate the linear deformation only, i.e. to solve a linear problem.

We assume that all inclusions in the composite are the same and are located far enough from each other. Therefore, interaction between inclusions can be neglected. Then the density of the elastic energy far from inclusions can be expressed as:

[TABLE]

where $u^{\infty}_{i\!j}\neq\overline{u}_{i\!j}$ is the distortion tensor far from inclusions. Therefore, we can write the average density of potential energy as

[TABLE]

One can see that ${\Pi(\vec{x})-\Pi^{\infty}}$ goes to zero far from inclusions. Therefore, the integral in Eq. (27) can be calculated using the contribution from each inclusion independently:

[TABLE]

where $\Pi^{1}(\vec{x})$ is the contribution to the potential energy density from one inclusion:

[TABLE]

with $u^{1}_{i\!j}(\vec{x})$ being the distortion around the inclusion $\omega_{1}$ and inside it. The averaging in Eq. (28) is performed using the following homogenization procedure:

[TABLE]

where $N$ is the number of inclusions and $c$ is the volume fraction of inclusions. The first integral in Eq. (30) is performed over an inclusion volume $\omega_{1}$. The second integral is performed over the surrounding matrix. The coefficient $(1-c)$ takes into account that the surrounding matrix is partially occupied by other inclusions. It is the simplest procedure, which takes into account the presence of neighbor inclusions. As we will see in Section 6, it gives accurate effective moduli in a wide range of the volume fraction $c$.

Similarly, the average distortion tensor $\overline{u}_{i\!j}$ can be written in the form

[TABLE]

Thus, we can consider inclusions separately from each other and find $u_{i\!j}^{1}(\vec{x})$ and $\Pi^{1}(\vec{x})$. Then, the effective elastic moduli can be found using Eqs. (28), (31) and (20).

4 Eshelby’s theory

A solution to the linear problem is known from the theory developed by Eshelby (Eshelby (1957); Colombo and Giordano (2011)). The strain field of a matrix with a single inclusion $\omega_{1}$ can be expressed using the “eigen-strain” parameter $u_{i\!j}^{*}$:

[TABLE]

where

[TABLE]

Tensor $\hat{L}^{\infty}$ is defined as

[TABLE]

The Eshelby’s tensor is defined as (Eshelby (1957))

[TABLE]

where harmonic $\Phi(\vec{x})$ and bi-harmonic $\Psi(\vec{x})$ potentials have the form

[TABLE]

For elliptic inclusions, the Eshelby’s tensor is constant inside inclusions: ${\hat{S}(\vec{x})=\hat{S}^{0}}$ for $\vec{x}\in\omega_{1}$. As a result, the potential energy density from a single particle $\Pi^{1}(\vec{x})$ has the form

[TABLE]

The Eshelby’s theory also allows to find the distortion tensor far from the inclusion $\hat{u}^{\infty}$:

[TABLE]

Thus, the density of potential energy far from inclusion has a form

[TABLE]

5 Effective elastic moduli of a composite with spherical inclusions

Substituting (38) and (40) into (28), we obtain the average deformation of the composite expressed through the eigen-strain parameter $\hat{u}^{*}$:

[TABLE]

where

[TABLE]

On the other hand, by substitution of Eqs. (39) and (32) into Eq. (31) we can also obtain the average deformation of the composite expressed through the eigen-strain parameter $\hat{u}^{*}$:

[TABLE]

where

[TABLE]

Then, according to (20), we can express the average density of potential energy by using the tensor $\hat{L}^{*}$ and the eigen-strain parameter $\hat{u}^{*}$:

[TABLE]

Comparing equations (41) and (48) for the average density of potential energy of the composite, we obtain two systems of algebraic equations for elements of the effective linear and nonlinear stiffness tensors $\hat{C}^{\rm eff}$ and $\hat{N}^{\rm eff}$. The system of equations for the linear moduli has a form:

[TABLE]

The system of equations for the nonlinear moduli has a similar form:

[TABLE]

The harmonic and bi-harmonic potentials have a simple form for spherical inclusions of radius $R$:

[TABLE]

Using the computer algebra system, we make the homogenization procedure in Eqs. (42), (45), (47) and solve systems of equations (49) and (50). As a result, we get effective linear and nonlinear moduli of a composite with spherical inclusions. The linear effective elastic moduli can be expressed in terms of elastic moduli of the matrix and volume fraction of inclusions $c=4\pi NR^{3}/3V$:

[TABLE]

This result coincides with Eqs. (15) and (16) in (Giordano et al. (2008)). Nonlinear effective elastic moduli of the composite take the form:

[TABLE]

Here $\mathbf{P}^{0}$ and $\mathbf{P}^{1}$ are $3\times 3$ matrices, which depend on the linear moduli of the matrix $K_{0},\mu_{0}$ and those of inclusions $K_{1},\mu_{1}$. In general, these matrices also depend on the volume fraction $c$. Nonlinear moduli $l_{g}$, $m_{g}$, $n_{g}$ provide the “geometric” nonlinear contribution which exists even if all Murnaghan moduli of the matrix and inclusions are zero. They also depend only on linear moduli $K_{0},\mu_{0},K_{1},\mu_{1}$ and the volume fraction $c$. The analytical expressions for $\mathbf{P}^{0}$, $\mathbf{P}^{1}$, $l_{g}$, $m_{g}$, $n_{g}$ are presented in A.

The same procedure can be applied to obtain effective nonlinear moduli for composites with ellipsoidal inclusions (aligned or randomly oriented). However, it involves more complicated integrals in the homogenization procedure (30) in Eqs. (42), (45), (47). For randomly oriented ellipsoids we expect the same form of the effective moduli (55) with coefficients depending also on the ellipsoidal aspect ratio. The proposed method can also be used to take into account spatial correlations of the inclusion positions. For example, one can consider so-called ellipsoidal microstructure with ellipsoidal symmetry of spatial correlations between the inclusions positions (Giordano (2017)). It requires an accurate modification of the homogenization procedure (30), which should be the subject of further work.

6 Numerical verification

In order to verify the obtained equations, we performed numerical calculations using the finite element method. Note that to obtain the most reliable results, it is important to verify a general case of the final equation (55). Observe that in the case of relatively rigid inclusions (e.g. silica inclusions in a polymeric matrix) equations can be additionally simplified. Therefore, we have chosen a pair of materials, with elastic moduli of the same order of magnitude. We use two materials with known linear and nonlinear moduli with different Poisson ratios: polycarbonate matrix and polystyrene inclusions. Elastic moduli of both materials are given in Table 1. We emphasize that the final result (55) can be applied for any pair of materials.

Let us consider a cubic sample $L\times L\times L$ with a small inclusion of a radius $R$ located in the center (Fig. 2). We apply a number of different trial distortions $u_{\smash{ij}}^{0}$ to obtain all elements of linear and nonlinear stiffness tensors $\hat{C}^{\rm eff}$ and $\hat{N}^{\rm eff}$ ($3^{4}+3^{6}=810$ in total). In order to reduce the influence of boundaries, we use periodic boundary conditions. In this case, we find the displacement in the form

[TABLE]

with $\tilde{u}_{i}(\vec{x})$ being the periodic function with the period $L$ in $x,y,z$ directions. The sample with one inclusion and periodic boundary conditions effectively acts as an infinite composite with inclusions placed at nodes of a simple cubic lattice.

We use trial distortions of the form $u^{0}_{\smash{ij}}=\epsilon k_{\smash{ij}}$ where $\epsilon=10^{\smash{-5}}$ and $k_{\smash{ij}}$ are integer numbers, satisfying the condition $\sum_{ij}k_{\smash{ij}}^{\smash{2}}\leq 2$. There are 163 trial distortions in total. For each distortion, we calculate all 9 elements of the average stress tensor $\hat{\overline{\sigma}}$, which have the following relation with the trial distortion $\hat{u}^{0}$:

[TABLE]

As a result, we obtain $9\cdot 163=1467$ equations to find 810 unknown elements of $\hat{C}^{\rm eff}$ and $\hat{N}^{\rm eff}$. Since the number of equations is bigger than the number of unknown parameters, we use the least squares method to find the best solution. After the solution was found, we have checked that all equations were valid with a high enough precision. We deliberately do not rely on various symmetry properties of $\hat{C}^{\rm eff}$ and $\hat{N}^{\rm eff}$ to additionally verify the computations.

After all elements of the stiffness tensors $C^{\rm eff}$ and $N^{\rm eff}$ have been obtained, we could find the effective elastic moduli of the composite, $\lambda_{\rm eff},\mu_{\rm eff},l_{\rm eff},m_{\rm eff},n_{\rm eff}$, using Eqs. (6) and (14). Once again, we used the least squares method to find the best values for the effective elastic moduli (5 values from 810 equations). For small-size inclusions the sample is effectively isotropic. For moderate values of the radius $R$, there is a small cubic anisotropy in the macroscopic behavior of the studied composite. However, the least squares method effectively takes the isotropic component of effective stiffness tensors.

We used Gmsh 2.10.1 (Geuzaine and Remacle (2009)) to prepare the hexagonal mesh (Fig. 2). The geometry of the mesh was chosen to maximize the precision of the results for given computation time. Then we used FEniCS 2018.2.0 to apply the finite element method (Alnæs et al. (2015)). It is open-source software, which allows to solve arbitrary variational problems, including the hyperelastic problems with arbitrary constituent relations.

The final equation (55) is precise for a small volume fraction of inclusions. In order to check this result, we consider a composite with a small inclusion volume fraction $c=0.4\%$, which corresponds to the inclusion radius $R\approx L/10$ (Fig. 2). We calculate effective moduli using several meshes with different mesh density. The basic mesh consists of $N_{0}=4736$ finite elements and is presented in Fig. 2. Then we iteratively double the number of finite elements in each direction to obtain two refined meshes with $N_{k}=8^{k}N_{0}$ finite elements for $k=1,2$. For each number of finite elements $N_{k}$, we obtain effective moduli $K_{\rm eff}^{(k)},\mu_{\rm eff}^{(k)},l_{\rm eff}^{(k)},m_{\rm eff}^{(k)},n_{\rm eff}^{(k)}$. The precision of the obtained moduli is scaled quadratically with the size of finite elements: $K_{\rm eff}^{(1)}-K_{\rm eff}^{(0)}\approx 4(K_{\rm eff}^{(2)}-K_{\rm eff}^{(1)})$ and the same for other moduli. As a result, we can extrapolate moduli for $N_{k}\to\infty$ as $K_{\rm eff}^{(\infty)}=(4K_{\rm eff}^{(2)}-K_{\rm eff}^{(1)})/3$. The results are shown in Table 2 as relative moduli $\Delta K/c=(K_{\rm eff}-K_{0})/c$, $\Delta\mu/c=(\mu_{\rm eff}-\mu_{0})/c$, etc, which mean the influence of inclusions on the effective moduli of the composite. Table 2 shows that the extrapolated values of the relative moduli are very close to the theoretical values: all of them coincide with 3–5 decimal places. Therefore, we have validated our main result (55) for a small volume fraction of inclusions.

In order to study variable elastic contrast of the composite, we multiply all elastic moduli of the inclusion (linear and nonlinear) by a coefficient $\alpha$. In this case matrix and inclusions have elastic moduli $K_{0},\mu_{0},l_{0},m_{0},n_{0}$ and $\alpha K_{1},\alpha\mu_{1},\alpha l_{1},\alpha m_{1},\alpha n_{1}$ respectively. For a small value of the parameter $\alpha\ll 1$, we have relatively soft inclusions. In the opposite case, $\alpha\gg 1$, inclusions are rigid. In both limiting cases, elastic moduli of inclusions do not affect the effective moduli of the composite.

For different values of the parameter $\alpha$, we apply the same finite element method with extrapolation procedure for the same volume fraction $c=0.4\%$. Figure 3 shows that the final equation (55) coincide with the numerical results for any value of the parameter $\alpha$. The point $\alpha=1$ corresponds to the values presented in Table 2. For soft inclusions ($\alpha\ll 1$) and rigid inclusions ($\alpha\gg 1$), one can see an increased influence of inclusions on effective elastic moduli.

The final equation (55) can be applied for the finite volume fraction of inclusions $c$. In this equation, we do not take into account a precise form of the interaction between inclusions. Therefore, the final equation (55) is approximate for finite values of $c$. In order to check the precision of the obtained effective moduli, we use the same finite element method for a variable inclusion radius $R$ using the mesh with a fixed number of finite elements $N_{1}=37888$. The effective moduli are shown in Fig. 4 for different volume fraction $c$ and different elastic contrasts: $\alpha=0.01$ (soft inclusions), $\alpha=1$ (polystyrene inclusions), and $\alpha=100$ (rigid inclusions). Despite the fact that the interaction between inclusions was not calculated precisely, the theory agrees well with the numerical results. The deviation of the effective moduli rarely exceeds 1%. The cases $\alpha=0.01$ and $\alpha=100$ demonstrate a noticeable nonlinear dependence on the effective moduli on the volume fraction $c$, which also coincides with the theory. There is only one visible deviation between the theory and the numerical calculation, which is the modulus $l_{\rm eff}$ for soft inclusions ($\alpha=0.01$) for $c\geq 0.1$. However, this modulus has the strongest sensitivity on the volume fraction $c$. This deviation can also be caused by the numerical calculation which considers one inclusion in a periodic cell. It effectively acts as an infinite composite with inclusions placed in a simple cubic lattice. The precise calculation of effective moduli of a large volume with randomly located inclusions requires much more computational resources. Accounting of the mutual arrangement of inclusions is the subject of future studies for both theory and numerical calculations.

7 Conclusions

In this paper we have derived analytical equations governing the effective nonlinear elastic moduli of a composite with a nonlinear matrix and nonlinear elastic spherical inclusions. The effective nonlinear moduli $l_{\rm eff},m_{\rm eff},n_{\rm eff}$ depend linearly on the nonlinear moduli of the constituents $l_{0},m_{0},n_{0},l_{1},m_{1},n_{1}$. At the same time, the dependence on the linear moduli $K_{0},\mu_{0},K_{1},\mu_{1}$ has a complicated nonlinear form. We have found the “geometric” contribution $l_{g},m_{g},n_{g}$ to the effective nonlinear moduli, which is due to the inhomogeneous structure of the composite and is defined solely by the linear moduli. The proposed methodology can be generalized to a case of more complex inclusions, such as e.g. those of ellipsoidal shape.

The obtained expressions have been verified using the finite element method for various elastic contrast of inclusions. For a small volume fraction of inclusions $c$, the obtained results were verified with 3–5 significant digits. For finite values of $c$, the analytical expressions agree well with the numerical data and the deviation of the effective moduli rarely exceeds 1%.

Despite a rather complicated analytical form, final equations can be easily used for computation of effective elastic moduli for given elastic moduli of the matrix and inclusions. Various examples of this procedure can be found on GitHub (Beltukov and Semenov (2019)) in the form of Wolfram Mathematica function, Python code, or an interactive XLS table.

The financial support from Russian Science Foundation (RSF), under the grant no. 17-72-20201, is gratefully acknowledged.

Appendix A Coefficients of the final equation

In this appendix we provide the matrices $\mathbf{P}^{0}$ and $\mathbf{P}^{1}$ and geometric nonlinearity modules $l_{g}$, $m_{g}$, $n_{g}$, which are included in the final equation (55). In order to shorten the expressions, we introduce the following notations:

[TABLE]

Note that $e$ and $f$ are denominators in the linear effective moduli (54) and (53) respectively. Using the above notation, elements of the matrices $\mathbf{P}^{0}$ and $\mathbf{P}^{1}$ can be written in the following form

[TABLE]

The geometrical contribution to nonlinear moduli has the following form:

[TABLE]