1 Introduction

Electroactive continua (see e.g. [31, 55] and references therein for a sufficiently comprehensive general theory of the continuum mechanics of electromagnetic solids) are deformable solids that respond by changing shape and size when subjected to the application of an electric field. Materials characterised by such properties have been increasingly investigated in the past decades due to their applicability to real-world scenarios of interest. Relevant contexts mostly involving electroactive (dielectric) polymers include, but are not limited to, bio-inspired, biomedical applications, such as artificial muscles [5], tunable optical lenses [13], and sustainable battery materials [52], as well as actuator technology [56].

The appearance of strains caused by the application of an electric field can occur via several mechanisms, which are different in nature. Mechanical strains linearly proportional to an applied electric field can develop in an electric continuum medium as a consequence of the so called inverse piezoelectric effect, whose direct counterpart, which is the celebrated direct piezoelectric effect, entails the generation of an electric polarisation via mechanical strains. Piezoelectric materials are widely used when the conversion of mechanical into electric energy, and vice versa, is desired, as for electromechanical transducers [1, 65].

Electroactive materials can deform via a second-order mechanism, known as electrostriction, which, as opposed to the inverse piezoelectricity, is quadratic in the electric field and therefore gives rise to strains and stresses that do not change when the electric field is reversed.

The latter phenomenon plays a major role in driving deformations in electromechanical continua for sufficiently high electric fields, and when piezoelectric effects are not relevant, i.e. when the material possesses a centre of symmetry. These materials do not change the state of polarisation when subjected to purely mechanical loading, i.e. in this sense, pure converse electrostriction does not appear in this case due to symmetry [31]. However, once an applied electric field has changed the polarisation state, the mechanical strains that arise due to electrostriction are technically affecting in turn the polarisation, as observed in [42], and exploited in [26] to measure the electrostrictive constants via an extremely sensitive instrumentation. There also exist higher-order electro-elastic effects and ferroelectric materials, the latter being characterised by a permanent change in their polarisation state even when the electric field is removed [31].

In this work, we focus on purely electrostrictive elastic composites (see e.g. [32, 49]), which are often encountered in bio-inspired materials and electronic devices, thus neglecting both piezoelectric and ferroelectric effects, as well as converse electrostriction and higher-order electro-elastic phenomena.

We aim at describing the effective behaviour of electrostrictive composites in a linear elastic regime, which can be relevant for applications, for example, related to enhancing the sensitivity of solid-state capacitance sensors, see e.g. [71].

In general, investigation of elastic composites (see, for example, [16, 39, 41, 51] among many others) and their response based on the individual composites’ constituents, as well as their geometrical arrangement and interactions, can be used both to optimise the performance of artificial media, and to achieve a more thorough understanding of specific physical phenomena (see e.g. [63] in the context of aged bone modelling).

As a matter of fact, it is often either non-trivial or virtually impossible to describe the behaviour of a whole composite (i.e. on the scale characterising its size) by resolving every detail of the composite’s micro-structure, so that the multi-scale, possibly hierarchical nature of these materials is typically addressed by means of suitable homogenisation techniques.

The latter are normally designed to provide a macro-scale description of the physical system at hand, which should (a) be conceived in the light of possible validation against experiments that take place on the coarse scale characterising the whole composite; (b) encode information concerning the micro-scale constituents in the composite and their mutual interplay; (c) be computationally feasible.

Although there exists a large variety of homogenisation approaches, these can be broadly summarised in two groups, i.e. average field techniques (these include both RVE and Eshelby based techniques) and the asymptotic homogenisation technique, see e.g. the review [38] and the work [61] for a theoretical and a computational comparison between them, respectively.

The representative volume element (RVE) approach relies on considering a volume element which is assumed to be (statistically and energetically) representative of the whole composite, and simultaneously sufficiently small such that the effective material properties can be computed at a reduced computational cost. A remarkable example of application of the latter technique to electro-sensitive elastomers can be found in [14, 43], where the authors apply the so called variational and computational homogenisation techniques, respectively, to nonlinear electro-sensitive elastomers by considering the energy functionals associated with a representative volume element.Footnote 1

Widely exploited average field techniques also involve the results by Eshelby [25], where the RVE is identified with an infinite medium characterised by uniform strain condition at infinity, and the individual constituents are represented as ellipsoidal inclusions. Relevant examples of Eshelby-based techniques, such as the Mori–Tanaka [53] and the self-consistent [36] schemes, applied to electrostrictive composites can be found in [8, 23, 48], respectively.

The alternative Asymptotic Homogenisation technique (see e.g. [3, 4, 37, 50, 57, 70]) can be used to determine the effective balance equations for composites exploiting the sharp length scale separation which exists in the system between the micro- and macro-scales.

Spatial variations are therefore decoupled, and a multiple scale expansion of the fields is performed to obtain differential relationships which are enforced, possibly, but not necessarily, under the assumption of micro-scale periodicity, to derive the effective balance equations for the medium.

The information concerning the fine scale properties and arrangement of the composite sub-phases is then embedded in the homogenised coefficients, which are to be computed by means of closed-form analytical formulae containing auxiliary quantities to be (often numerically) computed by solving elastic-type partial differential equations (PDEs).

The latter operation typically happens in the context of local periodicity, such that fine-scale geometrical and functional complexities are to be considered only on a periodic cell, thus enhancing computational feasibility.

There exists a number of paradigmatic works dealing with homogenised modelling and computation of the effective permittivity and piezoelectric coefficients via asymptotic homogenisation, see for example [9, 15] for laminated and fibre-reinforced piezoelectric composites, respectively.

Here we focus on establishing a homogenised modelling framework which is capable of providing the effective governing equations (and related effective coefficients) for an electrostrictive composite subjected to an applied given electric field. As we aim at retaining precise information concerning the composite micro-structure, we then embrace the asymptotic homogenisation technique to achieve our goal.

In [47], the Authors apply the asymptotic homogenisation technique to electrostrictive non-piezoelectric composites characterised by a regime similar to the one embraced in the current work, namely, small deformations and moderate electric fields. However, the Authors in [47] adopt a different standpoint as they consider a one-way coupled problem in which the electric potential is seen as a (multiscale) variable solving a generalised Poisson’s problem driven by rapidly oscillating charges.

Here we follow a different approach via considering the electric field as given, and addressing the asymptotic homogenisation of a linear elastic composite whose deformations are driven by the divergence of a tensor, referred to as generalised Maxwell stress tensor and supplied from the outset. As such, we provide a novel formulation in the context of asymptotic homogenisation by following an approach consistent with [64], where the Authors deduce the effective balance equations for a material composite subjected to a continuous volume force given by the Helmholtz decomposition. The two formulations are both addressing homogenisation of a locally unbounded problem and, for completeness, are shown to coincide under a consistent set of simplifying assumptions. In particular, the framework in [64] in absence of a vector potential can be recovered if the electric force is identified with the gradient of the potential introduced in [64] and if the electric properties are assumed to be continuous across the interface of the composite’s sub-phases. The formulation reported in [47] can be likewise considered as a particular case of our model for an appropriate specific choice of the imposed electric field and in absence of free charges.

The remainder of this work is organised as follows. In Sect. 2, we set the notation of the work and we illustrate the multi-scale balance equations of an electrostrictive composite subjected to the action of an imposed electric field. In Sect. 3, we introduce the hypothesis of sharp separation between the macro- and micro-scale, and we show the main mathematical assumptions and tools of the asymptotic homogenisation technique employed throughout the work. In Sect. 4, with reference to some specific choices of the elementary cell’s configuration, we review some topological and geometrical properties of the micro-structure of the considered composite. In Sect. 5, we derive the main result, i.e. the effective equilibrium equations for the composite under study. In Sect. 6 we specialise our model by focusing on two relevant particular cases. In Sect. 6.1, we investigate the case in which the imposed electric field is assumed to be uniform at the micro-scale. In Sect. 6.2, we assume a multiplicative decomposition of the electric field into components which are purely varying on the micro- and macro-scales. We further highlight under which circumstances this assumption leads to a reduced computational complexity of the model whilst retaining spatial variations of the imposed electric field on both scales. In Sect. 7, we compare our framework with other approaches to electrostriction available in the literature. Finally, in Sect. 8, we highlight the main results of the work and we outline some future developments of the modelling framework developed in this manuscript.

2 Kinematics and dynamics of electro-sensitive composites

Let us denote by \({\mathscr {C}}\) a material composite consisting of a host medium (or host phase), referred to as the matrix, in which a family of disjoint (solid) sub-phases is embedded. In this work, we specialise our research to investigate the mechanical behaviour of linear elastic electro-sensitive composites, for which both the matrix and the sub-phases behave like linear elastic media, whose macroscopic stress response is susceptible to the presence of an electric field [22, 75] (see the discussion after Eq. (64) below). Note that it is in the sense of this susceptibility of the stress that we speak of “electro-sensitive composites”.

Let \({\mathscr {C}}_{0}\) be the placement of the composite \({\mathscr {C}}\) in the physical space \({\mathscr {S}}\), corresponding, in this work, to the three dimensional Euclidean space. Furthermore, let us introduce two open subsets of \({\mathscr {C}}_{0}\), \({\mathscr {C}}_{\mathrm {M}}\) and \({\mathscr {C}}_{\mathrm {S}}\), representing the matrix and the sub-phases, respectively. In particular, \({\mathscr {C}}_{\mathrm {M}}\) and \({\mathscr {C}}_{\mathrm {S}}\) define a partition of \({\mathscr {C}}_{0}\), such that \({\bar{\mathscr {C}}}_{0}={\bar{\mathscr {C}}}_{\mathrm {M}}\cup {\bar{\mathscr {C}}}_{\mathrm {S}}\), \({\bar{\mathscr {C}}}_{\mathrm {M}}\cap {\mathscr {C}}_{\mathrm {S}}=\emptyset \) and \({\mathscr {C}}_{\mathrm {M}}\cap {\bar{\mathscr {C}}}_{\mathrm {S}}=\emptyset \), with the bar over a set denoting its topological closure. The (sharp) interface separating \({\mathscr {C}}_{\mathrm {M}}\) and \({\mathscr {C}}_{\mathrm {S}}\) is denoted by \(\Gamma _{0}\). Before going further, we mention that, using the jargon introduced in previous works [67,68,69], we distinguish among three different types of sub-phases, i.e. inclusions, fibres and strata, as will be clarified in Sect. 4. For future use, we introduce the three-dimensional inner product space \({\mathcal {V}}\), which represents the translational space associated with \({\mathscr {S}}\), and \(\mathrm {Lin}({\mathcal {V}})\), denoting the space of all linear endomorphisms from \({\mathcal {V}}\) into itself. In fact, the elements of \(\mathrm {Lin}({\mathcal {V}})\) are second order tensor fields [54].

Finally, we define a system of Cartesian coordinates on \({\mathscr {C}}_{0}\), namely, we associate a triple of real numbers \(x\equiv (x_{1},x_{2},x_{3})\) with each point \(x_{0}\in {\mathscr {C}}_{0}\), thereby representing the coordinates of \(x_{0}\) with respect to a Cartesian system of reference. We refer to x as spatial physical coordinates associated with the point \(x_{0}\) [3, 60]. From now on, each quantity of interest defined in \({\mathscr {C}}_{0}\), and with values in a suitable functional space, can be expressed as a function depending on the physical coordinates [24]. Hence, for the sake of a lighter notation, throughout the work, we use the same symbol both for the functional expression of a given quantity defined in \({\mathscr {C}}_{0}\) and for its counterpart rephrased in terms of the spatial physical coordinates [24].

Remark 1

(Component-wise representation of vectors and tensors) In this remark, we fuss over some notational conventions used to express vector and tensor quantities in components.

By adopting Cartesian coordinates and introducing the orthonormal vector basis \(\{{\varvec{i}}_{a}\}_{a=1}^{3}\subset {\mathcal {V}}\), the notation for the component-wise representation of vector and tensor quantities is as follows: for a vector field \({\varvec{v}}\in {\mathcal {V}}\) and a second-order tensor field \({\varvec{A}}\in \mathrm {Lin}\left( {\mathcal {V}} \right) \), we write

$$\begin{aligned}&{\varvec{v}}=v_{a}{\varvec{i}}_{a}, \end{aligned}$$
(1)
$$\begin{aligned}&{\varvec{A}}=A_{ab}\,{\varvec{i}}_{a}\otimes {\varvec{i}}_{b}, \end{aligned}$$
(2)

where Einstein’s convention of summation over repeated indices is employed. Analogously, a third-order tensor field \({\varvec{\beta }}\) and a fourth-order tensor field \({\mathbb {D}}\) admit the following representations in components

$$\begin{aligned}&{\varvec{\beta }}=\beta _{abc}\,{\varvec{i}}_{a}\otimes {\varvec{i}}_{b}\otimes {\varvec{i}}_{c},\end{aligned}$$
(3)
$$\begin{aligned}&{\mathbb {D}}={\mathbb {D}}_{abcd}\,{\varvec{i}}_{a}\otimes {\varvec{i}}_{b}\otimes {\varvec{i}}_{c}\otimes {\varvec{i}}_{d}. \end{aligned}$$
(4)

Note that, throughout this work, third-order tensors will be viewed as linear maps transforming second-order tensors into vectors, whereas fourth-order tensors will be defined as linear maps from \(\mathrm {Lin}({\mathcal {V}})\) into itself. Accordingly, we write \({\varvec{\beta }}:\mathrm {Lin}({\mathcal {V}})\rightarrow {\mathcal {V}}\) and \({\mathbb {D}}:\mathrm {Lin}({\mathcal {V}})\rightarrow \mathrm {Lin}({\mathcal {V}})\). In particular, when needed, we shall use the notation \({\varvec{\beta }}\in \mathrm {Lin}(\mathrm {Lin}({\mathcal {V}}),{\mathcal {V}})\) for third-order tensors and the notation \({\mathbb {D}}\in \mathrm {Lin}(\mathrm {Lin}({\mathcal {V}}))\) for fourth-order tensors. Note that, to represent endomorphisms, we use here the notation \(\mathrm {Lin}({\mathcal {V}})\) and \(\mathrm {Lin}(\mathrm {Lin}({\mathcal {V}}))\) in lieu of the less compact one \(\mathrm {Lin}({\mathcal {V}},{\mathcal {V}})\) and \(\mathrm {Lin}(\mathrm {Lin}({\mathcal {V}}),\mathrm {Lin}({\mathcal {V}}))\). Moreover, given \({\varvec{A}}\in \mathrm {Lin}({\mathcal {V}})\), \({\varvec{\beta }}\in \mathrm {Lin}(\mathrm {Lin}({\mathcal {V}}),{\mathcal {V}})\) and \({\mathbb {A}},{\mathbb {D}}\in \mathrm {Lin}(\mathrm {Lin}({\mathcal {V}}))\), for denoting double contractions, we write

$$\begin{aligned}&\left[ {\varvec{\beta }}:{\varvec{A}}\right] _{a}=\beta _{abc}\,A_{bc}, \end{aligned}$$
(5a)
$$\begin{aligned}&\left[ {\mathbb {D}}:{\varvec{A}}\right] _{ab}={\mathbb {D}}_{abcd}\,A_{cd}, \end{aligned}$$
(5b)
$$\begin{aligned}&\left[ {\mathbb {D}}:{\mathbb {A}}\right] _{abmn}={\mathbb {D}}_{abcd}\,A_{cdmn}. \end{aligned}$$
(5c)

We also define the transposed fourth-order tensor \({\mathbb {D}}^{\mathrm {T}}\) through the relation \({\varvec{A}}:{\mathbb {D}}:{\varvec{B}}={\varvec{B}}:{\mathbb {D}}^{\mathrm {T}}:{\varvec{A}}\), for all \({\varvec{A}},{\varvec{B}}\in \mathrm {Lin}({\mathcal {V}})\). In index notation, we have \([{\mathbb {D}}^{\mathrm {T}}]_{abcd}={\mathbb {D}}_{cdab}\). Finally, when a tensor quantity, such as a fourth-order tensor, is explicitly associated with the subset \({\mathscr {C}}_{\alpha }\) of \({\mathscr {C}}_{0}\), with \(\alpha = \mathrm {M},\mathrm {S}\), its component-wise notation will be expressed as \({\mathbb {D}}_{\alpha }=[{\mathbb {D}}_{\alpha }]_{abcd}\,{\varvec{i}}_{a}\otimes {\varvec{i}}_{b}\otimes {\varvec{i}}_{c}\otimes {\varvec{i}}_{d}\), i.e. by using square brackets for \({\mathbb {D}}_{\alpha }\). When necessary, such notation will be adopted also for vectors and second- and third-order tensors.

In light of the infinitesimal theory of continuum mechanics [28, 34], we describe the kinematics of the two phases constituting \({\mathscr {C}}\) in terms of the displacement fields

$$\begin{aligned} {\varvec{u}}_{\alpha }:={\varvec{u}}|_{\Omega _{\alpha }}:\Omega _{\alpha }\rightarrow {\mathcal {V}}, \quad \quad \alpha \in \{\mathrm {M},\mathrm {S}\}, \end{aligned}$$
(6)

and the displacement gradient fields

$$\begin{aligned} \nabla {\varvec{u}}_{\alpha }:\Omega _{\alpha }\rightarrow \mathrm {Lin}({\mathcal {V}}), \quad \quad \alpha \in \{\mathrm {M},\mathrm {S}\}. \end{aligned}$$
(7)

In (6), \({\varvec{u}}:{\mathscr {C}}_{0}\rightarrow {\mathcal {V}}\) is the displacement field of the composite material under study, so that \({\varvec{u}}_{\mathrm {M}}\) and \({\varvec{u}}_{\mathrm {S}}\) are the restrictions of \({\varvec{u}}\) to the matrix and to the family of sub-phases, respectively. Given \({\varvec{u}}_{\alpha }\), with \(\alpha \in \{\mathrm {M}, \mathrm {S}\}\), we introduce the infinitesimal Green–Lagrange strain tensor

$$\begin{aligned} {\varvec{\xi }}\left( {\varvec{u}}_\alpha \right) =\dfrac{\nabla {\varvec{u}}_{\alpha }+\left( \nabla {\varvec{u}}_{\alpha }\right) ^{\mathrm {T}}}{2}\in \mathrm {Sym}({\mathcal {V}}), \quad \quad \alpha \in \{\mathrm {M},\mathrm {S}\}, \end{aligned}$$
(8)

with \(\mathrm {Sym}({\mathcal {V}})\) being the subspace of \(\mathrm {Lin}({\mathcal {V}})\) consisting of all symmetric second order tensors. Note that \({\varvec{\xi }}\left( \,\cdot \,\right) \) is a tensor-valued differential operator that associates the vector field \({\varvec{u}}\) with the second-order tensor given by the symmetrised gradient of \({\varvec{u}}\).

We consider the case in which the composite \({\mathscr {C}}\) is subjected to the action of an imposed electric field \({\varvec{E}}:{\mathscr {C}}\rightarrow {\mathcal {V}}\), and we focus on the ideal situation in which the effects of possible magnetic fields can be disregarded. Within the picture described so far, and by neglecting inertial and mechanical body forces, the equations describing the point-wise equilibrium for the material composite \({\mathscr {C}}\) read

$$\begin{aligned} \nabla \,\cdot {\varvec{\sigma }}^{(\mathrm {t})}_{\mathrm {M}}&={\varvec{0}},&\mathrm {in}~ {\mathscr {C}}_{\mathrm {M}}, \end{aligned}$$
(9a)
$$\begin{aligned} \nabla \,\cdot {\varvec{\sigma }}^{(\mathrm {t})}_{\mathrm {S}}&={\varvec{0}},&\mathrm {in}~{\mathscr {C}}_{\mathrm {S}}, \end{aligned}$$
(9b)
$$\begin{aligned} {\varvec{\sigma }}^{(\mathrm {t})}_{\mathrm {M}}{\varvec{n}}_{0}&={\varvec{\sigma }}^{(\mathrm {t})}_{\mathrm {S}}{\varvec{n}}_{0},&\mathrm {on}~ \Gamma _{0}, \end{aligned}$$
(9c)
$$\begin{aligned} {\varvec{u}}_{\mathrm {M}}&={\varvec{u}}_{\mathrm {S}},&\mathrm {on}~ \Gamma _{0}, \end{aligned}$$
(9d)

with \({\varvec{n}}_{0}\) being the unit vector normal to the interface \(\Gamma _{0}\), and Eq. (9d) prescribing the continuity of the displacements on \(\Gamma _{0}\). In addition, suitable boundary conditions have to be specified on the boundary \(\partial {\mathscr {C}}_{0}\) of \({\mathscr {C}}_{0}\). In (9a)–(9c), and with \(\alpha \in \{\mathrm {M}, \mathrm {S}\}\), \({\varvec{\sigma }}^{(\mathrm {t})}_{\alpha }:{\mathscr {C}}_{0}\rightarrow \mathrm {Lin}({\mathcal {V}})\) is the total, generalised Cauchy stress tensor of the phase \({\mathscr {C}}_{\alpha }\), defined as the work-conjugate of the displacement gradient \(\nabla \,{\varvec{u}}_{\alpha }\).

Within the present framework, \({\varvec{\sigma }}^{(\mathrm {t})}_{\alpha }\) is the sum of two terms [22], i.e.

$$\begin{aligned} {\varvec{\sigma }}^{(\mathrm {t})}_{\alpha }={\varvec{\sigma }}^{(\mathrm {m})}_{\alpha }+{\varvec{\sigma }}^{(\mathrm {e})}_{\alpha }, \quad \alpha \in \{\mathrm {M},\mathrm {S}\}. \end{aligned}$$
(10)

The first summand, \({\varvec{\sigma }}^{(\mathrm {m})}_{\alpha }\), is the “classical” Cauchy stress tensor of linear elasticity, given by

$$\begin{aligned} {\varvec{\sigma }}^{(\mathrm {m})}_{\alpha }:={\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha }={\mathbb {C}}_{\alpha }:\varvec{{\varvec{\xi }}}\left( {\varvec{u}}_{\alpha }\right) , \quad \alpha =\{\mathrm {M},\mathrm {S}\}, \end{aligned}$$
(11)

and endowed with the symmetry property \({\varvec{\sigma }}^{(\mathrm {m})}_{\alpha }=({\varvec{\sigma }}^{(\mathrm {m})}_{\alpha })^{\mathrm {T}}\), and \({\mathbb {C}}_{\alpha }\in \mathrm {Lin}(\mathrm {Lin}({\mathcal {V}}))\) is the fourth-order elasticity tensor of \({\mathscr {C}}_{\alpha }\), which enjoys the major and the minor symmetries,

$$\begin{aligned} \left[ {\mathbb {C}}_{\alpha }\right] _{ijkl}=\left[ {\mathbb {C}}_{\alpha }\right] _{klij}, \quad \left[ {\mathbb {C}}_{\alpha }\right] _{ijkl}=\left[ {\mathbb {C}}_{\alpha }\right] _{jikl}, \quad \left[ {\mathbb {C}}_{\alpha }\right] _{ijkl}=\left[ {\mathbb {C}}_{\alpha }\right] _{ijlk}, \quad \quad \alpha \in \{\mathrm {M},\mathrm {S}\}. \end{aligned}$$
(12)

The second term, \({\varvec{\sigma }}^{(\mathrm {e})}_{\alpha }\), is a given tensor, here identified with the generalised Maxwell stress tensor, which, in the absence of magnetic contributions, can be proven to take on the expression (see Appendix A)

$$\begin{aligned} {\varvec{\sigma }}^{(\mathrm {e})}_{\alpha }=-\tfrac{1}{2}{\mathbb {B}}_{\alpha }^{\mathrm {T}}:({\varvec{E}}\otimes {\varvec{E}}), \end{aligned}$$
(13)

where the fourth-order tensor \({\mathbb {B}}_{\alpha }\) is referred to as electrostriction tensor and is required to enjoy the pair (minor) symmetries \([{\mathbb {B}}_{\alpha }]_{abcd}=[{\mathbb {B}}_{\alpha }]_{bacd}=[{\mathbb {B}}_{\alpha }]_{abdc}\). Such symmetry properties guarantee the symmetry of the generalised Maxwell stress tensor, \({\varvec{\sigma }}^{(\mathrm {e})}_{\alpha }\). The tensor \({\mathbb {B}}_{\alpha }\) introduced here is rather general, but it can be specialised to the case of interest according to the material symmetries of the \(\alpha \)th phase. For instance, if we assume the matrix and the sub-phases to be electrically isotropic, Eq. (13) can be recast as follows

$$\begin{aligned} {\varvec{\sigma }}^{(\mathrm {e})}_{\alpha }=(\epsilon _{0}+\epsilon _{0\alpha }-\epsilon _{2\alpha })\,{\varvec{E}}\otimes {\varvec{E}}-\tfrac{1}{2}\,(\epsilon _{0}+\epsilon _{0\alpha }-\epsilon _{1\alpha }-3\epsilon _{2\alpha }) \left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}={\mathbb {T}}_{\alpha }:({\varvec{E}}\otimes {\varvec{E}}), \end{aligned}$$
(14)

where \({\varvec{I}} \in \mathrm {Lin}({\mathcal {V}})\) is the second-order identity tensor, \(\epsilon _{0\alpha }, \epsilon _{1\alpha } \) and \(\epsilon _{2\alpha }\) are the electric permittivities of \({\mathscr {C}}_{\alpha }\), \(\epsilon _{0}\) is the permittivity of the vacuum and the fourth-order tensor \({\mathbb {T}}_{\alpha }\) is defined in Appendix A. In particular, \({\mathbb {T}}_{\alpha }\) is a fully symmetric fourth-order tensor and the identification \({\mathbb {B}}_{\alpha }={\mathbb {B}}_{\alpha }^{\mathrm {T}}=-2{\mathbb {T}}_{\alpha }\) holds true. The Reader is referred to Appendix A also for a derivation of the balance equations (9) and of the additive splitting of the total stress tensor (10). We remark that the superscripts “(m)” and “(e)”, given to \({\varvec{\sigma }}_{\alpha }^{\mathrm {(m)}}\) and \({\varvec{\sigma }}_{\alpha }^{\mathrm {(e)}}\), stand for “mechanical” and “electric” part of the total stress tensor, \({\varvec{\sigma }}_{\alpha }^{\mathrm {(t)}}\), respectively. We highlight that the expression (14) of the generalised Maxwell stress tensor descends directly from the coupling between the mechanical and the electrical behaviour of the composite. As explained in Appendix A, this coupling has been established by means of a suitable Lagrangian density function, obtained by taking inspiration from [22]. In particular, for addressing the problem at hand, we have considered a linearised version of the constitutive framework outlined in [22]. In doing this, we have linearised a Lagrangian density function featuring an isotropic fourth-order permittivity tensor, which complies with the material symmetries of the composite studied in this work. Finally, we remark that we have referred to \({\varvec{\sigma }}^{\mathrm {(e)}}\) as “generalised Maxwell stress tensor” because the “true” Maxwell stress tensor used in [22] features solely the dielectric constant of the vacuum, \(\varepsilon _{0}\).

With the introduction of \({\varvec{\sigma }}^{(\mathrm {m})}_{\alpha }\) and of \({\varvec{\sigma }}^{(\mathrm {e})}_{\alpha }\), the system of equations (9) can be recast as follows

$$\begin{aligned} \nabla \,\cdot \left( {\varvec{\sigma }}^{(\mathrm {m})}_{\mathrm {M}} + {\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}\right)&={\varvec{0}},&\quad \mathrm {in}~{\mathscr {C}}_{\mathrm {M}}, \end{aligned}$$
(15a)
$$\begin{aligned} \nabla \,\cdot \left( {\varvec{\sigma }}^{(\mathrm {m})}_{\mathrm {S}} + {\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}\right)&={\varvec{0}},&\quad \mathrm {in}~{\mathscr {C}}_{\mathrm {S}}, \end{aligned}$$
(15b)
$$\begin{aligned} \left( {\varvec{\sigma }}^{(\mathrm {m})}_{\mathrm {M}} + {\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}\right) {\varvec{n}}_{0}&=\left( {\varvec{\sigma }}^{(\mathrm {m})}_{\mathrm {S}} + {\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}\right) {\varvec{n}}_{0},&\quad \mathrm {on}~\Gamma _{0}, \end{aligned}$$
(15c)
$$\begin{aligned} {\varvec{u}}_{\mathrm {M}}&={\varvec{u}}_{\mathrm {S}},&\quad \mathrm {on}~\Gamma _{0}. \end{aligned}$$
(15d)

Remark 2

With reference to Eq. (14), we remark that the “classical” case of isotropic material characterised by one electric permittivity only can be retrieved by setting \(\epsilon _{1\alpha }\) and \(\epsilon _{2\alpha }\) equal to zero. By doing so, the permittivity tensor becomes proportional to the identity tensor, as expected for linear materials or for non-linear models of “unconditionally isotropic” materials (to use the jargon of [2]). However, our model is obtained by linearising a non-linear constitutive expression that relates, through the Representation Theorem of tensor-valued functions, the dielectric tensor of the material, denoted by \({\hat{\varvec{\epsilon }}}_{\alpha }^{\mathrm {(c)}}\), with the deformation (see Appendix A). Although such a constitutive description characterises an isotropic material, it accounts for the “variability” of the material’s isotropy with the deformation and, consequently, it requires the presence of suitable coefficients that multiply the left Cauchy–Green tensor, i.e., \({\varvec{b}}_{\alpha }\), and \({\varvec{b}}^{2}_{\alpha }\) in the tensor representation of \({\hat{\varvec{\epsilon }}}_{\alpha }^{\mathrm {(c)}}\) (see Appendix A). Such coefficients are, in fact, \(\epsilon _{1\alpha }\) and \(\epsilon _{2\alpha }\) and they themselves, or combinations of them, should be determined experimentally. In the linear case, however, it can be inferred from Eq. (14) that it is possible to introduce the auxiliary coefficients \({\tilde{\epsilon }}_{1\alpha }=\epsilon _{0\alpha }-\epsilon _{2\alpha }\) and \({\tilde{\epsilon }}_{2\alpha }=\epsilon _{0\alpha }-\epsilon _{1\alpha }-3\epsilon _{2\alpha }\) and write

$$\begin{aligned} {\varvec{\sigma }}_{\alpha }^{\mathrm {(e)}}&=(\epsilon _{0}+{\tilde{\epsilon }}_{1\alpha }){\varvec{E}}\otimes {\varvec{E}}-\tfrac{1}{2}(\epsilon _{0}+{\tilde{\epsilon }}_{2\alpha })\Vert {\varvec{E}} \Vert ^{2}{\varvec{I}}\nonumber \\&=(\epsilon _{0}+{\tilde{\epsilon }}_{1\alpha }){\mathbb {I}}:{\varvec{E}}\otimes {\varvec{E}}-\tfrac{3}{2}(\epsilon _{0}+{\tilde{\epsilon }}_{2\alpha }){\mathbb {K}}:{\varvec{E}}\otimes {\varvec{E}}\nonumber \\&=[(\epsilon _{0}+{\tilde{\epsilon }}_{1\alpha }){\mathbb {I}}-\tfrac{3}{2}(\epsilon _{0}+{\tilde{\epsilon }}_{2\alpha }){\mathbb {K}}]:{\varvec{E}}\otimes {\varvec{E}}, \end{aligned}$$
(16)

where the fourth-order tensors \({\mathbb {I}}\) and \({\mathbb {K}}\) are defined in Appendix A. This way, and by introducing the fourth-order tensor \({\mathbb {M}}={\mathbb {I}}-{\mathbb {K}}\) [73], which extracts the deviatoric part of a symmetric second-order tensor, \({\varvec{\sigma }}_{\alpha }^{\mathrm {(e)}}\) is related to \({\varvec{E}}\otimes {\varvec{E}}\) through the isotropic fourth-order tensor

$$\begin{aligned} {\mathbb {T}}_{\alpha }&=(\epsilon _{0}+{\tilde{\epsilon }}_{1\alpha }){\mathbb {I}}-\tfrac{3}{2}(\epsilon _{0}+{\tilde{\epsilon }}_{2\alpha }){\mathbb {K}}\nonumber \\&=\left( -\tfrac{1}{2}\epsilon _{0}+{\tilde{\epsilon }}_{1\alpha }-\tfrac{3}{2}{\tilde{\epsilon }}_{2\alpha }\right) {\mathbb {K}}+(\epsilon _{0}+{\tilde{\epsilon }}_{1\alpha }){\mathbb {M}}, \end{aligned}$$
(17)

which, in general, apart from the permittivity of the vacuum, \(\epsilon _{0}\), requires the two independent coefficients \({\tilde{\epsilon }}_{1\alpha }\) and \({\tilde{\epsilon }}_{2\alpha }\) to be represented. Of course, if we set again \(\epsilon _{1\alpha }=0\) and \(\epsilon _{2\alpha }=0\), the resulting expression for \({\mathbb {T}}_{\alpha }\) becomes

$$\begin{aligned} {\mathbb {T}}_{\alpha }&=(\epsilon _{0}+\epsilon _{0\alpha }){\mathbb {I}}-\tfrac{3}{2}(\epsilon _{0}+\epsilon _{0\alpha }){\mathbb {K}}\nonumber \\&=\left( -\tfrac{1}{2}\epsilon _{0}-\tfrac{1}{2}\epsilon _{0\alpha }\right) {\mathbb {K}}+(\epsilon _{0}+\epsilon _{0\alpha }){\mathbb {M}}, \end{aligned}$$
(18)

which amounts to “rescale” the vacuum by means of the material’s permittivity \(\epsilon _{0\alpha }\). Clearly, this is a limit case, since an isotropic fourth-order tensor, like \({\mathbb {T}}_{\alpha }\), necessitates, in general, two coefficients to be represented. We conclude this remark by mentioning that, following [73], the expressions of \({\mathbb {T}}_{\alpha }\) featuring in the far right-hand-side of (17) and (18) are introduced because they provide the representations of isotropic fourth-order tensors in their natural tensor basis \(\{{\mathbb {K}},{\mathbb {M}}\}\) consisting of \({\mathbb {K}}\), which extracts the spherical (or volumetric) part of a symmetric second-order tensor, and \({\mathbb {M}}\), which extracts the isochoric part of a symmetric second-order tensor (for details, the Reader is referred to [28,29,30, 73]).

3 Multiscale modelling

In the following, the system of equations (15a)–(15d) is studied by means of a technique known as asymptotic homogenisation, which allows to obtain an effective description of the electro-mechanical properties of the composite as a whole. This can be achieved by accounting for the distribution of the phases \({\mathscr {C}}_{\mathrm {M}}\) and \({\mathscr {C}}_{\mathrm {S}}\) within the composite, specifying the elasticity and electric permittivity tensors of \({\mathscr {C}}_{\mathrm {M}}\) and \({\mathscr {C}}_{\mathrm {S}}\), resolving the interactions exchanged by \({\mathscr {C}}_{\mathrm {M}}\) and \({\mathscr {C}}_{\mathrm {S}}\) across \(\Gamma _{0}\), and evaluating \({\varvec{u}}_{\mathrm {M}}\) and \({\varvec{u}}_{\mathrm {S}}\) in a scale-dependent manner. The distribution of the phases in the composite constituents is strongly related to the geometric configuration of the micro-structure. In this regard, we will sketch out the way in which the geometry and the topology of the micro-structure comes into play in the following Sect. 4.

The employment of the asymptotic homogenisation technique requires a sharp separation of the length scales that characterise the hierarchical levels with which the composite’s internal structure is organised. For the problem at hand, we stipulate the existence of two relevant length scales. More specifically, we denote by \(\ell \) the representative length of the internal structure of the composite \({\mathscr {C}}\), which characterises, for instance, the distance between two neighbouring sub-phases, and we let L describe the characteristic size of the domain occupied by the composite as a whole. The condition for which the two considered scales are well separated can be formalised by requiring that the ratio between the fine length scale, \(\ell \), and the coarse one, L, is much smaller than unity, i.e.

$$\begin{aligned} \varepsilon :=\dfrac{\ell }{L}\ll 1. \end{aligned}$$
(19)

The strictly positive, dimensionless quantity \(\varepsilon \) is called scaling parameter and is representative of the two-scale heterogeneity of the composite. We recall that, as is customary in several works based on the asymptotic homogenisation technique, the definition (19) of \(\varepsilon \) applies to the case in which L and \(\ell \) can be regarded as constant parameters of the composite’s structural organisation. In a more general setting, however, when mechanical processes other than those studied in this work are accounted for, it could be necessary to specify \(\varepsilon \) as a function of the material points of the composite and, in principle, also of time. As noticed in [66], this may occur, for instance, when the composite under study exhibits large deformations or when it undergoes remodelling or growth. Nevertheless, if the scaling parameter is bounded from above by a constant upper bound, and if such bound is much smaller than unity, the asymptotic homogenisation technique can still be adopted.

For the class of problems addressed in this work, it is important to establish a direct relationship between a given physical quantity F(x) and the composite’s characteristic length scales L and \(\ell \). This can be achieved by rewriting F(x) as \(F(x)={\check{F}}(x;L,\ell )\), with L and \(\ell \) regarded as parameters. In fact, this way of writing is very general and applies also when the ratio \(\ell /L\) is not much smaller than unity. However, as clearly shown in [60], a particularly relevant case occurs when the dependence of \({\check{F}}\) on L and \(\ell \) is such that the quantity of interest is, in fact, a function of the ratios x/L and \(x/\ell \), i.e. when it holds true that \(F(x)={\check{F}}(x;L,\ell )={\tilde{F}}(x/L,x/\ell )\), thereby making it explicit that the considered physical quantity depends on two different resolutions of the spatial variable x. To let \(\varepsilon \) feature explicitly, it suffices to write

$$\begin{aligned} F(x)={\check{F}}(x;L,\ell )={\tilde{F}}\left( \frac{x}{L},\frac{x}{\ell }\right) ={\tilde{F}}\left( \frac{x}{L},\frac{x}{L\frac{\ell }{L}}\right) ={\tilde{F}}\left( {\tilde{x}},\frac{{\tilde{x}}}{\varepsilon }\right) , \end{aligned}$$
(20)

where \({\tilde{x}}=x/L\) is the spatial variable made dimensionless, and the inhomogeneities arising from the micro-structure are kept track of through the “zoom lens” provided by \(\varepsilon \). Finally, by introducing an auxiliary function \(\phi _{\varepsilon }\) defined by

$$\begin{aligned} \phi _{\varepsilon }({\tilde{x}})={\tilde{x}}/\varepsilon \equiv {\tilde{y}}, \end{aligned}$$
(21)

Eq. (20) becomes

$$\begin{aligned} F(x)={\tilde{F}}\left( {\tilde{x}},\frac{{\tilde{x}}}{\varepsilon }\right) ={\tilde{F}}({\tilde{x}},\phi _{\varepsilon }({\tilde{x}}))={\tilde{F}}({\tilde{x}},{\tilde{y}}), \end{aligned}$$
(22)

so that F is reformulated as a function \({\tilde{F}}\) of two variables, \({\tilde{x}}\) and \({\tilde{y}}\), that are independent of each other. In fact, the two arguments of \({\tilde{F}}\) are representative of two different length scales, and the advantage of the formulation (22) is that it highlights the two-scale nature of the considered physical quantity. Indeed, \({\tilde{x}}\) is referred to as macro-scale or slow variable, because it describes the variations of \({\tilde{F}}\) resolved at the scale L, whereas \({\tilde{y}}\) is said to be the micro-scale or fast variable, as it captures the high-frequency variations of \({\tilde{F}}\), associated with \(\ell \). The two-scale representation (20) should be performed under suitable hypotheses of regularity of the function \({\tilde{F}}\), such as local boundedness with respect to the second slot, i.e. the one “occupied” by the micro-scale variable. This question and the role played by this kind of regularity assumption are analysed, for example, in [60]. In spite of their importance, such mathematical aspects are not discussed in this work. However, the Reader interested in these technical issues is referred to the literature of asymptotic homogenisation [4, 17, 60].

For the problem under study, we assume that all the material properties considered in the following can be expressed ab initio as functions on the rescaled spatial variable \({\tilde{x}}=x/L\) and, by virtue of Eq. (22), we enforce the rewriting

$$\begin{aligned} {\mathbb {C}}_{\alpha }({\tilde{x}})&={\tilde{\mathbb {C}}}_{\alpha }\left( {\tilde{x}},{\tilde{y}}\right) , \quad \alpha =\{\mathrm {M},\mathrm {S}\}, \end{aligned}$$
(23a)
$$\begin{aligned} \epsilon _{\alpha }({\tilde{x}})&={\tilde{\epsilon }}_{\alpha }\left( {\tilde{x}},{\tilde{y}}\right) , \quad \,\,\alpha =\{\mathrm {M},\mathrm {S}\}, \end{aligned}$$
(23b)
$$\begin{aligned} {\varvec{E}}({\tilde{x}})&={\tilde{\varvec{E}}}\,({\tilde{x}},{\tilde{y}}). \end{aligned}$$
(23c)

The same representation is applied to the displacements of the two constituents of the composite. This time, however, we write

$$\begin{aligned} {\varvec{u}}_{\alpha }({\tilde{x}})&=L\,{\tilde{\varvec{u}}}_{\alpha }\left( {\tilde{x}},{\tilde{y}}\right) , \quad \alpha =\{\mathrm {M},\mathrm {S}\}, \end{aligned}$$
(24)

which amounts to imposing that \({\varvec{u}}_{\alpha }\) scales with L, with \({\tilde{\varvec{u}}}_{\alpha }\) being a dimensionless displacement.

Moreover, upon using the identities \(f({\tilde{x}})={\tilde{f}}({\tilde{x}},\phi _{\varepsilon }({\tilde{x}}))={\tilde{f}}({\tilde{x}},{\tilde{y}})\), and applying the chain rule of differentiation of composed maps, we can write the gradient of f as

$$\begin{aligned} \nabla _{{\tilde{x}}}f=\nabla _{{\tilde{x}}}{\tilde{f}}+\varepsilon ^{-1}\nabla _{{\tilde{y}}}{\tilde{f}}, \end{aligned}$$
(25)

where \(\nabla _{{\tilde{x}}}\) and \(\nabla _{{\tilde{y}}}\) denote the gradients of \({\tilde{f}}\) in the slow variable, \({\tilde{x}}\), and in the fast variable, \({\tilde{y}}\), respectively. For instance, the infinitesimal Green–Lagrange strain tensor, defined in Eq. (8), reads

$$\begin{aligned} {\varvec{\xi }}\left( {\varvec{u}}_{\alpha }\right) =L^{-1}{\varvec{\xi }}_{{\tilde{x}}}\left( {\varvec{u}}_{\alpha }\right) ={\varvec{\xi }}_{{\tilde{x}}}\left( {\tilde{\varvec{u}}}_{\alpha }\right) +\varepsilon ^{-1}{\varvec{\xi }}_{{\tilde{y}}}\left( {\tilde{\varvec{u}}}_{\alpha }\right) , \quad \alpha \in \{\mathrm {M},\mathrm {S}\}, \end{aligned}$$
(26)

with

$$\begin{aligned} {\varvec{\xi }}_{{\tilde{x}}}({\tilde{\varvec{u}}}_{\alpha })=\dfrac{\nabla _{{\tilde{x}}}{\tilde{\varvec{u}}}_{\alpha }+\left( \nabla _{{\tilde{x}}}{\tilde{\varvec{u}}}_{\alpha }\right) ^{\mathrm {T}}}{2}, \quad {\varvec{\xi }}_{{\tilde{y}}}({\tilde{\varvec{u}}}_{\alpha })=\dfrac{\nabla _{{\tilde{y}}}{\tilde{\varvec{u}}}_{\alpha }+\left( \nabla _{{\tilde{y}}}{\tilde{\varvec{u}}}_{\alpha }\right) ^{\mathrm {T}}}{2}, \quad \alpha \in \{\mathrm {M},\mathrm {S}\}. \end{aligned}$$
(27)

Note that, whereas \({\varvec{\xi }}({\varvec{u}}_{\alpha })\) is by definition dimensionless, \({\varvec{\xi }}_{{\tilde{x}}}\left( {\varvec{u}}_{\alpha }\right) \) has dimensions of length, while \({\varvec{\xi }}_{{\tilde{x}}}\left( {\tilde{\varvec{u}}}_{\alpha }\right) \) and \({\varvec{\xi }}_{{\tilde{y}}}\left( {\tilde{\varvec{u}}}_{\alpha }\right) \) are dimensionless by virtue of (24).

Remark 3

(The role of locally unbounded generalised forces)

We recall that a vector-valued function \({\varvec{f}}(x)={\tilde{\varvec{f}}}({\tilde{x}},\frac{{\tilde{x}}}{\varepsilon })={\tilde{\varvec{f}}}({\tilde{x}},{\tilde{y}})\) is said to be locally bounded if, for every \({\tilde{x}}\), its two-scale representation satisfies

$$\begin{aligned} \lim _{\left\Vert {\tilde{y}}\right\Vert \rightarrow +\infty }\left\Vert {\tilde{\varvec{f}}}({\tilde{x}},{\tilde{y}})\right\Vert <+\infty \end{aligned}$$
(28)

(see also [64, 70]). By noticing that the macro-scale variable x spans a bounded domain of \({\mathbb {R}}^{3}\), that is the representation of a bounded domain of the three-dimensional Euclidean space \({\mathscr {S}}\), and by virtue of Eq. (21), the definition of locally bounded functions given in (28) can be recast in the following equivalent way

$$\begin{aligned} \lim _{\varepsilon \rightarrow 0}\left\Vert {\tilde{\varvec{f}}}\left( {\tilde{x}},\frac{{\tilde{x}}}{\varepsilon }\right) \right\Vert <+\infty . \end{aligned}$$
(29)

Whenever the vector-valued function \({\varvec{f}}\) does not comply with Eq. (28), or equivalently with Eq. (29), we say that \({\varvec{f}}\) is locally unbounded.

By focusing on the system of PDEs (9), together with the expression for the total stress tensor (10), we can identify such a vector-valued function with the divergence of the given generalised Maxwell tensor \({\varvec{\sigma }}_{\alpha }^{(\mathrm {e})}\), i.e.

$$\begin{aligned} \tilde{\varvec{f}}_{\alpha }\left( {\tilde{x}},\frac{{\tilde{x}}}{\varepsilon }\right) :=\nabla _{{\tilde{x}}} \cdot {\varvec{\sigma }}_{\alpha }^{(\mathrm {e})}. \end{aligned}$$
(30)

For the case studied in [64] and in this work, Eq. (30) yields a locally unbounded force and the differential problem given by (910) is to be considered as being driven by a locally unbounded generalised force \(\displaystyle \tilde{\varvec{f}}_{\alpha }\left( {\tilde{x}},\frac{{\tilde{x}}}{\varepsilon }\right) \)Footnote 2, which is not fulfilling condition (29). The systems of PDEs (15) is nonetheless shown to be homogenisable [64] by embracing micro-scale periodicity, as done in the remainder of this work.

In spite of the formalism adopted so far, following [3, 60], from here on we simplify the notation by dropping the tilde, thereby indicating \({\tilde{x}}\) and \({\tilde{y}}\) with x and y, respectively, \({\tilde{\varvec{u}}}_{\alpha }\) with \({\varvec{u}}_{\alpha }\), and \({\tilde{\mathbb {C}}}_{\alpha }\), \({\tilde{\epsilon }}_{\alpha }\) and \({\tilde{\varvec{E}}}_{\alpha }\) with \({\mathbb {C}}_{\alpha }\), \(\epsilon _{\alpha }\) and \({\varvec{E}}_{\alpha }\), respectively. Furthermore, for the two-scale representation of the displacement fields (24), we admit the existence of a formal (two-scale) expansion in power series of \(\varepsilon \), namely,

$$\begin{aligned} {\varvec{u}}_{{\alpha }}(x,y)\equiv {\varvec{u}}_{{\alpha }}^{\varepsilon }(x,y)=\sum _{n=0}^{+\infty }{\varvec{u}}_{{\alpha }}^{(n)}(x,y)\,\varepsilon ^n, \quad \alpha \in \{\mathrm {M},\mathrm {S}\}. \end{aligned}$$
(31)

In this work, we do not fuss over the technicalities and the regularity hypothesis needed to justify in a rigorous way the method of power series expansion (31). The Reader, however, can take a deeper look into such issues by referring to the huge literature in the field (see, for example, [4, 17] and references therein).

4 Geometry and topology of the micro-structure

We assume that all the fields of interest, i.e. \({\varvec{u}}_{\alpha }\), with \(\alpha \in \{\mathrm {M}, \mathrm {S}\}\), and the material parameters characterising the matrix and the sub-phases of the composite are periodic with respect to the micro-scale variable y. From a technical point of view, the micro-scale periodicity implies that y can be assumed to span an open subset of \({\mathbb {R}}^3\), denoted with \(\Omega \) and referred to as the elementary cell. Analogously to what we have previously done, we let \(\Omega _{\mathrm {M}}\) and \(\Omega _{\mathrm {S}}\) be the portions of \(\Omega \) occupied by the matrix and the sub-phases, respectively. Here, both sets are taken to be open and, depending on the topology of the elementary cell, there can be cases in which \(\Omega _{\mathrm {M}}\) and \(\Omega _{\mathrm {S}}\) are either both connected or both disconnected, cases in which \(\Omega _{\mathrm {M}}\) is disconnected and \(\Omega _{\mathrm {S}}\) is connected, and cases in which \(\Omega _{\mathrm {M}}\) is connected and \(\Omega _{\mathrm {S}}\) is disconnected. At any rate, the interface between \(\Omega _{\mathrm {M}}\) and \(\Omega _{\mathrm {S}}\) is the surface \(\Gamma \) identified with the intersections of the boundaries of \(\Omega _{\mathrm {M}}\) and \(\Omega _{\mathrm {S}}\), i.e. \(\Gamma = \partial \Omega _{\mathrm {M}}\cap \partial \Omega _{\mathrm {S}}\).

In the forthcoming sections, we shall concentrate exclusively on the three cell configurations described in the sequel, which refer to the cases in which the sub-phases contained in the cell are either an inclusion or a fibre or a stratum. In all these situations, our results apply because the elementary cell can be decomposed in regular sub-domains with Lipschitz boundaries.

Inclusion-like sub-phases A sub-phase is said to be an inclusion when it is completely embedded in \(\Omega \). More specifically, we assume that \(\Omega _{\mathrm {M}}\) and \(\Omega _{\mathrm {S}}\) are both connected and, in particular, that the inclusion is contained in the elementary cell in a way such that the closure of \(\Omega _{\mathrm {S}}\) is contained in \(\Omega \), i.e. \({\bar{\Omega }}_{\mathrm {S}}\subset \Omega \). Accordingly, the interface \(\Gamma \) between \(\Omega _{\mathrm {M}}\) and \(\Omega _{\mathrm {S}}\) consists of sole inner points of \(\Omega \). Moreover, the boundary of \(\Omega _{\mathrm {M}}\) is given by \(\partial \Omega _{\mathrm {M}} = \partial _{\mathrm {o}} \Omega _{\mathrm {M}}\cup \partial _{\mathrm {i}}\Omega _{\mathrm {M}}\), where \(\partial _{\mathrm {o}}\Omega _{\mathrm {M}}\) and \(\partial _{\mathrm {i}}\Omega _{\mathrm {M}}\) are the outer and the inner boundaries of \(\Omega _{\mathrm {M}}\), respectively, and are such that \(\partial _{\mathrm {o}} \Omega _{\mathrm {M}}\cap \partial _{\mathrm {i}}\Omega _{\mathrm {M}}=\emptyset \), \(\partial _{\mathrm {o}}\Omega _{\mathrm {M}}\equiv \partial \Omega \), and \(\partial _{\mathrm {i}}\Omega _{\mathrm {M}}=\Gamma \). It follows from this topological description of the elementary cell that the closure of \(\Omega \) is given by \({\bar{\Omega }} ={\bar{\Omega }}_{\mathrm {M}}\cup \Omega _{\mathrm {S}}\), with \({\bar{\Omega }}_{\mathrm {M}}\cap \Omega _{\mathrm {S}}=\emptyset \) and \(\Omega _{\mathrm {M}}\cap {\bar{\Omega }}_{\mathrm {S}}=\emptyset \), and that the interface reads \(\Gamma =\partial \Omega _{\mathrm {M}}\cap \partial \Omega _{\mathrm {S}}={\bar{\Omega }}_{\mathrm {M}}\cap {\bar{\Omega }}_{\mathrm {S}}\) (see the left panel of Fig. 1).

Fibre-like sub-phases A sub-phase is said to be a fibre when \(\Omega _{\mathrm {S}}\) is contained in \(\Omega \) (recall that \(\Omega _{\mathrm {S}}\) is an open set and it, thus, coincides with the set of its internal points) in such a way that \(\Omega _{\mathrm {M}}\) remains connected and the intersection \({\bar{\Omega }}_{\mathrm {S}}\cap \partial \Omega \) is non-empty and respectful of the periodicity of the composite. From a more geometric perspective, we can say that a sub-phase is a fibre if its boundary, \(\partial \Omega _{\mathrm {S}}\), features two disjoint sets, denoted by \((\partial \Omega _{\mathrm {S}})_{\mathrm {en}}\) and \((\partial \Omega _{\mathrm {S}})_{\mathrm {ex}}\), such that \({\bar{\Omega }}_{\mathrm {S}}\cap \partial \Omega =(\partial \Omega _{\mathrm {S}})_{\mathrm {en}}\sqcup (\partial \Omega _{\mathrm {S}})_{\mathrm {ex}}\) (see the right panel of Fig. 1 and the right panel of Fig. 2). Here, the subscripts “en” and “ex” indicate the portions of the fibre’s boundary corresponding to the “entry side” and the “exit side” of the elementary cell.

Stratum-like sub-phases A sub-phase is said to be a stratum when \(\Omega _{\mathrm {S}}\) is contained in \(\Omega \) in such a way that \(\Omega _{\mathrm {M}}\) is a disconnected set and the periodicity of the composite is maintained. More precisely, we consider an elementary cell in which the portion occupied by the matrix, \(\Omega _{\mathrm {M}}\), is the union of a (natural) number P of sets, \(P\ge 2\), separated by \(\Omega _{\mathrm {S}}\). Thus, \(\Omega _{\mathrm {M}}\) can be written as \(\Omega _{\mathrm {M}}=\cup _{q=1}^{P}\Omega _{q\mathrm {M}}\) and its boundary, \(\partial \Omega _{\mathrm {M}}\), is given by \(\partial \Omega _{\mathrm {M}}=\cup _{q=1}^{P}\partial \Omega _{q\mathrm {M}}\), whereas the interface between the stratum and the matrix is identified with \(\Gamma = \partial \Omega _{\mathrm {M}}\cap \partial \Omega _{\mathrm {S}}\). Note that, in the just depicted situation, it holds that \(\Gamma \ne {\bar{\Omega }}_{\mathrm {M}}\cap {\bar{\Omega }}_{\mathrm {S}}\) because the set \({\bar{\Omega }}_{\mathrm {M}}\cap {\bar{\Omega }}_{\mathrm {S}}\) contains portions of the boundary of \(\Omega \) that are not contained in the interface, these portions being given by the inner points of \(\partial \Omega \cap \partial \Omega _{\mathrm {S}}\) (see the right panel of Fig. 1). We emphasise that the distinction between a fibre and a stratum is characterised by the fact that a stratum makes the matrix \(\Omega _{\mathrm {M}}\) disconnected, whereas a fibre does not break the connectedness of \(\Omega _{\mathrm {M}}\). Therefore, what distinguishes a fibre from a stratum is a topological, rather than geometrical, fact. Indeed, considering for instance a sub-phase shaped as a thin parallelepiped, the sub-phase is a fibre even if the length of the long edges of its transversal section is comparable with the length of the edges of its longitudinal section, because such configuration makes the matrix connected (see Fig. 2).

Fig. 1
figure 1

Two examples of elementary cell, representative of the micro-structural topological configurations discussed in this work. Left panel: Elementary cell containing an inclusion. Right panel: Elementary cell containing a fibre

To obtain a macro-scale description of the composite, it is useful to define the following averages for a generic quantity of interest, F:

$$\begin{aligned} \left\langle F \right\rangle (x):=\dfrac{1}{|\Omega |}\int \limits _{\Omega } F(x,y)\,\mathrm {d}y, \quad \quad \left\langle F \right\rangle _{\Omega _{\alpha }}(x):=\dfrac{1}{|\Omega |}\int \limits _{\Omega _{\alpha }} F(x,y)\,\mathrm {d}y, \end{aligned}$$
(32)

where \(|\Omega |\) is the size of \(\Omega \). Note that \(\left\langle F \right\rangle \) is the average of F over the whole elementary cell, \(\Omega \), whereas \(\left\langle F \right\rangle _{\Omega _{\alpha }}\) is the average over the portion \(\Omega _{\alpha }\) of \(\Omega \).

Fig. 2
figure 2

Two examples of elementary cell, representative of the micro-structural topological configurations discussed in this work. Left panel: Elementary cell containing a stratum. Right panel: Elementary cell containing a fibre

In light of the hypothesis of y-periodicity, a remark is in order. This hypothesis and, consequently, the definitions (32) rely on the fact that the elementary cell is supposed to be representative of the local structure of the composite. In other words, the topological and geometric information on the composite’s micro-structure is kept track of, even though only a portion of it, i.e. the elementary cell, is observed. In general, the choice of the elementary cell may vary with respect to x, and hence, the integrals defined in (32) may depend on the way in which the elementary cell is chosen, see e.g. [10, 21, 37, 58, 59].

In the following, for the sake of simplicity, we adopt a procedure followed also in other works (see, for example, [64]), in which the asymptotic homogenisation technique is employed under the assumption of macroscopic uniformity, which permits to select the elementary cell once for all, i.e. independently of x. The main consequence of this assumption is that the “macro-scale divergence operator” and the average operator over the elementary cell \(\Omega \), defined in (32), commute, i.e.

$$\begin{aligned} \nabla _{x}\cdot \left\langle {\varvec{v}}\right\rangle =\left\langle \nabla _{x}\cdot {\varvec{v}}\right\rangle , \end{aligned}$$
(33)

for every suitable vector field \({\varvec{v}}\in {\mathcal {V}}\). Indeed, the macroscopic uniformity hypothesis permits to write

$$\begin{aligned} \left[ \nabla _{x}\cdot \int \limits _{\Omega } {\varvec{v}}(\cdot ,y)\,\mathrm {d}y\right] (x) =\int \limits _{\Omega } \left( \nabla _{x}\cdot {\varvec{v}}\right) (x,y)\,\mathrm {d}y. \end{aligned}$$
(34)

Hereafter, for the sake of a lighter notation, we will omit the explicit indication of the independent variables and of the symbol of the differential in the writing of volume and surface integrals.

5 The macro-scale model

By making use of the asymptotic homogenisation technique introduced so far, and by enforcing the technical assumptions discussed in Sects. 3 and 4, we rephrase the boundary value problem with interface conditions (15c) in the following two-scale fashion,

$$\begin{aligned} \left( \nabla _x+\varepsilon ^{-1}\nabla _y \right) \,\cdot \left( {\varvec{\sigma }}^{(\mathrm {m})}_{\mathrm {M}} + {\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}\right)&={\varvec{0}},&\mathrm {in}~{\Omega }_{\mathrm {M}}, \end{aligned}$$
(35a)
$$\begin{aligned} \left( \nabla _x+\varepsilon ^{-1}\nabla _y \right) \,\cdot \left( {\varvec{\sigma }}^{(\mathrm {m})}_{\mathrm {S}} + {\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}\right)&={\varvec{0}},&\mathrm {in}~{\Omega }_{\mathrm {S}}, \end{aligned}$$
(35b)
$$\begin{aligned} \left( {\varvec{\sigma }}^{(\mathrm {m})}_{\mathrm {M}} + {\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}\right) {\varvec{n}}&=\left( {\varvec{\sigma }}^{(\mathrm {m})}_{\mathrm {S}} + {\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}\right) {\varvec{n}},&\mathrm {on}~\Gamma , \end{aligned}$$
(35c)
$$\begin{aligned} {\varvec{u}}^{\varepsilon }_{\mathrm {M}}&={\varvec{u}}^{\varepsilon }_{\mathrm {S}},&\mathrm {on}~\Gamma , \end{aligned}$$
(35d)

with \({\varvec{n}}\) being the unit vector normal to the interface \(\Gamma \). In particular, we assume that \({\varvec{n}}\) points from the matrix to the sub-phases. It is worth noticing that, with the procedure outlined in the previous sections, we have transformed the problem (15) in such a way that the unknowns are not the two displacement fields \({\varvec{u}}_{\mathrm {M}}\) and \({\varvec{u}}_{\mathrm {S}}\), but the coefficients of the power series expansion (31), which are, in principle, infinite. Moreover, we have highlighted how the fields of the theory depend on the micro-scale variable y. Although such a formulation seems to complicate the solution of the problem, this is, in fact, not the case. Indeed, since our purpose here is to derive an effective description of the composite under study, we will focus only on the leading orders of the power series (31), \({\varvec{u}}_{{\alpha }}^{(0)}\), with \(\alpha \in \{\mathrm {M}, \mathrm {S}\}\). To this end, we will derive a boundary value problem for each order of the power series that is of interest for the present study. Such a boundary value problem will encode all the pieces of information arising from the fine scale, expressed in terms of heterogeneities, geometry and coupling between the mechanical behaviour of the two phases of the composite \({\mathscr {C}}\), with the latter manifesting themselves through suitable effective coefficients [4].

Now, let us proceed by substituting the definition (27) and the constitutive equation (11) into (35). By multiplying the differential conditions (35a), and (35b) by \(\varepsilon ^{2}\) and (35c) by \(\varepsilon \), we obtain the following system of partial differential equations for the linear elastic electro-active composite \({\mathscr {C}}\)

$$\begin{aligned} \left( \varepsilon ^2\nabla _x+\varepsilon \nabla _y \right) \,\cdot \left( {\mathbb {C}}_\mathrm {M}:{\varvec{\xi }}_{x}\left( {\varvec{u}}^{\varepsilon }_{\mathrm {M}}\right) +\varepsilon ^{-1}{\mathbb {C}}_\mathrm {M}:{\varvec{\xi }}_{y}\left( {\varvec{u}}^{\varepsilon }_{\mathrm {M}}\right) +{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}\right)&={\varvec{0}},&\mathrm {in}~{\Omega }_{\mathrm {M}}, \end{aligned}$$
(36a)
$$\begin{aligned} \left( \varepsilon ^2\nabla _x+\varepsilon \nabla _y \right) \,\cdot \left( {\mathbb {C}}_\mathrm {S}:{\varvec{\xi }}_{x}\left( {\varvec{u}}^{\varepsilon }_{\mathrm {S}}\right) +\varepsilon ^{-1}{\mathbb {C}}_\mathrm {S}:{\varvec{\xi }}_{y}\left( {\varvec{u}}^{\varepsilon }_{\mathrm {S}}\right) +{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}\right)&={\varvec{0}},&\mathrm {in}~{\Omega }_{\mathrm {S}}, \end{aligned}$$
(36b)
$$\begin{aligned} \left\{ {\mathbb {C}}_{\mathrm {M}}:\left[ \varepsilon \,{\varvec{\xi }}_{x}\left( {\varvec{u}}^{\varepsilon }_{\mathrm {M}}\right) +{\varvec{\xi }}_{y}\left( {\varvec{u}}^{\varepsilon }_{\mathrm {M}}\right) \right] +\varepsilon {\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}-{\mathbb {C}}_{\mathrm {S}}:\left[ \varepsilon \,{\varvec{\xi }}_{x}\left( {\varvec{u}}^{\varepsilon }_{\mathrm {S}}\right) +{\varvec{\xi }}_{y}\left( {\varvec{u}}^{\varepsilon }_{\mathrm {S}}\right) \right] -\varepsilon {\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}\right\} {\varvec{n}}&={\varvec{0}},&\mathrm {on}~\Gamma , \end{aligned}$$
(36c)
$$\begin{aligned} {\varvec{u}}^{\varepsilon }_{\mathrm {M}}&={\varvec{u}}^{\varepsilon }_{\mathrm {S}},&\mathrm {on}~\Gamma . \end{aligned}$$
(36d)

We now proceed by equating the coefficients of the same powers of \(\varepsilon \). Hence, we obtain a set of differential conditions at various orders of \(\varepsilon \), holding in the periodic cell \(\Omega \), and parametrised by the macro-scale variable x.

5.1 First cell problem

By equating the coefficients of \(\varepsilon ^{0}\) in (36), it holds that

$$\begin{aligned} \nabla _y \,\cdot \left[ {\mathbb {C}}_{\mathrm {M}}:{\varvec{\xi }}_y\left( {\varvec{u}}^{(0)}_{\mathrm {M}}\right) \right]&={\varvec{0}},&\quad \mathrm {in}~\Omega _{\mathrm {M}}, \end{aligned}$$
(37a)
$$\begin{aligned} \nabla _y \,\cdot \left[ {\mathbb {C}}_{\mathrm {S}}:{\varvec{\xi }}_y\left( {\varvec{u}}^{(0)}_{\mathrm {S}}\right) \right]&={\varvec{0}},&\quad \mathrm {in}~\Omega _{\mathrm {S}}, \end{aligned}$$
(37b)
$$\begin{aligned} \left[ {\mathbb {C}}_\mathrm {M}:{\varvec{\xi }}_y\left( {\varvec{u}}^{(0)}_\mathrm {M}\right) -{\mathbb {C}}_\mathrm {S}:{\varvec{\xi }}_y\left( {\varvec{u}}^{(0)}_\mathrm {S}\right) \right] {\varvec{n}}&={\varvec{0}},&\quad \mathrm {on}~\Gamma , \end{aligned}$$
(37c)
$$\begin{aligned} {\varvec{u}}^{(0)}_\mathrm {M}&={\varvec{u}}^{(0)}_\mathrm {S},&\quad \mathrm {on}~\Gamma . \end{aligned}$$
(37d)

This boundary value problem (37), defined in the periodic cell \(\Omega \), is a linear-elastic type problem with zero source terms in (37a) and (37b). It is equipped with the stress jump condition (37c) and with the continuity condition on the zero-order displacements (37d), both holding on \(\Gamma \).

It has been proven that, for differential problems of the kind (37a)–(37d), the only periodic solutions are constant functions with respect to the macro-scale variable y (see [4, 17, 62]). Consequently, the leading order terms of the asymptotic expansions (31), \({\varvec{u}}_{\mathrm {M}}^{(0)}\) and \({\varvec{u}}_{\mathrm {S}}^{(0)}\), depend solely on the macro-scale variable x, thereby leading to the following writing

$$\begin{aligned} {\varvec{u}}^{(0)}_{\mathrm {M}}(x,y)={\varvec{u}}_{0\mathrm {M}}(x), \quad {\varvec{u}}^{(0)}_{\mathrm {S}}(x,y)={\varvec{u}}_{0\mathrm {S}}(x). \end{aligned}$$
(38)

In addition, the functions \({\varvec{u}}_{0\mathrm {M}}\) and \({\varvec{u}}_{0\mathrm {S}}\) satisfy identically (37c) and, by virtue of (37d), they coincide with each other, i.e.

$$\begin{aligned} {\varvec{u}}_{0\mathrm {M}}(x)={\varvec{u}}_{0\mathrm {S}}(x):={\varvec{u}}_{0}(x)\,. \end{aligned}$$
(39)

5.2 Second cell problem

Let us equate the coefficients of \(\varepsilon ^1\) in (36) and, by making use of (39), we are able to write the second cell problem in the form

$$\begin{aligned} \nabla _y \,\cdot \left[ {\mathbb {C}}_{\mathrm {M}}:{\varvec{\xi }}_y\left( {\varvec{u}}^{(1)}_{\mathrm {M}}\right) +{\mathbb {C}}_{\mathrm {M}}:{\varvec{\xi }}_x\left( {\varvec{u}}_{0}\right) +{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}\right]&={\varvec{0}} ,&\quad \mathrm {in}~\Omega _{\mathrm {M}}, \end{aligned}$$
(40a)
$$\begin{aligned} \nabla _y \,\cdot \left[ {\mathbb {C}}_{\mathrm {S}}:{\varvec{\xi }}_y\left( {\varvec{u}}^{(1)}_{\mathrm {S}}\right) +{\mathbb {C}}_{\mathrm {S}}:{\varvec{\xi }}_x\left( {\varvec{u}}_{0}\right) +{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}\right]&={\varvec{0}},&\quad \mathrm {in}~\Omega _{\mathrm {S}}, \end{aligned}$$
(40b)
$$\begin{aligned} \left[ {\mathbb {C}}_\mathrm {M}:{\varvec{\xi }}_y\left( {\varvec{u}}^{(1)}_\mathrm {M}\right) -{\mathbb {C}}_\mathrm {S}:{\varvec{\xi }}_y\left( {\varvec{u}}^{(1)}_\mathrm {S}\right) +\left( {\mathbb {C}}_\mathrm {M}-{\mathbb {C}}_\mathrm {S}\right) :{\varvec{\xi }}_x\left( {\varvec{u}}_{0}\right) +\left( {\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}-{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}\right) \right] {\varvec{n}}&={\varvec{0}},&\quad \mathrm {on}~\Gamma , \end{aligned}$$
(40c)
$$\begin{aligned} {\varvec{u}}^{(1)}_\mathrm {M}&={\varvec{u}}^{(1)}_\mathrm {S},&\quad \mathrm {on}~\Gamma . \end{aligned}$$
(40d)

In spite of the presence of the generalised Maxwell stress tensors \({\varvec{\sigma }}^{\mathrm {(e)}}_{\mathrm {M}}\) and \({\varvec{\sigma }}^{\mathrm {(e)}}_{\mathrm {S}}\), which constitute the main novelty of this work, the cell-problem (40a)–(40d) naturally yields the integral identity

$$\begin{aligned} \int \limits _{\partial \Omega _{\mathrm {M}}}\varvec{\Upsilon }_{\mathrm {M}}\,{\varvec{n}}+\int \limits _{\partial \Omega _{\mathrm {S}}}\varvec{\Upsilon }_{\mathrm {S}}\,{\varvec{n}}=\int \limits _{\Gamma }\left[\varvec{\Upsilon }_{\mathrm {M}}-\varvec{\Upsilon }_{\mathrm {S}}\right]\,{\varvec{n}}={\varvec{0}}, \end{aligned}$$
(41)

where \(\varvec{\Upsilon }_{\mathrm {M}}\) and \(\varvec{\Upsilon }_{\mathrm {S}}\) are abbreviations for the sums in the square brackets of Eqs. (40a) and (40b), respectively, and the difference \(\varvec{\Upsilon }_{\mathrm {M}}-\varvec{\Upsilon }_{\mathrm {S}}\) is thus the combination of the terms within the square brackets of (40c). The identity (41) is obtained by integrating Eqs. (40a) and (40b) over \(\Omega _{\mathrm {M}}\) and \(\Omega _{\mathrm {S}}\), respectively, applying Gauss’ Theorem and invoking the y-periodicity of the involved quantities. Clearly, the integrand \(\left[\varvec{\Upsilon }_{\mathrm {M}}-\varvec{\Upsilon }_{\mathrm {S}}\right]\,{\varvec{n}}\) coincides with the left-hand-side of (40c). Moreover, the fact that the sum of the surface integrals on the left-hand-side of (41) is null is consistent with the absence of external body forces of order \(\varepsilon ^{1}\) acting on the composite.

Before moving on with the solution to the problem (40), we need to ensure its well-posedness. As is customary in classical asymptotic homogenisation, this requires a compatibility, or solvability, condition for the specific differential problem under study (see [4, 64]). We will discuss about this aspect in the following Remark 4.

Remark 4

(An a priori given compatibility condition) If external forces appear on the right-hand sides of (40a) and (40b), the theory of asymptotic homogenisation requires the determination of a compatibility, or solvability, condition, for such forces [4]. In fact, their cell average must be null [4]. Hence, by indicating with \({\varvec{f}}_{ \,y}\) the term of order \(\varepsilon ^{1}\) with which such forces feature in the cell problem, the condition \(\int \limits _{\Omega }{\varvec{f}}_{ \,y}={\varvec{0}}\) has to be met [4]. Penta et al. [64] obtained this condition for the case of external forces of electric type and, more specifically, upon setting \({\varvec{f}}_{ \,y}=\nabla _{ \,y}\phi + \nabla _{ \,y}\times {\varvec{A}}\), where \(\phi \) and \({\varvec{A}}\) represent a scalar field and a vector field, which remind of the scalar potential and of the vector potential used in Electromagnetism, respectively. In the present study, however, the situation is substantially different. Indeed, although Eqs. (40a) and (40b) present the divergences of the electric part of the Cauchy stresses, i.e. \(\nabla _{ \,y}\cdot {\varvec{\sigma }}_{\mathrm {M}}^{\mathrm {(e)}}\) and \(\nabla _{ \,y}\cdot {\varvec{\sigma }}_{\mathrm {S}}^{\mathrm {(e)}}\), as inputs for the displacements \({\varvec{u}}_{\mathrm {M}}^{(1)}\) and \({\varvec{u}}_{\mathrm {S}}^{(1)}\), these divergences are not external forces. Rather, they stem from the electromechanical coupling discussed in the model, which results in the additive split of the overall stresses \({\varvec{\sigma }}_{\alpha }^{\mathrm {(t)}}={\varvec{\sigma }}_{\alpha }^{\mathrm {(m)}}+{\varvec{\sigma }}_{\alpha }^{\mathrm {(e)}}\) (see Appendix A). Moreover, although in this work the electric field is prescribed from the outset, the electric-type forces that it generates rely totally on the constitutive framework outlined in Sect. 2 and Appendix A. For this reason, and since Eqs. (40a) and (40b) express the static equilibrium of the composite at the order \(\varepsilon ^{1}\), we collected them on the left-hand side of (40a) and (40b). Consistently, the cell problem (40) turns out to be an elastic-type problem in the unknowns \({\varvec{u}}_{\alpha }^{(1)}\), \(\alpha \in \{\mathrm {M},\mathrm {S}\}\), with vanishing source terms, and equipped with a continuity condition for the displacement and a vanishing jump condition for the quantity \([\varvec{\Upsilon }_{\mathrm {M}}-\varvec{\Upsilon }_{\mathrm {S}}]{\varvec{n}}\) on the interface \(\Gamma \). Hence, by employing some classical results of asymptotic homogenisation, the existence of a unique solution is ensured and no further compatibility condition in required. Indeed, by applying Gauss’ Theorem and invoking the y-periodicity of the quantities involved in the expressions of order \(\varepsilon ^{1}\) of the overall Cauchy stresses, Eq. (41) is obtained naturally, i.e. without the need of looking for further compatibility conditions. In this perspective, we speak of Eq. (41) in terms of “an a priori given compatibility condition”.

According to the compatibility condition discussed in Remark 4, classical results of asymptotic homogenisation ensure that the boundary value problem (40) admits a unique solution, defined up to a function depending solely on x [4, 17]. Moreover, by virtue of linearity, a standard procedure (see for example [60, 62, 64]) allows to express \({\varvec{u}}_{\alpha }^{(1)}\), with \(\alpha \in \{\mathrm {M}, \mathrm {S}\}\), through the ansatz

$$\begin{aligned}&{\varvec{u}}^{(1)}_{\alpha }=\varvec{\chi }_{\alpha }:{\varvec{\xi }}_{x}\left( {\varvec{u}}_{0}\right) +{\varvec{\omega }}_{\alpha },&\alpha \in \{\mathrm {M},\mathrm {S}\}, \end{aligned}$$
(42a)
$$\begin{aligned}&[{\varvec{u}}^{(1)}_{\alpha }]_{a}=[\varvec{\chi }_{\alpha }]_{abc}[{\varvec{\xi }}_{x}\left( {\varvec{u}}_{0}\right) ]_{bc}+[{\varvec{\omega }}_{\alpha }]_{a},&\alpha \in \{\mathrm {M},\mathrm {S}\}, \end{aligned}$$
(42b)

where \(\varvec{\chi }_{\alpha }\in \mathrm {Lin}(\mathrm {Lin}({\mathcal {V}}),{\mathcal {V}})\) and \({\varvec{\omega }}_{\alpha }\in {\mathcal {V}}\) are a third-order tensor field and a vector field, respectively, that, in general, depend both on the macro-scale variable, x, and on the micro-scale variable, y. In addition, for consistency with the hypothesis of y-periodicity of \({\varvec{u}}^{(1)}_{\alpha }\), also \(\varvec{\chi }_{\alpha }\) and \({\varvec{\omega }}_{\alpha }\) must be y-periodic. Note that, in general, the tensor field \(\varvec{\chi }_{\alpha }\) possesses 27 independent components. However, since in (42) \(\varvec{\chi }_{\alpha }\) is double-contracted with \({\varvec{\xi }}_{x}({\varvec{u}}_{0})\), which is a symmetric second-order tensor field, only 18 independent components of \(\varvec{\chi }_{\alpha }\) are needed and these are given by \([\varvec{\chi }_{\alpha }]_{a[bc]}=\tfrac{1}{2}([\varvec{\chi }_{\alpha }]_{abc}+[\varvec{\chi }_{\alpha }]_{acb})\).

By substituting (42) into (40), two auxiliary differential problems are formulated, one for \(\varvec{\chi }_{\alpha }\) and one for \({\varvec{\omega }}_{\alpha }\). To make the notation as compact as possible, we introduce the auxiliary fourth-order tensor \({\mathbb {X}}_{y}(\varvec{\chi }_{\alpha })\), whose components are defined as

$$\begin{aligned} \left[ {\mathbb {X}}_{y}(\varvec{\chi }_{\alpha })\right] _{abcd}:=\dfrac{1}{2}\left( \dfrac{\partial \left[ \varvec{\chi }_{\alpha }\right] _{acd}}{\partial y_{b}}+\dfrac{\partial \left[ \varvec{\chi }_{\alpha }\right] _{bcd}}{\partial y_{a}}\right) ,\quad \alpha \in \{\mathrm {M}, \mathrm {S}\}, \end{aligned}$$
(43)

with the first two indices, a and b, being determined by the first index of \(\varvec{\chi }_{\alpha }\) and the direction along which the partial derivatives are computed, respectively. By construction, \({\mathbb {X}}_{y}(\varvec{\chi }_{\alpha })\) is symmetric in its first pair of indices, so that the equality \(\left[ {\mathbb {X}}_{y}(\varvec{\chi }_{\alpha })\right] _{abcd}=\left[ {\mathbb {X}}_{y}(\varvec{\chi }_{\alpha })\right] _{bacd}\) holds true. We notice that Eq. (43) defines \({\mathbb {X}}_{y}(\varvec{\chi }_{\alpha })\) as the result of the action of the operator \({\mathbb {X}}_{y}(\,\cdot \,)\) onto \(\varvec{\chi }_{\alpha }\), thereby extending the definition of \({\varvec{\xi }}(\,\cdot \,)\), given in (8), to tensor fields of the third order.

Finally, by using (43), the boundary-value problem for \(\varvec{\chi }_{\alpha }\) reads

$$\begin{aligned} \dfrac{\partial \left( \left[ {\mathbb {C}}_{\mathrm {M}}:{\mathbb {X}}_{y}(\varvec{\chi }_{\mathrm {M}})\right] _{abcd}+\left[ {\mathbb {C}}_{\mathrm {M}}\right] _{abcd}\right) }{\partial y_{b}}&=0,&\quad \mathrm {in}~\Omega _{\mathrm {M}}, \end{aligned}$$
(44a)
$$\begin{aligned} \dfrac{\partial \left( \left[ {\mathbb {C}}_{\mathrm {S}}:{\mathbb {X}}_{y}(\varvec{\chi }_{\mathrm {S}})\right] _{abcd}+\left[ {\mathbb {C}}_{\mathrm {S}}\right] _{abcd}\right) }{\partial y_{b}}&=0,&\quad \mathrm {in}~\Omega _{\mathrm {S}}, \end{aligned}$$
(44b)
$$\begin{aligned} \left[ {\mathbb {C}}_{\mathrm {M}}:{\mathbb {X}}_{y}(\varvec{\chi }_{\mathrm {M}})-{\mathbb {C}}_{\mathrm {S}}:{\mathbb {X}}_{y}(\varvec{\chi }_{\mathrm {S}})+{\mathbb {C}}_\mathrm {M}-{\mathbb {C}}_\mathrm {S}\right] _{abcd}n_b&=0,&\!\!\!\quad \mathrm {on}~\Gamma , \end{aligned}$$
(44c)
$$\begin{aligned} \left[ \varvec{\chi }_{\mathrm {M}}\right] _{acd}&=\left[ \varvec{\chi }_{\mathrm {S}}\right] _{acd},&\quad \!\!\! \mathrm {on}~\Gamma , \end{aligned}$$
(44d)

where each equation has to hold for all the values taken by the indices a, c and d. Analogously, the boundary value problem for \({\varvec{\omega }}_{\alpha }\) is given by

$$\begin{aligned} \dfrac{\partial \left( [{\mathbb {C}}_{\mathrm {M}}:{\varvec{\xi }}_y\left( {\varvec{\omega }}_{\mathrm {M}}\right) ]_{ab}+[{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}]_{ab}\right) }{\partial y_{b}}&=0,&\quad \mathrm {in}~\Omega _{\mathrm {M}}, \end{aligned}$$
(45a)
$$\begin{aligned} \dfrac{\partial \left( [{\mathbb {C}}_{\mathrm {S}}:{\varvec{\xi }}_y\left( {\varvec{\omega }}_{\mathrm {S}}\right) ]_{ab}+[{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}]_{ab}\right) }{\partial y_{b}}&=0,&\quad \mathrm {in}~\Omega _{\mathrm {S}}, \end{aligned}$$
(45b)
$$\begin{aligned} \left[ {\mathbb {C}}_{\mathrm {M}}:{\varvec{\xi }}_y\left( {\varvec{\omega }}_{\mathrm {M}}\right) -{\mathbb {C}}_{\mathrm {S}}:{\varvec{\xi }}_y\left( {\varvec{\omega }}_{\mathrm {S}}\right) + {\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}-{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}\right] _{ab}n_b&=0,&\quad \mathrm {on}~\Gamma , \end{aligned}$$
(45c)
$$\begin{aligned} \left[ {\varvec{\omega }}_{\mathrm {M}}\right] _a&=\left[ {\varvec{\omega }}_{\mathrm {S}}\right] _a,&\quad \mathrm {on}~\Gamma , \end{aligned}$$
(45d)

for all values taken by the index a. The problems (44) and (45) are completed with periodic boundary conditions on \(\partial \Omega \) and, for both of them, it is guaranteed the existence of a unique solution, defined up to a y-constant function [4, 17].

Remark 5

(Role of material properties in the computation of \(\varvec{\chi }_{\alpha }\) and \({\varvec{\omega }}_{\alpha }\)) Let us consider the weak form associated with the boundary value problem (44), i.e.

$$\begin{aligned}&\int \limits _{\Omega _{\mathrm {M}}}\left\{ \left[ {\mathbb {C}}_{\mathrm {M}}:{\mathbb {X}}_{y}(\varvec{\chi }_{\mathrm {M}})+{\mathbb {C}}_{\mathrm {M}}\right] _{abcd}\left[ {\mathbb {X}}_{y}(\varvec{\eta }_{\mathrm {M}})\right] _{abcd}\right\} +\int \limits _{\Omega _{\mathrm {S}}}\left\{ \left[ {\mathbb {C}}_{\mathrm {S}}:{\mathbb {X}}_{y}(\varvec{\chi }_{\mathrm {S}})+{\mathbb {C}}_{\mathrm {S}}\right] _{abcd} \left[ {\mathbb {X}}_{y}(\varvec{\eta }_{\mathrm {S}})\right] _{abcd}\right\} \nonumber \\&\quad =\int \limits _{\Omega _{\mathrm {M}}}\left[ {\mathbb {C}}_{\mathrm {M}}:{\mathbb {X}}_{y}(\varvec{\chi }_{\mathrm {M}})\right] _{abcd}\left[ {\mathbb {X}}_{y}(\varvec{\eta }_{\mathrm {M}})\right] _{abcd}+\int \limits _{\Omega _{\mathrm {S}}}\left[ {\mathbb {C}}_{\mathrm {S}}:{\mathbb {X}}_{y}(\varvec{\chi }_{\mathrm {S}})\right] _{abcd}\left[ {\mathbb {X}}_{y}(\varvec{\eta }_{\mathrm {S}})\right] _{abcd}\nonumber \\&\quad \quad +\,\int \limits _{\Omega _{\mathrm {M}}}\left[ {\mathbb {C}}_{\mathrm {M}}\right] _{abcd}\left[ {\mathbb {X}}_{y}(\varvec{\eta }_{\mathrm {M}})\right] _{abcd}+\int \limits _{\Omega _{\mathrm {S}}}\left[ {\mathbb {C}}_{\mathrm {S}}\right] _{abcd}\left[ {\mathbb {X}}_{y}(\varvec{\eta }_{\mathrm {S}})\right] _{abcd}=0, \end{aligned}$$
(46)

obtained by introducing suitable tensor-valued test functions \(\varvec{\eta }_{\mathrm {M}}\) and \(\varvec{\eta }_{\mathrm {S}}\), each of which representing a third-order, y-periodic tensor field. We notice that, in general, \({\mathbb {C}}_{\mathrm {M}}\) and \({\mathbb {C}}_{\mathrm {S}}\), besides depending on the material to which they refer, can vary within the matrix and the sub-phases, respectively. Such variability, as shown by (46), modulates the dependence of \(\varvec{\chi }_{\mathrm {M}}\) and \(\varvec{\chi }_{\mathrm {S}}\) on the fine-scale variable y. To provide a deeper insight in the way in which the elastic coefficients contribute to determine \(\varvec{\chi }_{\mathrm {M}}\) and \(\varvec{\chi }_{\mathrm {S}}\), let us assume the matrix and the sub-phases to be homogeneous, so that their elastic coefficients are phase-wise constant, and Eq. (46) rewrites as

$$\begin{aligned}&\int \limits _{\Omega _{\mathrm {M}}}\left[ {\mathbb {C}}_{\mathrm {M}}:{\mathbb {X}}_{y}(\varvec{\chi }_{\mathrm {M}})\right] _{abcd}\left[ {\mathbb {X}}_{y}(\varvec{\eta }_{\mathrm {M}})\right] _{abcd}+\int \limits _{\Omega _{\mathrm {S}}}\left[ {\mathbb {C}}_{\mathrm {S}}:{\mathbb {X}}_{y}(\varvec{\chi }_{\mathrm {S}})\right] _{abcd} \left[ {\mathbb {X}}_{y}(\varvec{\eta }_{\mathrm {S}})\right] _{abcd}\nonumber \\&\quad +\int \limits _{\Gamma }\left( \left[ {\mathbb {C}}_{\mathrm {M}}-{\mathbb {C}}_{\mathrm {S}}\right] _{abcd}n_{b}\right) \eta _{acd}=0, \end{aligned}$$
(47)

where \(\varvec{\eta }\) is defined as the common restriction of \(\varvec{\eta }_{\mathrm {M}}\) and \(\varvec{\eta }_{\mathrm {S}}\) to the interface \(\Gamma \), since \(\varvec{\eta }_{\mathrm {M}}\) and \(\varvec{\eta }_{\mathrm {S}}\) coincide on \(\Gamma \) by virtue of (44d). By looking at (47), we note that, if \({\mathbb {C}}_{\mathrm {M}}\) were equal to \({\mathbb {C}}_{\mathrm {S}}\), the integral over \(\Gamma \) would vanish identically, and the solutions to (47) and (44d) would be all those functions, and only those functions, \(\varvec{\chi }_{0}\) that are y-independent and satisfy the equalities \(\varvec{\chi }_{0}(x)=\varvec{\chi }_{\mathrm {M}}(x)=\varvec{\chi }_{\mathrm {S}}(x)\). However, since the inequality \({\mathbb {C}}_{\mathrm {M}}\ne {\mathbb {C}}_{\mathrm {S}}\) holds true throughout this work, the jump of the coefficients \({\mathbb {C}}_{\mathrm {M}}\) and \({\mathbb {C}}_{\mathrm {S}}\) across the interface \(\Gamma \) gives rise to a “surface generalised force”, with components \(\left[ {\mathbb {C}}_{\mathrm {M}}-{\mathbb {C}}_{\mathrm {S}}\right] _{abcd}n_{b}\), which expends work on \(\varvec{\eta }\) and generates non-trivial solutions for \(\varvec{\chi }_{\mathrm {M}}\) and \(\varvec{\chi }_{\mathrm {S}}\).

Analogously to (46), the weak form of the boundary value problem (45) reads

$$\begin{aligned}&\int \limits _{\Omega _{\mathrm {M}}}\big [{\mathbb {C}}_{\mathrm {M}}:{\varvec{\xi }}_y\left( {\varvec{\omega }}_{\mathrm {M}}\right) +{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}\big ]_{ab}\big [{\varvec{\xi }}_y\left( \varvec{\vartheta }_\mathrm {M}\right) \big ]_{ab}+\int \limits _{\Omega _{\mathrm {S}}}\big [{\mathbb {C}}_{\mathrm {S}}:{\varvec{\xi }}_y\left( {\varvec{\omega }}_{\mathrm {S}}\right) +{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}\big ]_{ab}\big [{\varvec{\xi }}_y\left( \varvec{\vartheta }_\mathrm {S}\right) \big ]_{ab}\nonumber \\&\quad =\int \limits _{\Omega _{\mathrm {M}}}\left[ {\mathbb {C}}_{\mathrm {M}}:{\varvec{\xi }}_y\left( {\varvec{\omega }}_{\mathrm {M}}\right) \right] _{ab}\left[ {\varvec{\xi }}_y\left( \varvec{\vartheta }_\mathrm {M}\right) \right] _{ab}+\int \limits _{\Omega _{\mathrm {S}}}\left[ {\mathbb {C}}_{\mathrm {S}}:{\varvec{\xi }}_y\left( {\varvec{\omega }}_{\mathrm {S}}\right) \right] _{ab}\left[ {\varvec{\xi }}_y\left( \varvec{\vartheta }_\mathrm {S}\right) \right] _{ab}\nonumber \\&\quad \quad +\int \limits _{\Omega _{\mathrm {M}}}\big [{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}\big ]_{ab}\big [{\varvec{\xi }}_y\left( \varvec{\vartheta }_\mathrm {M}\right) \big ]_{ab}+\int \limits _{\Omega _{\mathrm {S}}}\big [{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}\big ]_{ab}\big [{\varvec{\xi }}_y\left( \varvec{\vartheta }_\mathrm {S}\right) \big ]_{ab}=0, \end{aligned}$$
(48)

with \(\varvec{\vartheta }_{\mathrm {M}}\) and \(\varvec{\vartheta }_{\mathrm {S}}\) being suitable vector-valued y-periodic test functions associated with \({\varvec{\omega }}_{\mathrm {M}}\) and \({\varvec{\omega }}_{\mathrm {S}}\), respectively, and complying with the restriction \(\varvec{\vartheta }_{\mathrm {M}}=\varvec{\vartheta }_{\mathrm {S}}=\varvec{\vartheta }\) on \(\Gamma \). In this case, the fields \({\varvec{\omega }}_{\mathrm {M}}\) and \({\varvec{\omega }}_{\mathrm {S}}\) are modulated by the variability of both the elasticity tensor and the generalised Maxwell stress tensor in each phase. As above, however, an insightful simplification occurs when \({\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}\) and \({\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}\) are phase-wise constant, which allows to rewrite Eq. (48) as

$$\begin{aligned}&\int \limits _{\Omega _{\mathrm {M}}}\left[ {\mathbb {C}}_{\mathrm {M}}:{\varvec{\xi }}_y\left( {\varvec{\omega }}_{\mathrm {M}}\right) \right] _{ab}\left[ {\varvec{\xi }}_y\left( \varvec{\vartheta }_\mathrm {M}\right) \right] _{ab}+\int \limits _{\Omega _{\mathrm {S}}}\left[ {\mathbb {C}}_{\mathrm {S}}:{\varvec{\xi }}_y\left( {\varvec{\omega }}_{\mathrm {S}}\right) \right] _{ab}\left[ {\varvec{\xi }}_y\left( \varvec{\vartheta }_\mathrm {S}\right) \right] _{ab}\nonumber \\&\quad +\int \limits _{\Gamma }\left( [{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}-{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}]_{ab}n_{b}\right) \vartheta _{a}=0. \end{aligned}$$
(49)

Hence, although the matrix and the sub-phases are electrically homogeneous (i.e. their electric permittivities and, thus, the corresponding generalised Maxwell stress tensors are constant in \(\Omega _{\mathrm {M}}\) and \(\Omega _{\mathrm {S}}\)), the heterogeneity condition \({\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}\ne {\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}\) results into the contact force \([{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}-{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}]{\varvec{n}}\), which performs virtual work on \(\varvec{\vartheta }\) and yields non-trivial solutions for \({\varvec{\omega }}_{\mathrm {M}}\) and \({\varvec{\omega }}_{\mathrm {S}}\). We further emphasise that, for the problem at hand, the contact force \([{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}-{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}]{\varvec{n}}\) originates from the jump of the permittivities at the interface \(\Gamma \). Indeed, as anticipated in Eq. (14), the relation \({\varvec{\sigma }}^{(\mathrm {e})}_{\alpha }={\mathbb {T}}_{\alpha }:({\varvec{E}}\otimes {\varvec{E}})\) applies (see also Appendix A), so that the contact force reads

$$\begin{aligned}{}[{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}-{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}]{\varvec{n}}=\{[{\mathbb {T}}_{\mathrm {M}}-{\mathbb {T}}_{\mathrm {S}}]:({\varvec{E}}\otimes {\varvec{E}})\}{\varvec{n}} \; \Rightarrow \; [{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}-{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}]_{ab}n_{b}=\{[{\mathbb {T}}_{\mathrm {M}}-{\mathbb {T}}_{\mathrm {S}}]_{abcd}E_{c}E_{d}\}n_{b}. \end{aligned}$$
(50)

The terms \([{\mathbb {C}}_{\alpha }]_{abcd}[{\mathbb {X}}_{y}(\varvec{\eta }_{\alpha })]_{abcd}\) and \([{\varvec{\sigma }}_{\alpha }^{\mathrm {(e)}}]_{ab}[{\varvec{\xi }}_{y}(\varvec{\vartheta }_{\alpha })]_{ab}\) featuring in the weak forms (46) and (48) play a specific role in the determination of \(\varvec{\chi }_{\alpha }\) and \({\varvec{\omega }}_{\alpha }\). This role becomes evident when these weak forms are discretised according to a suitable numerical scheme, like the finite element method. Indeed, when this is done, two linear systems are obtained, in which the unknowns are the values taken by \(\varvec{\chi }_{\alpha }\) and \({\varvec{\omega }}_{\alpha }\) on a discrete set of nodes, while the source terms are given by the discretisation of \([{\mathbb {C}}_{\alpha }]_{abcd}[{\mathbb {X}}_{y}(\varvec{\eta }_{\alpha })]_{abcd}\) and \(\big [{\varvec{\sigma }}_{\alpha }^{\mathrm {(e)}}\big ]_{ab}\big [{\varvec{\xi }}_{y}(\varvec{\vartheta }_{\alpha })\big ]_{ab}\). In this sense, the variability of the material coefficients gives rise to generalised internal “forces”, which result into the source terms for the discretised versions of the two boundary value problems for \(\varvec{\chi }_{\alpha }\) and \({\varvec{\omega }}_{\alpha }\), respectively.

We conclude this section by proving that, in addition to the left minor symmetry discussed above, the fourth-order tensor \({\mathbb {X}}_{y}(\varvec{\chi }_{\alpha })\) is also symmetric with respect its second pair of indices, i.e. \([{\mathbb {X}}_{y}(\varvec{\chi }_{\alpha })]_{abcd}=[{\mathbb {X}}_{y}(\varvec{\chi }_{\alpha })]_{abdc}\). To deduce this property, let us rewrite the differential problem (44) in the following way,

$$\begin{aligned} \dfrac{\partial }{\partial y_{b}}\left( \left[ {\mathbb {C}}_{\mathrm {M}}\right] _{abmn}\dfrac{\partial \left[ \varvec{\chi }_{\mathrm {M}}\right] _{mcd}}{\partial y_{n}}+\left[ {\mathbb {C}}_{\mathrm {M}}\right] _{abcd}\right)&=0,&\quad \mathrm {in}~\Omega _{\mathrm {M}}, \end{aligned}$$
(51a)
$$\begin{aligned} \dfrac{\partial }{\partial y_{b}}\left( \left[ {\mathbb {C}}_{\mathrm {S}}\right] _{abmn}\dfrac{\partial \left[ \varvec{\chi }_{\mathrm {S}}\right] _{mcd}}{\partial y_{n}}+\left[ {\mathbb {C}}_{\mathrm {S}}\right] _{abcd}\right)&=0,&\quad \mathrm {in}~\Omega _{\mathrm {S}}, \end{aligned}$$
(51b)
$$\begin{aligned} \left( \left[ {\mathbb {C}}_{\mathrm {M}}\right] _{abmn}\dfrac{\partial \left[ \varvec{\chi }_{\mathrm {M}}\right] _{mcd}}{\partial y_{n}}-\left[ {\mathbb {C}}_{\mathrm {S}}\right] _{abmn}\dfrac{\partial \left[ \varvec{\chi }_{\mathrm {S}}\right] _{mcd}}{\partial y_{n}}\right) n_b+\left[ {\mathbb {C}}_\mathrm {M}-{\mathbb {C}}_\mathrm {S}\right] _{abcd}n_b&=0,&\quad \mathrm {on}~\Gamma , \end{aligned}$$
(51c)
$$\begin{aligned} \left[ \varvec{\chi }_{\mathrm {M}}\right] _{abc}&=\left[ \varvec{\chi }_{\mathrm {S}}\right] _{abc},&\quad \mathrm {on}~\Gamma , \end{aligned}$$
(51d)

where we have taken advantage of the right minor symmetry of the fourth-order elasticity tensor \({\mathbb {C}}_{\alpha }\). We note that, because of the right minor symmetry of \({\mathbb {C}}_{\alpha }\), the components \(\left[ \varvec{\chi }_{\alpha }\right] _{mcd}\) and \(\left[ \varvec{\chi }_{\alpha }\right] _{mdc}\) solve the same differential problem. This fact, together with the hypothesis of existence and uniqueness of the solution to the problem (44), implies the equality \(\left[ \varvec{\chi }_{\alpha }\right] _{acd}=\left[ \varvec{\chi }_{\alpha }\right] _{adc}\). Hence, the third-order tensor \(\varvec{\chi }_{\alpha }\) is symmetric in its last two indices, i.e. \([\varvec{\chi }_{\alpha }]_{acd}=[\varvec{\chi }_{\alpha }]_{a[cd]},\) which yields the left minor symmetry of \({\mathbb {X}}_{y}(\varvec{\chi }_{\alpha })\), i.e.

$$\begin{aligned} \left[ {\mathbb {X}}_{y}(\varvec{\chi }_{\alpha })\right] _{abcd}=\dfrac{1}{2}\left( \dfrac{\partial \left[ \varvec{\chi }_{\alpha }\right] _{acd}}{\partial y_{b}}+\dfrac{\partial \left[ \varvec{\chi }_{\alpha }\right] _{bcd}}{\partial y_{a}}\right) =\dfrac{1}{2}\left( \dfrac{\partial \left[ \varvec{\chi }_{\alpha }\right] _{adc}}{\partial y_{b}}+\dfrac{\partial \left[ \varvec{\chi }_{\alpha }\right] _{bdc}}{\partial y_{a}}\right) =\left[ {\mathbb {X}}_{y}(\varvec{\chi }_{\alpha })\right] _{abdc}, \end{aligned}$$
(52)

for \(\alpha \in \{\mathrm {M}, \mathrm {S}\}\). In general, nothing can be foreseen about the major symmetry of \({\mathbb {X}}_{y}(\varvec{\chi }_{\alpha })\).

By summarising the previous discussion, the differential problem (44) corresponds to a collection of three linear elastic problems, each of which individuated by the index a. Once a is fixed, by virtue of the symmetry properties of \({\mathbb {C}}_{\alpha }\), we obtain six independent scalar differential problems. Hence, we have a total number of eighteen scalar differential problems corresponding to the eighteen independent components of \(\varvec{\chi }_\alpha \). Also the boundary value problem (45) is formally of elastic-type, and it corresponds to three scalar independent differential equations for the three components of the vector field \({\varvec{\omega }}_\alpha \). For future use we highlight that the following identity holds true

$$\begin{aligned} {\varvec{\xi }}_y({\varvec{u}}^{(1)}_{\alpha })={\mathbb {X}}_{y}(\varvec{\chi }_{\alpha }):{\varvec{\xi }}_x\left( {\varvec{u}}_{0}\right) +{\varvec{\xi }}_y\left( {\varvec{\omega }}_{\alpha }\right) , \quad \alpha \in \{\mathrm {M}, \mathrm {S}\}. \end{aligned}$$
(53)

5.3 Homogenised problem

In this section, we show how it is possible to achieve two equivalent, yet conceptually different, effective formulations of the electro-mechanical model of the linear active composite under study. For this purpose, we start from the conditions

$$\begin{aligned} \nabla _x \cdot \left\{ {\mathbb {C}}_{\mathrm {M}} \,: \,\left[ \varvec{\xi }_x\left( {\varvec{u}}_{0}\right) \,+ \,\varvec{\xi }_y\left( {\varvec{u}}^{(1)}_{\mathrm {M}}\right) \right] \,+ \,{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {M}}\right\} \,+ \,\nabla _y \cdot \left\{ {\mathbb {C}}_{\mathrm {M}} \,: \,\left[ \varvec{\xi }_x\left( {\varvec{u}}^{(1)}_{\mathrm {M}}\right) \,+ \,\varvec{\xi }_y\left( {\varvec{u}}^{(2)}_{\mathrm {M}}\right) \right] \right\}&={\varvec{0}},&\mathrm {in}~\Omega _{\mathrm {M}}, \end{aligned}$$
(54a)
$$\begin{aligned} \nabla _x \cdot \left\{ {\mathbb {C}}_{\mathrm {S}} \,: \,\left[ \varvec{\xi }_x\left( {\varvec{u}}_{0}\right) \,+ \,\varvec{\xi }_y\left( {\varvec{u}}^{(1)}_{\mathrm {S}}\right) \right] \,+ \,{\varvec{\sigma }}^{(\mathrm {e})}_{\mathrm {S}}\right\} \,+ \,\nabla _y \cdot \left\{ {\mathbb {C}}_{\mathrm {S}} \,: \,\left[ \varvec{\xi }_x\left( {\varvec{u}}^{(1)}_{\mathrm {S}}\right) \,+ \,\varvec{\xi }_y\left( {\varvec{u}}^{(2)}_{\mathrm {S}}\right) \right] \right\}&={\varvec{0}},&\mathrm {in}~\Omega _{\mathrm {S}}, \end{aligned}$$
(54b)
$$\begin{aligned} \left\{ {\mathbb {C}}_\mathrm {M}:\left[ \varvec{\xi }_x\left( {\varvec{u}}^{(1)}_\mathrm {M}\right) +\varvec{\xi }_y\left( {\varvec{u}}^{(2)}_\mathrm {M}\right) \right] -{\mathbb {C}}_\mathrm {S}:\left[ \varvec{\xi }_x\left( {\varvec{u}}^{(1)}_\mathrm {S}\right) +\varvec{\xi }_y\left( {\varvec{u}}^{(2)}_\mathrm {S}\right) \right] \right\} {\varvec{n}}&={\varvec{0}},&\mathrm {on}~\Gamma , \end{aligned}$$
(54c)
$$\begin{aligned} {\varvec{u}}^{(2)}_\mathrm {M}&={\varvec{u}}^{(2)}_\mathrm {S},&\mathrm {on}~\Gamma , \end{aligned}$$
(54d)

which we obtain by matching the terms multiplied by \(\varepsilon ^2\) in Eq. (36) and employing the result \({\varvec{u}}_{\mathrm {M}}^{(0)}(x)={\varvec{u}}_{\mathrm {S}}^{(0)}(x)={\varvec{u}}_{0}(x)\).

Secondly, by averaging Eqs. (54a) and (54b) by means of the operator \(\langle \,\cdot \,\rangle _{\Omega _{\alpha }}\), \(\alpha \in \{\mathrm {M}, \mathrm {S}\}\), defined in the second expression of (32), we write

$$\begin{aligned}&\left\langle \nabla _x\cdot \left\{ {\mathbb {C}}_{\alpha } \,: \,\left[ \varvec{\xi }_x\left( {\varvec{u}}_{0}\right) \,+ \,\varvec{\xi }_y\left( {\varvec{u}}^{(1)}_{\alpha }\right) \right] \,+ \,{\varvec{\sigma }}^{(\mathrm {e})}_{\alpha }\right\} \right\rangle _{\Omega _{\alpha }} \,+ \,\left\langle \nabla _y \cdot \left\{ {\mathbb {C}}_{\alpha } \,: \,\left[ \varvec{\xi }_x\left( {\varvec{u}}^{(1)}_{\alpha }\right) \,+ \,\varvec{\xi }_y\left( {\varvec{u}}^{(2)}_{\alpha }\right) \right] \right\} \right\rangle _{\Omega _{\alpha }}={\varvec{0}}. \end{aligned}$$
(55)

Furthermore, we recast the arguments of the first angular brackets of Eq. (55) in the form

$$\begin{aligned} \left\langle \nabla _x \,\cdot \left[ {\mathbb {C}}_{\alpha }:{\varvec{\xi }}_x\left( {\varvec{u}}_{0}\right) \right] \right\rangle _{\Omega _{\alpha }}&=\nabla _x \,\cdot \left[ \left\langle {\mathbb {C}}_{\alpha }\right\rangle _{\Omega _{\alpha }}:{\varvec{\xi }}_x\left( {\varvec{u}}_{0}\right) \right] , \end{aligned}$$
(56a)
$$\begin{aligned} \left\langle \nabla _x \,\cdot \left[ {\mathbb {C}}_{\alpha }:{\varvec{\xi }}_y\left( {\varvec{u}}^{(1)}_{\alpha }\right) \right] \right\rangle _{\Omega _{\alpha }}&=\nabla _x \,\cdot \left[ \left\langle {\mathbb {C}}_{\alpha }:{\mathbb {X}}_{y}\left( \varvec{\chi }_{\alpha }\right) \right\rangle _{\Omega _{\alpha }}:{\varvec{\xi }}_x\left( {\varvec{u}}_{0}\right) \right] + \nabla _x \,\cdot \left[ \left\langle {\mathbb {C}}_{\alpha }:{\varvec{\xi }}_y\left( {\varvec{\omega }}_{\alpha }\right) \right\rangle _{\Omega _{\alpha }}\right] , \end{aligned}$$
(56b)
$$\begin{aligned} \left\langle \nabla _x \,\cdot {\varvec{\sigma }}^{(\mathrm {e})}_{\alpha }\right\rangle _{\Omega _{\alpha }}&=\nabla _x \,\cdot \langle {\varvec{\sigma }}^{(\mathrm {e})}_{\alpha }\rangle _{\Omega _{\alpha }}, \end{aligned}$$
(56c)

where we used the hypothesis of macroscopic uniformity reported in Eqs. (33) and (34), the expression of the symmetrised y-gradient of \({\varvec{u}}_{\alpha }^{(1)}\) in Eq. (53), and the fact that \({\varvec{u}}_{0}\) depends on x only. Now, looking at the arguments of the second angular brackets of Eq. (55), we notice that the topological and geometrical properties of the periodic cell (see Sect. 4) and the y-periodicity of \({\varvec{u}}_{\alpha }^{(1)}\), \({\varvec{u}}_{\alpha }^{(2)}\) and \({\mathbb {C}}_{\alpha }\) imply

$$\begin{aligned} \left\langle \nabla _y \,\cdot \left[ {\mathbb {C}}_{\alpha }:{\varvec{\xi }}_x ({\varvec{u}}^{(1)}_{\alpha }) \right] \right\rangle _{\Omega _{\alpha }}&\,= \,\dfrac{1}{|\Omega |}\int \limits _{\Omega _{\alpha }} \,\nabla _y \cdot \left[ {\mathbb {C}}_{\alpha }:{\varvec{\xi }}_x ({\varvec{u}}^{(1)}_{\alpha }) \right] \nonumber \\&=\dfrac{1}{|\Omega |}\int \limits _{\partial \Omega _{\alpha }}\left[ {\mathbb {C}}_{\alpha }:{\varvec{\xi }}_x ({\varvec{u}}^{(1)}_{\alpha }) \right] {\varvec{n}}_{\alpha }\,\mathrm {d}S=\dfrac{1}{|\Omega |}\int \limits _{\Gamma }\left[ {\mathbb {C}}_{\alpha }:{\varvec{\xi }}_x ({\varvec{u}}^{(1)}_{\alpha }) \right] {\varvec{n}}_{\alpha }\,\mathrm {d}S, \end{aligned}$$
(57a)
$$\begin{aligned} \left\langle \nabla _y \,\cdot \left[ {\mathbb {C}}_{\alpha }:{\varvec{\xi }}_y ({\varvec{u}}^{(2)}_{\alpha }) \right] \right\rangle _{\Omega _{\alpha }}&\,= \,\dfrac{1}{|\Omega |}\int \limits _{\Omega _{\alpha }} \,\nabla _y \cdot \left[ {\mathbb {C}}_{\alpha }:{\varvec{\xi }}_y ({\varvec{u}}^{(2)}_{\alpha }) \right] \nonumber \\&=\dfrac{1}{|\Omega |}\int \limits _{\partial \Omega _{\alpha }}\left[ {\mathbb {C}}_{\alpha }:{\varvec{\xi }}_y ({\varvec{u}}^{(2)}_{\alpha }) \right] {\varvec{n}}_{\alpha }\,\mathrm {d}S=\dfrac{1}{|\Omega |}\int \limits _{\Gamma }\left[ {\mathbb {C}}_{\alpha }:{\varvec{\xi }}_y ({\varvec{u}}^{(2)}_{\alpha }) \right] {\varvec{n}}_{\alpha }\,\mathrm {d}S, \end{aligned}$$
(57b)

where \({\varvec{n}}_{\alpha }={\varvec{n}}\) for \(\alpha =\mathrm {M}\) and \({\varvec{n}}_{\alpha }=-{\varvec{n}}\) for \(\alpha =\mathrm {S}\). In particular, for \(\Omega _{\alpha }=\Omega _{\mathrm {S}}\) the surface integrals over \(\partial \Omega _{\mathrm {S}}\) reduce to surface integrals over \(\Gamma \) because \(\Gamma \) coincides with \(\partial \Omega _{\mathrm {S}}\) in the case of inclusions and because of the y-periodicity of the integrands on \(\partial \Omega _{\mathrm {S}}\cap \partial \Omega \) in any other case considered in this work. For \(\Omega _{\alpha }=\Omega _{\mathrm {M}}\), an analogous result applies, with the integrands being this time y-periodic on \(\partial \Omega _{\mathrm {M}}\cap \partial \Omega \). Note that this analogy holds true for the class of cells chosen in this work, although the sets \(\partial \Omega _{\mathrm {S}}\cap \partial \Omega \) and \(\partial \Omega _{\mathrm {M}}\cap \partial \Omega \) can be topologically different from each other (in particular, \(\partial \Omega _{\mathrm {S}}\cap \partial \Omega \) is the empty set in the case of inclusions).

Finally, by substituting (56) and (57) in (55), we rewrite Eq. (55) as

$$\begin{aligned} \nabla _x \,\cdot \left[ {\mathbb {C}}_{\alpha ,\mathrm {eff}}:{\varvec{\xi }}_x\left( {\varvec{u}}_{0}\right) +\varvec{\tau }_{\alpha }\right] +\dfrac{1}{|\Omega |}\int \limits _{\Gamma }{\mathbb {C}}_{\alpha }:\left[ {\varvec{\xi }}_x ({\varvec{u}}^{(1)}_{\alpha })+{\varvec{\xi }}_y ({\varvec{u}}^{(2)}_{\alpha }) \right] {\varvec{n}}_{\alpha }\,\mathrm {d}S={\varvec{0}}, \end{aligned}$$
(58)

where we introduced the auxiliary notation

$$\begin{aligned} {\mathbb {C}}_{\alpha ,\mathrm {eff}}:=\langle {\mathbb {C}}_{\alpha } + {\mathbb {C}}_{\alpha }:{\mathbb {X}}_{y}(\varvec{\chi }_{\alpha })\rangle _{\Omega _{\alpha }}, \quad \varvec{\tau }_{\alpha }:=\langle {\varvec{\sigma }}_{\alpha }^{\mathrm {(e)}}+{\mathbb {C}}_{\alpha }:{\varvec{\xi }}_{y}({\varvec{\omega }}_{\alpha })\rangle _{\Omega _{\alpha }},\quad \alpha \in \{\mathrm {M}, \mathrm {S}\}. \end{aligned}$$
(59)

At this stage, the surface integrals over \(\Gamma \) can be eliminated by summing up Eq. (58) over \(\alpha \) and invoking (54c). This, in fact, yields the homogenised problem for the leading order term \({\varvec{u}}_{0}\), i.e.

$$\begin{aligned}&\nabla _x \,\cdot \left[ {\mathbb {C}}_{\mathrm {eff}}:{\varvec{\xi }}_x\left( {\varvec{u}}_{0}\right) +\varvec{\tau }\right] ={\varvec{0}}, \end{aligned}$$
(60)

where we have set

$$\begin{aligned} {\mathbb {C}}_{\mathrm {eff}}&:={\mathbb {C}}_{\mathrm {M,eff}}+{\mathbb {C}}_{\mathrm {S,eff}}=\langle {\mathbb {C}}_{\mathrm {M}} + {\mathbb {C}}_{\mathrm {M}}:{\mathbb {X}}_{y}(\varvec{\chi }_{\mathrm {M}})\rangle _{\Omega _{\mathrm {M}}}+\langle {\mathbb {C}}_{\mathrm {S}} + {\mathbb {C}}_{\mathrm {S}}:{\mathbb {X}}_{y}(\varvec{\chi }_{\mathrm {S}})\rangle _{\Omega _{\mathrm {S}}}, \end{aligned}$$
(61a)
$$\begin{aligned} \varvec{\tau }&:=\varvec{\tau }_{\mathrm {M}}+\varvec{\tau }_{\mathrm {S}}=\langle {\mathbb {C}}_{\mathrm {M}}:{\varvec{\xi }}_{y}({\varvec{\omega }}_{\mathrm {M}})+{\varvec{\sigma }}_{\mathrm {M}}^{\mathrm {(e)}}\rangle _{\Omega _{\mathrm {M}}}+\langle {\mathbb {C}}_{\mathrm {S}}:{\varvec{\xi }}_{y}({\varvec{\omega }}_{\mathrm {S}})+{\varvec{\sigma }}_{\mathrm {S}}^{\mathrm {(e)}}\rangle _{\Omega _{\mathrm {S}}}. \end{aligned}$$
(61b)

The effective elasticity tensor \({\mathbb {C}}_{\mathrm {eff}}\) is the sum of the two effective elasticity tensors \({\mathbb {C}}_{\mathrm {M,eff}}\) and \({\mathbb {C}}_{\mathrm {S,eff}}\), and is determined by “mixing” the intrinsic elastic properties of each phase with the pieces of micro-structural information enclosed in \({\mathbb {X}}_{y}(\varvec{\chi }_{\mathrm {M}})\) and \({\mathbb {X}}_{y}(\varvec{\chi }_{\mathrm {S}})\). In addition, \(\varvec{\tau }\) plays the role of an “extra” stress tensor, which adds itself to the homogenised version of the elastic Cauchy stress tensor \({\mathbb {C}}_{\mathrm {eff}}:{\varvec{\xi }}_x({\varvec{u}}_{0})\). More specifically, \(\varvec{\tau }\) describes the electro-mechanical coupling that originates, through the auxiliary fields \({\varvec{\omega }}_{\mathrm {M}}\) and \({\varvec{\omega }}_{\mathrm {S}}\), from the multi-scale nature of the composite under study.

We conclude this section noticing that it is also possible to provide an alternative formulation of the problem at hand. Indeed, by suitably rearranging the terms under divergence in Eq. (60), it can be shown that Eq. (60) becomes

$$\begin{aligned} \nabla _{x}\cdot \left[ {\varvec{\sigma }}_{\mathrm {eff}}^{\mathrm {(m)}} + {\varvec{\sigma }}_{\mathrm {eff}}^{\mathrm {(e)}}\right] ={\varvec{0}}, \end{aligned}$$
(62)

where we defined the effective mechanical stress tensor

$$\begin{aligned} {\varvec{\sigma }}_{\mathrm {eff}}^{\mathrm {(m)}} := {\mathbb {C}}_{\mathrm {eff}}:{\varvec{\xi }}_{x}({\varvec{u}}_{0})\,+\,\langle {\mathbb {C}}_{\mathrm {M}}:{\varvec{\xi }}_{y}({\varvec{\omega }}_{\mathrm {M}})\rangle _{\Omega _{\mathrm {M}}}\,+\,\langle {\mathbb {C}}_{\mathrm {S}}:{\varvec{\xi }}_{y}({\varvec{\omega }}_{\mathrm {S}})\rangle _{\Omega _{\mathrm {S}}}, \end{aligned}$$
(63)

and the effective generalised Maxwell stress tensor

$$\begin{aligned} {\varvec{\sigma }}^{\mathrm {(e)}}_{\mathrm {eff}}:=\langle {\varvec{\sigma }}_{\mathrm {M}}^{\mathrm {(e)}}\rangle _{\Omega _{\mathrm {M}}}+\langle {\varvec{\sigma }}_{\mathrm {S}}^{\mathrm {(e)}}\rangle _{\Omega _{\mathrm {S}}}. \end{aligned}$$
(64)

Clearly, either Eq. (60) or Eq. (62) has to be completed by boundary conditions to be specified for the only remaining unknown, i.e. \({\varvec{u}}_{0}\), on the boundary of the composite medium as a whole. Indeed, we recall that, within the present framework, the electric field is given from the outset, thereby completely defining \({\varvec{\sigma }}_{\mathrm {eff}}^{\mathrm {(e)}}\), the fields \({\varvec{\omega }}_{\mathrm {M}}\) and \({\varvec{\omega }}_{\mathrm {S}}\) are computed separately by solving the auxiliary problems (45a)–(45d), and \({\mathbb {C}}_{\mathrm {eff}}\) is determined by calculating \(\varvec{\chi }_{\mathrm {M}}\) and \(\varvec{\chi }_{\mathrm {S}}\).

It should be emphasised that the effective equation (62) has the same structure as (15a) or (15b), which represent the micro-scale force balances in \({\mathscr {C}}_{\mathrm {M}}\) and \({\mathscr {C}}_{\mathrm {S}}\), respectively. In spite of this “invariance”, the micro-scale formulation and the effective one feature an essential difference that manifests itself through the constitutive expression of \({\varvec{\sigma }}_{\mathrm {eff}}^{\mathrm {(m)}}\). The constitutive representation of this tensor, indeed, does not reduce to the double contraction of the effective elasticity tensor, \({\mathbb {C}}_{\mathrm {eff}}\), with \({\varvec{\xi }}_{x}({\varvec{u}}_{0})\). Rather, \({\varvec{\sigma }}_{\mathrm {eff}}^{\mathrm {(m)}}\) requires the contributions \(\langle {\mathbb {C}}_{\mathrm {M}}:{\varvec{\xi }}_{y}({\varvec{\omega }}_{\mathrm {M}})\rangle _{\Omega _{\mathrm {M}}}\) and \(\langle {\mathbb {C}}_{\mathrm {S}}:{\varvec{\xi }}_{y}({\varvec{\omega }}_{\mathrm {S}})\rangle _{\Omega _{\mathrm {S}}}\), which upscale the micro-scale heterogeneities related to the geometry and to the material coefficients of the composite, and “measure” the susceptibility of the stress response to the imposed electric field. Indeed, we recall that \({\varvec{\omega }}_{\mathrm {M}}\) and \({\varvec{\omega }}_{\mathrm {S}}\) are modulated by the electric field through \({\varvec{\sigma }}_{\mathrm {M}}^{\mathrm {(e)}}\) and \({\varvec{\sigma }}_{\mathrm {S}}^{\mathrm {(e)}}\), as prescribed by the problem (45).

A last comment pertains to the role played by \(\langle {\mathbb {C}}_{\mathrm {M}}:{\varvec{\xi }}_{y}({\varvec{\omega }}_{\mathrm {M}})\rangle _{\Omega _{\mathrm {M}}}\) and \(\langle {\mathbb {C}}_{\mathrm {S}}:{\varvec{\xi }}_{y}({\varvec{\omega }}_{\mathrm {S}})\rangle _{\Omega _{\mathrm {S}}}\) in the homogenised problem, which we formulate here with reference to a homogenised version of the original composite [60, 62, 64, 66], hereafter denoted by \({\mathscr {C}}_{\mathrm {H}}\). For this purpose, we assign the homogenised boundary value problem in the following form

$$\begin{aligned}&\nabla _{x}\cdot [{\varvec{\sigma }}^{\mathrm {(m)}}_{\mathrm {eff}}+{\varvec{\sigma }}^{\mathrm {(e)}}_{\mathrm {eff}}]={\varvec{0}},&\text {in} \, {\mathscr {C}}_{\mathrm {H}},\end{aligned}$$
(65a)
$$\begin{aligned}&{\varvec{u}}_{0}(x)={\varvec{u}}_{\mathrm {p}}(x),&\text {on} \, \partial _{\mathrm {D}}{\mathscr {C}}_{\mathrm {H}}, \end{aligned}$$
(65b)
$$\begin{aligned}&[{\varvec{\sigma }}^{\mathrm {(m)}}_{\mathrm {eff}}+{\varvec{\sigma }}^{\mathrm {(e)}}_{\mathrm {eff}}]{\varvec{n}}={\varvec{\sigma }}^{\mathrm {(e)}}_{\mathrm {ext}}{\varvec{n}},&\text {on} \, \partial _{\mathrm {N}}{\mathscr {C}}_{\mathrm {H}}, \end{aligned}$$
(65c)

where \({\varvec{n}}\) is the field of unit vectors normal to \(\partial {\mathscr {C}}_{\mathrm {H}}\), \({\varvec{\sigma }}^{\mathrm {(e)}}_{\mathrm {ext}}\) denotes here the generalised Maxwell stress tensor generated by the external electric field applied to \({\mathscr {C}}_{\mathrm {H}}\), and \({\varvec{u}}_{\mathrm {p}}(x)\) is the displacement prescribed on the Dirichlet boundary of \({\mathscr {C}}_{\mathrm {H}}\). Note that, while Eq. (65a) is the same as (62), the boundary conditions (65b) and (65c) prescribe a displacement on the Dirichlet boundary of \({\mathscr {C}}_{\mathrm {H}}\), i.e. \(\partial _{\mathrm {D}}{\mathscr {C}}_{\mathrm {H}}\), and a balance of contact forces on the Neumann boundary of \({\mathscr {C}}_{\mathrm {H}}\), i.e. \(\partial _{\mathrm {N}}{\mathscr {C}}_{\mathrm {H}}\). We remark that the homogenised composite, \({\mathscr {C}}_{\mathrm {H}}\), is an approximation of the original one and differs from this one in that its elastic properties are represented by the effective elasticity tensor, \({\mathbb {C}}_{\mathrm {eff}}\), and its electric response is described by the effective generalised Maxwell stress tensor \({\varvec{\sigma }}^{\mathrm {(e)}}_{\mathrm {eff}}\). In addition to this consideration, for the purposes of this work it is important to look at Eq. (65c), which considers the very instructive case in which, even in the absence of contact forces of mechanical type imposed on \(\partial _{\mathrm {N}}{\mathscr {C}}_{\mathrm {H}}\), the non-trivial traction \({\varvec{\sigma }}^{\mathrm {(e)}}_{\mathrm {ext}}{\varvec{n}}\) exists on this boundary [72]. Such traction, in fact, is the “contact force” resulting from the exposure of the composite to the electric field.

In our opinion, a deeper understanding of the problem (65) can be achieved by employing its weak formulation, which, upon introducing suitable test functions \({\varvec{v}}_{0}\), reads

$$\begin{aligned} \int \limits _{{\mathscr {C}}_{\mathrm {H}}}{\varvec{\xi }}_{x}({\varvec{v}}_{0}):{\mathbb {C}}_{\mathrm {eff}}:{\varvec{\xi }}_{x}({\varvec{u}}_{0}) =&- \int \limits _{{\mathscr {C}}_{\mathrm {H}}}\nabla _{x}{\varvec{v}}_{0}:\left[ \langle {\mathbb {C}}_{\mathrm {M}}:{\varvec{\xi }}_{y}({\varvec{\omega }}_{\mathrm {M}})\rangle _{\Omega _{\mathrm {M}}}+\langle {\mathbb {C}}_{\mathrm {S}}:{\varvec{\xi }}_{y}({\varvec{\omega }}_{\mathrm {S}})\rangle _{\Omega _{\mathrm {S}}}\right] \nonumber \\&-\int \limits _{{\mathscr {C}}_{\mathrm {H}}}\nabla _{x}{\varvec{v}}_{0}:{\varvec{\sigma }}_{\mathrm {eff}}^{\mathrm {(e)}} + \int \limits _{\partial _{\mathrm {N}}{\mathscr {C}}_{\mathrm {H}}}({\varvec{\sigma }}_{\mathrm {ext}}^{\mathrm {(e)}}{\varvec{n}})\cdot {\varvec{v}}_{0}. \end{aligned}$$
(66)

The terms on the right-hand-side of Eq. (66) can be interpreted as the virtual work done by generalised forces “external” to the macro-scale mechanical system described by \({\varvec{u}}_{0}\). More specifically, the first integral on the right-hand-side of (66) stems from the micro-structural descriptors \({\varvec{\omega }}_{\mathrm {M}}\) and \({\varvec{\omega }}_{\mathrm {S}}\) through the stress-like quantities \(\langle {\mathbb {C}}_{\mathrm {M}}:{\varvec{\xi }}_{y}({\varvec{\omega }}_{\mathrm {M}})\rangle _{\Omega _{\mathrm {M}}}\) and \(\langle {\mathbb {C}}_{\mathrm {S}}:{\varvec{\xi }}_{y}({\varvec{\omega }}_{\mathrm {S}})\rangle _{\Omega _{\mathrm {S}}}\), which are completely determined by the micro-scale problems (45), and are, in this regard, external to \({\varvec{u}}_{0}\). The last two integrals of (66), instead, are consequences of the electric field through the effective generalised Maxwell stress tensor, \({\varvec{\sigma }}_{\mathrm {eff}}^{\mathrm {(e)}}\), and the electric-type tractions \({\varvec{\sigma }}_{\mathrm {ext}}^{\mathrm {(e)}}{\varvec{n}}\). In particular, the structure of these last integrals permits to rephrase Eq. (66) in the suggestive form

$$\begin{aligned} \int \limits _{{\mathscr {C}}_{\mathrm {H}}}{\varvec{\xi }}_{x}({\varvec{v}}_{0}):{\mathbb {C}}_{\mathrm {eff}}:{\varvec{\xi }}_{x}({\varvec{u}}_{0}) =&- \int \limits _{{\mathscr {C}}_{\mathrm {H}}}\nabla _{x}{\varvec{v}}_{0}:\left[ \langle {\mathbb {C}}_{\mathrm {M}}:{\varvec{\xi }}_{y}({\varvec{\omega }}_{\mathrm {M}})\rangle _{\Omega _{\mathrm {M}}}+\langle {\mathbb {C}}_{\mathrm {S}}:{\varvec{\xi }}_{y}({\varvec{\omega }}_{\mathrm {S}})\rangle _{\Omega _{\mathrm {S}}}\right] \nonumber \\&+\int \limits _{{\mathscr {C}}_{\mathrm {H}}}(\nabla _{x}\cdot {\varvec{\sigma }}_{\mathrm {eff}}^{\mathrm {(e)}})\cdot {\varvec{v}}_{0} + \int \limits _{\partial _{\mathrm {N}}{\mathscr {C}}_{\mathrm {H}}}\llbracket {\varvec{\sigma }}_{\mathrm {ext}}^{\mathrm {(e)}}-{\varvec{\sigma }}_{\mathrm {eff}}^{\mathrm {(e)}} \rrbracket {\varvec{n}}\cdot {\varvec{v}}_{0}, \end{aligned}$$
(67)

where it is possible to identify the bulk force \(\nabla _{x}\cdot {\varvec{\sigma }}_{\mathrm {eff}}^{\mathrm {(e)}}\) and the surface force \(\llbracket {\varvec{\sigma }}_{\mathrm {ext}}^{\mathrm {(e)}}-{\varvec{\sigma }}_{\mathrm {eff}}^{\mathrm {(e)}} \rrbracket {\varvec{n}}\) triggered by the jump between the external and the effective generalised Maxwell stress tensors.

6 Particular cases

In this section we focus on relevant particular cases by assuming specific functional forms for the imposed electric field. In fact, we present how the model specialises when assuming that the electric field is either uniform on the micro-scale or multiplicatively decomposed into purely micro- and macro-scale components. The former case is frequently encountered in the context of experimental measurements and represents the simplest possible case also from a computational standpoint. The latter one is indeed a more general case, where we show how to reduce the computational complexity of the model by nonetheless retaining both micro- and macro-spatial variations of the electric field.

6.1 An imposed electric field which is uniform on the micro-scale

We now specialise the framework outlined in the previous sections to the case in which the applied electric field is independent on the micro-scale variable, y. In other words, we assume that the variability of the imposed electric field manifests itself solely at the scale of the composite as a whole, rather than at the micro-scale.

To the best of our knowledge, several experimental protocols are used to test the electro-mechanical response of active composites employed in different industrial contexts, and such procedures are often based on the use of electric fields depending on the macro-scale variable only.

By neglecting the micro-scale variation of the electric field, we limit ourselves to the case in which the imposed electric field varies inside the composite, but remains constant within the periodic cell. In light of this assumption, and looking at Eq. (14), we prescribe the generalised Maxwell stress tensor to depend on x and y in the following way:

$$\begin{aligned} {\varvec{\sigma }}^{(\mathrm {e})}_{\alpha }(x,y)={\mathbb {T}}_{\alpha }(x,y):({\varvec{E}}(x)\otimes {\varvec{E}}(x)), \quad \alpha \in \{\mathrm {M}, \mathrm {S}\}, \end{aligned}$$
(68)

where, as specified in Sect. 2, \({\mathbb {T}}_{\alpha }\) is a y-periodic symmetric fourth-order tensor whose components are functions of the electric permittivity of the \(\alpha \)-phase, \(\alpha \in \{\mathrm {M}, \mathrm {S}\}\) (see also Appendix A).

By taking inspiration from the procedure followed in the previous sections, we propose the following ansatz for the vector field \({\varvec{\omega }}_\alpha \):

$$\begin{aligned}&{\varvec{\omega }}_{\alpha }(x,y)={\varvec{\beta }}_{\alpha }(x,y):({\varvec{E}}(x)\otimes {\varvec{E}}(x)),&\alpha \in \{\mathrm {M},\mathrm {S}\}, \end{aligned}$$
(69a)
$$\begin{aligned}&[{\varvec{\omega }}_{\alpha }(x,y)]_{a}=[{\varvec{\beta }}_{\alpha }(x,y)]_{abc}\,E_{b}(x)\,E_{c}(x),&\alpha \in \{\mathrm {M},\mathrm {S}\}, \end{aligned}$$
(69b)

where \({\varvec{\beta }}_{\alpha }(x,\,\cdot \,)\) is a y-periodic third-order tensor field. Moreover, analogously to (43) and (53), we introduce the fourth-order tensor field \({\mathbb {X}}_y({\varvec{\beta }}_{\alpha })\), whose component expression is formally identical to the one defining \({\mathbb {X}}_y(\varvec{\chi }_\alpha )\) in Eq. (43), i.e.

$$\begin{aligned} \left[ {\mathbb {X}}_{y}({\varvec{\beta }}_{\alpha })\right] _{abcd}:=\dfrac{1}{2}\left( \dfrac{\partial \left[ {\varvec{\beta }}_{\alpha }\right] _{acd}}{\partial y_{b}}+\dfrac{\partial \left[ {\varvec{\beta }}_{\alpha }\right] _{bcd}}{\partial y_{a}}\right) ,\quad \alpha \in \{\mathrm {M}, \mathrm {S}\}. \end{aligned}$$
(70)

By employing Eq. (70), we can prove that the symmetrised micro-scale gradient of \({\varvec{\omega }}\), \({\varvec{\xi }}_y({\varvec{\omega }}_{\alpha })\), is given by

$$\begin{aligned}&{\varvec{\xi }}_y({\varvec{\omega }}_{\alpha })={\mathbb {X}}_{y}({\varvec{\beta }}_{\alpha }):({\varvec{E}}\otimes {\varvec{E}}), \end{aligned}$$
(71a)
$$\begin{aligned}&[{\varvec{\xi }}_{y}({\varvec{\omega }}_{\alpha })]_{mn}=[{\mathbb {X}}_{y}({\varvec{\beta }}_{\alpha })]_{mncd}E_{c}E_{d}, \end{aligned}$$
(71b)

where we have used the fact that the vector field \({\varvec{E}}\) is independent on the micro-scale variable y. By substituting Eqs. (71) in (45), we obtain the following differential conditions for \({\varvec{\beta }}_{\alpha }\):

$$\begin{aligned} \dfrac{\partial \left( [{\mathbb {C}}_{\mathrm {M}}:{\mathbb {X}}_y\left( {\varvec{\beta }}_{\mathrm {M}}\right) ]_{abcd}+[{\mathbb {T}}_\mathrm {M}]_{abcd}\right) }{\partial y_{b}}&=0,&\quad \mathrm {in}~\Omega _{\mathrm {M}}, \end{aligned}$$
(72a)
$$\begin{aligned} \dfrac{\partial \left( [{\mathbb {C}}_{\mathrm {S}}:{\mathbb {X}}_y\left( {\varvec{\beta }}_{\mathrm {S}}\right) ]_{abcd}+[{\mathbb {T}}_\mathrm {S}]_{abcd}\right) }{\partial y_{b}}&=0,&\quad \mathrm {in}~\Omega _{\mathrm {S}}, \end{aligned}$$
(72b)
$$\begin{aligned} \left[ {\mathbb {C}}_{\mathrm {M}}:{\mathbb {X}}_y\left( {\varvec{\beta }}_{\mathrm {M}}\right) -{\mathbb {C}}_{\mathrm {S}}:{\mathbb {X}}_y\left( {\varvec{\beta }}_{\mathrm {S}}\right) + {\mathbb {T}}_{\mathrm {M}}-{\mathbb {T}}_{\mathrm {S}}\right] _{abcd}n_b&=0,&\quad \mathrm {on}~\Gamma , \end{aligned}$$
(72c)
$$\begin{aligned} \left[ {\varvec{\beta }}_{\mathrm {M}}\right] _{acd}&=\left[ {\varvec{\beta }}_{\mathrm {S}}\right] _{acd},&\quad \mathrm {on}~ \Gamma . \end{aligned}$$
(72d)

We remark that the boundary value problem (72) is formally identical to (44). In particular, it consists of eighteen scalar partial differential equations for the eighteen independent components of \({\varvec{\beta }}_{\alpha }\). Indeed, by using the same arguments outlined in the previous section to prove the symmetry properties of \(\varvec{\chi }_{\alpha }\), and by virtue of the right minor symmetry of \({\mathbb {T}}_{\alpha }\), it is possible to prove that \({\varvec{\beta }}_{\alpha }\) has to be symmetric with respect to the last two indices.

Finally, by inserting (71) in (61b), we obtain that the average electro-mechanical contribution to the elastic stress tensor in (60) results to be

$$\begin{aligned} \varvec{\tau }={\mathbb {T}}_{\mathrm {eff}}:({\varvec{E}}\otimes {\varvec{E}}), \end{aligned}$$
(73)

where the fourth-order tensor \({\mathbb {T}}_{\mathrm {eff}}\) is given by

$$\begin{aligned} {\mathbb {T}}_{\mathrm {eff}}&:=\left\langle {\mathbb {T}}_\mathrm {M}+{\mathbb {C}}_{\mathrm {M}}:{\mathbb {X}}_{y}({\varvec{\beta }}_{\mathrm {M}})\right\rangle _{\Omega _{\mathrm {M}}}+\left\langle {\mathbb {T}}_\mathrm {S}+{\mathbb {C}}_{\mathrm {S}}:{\mathbb {X}}_{y}({\varvec{\beta }}_{\mathrm {S}})\right\rangle _{\Omega _{\mathrm {S}}} \end{aligned}$$
(74)

and contains both the micro-scale electric permittivities of the phases and the elastic contributions associated with \({\mathbb {C}}_{\mathrm {M}}\) and \({\mathbb {C}}_{\mathrm {S}}\). Finally it holds true that

$$\begin{aligned}&{\varvec{\sigma }}_{\mathrm {eff}}^{(\mathrm {m})}={\mathbb {C}}_{\mathrm {eff}}:{\varvec{\xi }}_x\left( {\varvec{u}}_{0}\right) +\left[ \langle \mathbb {C}_{\mathrm {M}}:\mathbb {X}_{y}(\varvec{\beta }_{\mathrm {M}})\rangle _{\Omega _{\mathrm {M}}}+\langle \mathbb {C}_{\mathrm {S}}:\mathbb {X}_{y}(\varvec{\beta }_{\mathrm {S}})\rangle _{\Omega _{\mathrm {S}}}\right] :({\varvec{E}}\otimes {\varvec{E}}), \end{aligned}$$
(75a)
$$\begin{aligned}&{\varvec{\sigma }}_{\mathrm {eff}}^{(\mathrm {e})}=\left( \left\langle {\mathbb {T}}_{\mathrm {M}}\right\rangle _{\Omega _{\mathrm {M}}}+\left\langle {\mathbb {T}}_{\mathrm {S}}\right\rangle _{\Omega _{\mathrm {S}}}\right) :({\varvec{E}}\otimes {\varvec{E}}). \end{aligned}$$
(75b)

6.2 Multiplicative decomposition of the electric field

The effective fourth-order elasticity tensor \({\mathbb {C}}_{\mathrm {eff}}\) in Eq. (61a), and the effective stress-like tensor \(\varvec{\tau }\) in Eq. (61b), which arise in the homogenised problem (60), are to be computed by solving the two cell problems (44) and (45). These two local problems are defined on the periodic cell \(\Omega \) and are parametrised by the macro-scale variable x. Such a dependence on the coarse scale, beyond the assumption of macroscopic uniformity, is retained in the elasticity tensors \({\mathbb {C}}_{\mathrm {M}}\) and \({\mathbb {C}}_{\mathrm {S}}\) and in the vector fields \({\varvec{\omega }}_{\mathrm {M}}\) and \({\varvec{\omega }}_{\mathrm {S}}\). In principle, this implies that (44) and (45) should be solved for each point x of the homogenised domain \({\mathscr {C}}_{\mathrm {H}}\) and, in general, this could require significant computational resources, in terms of calculation time and of memory storage.

In order to reduce the computational complexity of the two cell problems at hand, we make some simplifying hypotheses. By taking inspiration from [64], we assume that the fourth-order elasticity tensors, \({\mathbb {C}}_{\mathrm {M}}\) and \({\mathbb {C}}_{\mathrm {S}}\), and the fourth-order tensors accounting for the electric properties of the composite, \({\mathbb {T}}_{\mathrm {M}}\) and \({\mathbb {T}}_{\mathrm {S}}\), are x-constant, i.e. \({\mathbb {C}}_{\alpha }(x,y)={\mathbb {C}}_{\alpha }(y)\) and \({\mathbb {T}}_{\alpha }(x,y)={\mathbb {T}}_{\alpha }(y)\), for \(\alpha \in \{\mathrm {M}, \mathrm {S}\}\). Moreover, we propose that the components of the electric field \({\varvec{E}}\) can be multiplicatively decomposed into a locally varying, x-constant factor, and a macro-scale, y-constant, contribution. Hence, we prescribe that there exist functions

$$\begin{aligned} {\mathfrak {E}}_{(a)}:{\mathscr {C}}_{\mathrm {H}}\rightarrow {\mathbb {R}},\quad \quad {\mathfrak {e}}_{(a)}:\Omega \rightarrow {\mathbb {R}},\quad a=1,2,3, \end{aligned}$$
(76)

such that \(E_{a}(x,y)={\mathfrak {E}}_{(a)}(x){\mathfrak {e}}_{(a)}(y)\) and

$$\begin{aligned} {\varvec{E}}(x,y)={\mathfrak {E}}_{(a)}(x){\mathfrak {e}}_{(a)}(y)\,{\varvec{i}}_{a}={\mathfrak {E}}_{(1)}(x){\mathfrak {e}}_{(1)}(y)\,{\varvec{i}}_{1}+{\mathfrak {E}}_{(2)}(x){\mathfrak {e}}_{(2)}(y)\,{\varvec{i}}_{2}+{\mathfrak {E}}_{(3)}(x){\mathfrak {e}}_{(3)}(y)\,{\varvec{i}}_{3}, \end{aligned}$$
(77)

where the summation over the index a is omitted but understood. Enforcing (77), we can write the generalised Maxwell stress tensor as

$$\begin{aligned} {\varvec{\sigma }}_{\alpha }^{(\mathrm {e})}(x,y)={\mathfrak {E}}_{(a)}(x){\mathfrak {E}}_{(b)}(x){\mathfrak {e}}_{(a)}(y){\mathfrak {e}}_{(b)}(y)\,{\mathbb {T}}_{\alpha }(x,y):(\,{\varvec{i}}_{a}\otimes {\varvec{i}}_{b}), \end{aligned}$$
(78)

where we sum up over the indices a and b. With Eqs. (77) and (78) in mind, we introduce the following ansatz to the problem (45):

$$\begin{aligned} {\varvec{\xi }}_y({\varvec{\omega }}_{\alpha })={\mathfrak {E}}_{(a)}{\mathfrak {E}}_{(b)}\,{\mathbb {X}}_{y}({\varvec{\beta }}_{\alpha }^{ab}):({\varvec{i}}_{a}\otimes {\varvec{i}}_{b}). \end{aligned}$$
(79)

By substituting this ansatz, along with (78), into (45), we obtain six cell problems that do not retain any macro-scale dependency. In particular, by fixing the indices \(m,n=1,2,3\), the six cell problems are given by (no summation on the pair of indices m and n)

$$\begin{aligned} \dfrac{\partial \left( [{\mathbb {C}}_{\mathrm {M}}:{\mathbb {X}}_y\left( {\varvec{\beta }}_{\mathrm {M}}^{mn}\right) ]_{abmn}+[{\mathfrak {e}}_{(m)}{\mathfrak {e}}_{(n)}{\mathbb {T}}_\mathrm {M}]_{abmn}\right) }{\partial y_{b}}&=0,&\quad \mathrm {in}~\Omega _{\mathrm {M}}, \end{aligned}$$
(80a)
$$\begin{aligned} \dfrac{\partial \left( [{\mathbb {C}}_{\mathrm {S}}:{\mathbb {X}}_y\left( {\varvec{\beta }}_{\mathrm {S}}^{mn}\right) ]_{abmn}+[{\mathfrak {e}}_{(m)}{\mathfrak {e}}_{(n)} {\mathbb {T}}_\mathrm {S}]_{abmn}\right) }{\partial y_{b}}&=0,&\quad \mathrm {in}~\Omega _{\mathrm {S}}, \end{aligned}$$
(80b)
$$\begin{aligned} \left[ {\mathbb {C}}_{\mathrm {M}}:{\mathbb {X}}_y\left( {\varvec{\beta }}_{\mathrm {M}}^{mn}\right) -{\mathbb {C}}_{\mathrm {S}}:{\mathbb {X}}_y\left( {\varvec{\beta }}_{\mathrm {S}}^{mn}\right) \right] _{abmn}n_b&=\left[ {\mathfrak {e}}_{(m)}{\mathfrak {e}}_{(n)}({\mathbb {T}}_{\mathrm {S}}-{\mathbb {T}}_{\mathrm {M}})\right] _{abmn}n_b,&\quad \mathrm {on}~\Gamma , \end{aligned}$$
(80c)
$$\begin{aligned} \left[ {\varvec{\beta }}_{\mathrm {M}}^{mn}\right] _{acd}&=\left[ {\varvec{\beta }}_{\mathrm {S}}^{mn}\right] _{acd},&\quad \mathrm {on}~\Gamma . \end{aligned}$$
(80d)

We finally substitute (79) and (78) into (61b) to obtain the average electro-mechanical contribution to the elastic stress tensor in the homogenised problem (60). That is,

$$\begin{aligned} \varvec{\tau } \,= \,{\mathfrak {E}}_{(a)}{\mathfrak {E}}_{(b)} \,\left\{ \langle {\mathbb {C}}_{\mathrm {M}}:{\mathbb {X}}_{y}(\varvec{\beta }_{\mathrm {M}}^{ab})+{\mathfrak {e}}_{(a)}{\mathfrak {e}}_{(b)}{\mathbb {T}}_{\mathrm {M}}\rangle _{\Omega _{\mathrm {M}}}+\langle {\mathbb {C}}_{\mathrm {S}}:{\mathbb {X}}_{y}(\varvec{\beta }_{\mathrm {S}}^{ab})+{\mathfrak {e}}_{(a)}{\mathfrak {e}}_{(b)}{\mathbb {T}}_{\mathrm {S}}\rangle _{\Omega _{\mathrm {S}}}\right\} :({\varvec{i}}_{a}\otimes {\varvec{i}}_{b}). \end{aligned}$$
(81)

Also considering the alternative formulation for the homogenised problem, given by (62), we can write the effective mechanical stress tensor using the decomposition. That is,

$$\begin{aligned} {\varvec{\sigma }}_{\mathrm {eff}}^{\mathrm {(m)}}&= {\mathbb {C}}_{\mathrm {eff}}:{\varvec{\xi }}_{x}({\varvec{u}}_{0})+{\mathfrak {E}}_{(a)}{\mathfrak {E}}_{(b)} \,\left\{ \langle {\mathbb {C}}_{\mathrm {M}}:{\mathbb {X}}_{y}(\varvec{\beta }_{\mathrm {M}}^{ab})\rangle _{\Omega _{\mathrm {M}}}+\langle {\mathbb {C}}_{\mathrm {S}}:{\mathbb {X}}_{y}(\varvec{\beta }_{\mathrm {S}}^{ab})\rangle _{\Omega _{\mathrm {S}}}\right\} :({\varvec{i}}_{a}\otimes {\varvec{i}}_{b}). \end{aligned}$$
(82)

We can also write the effective generalised Maxwell stress tensor using the decomposition as

$$\begin{aligned} {\varvec{\sigma }}_{\mathrm {eff}}^{(\mathrm {e})}={\mathfrak {E}}_{(a)}{\mathfrak {E}}_{(b)}\left\{ \langle {\mathfrak {e}}_{(a)}{\mathfrak {e}}_{(b)}\,{\mathbb {T}}_{\mathrm {M}}\rangle _{\Omega _{\mathrm {M}}}+\langle {\mathfrak {e}}_{(a)}{\mathfrak {e}}_{(b)}\,{\mathbb {T}}_{\mathrm {S}}\rangle _{\Omega _{\mathrm {S}}}\right\} :(\,{\varvec{i}}_{a}\otimes {\varvec{i}}_{b}). \end{aligned}$$
(83)

Despite the introduction of a more complicated mathematical framework, the resulting decomposed model can be solved at a reduced computational cost when compared with the general model. The multiplicative decomposition of the electric field leads to the decoupling of the macro-scale and micro-scale problems. If we assume the elasticity tensors to be macroscopically uniform, then it is sufficient to solve the six elastic-type problems (44) and the six decomposed cell problems (80a)–(80d). The results from solving these problems can then be used to completely specify the homogenised problem (60).

We conclude this section by noticing that one crucial advantage of our formulation is the flexibility in the choice of the prescribed electric field. For instance, it is worth remarking that, in absence of free charges, the model proposed in [47] can be obtained as a particular case of our formulation (up to notation and rearranging terms) by assuming a functional form of the electric field as the leading (zeroth)-order electric field given in the latter work. In fact, in [47] the electric field is expressed in terms of the gradient of a scalar potential, which is in turn considered as a multiscale variable, which however (as in our case) does not depend on the elastic deformations. Although we could potentially prescribe an electric field which accounts a priori for the presence of a specific electric charge distribution, our modeling framework is designed for given, imposed electric fields. As such, whenever spatial variations of the latter are supposed to be influenced by the specific composite and/or its physical properties (as in the case of free charges or converse electrostriction), it could be more appropriate to proceed via considering the electric field as a multiscale variable, as done in [47]. In the next section we compare our formulation with different approaches to electrostriction by indeed assuming the absence of free charges.

7 Comparison with previous approaches to electrostriction

It is worth mentioning that, by employing the results summarised in the Appendix A, and specialising them for simplicity to the case of electrically isotropic materials, the total stress can be written as (cf. Eq. (132))

$$\begin{aligned} {\varvec{\sigma }}_{\alpha }^{\mathrm {(t)}}={\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha }+{\mathbb {T}}_{\alpha }:({\varvec{E}}\otimes {\varvec{E}}). \end{aligned}$$
(84)

Upon writing the fourth-order tensor \({\mathbb {T}}_{\alpha }\) in explicit form, i.e. as \({\mathbb {T}}_{\alpha }={\mathbb {E}}_{\alpha }+{\mathbb {E}}_{0}\) (cf. Eq. (129)), the total stress can be rearranged as

$$\begin{aligned} {\varvec{\sigma }}_{\alpha }^{\mathrm {(t)}}= \underbrace{{\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha } \,+ \,{\mathbb {E}}_{\alpha }:({\varvec{E}} \otimes {\varvec{E}})}_{:={\varvec{\sigma }}_{\alpha }^{\mathrm {(me)}}} +\underbrace{{\mathbb {E}}_{0}:({\varvec{E}}\otimes {\varvec{E}})}_{:={\varvec{\sigma }}^{\mathrm {(Max)}}}. \end{aligned}$$
(85)

Note that we have denoted the sum of the first two stress contributions by \({\varvec{\sigma }}_{\alpha }^{\mathrm {(me)}}\), where the superscript “\(\mathrm {(me)}\)” stands for mechano-electric stress tensor, because \({\varvec{\sigma }}_{\alpha }^{\mathrm {(me)}}\) features both the purely mechanical term, \({\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha }\), and the purely electric one, \({\mathbb {E}}_{\alpha }:({\varvec{E}}\otimes {\varvec{E}})\). Moreover, we have identified the last summand of Eq. (85) with the “true” Maxwell stress tensor, here denoted by \({\varvec{\sigma }}^{\mathrm {(Max)}}\).

7.1 Comparison with the Dorfmann & Ogden model [22]

To check for consistency with [22], we now subtract and add the tensor \({\varvec{E}}\otimes {\varvec{P}}_{\alpha }\) on the right-hand-side of Eq. (85), where \({\varvec{P}}_{\alpha }\) is given in Eq. (135b) of the Appendix A (for the sake of conciseness, we write \({\varvec{P}}_{\alpha }\) as \({\varvec{P}}_{\alpha }={\varvec{D}}_{\alpha }-\varepsilon _{0}{\varvec{E}}\) only in the term that is added). This yields

$$\begin{aligned} {\varvec{\sigma }}_{\alpha }^{\mathrm {(t)}}=&{\varvec{\sigma }}_{\alpha }^{\mathrm {(me)}}-{\varvec{E}}\otimes {\varvec{P}}_{\alpha }+{\varvec{E}}\otimes \left\{ {\varvec{D}}_{\alpha }-\epsilon _{0}{\varvec{E}}\right\} \,+{\mathbb {E}}_{0}:({\varvec{E}}\otimes {\varvec{E}}), \end{aligned}$$
(86)

and, by means of a direct calculation, involving the definition of \({\mathbb {E}}_{0}\) (cf. Eq. (125b)), we find

$$\begin{aligned} {\varvec{\sigma }}_{\alpha }^{\mathrm {(t)}}&={\varvec{\sigma }}_{\alpha }^{\mathrm {(me)}}-{\varvec{E}}\otimes {\varvec{P}}_{\alpha }+{\varvec{E}}\otimes {\varvec{D}}_{\alpha }-{\varvec{E}}\otimes \epsilon _{0}{\varvec{E}}\,+\epsilon _{0}{\varvec{E}}\otimes {\varvec{E}}-\tfrac{1}{2}\epsilon _{0}\Vert {\varvec{E}}\Vert ^{2}{\varvec{I}}\nonumber \\&=\left( {\varvec{\sigma }}_{\alpha }^{\mathrm {(me)}}-{\varvec{E}}\otimes {\varvec{P}}_{\alpha }\right) +\left( {\varvec{E}}\otimes {\varvec{D}}_{\alpha }-\tfrac{1}{2}\epsilon _{0}\Vert {\varvec{E}}\Vert ^{2}{\varvec{I}}\right) . \end{aligned}$$
(87)

Next, for comparison with Dorfmann & Ogden [22], and slightly adapting their notation to our framework, we set

$$\begin{aligned} {\varvec{\sigma }}_{\alpha }&:={\varvec{\sigma }}_{\alpha }^{\mathrm {(me)}}-{\varvec{E}}\otimes {\varvec{P}}_{\alpha }, \end{aligned}$$
(88a)
$$\begin{aligned} \varvec{\tau }_{\alpha }^{\mathrm {(e)}}&:={\varvec{E}}\otimes {\varvec{D}}_{\alpha }-\tfrac{1}{2}\epsilon _{0}\Vert {\varvec{E}}\Vert ^{2}{\varvec{I}}. \end{aligned}$$
(88b)

Then, by exploiting Maxwell’s equation \(\nabla \cdot {\varvec{D}}_{\alpha }=0\), for \(\alpha \in \{\mathrm {M},\mathrm {S}\}\), and the relation \({\varvec{D}}_{\alpha }={\varvec{P}}_{\alpha }+\epsilon _{0}{\varvec{E}}\), we notice that the divergence of \(\varvec{\tau }_{\alpha }^{\mathrm {(e)}}\) reads

$$\begin{aligned} \nabla \cdot \varvec{\tau }_{\alpha }^{\mathrm {(e)}}&=(\nabla {\varvec{E}}){\varvec{D}}_{\alpha }+{\varvec{E}}\,\nabla \cdot {\varvec{D}}_{\alpha }-\epsilon _{0}(\nabla {\varvec{E}})^{\mathrm {T}}{\varvec{E}} \nonumber \\&=(\nabla {\varvec{E}}){\varvec{P}}_{\alpha }+(\nabla {\varvec{E}})\epsilon _{0}{\varvec{E}}-\epsilon _{0}(\nabla {\varvec{E}})^{\mathrm {T}}{\varvec{E}}, \end{aligned}$$
(89)

and, since the last two terms have to cancel each other because of the constraint \(\nabla \times {\varvec{E}}={\varvec{0}}\) (we recall that magnetic effects are neglected throughout this work), we find

$$\begin{aligned} \nabla \cdot \varvec{\tau }_{\alpha }^{\mathrm {(e)}}&=(\nabla {\varvec{E}}){\varvec{P}}_{\alpha }, \end{aligned}$$
(90)

which suggests to identify \(\nabla \cdot \varvec{\tau }_{\alpha }^{\mathrm {(e)}}\) with the electric body force \({\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}:=\nabla \cdot \varvec{\tau }_{\alpha }^{\mathrm {(e)}}=(\nabla {\varvec{E}}){\varvec{P}}_{\alpha }\). Accordingly, the equilibrium equation \(\nabla \cdot {\varvec{\sigma }}_{\alpha }^{\mathrm {(t)}}={\varvec{0}}\) rewrites

$$\begin{aligned} \nabla \cdot {\varvec{\sigma }}_{\alpha }+{\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}={\varvec{0}}. \end{aligned}$$
(91)

7.2 Comparison with the model by Penta et al. [64]

Another comparison should be done with the results presented in [64], which are based on another definition of the electric body force \({\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}\). In [64], the overall approach is rather different from the one presented here, and the Authors define \({\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}\) as

$$\begin{aligned} {\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}:=\nabla \left( \tfrac{1}{2}\kappa _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}\right) =\nabla \cdot \left( \tfrac{1}{2}\kappa _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}\right) =\kappa _{\alpha }(\nabla {\varvec{E}})^{\mathrm {T}}{\varvec{E}}=\kappa _{\alpha }(\nabla {\varvec{E}}){\varvec{E}}, \end{aligned}$$
(92)

where we have slightly adapted the notation to our framework, we have assumed \(\kappa _{\alpha }\) to be constant and have used \((\nabla {\varvec{E}})^{\mathrm {T}}=\nabla {\varvec{E}}\).

In the limit of a very weak coupling between the mechanical displacement and the electric field, it is plausible to regard the electric displacement as independent of the displacement gradient, so that one may write \({\varvec{D}}_{\alpha }=(\epsilon _{0}+\epsilon _{\alpha }){\varvec{E}}\), which defines the polarisation vector as \({\varvec{P}}_{\alpha }=\epsilon _{\alpha }{\varvec{E}}\). Hence, by setting \(\kappa _{\alpha }=\epsilon _{\alpha }\), Eq. (92) becomes

$$\begin{aligned} {\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}=\nabla \cdot \left( \tfrac{1}{2}\epsilon _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}\right) =\epsilon _{\alpha }(\nabla {\varvec{E}}){\varvec{E}}=(\nabla {\varvec{E}}){\varvec{P}}_{\alpha }, \end{aligned}$$
(93)

and represents the body force known as electrostriction force. It is interesting to notice that, in the absence of free electric charges, it must hold true that \( \nabla \cdot {\varvec{D}}_{\alpha }=0\) and, if \(\epsilon _{0}\) and \(\epsilon _{\alpha }\) are constant, then \(\nabla \cdot {\varvec{E}}=0\) applies too. By assuming that these two conditions are fulfilled, and recalling Eq. (88b), \({\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}\) admits the equivalent expression

$$\begin{aligned} {\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}=\nabla \cdot \left( \tfrac{1}{2}\epsilon _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}\right) =(\nabla {\varvec{E}}){\varvec{P}}_{\alpha }=\nabla \cdot \left( {\varvec{E}}\otimes {\varvec{D}}_{\alpha }-\tfrac{1}{2}\epsilon _{0}\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}\right) \equiv \nabla \cdot \varvec{\tau }_{\alpha }^{\mathrm {(e)}}, \end{aligned}$$
(94)

which implies that the difference between \(\varvec{\tau }_{\alpha }^{\mathrm {(e)}}\) and \(\tfrac{1}{2}\epsilon _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}\) is a divergence-free second-order (stress) tensor, i.e.

$$\begin{aligned} \nabla \cdot \left[ \varvec{\tau }_{\alpha }^{\mathrm {(e)}}-\tfrac{1}{2}\epsilon _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}\right] \equiv \nabla \cdot \left[ \left( 1+\frac{\epsilon _{\alpha }}{\epsilon _{0}}\right) {\varvec{\sigma }}^{\mathrm {(Max)}}\right] =\left( 1+\frac{\epsilon _{\alpha }}{\epsilon _{0}}\right) \nabla \cdot {\varvec{\sigma }}^{\mathrm {(Max)}}={\varvec{0}}. \end{aligned}$$
(95)

To us, the physical meaning of this result is that, in the “gauge” \(\nabla \cdot {\varvec{\sigma }}^{\mathrm {(Max)}}={\varvec{0}}\), the tensor \(\tfrac{1}{2}\epsilon _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}\) is the “effective” part of \(\varvec{\tau }_{\alpha }^{\mathrm {(e)}}\), since both tensors have the same divergence and, thus, they both produce the same electrostriction force \({\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}\).

Variational deduction of the model by Penta et al. [64] Apparently, if Cauchy stress tensor is assumed to consist of the purely mechanical part only, i.e. \({\varvec{\sigma }}_{\alpha }={\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha }\), and if the electric field is regarded as given from the outset, provided it is both irrotational and solenoidal, then the force \({\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}\) can be viewed as an external force for the problem [64]

$$\begin{aligned} \nabla \cdot ({\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha })+{\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}={\varvec{0}}, \end{aligned}$$
(96)

in the unknown displacement field \({\varvec{u}}_{\alpha }\). Still, according to Eqs. (92)–(95), the force balance (96) can be rephrased as

$$\begin{aligned} \nabla \cdot ({\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha })+{\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}=\nabla \cdot \left( {\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha }+\tfrac{1}{2}\epsilon _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}\right) ={\varvec{0}}, \end{aligned}$$
(97)

and can be obtained by differentiating a suitable Lagrangian density function with respect to \(\nabla {\varvec{u}}_{\alpha }\). Indeed, a direct inspection shows that such a Lagrangian can be written as

$$\begin{aligned} {\hat{\mathcal {U}}}_{\alpha }(\nabla {\varvec{u}}_{\alpha })=-\tfrac{1}{2}\nabla {\varvec{u}}_{\alpha }:{\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha }-\tfrac{1}{2}\epsilon _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}:\nabla {\varvec{u}}_{\alpha }+h_{\alpha }, \end{aligned}$$
(98)

where \(h_{\alpha }\) is an arbitrary function independent of \(\nabla {\varvec{u}}_{\alpha }\), and that Eq. (97) becomes

$$\begin{aligned} -\nabla \cdot \left( \frac{\partial {\hat{\mathcal {U}}}_{\alpha }}{\partial \nabla {\varvec{u}}_{\alpha }}(\nabla {\varvec{u}}_{\alpha })\right) =\nabla \cdot \left( {\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha }+\tfrac{1}{2}\epsilon _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}\right) ={\varvec{0}}. \end{aligned}$$
(99)

Finally, we remark that, looking at Eq. (95), it holds true that \(\nabla \cdot \varvec{\tau }_{\alpha }^{\mathrm {(e)}}=\nabla \cdot \left( \tfrac{1}{2}\epsilon _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}\right) \) and, thus, the Euler–Lagrange equation (99) remains invariant if the Lagrangian \({\hat{\mathcal {U}}}_{\alpha }(\nabla {\varvec{u}}_{\alpha })\) is replaced with

$$\begin{aligned} {\hat{\mathcal {U}}}_{\alpha }^{\mathrm {(new)}}(\nabla {\varvec{u}}_{\alpha })= -\tfrac{1}{2}\nabla {\varvec{u}}_{\alpha }:{\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha }-\varvec{\tau }_{\alpha }^{\mathrm {(e)}}:\nabla {\varvec{u}}_{\alpha }+h_{\alpha }. \end{aligned}$$
(100)

Indeed, employing \({\hat{\mathcal {U}}}_{\alpha }^{\mathrm {(new)}}(\nabla {\varvec{u}}_{\alpha })\) in (99) and using the “gauge” condition (95) lead to

$$\begin{aligned} -\nabla \cdot \left( \frac{\partial {\hat{\mathcal {U}}}_{\alpha }^{\mathrm {(new)}}}{\partial \nabla {\varvec{u}}_{\alpha }}(\nabla {\varvec{u}}_{\alpha })\right)&=\nabla \cdot \left( {\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha }+\varvec{\tau }_{\alpha }^{\mathrm {(e)}}\right) \nonumber \\&=\nabla \cdot \left( {\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha }+\varvec{\tau }_{\alpha }^{\mathrm {(e)}}-\tfrac{1}{2}\epsilon _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}+\tfrac{1}{2}\epsilon _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}\right) \nonumber \\&=\nabla \cdot \left( {\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha }+\tfrac{1}{2}\epsilon _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}\right) +\underbrace{\nabla \cdot \left( \varvec{\tau }_{\alpha }^{\mathrm {(e)}}-\tfrac{1}{2}\epsilon _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}\right) }_{={\varvec{0}}}\nonumber \\&=\nabla \cdot \left( {\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha }+\tfrac{1}{2}\epsilon _{\alpha }\left\Vert {\varvec{E}}\right\Vert ^{2}{\varvec{I}}\right) ={\varvec{0}}. \end{aligned}$$
(101)

This result, in fact, is consistent with a property of Lagrangian density functions (see e.g. [35]), according to which \({\hat{\mathcal {U}}}_{\alpha }(\nabla {\varvec{u}}_{\alpha })\) and \({\hat{\mathcal {U}}}_{\alpha }^{\mathrm {(new)}}(\nabla {\varvec{u}}_{\alpha })\) lead to the same Euler–Lagrange equations, since they differ by the divergence of a vector-valued function depending only on \({\varvec{u}}_{\alpha }\) (and on the electric field, through \({\varvec{\sigma }}^{\mathrm {(Max)}}\)). Indeed, by exploiting this property in Eq. (100), one obtains

$$\begin{aligned} {\hat{\mathcal {U}}}_{\alpha }^{\mathrm {(new)}}(\nabla {\varvec{u}}_{\alpha })={\hat{\mathcal {U}}}_{\alpha }(\nabla {\varvec{u}}_{\alpha })-\nabla \cdot \left[ \left( 1+\frac{\epsilon _{\alpha }}{\epsilon _{0}}\right) {\varvec{\sigma }}^{\mathrm {(Max)T}}{\varvec{u}}_{\alpha }\right] , \end{aligned}$$
(102)

where we emphasise that the term between brackets on the right-hand-side of Eq. (102) has the structure of the Noether current generated by the Maxwell stress tensor (in fact, an energy-momentum tensor) [46].

The electrostriction force viewed as a “polygenic” force A completely different picture can be conceived if the electric field is understood as a variable of the problem, like the displacement, in spite of the fact that, under the hypotheses done in this section (i.e. \(\nabla \cdot {\varvec{D}}_{\alpha }=0\) and constant \(\epsilon _{0}\) and \(\epsilon _{\alpha }\)), it is a priori known to be solenoidal. Indeed, by looking at the force balance (96), with \({\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}\) regarded as an external force, no matter whether or not it is obtainable from a potential, and studying it in conjunction with Maxwell’s equation \(\nabla \cdot {\varvec{D}}_{\alpha }=0\), with \({\varvec{D}}_{\alpha }=(\epsilon _{0}+\epsilon _{\alpha }){\varvec{E}}\), we end up with the system of (decoupled) “equations of motion”

$$\begin{aligned}&\nabla \cdot \left( {\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha }\right) =-{\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}, \end{aligned}$$
(103a)
$$\begin{aligned}&\nabla \cdot [(\epsilon _{0}+\epsilon _{\alpha }){\varvec{E}}]=0. \end{aligned}$$
(103b)

In fact, by invoking the Lagrangian density function

$$\begin{aligned} {\hat{\mathcal {W}}}_{\alpha }(\nabla {\varvec{u}}_{\alpha },{\varvec{E}})=-\tfrac{1}{2}\nabla {\varvec{u}}_{\alpha }:{\mathbb {C}}_{\alpha }:\nabla {\varvec{u}}_{\alpha }+\tfrac{1}{2}(\epsilon _{0}+\epsilon _{\alpha })\left\Vert {\varvec{E}}\right\Vert ^{2}, \end{aligned}$$
(104)

Eqs. (103a) and (103b) become

$$\begin{aligned} -\nabla \cdot \left( \frac{\partial {\hat{\mathcal {W}}}_{\alpha }}{\partial \nabla {\varvec{u}}_{\alpha }}(\nabla {\varvec{u}}_{\alpha },{\varvec{E}})\right)&=-{\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}, \end{aligned}$$
(105a)
$$\begin{aligned} \nabla \cdot \left( \frac{\partial {\hat{\mathcal {W}}}_{\alpha }}{\partial {\varvec{E}}}(\nabla {\varvec{u}}_{\alpha },{\varvec{E}})\right)&=0, \end{aligned}$$
(105b)

provided we renounce to find a unique Lagrangian density function that is simultaneously compatible with the “equations of motion” and capable of returning \({\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}\). In other words, there exists no single scalar function from which the system of equations (103a) and (103b) descends. For this reason, by adopting the jargon of Lanczos [44], we refer to \({\varvec{f}}_{ \,\alpha }^{\mathrm {(e)}}\) as “polygenic force”.

8 Concluding remarks

In this work, we have studied the mechanical response of linear elastic active composites, i.e. composite materials characterised by a constitutive relationship in which the total stress can be decomposed into a purely linear elastic contribution and another component, which is assumed to be given. As our chief motivation has been the study of electrostrictive composites, we have specialised our formulation accordingly.

This class of active materials includes, among others, composites whose matrix is made of a polymer that encloses several dielectric sub-phases. Such materials have been investigated both from the theoretical and the experimental point of view [6, 7], and are largely employed in the field of biomechanics, as is the case of artificial muscles [49].

Within the framework sketched above, we have specialised our study to the case of electro-sensitive, linear elastic composites subjected to an imposed electric field. This setting, in fact, reflects various experimental procedures used to test, for instance, the electro-mechanical properties of materials in several industrial contexts.

Therefore, we have determined the effective electro-mechanical properties by having recourse to the theory of asymptotic homogenisation.

The starting point of our work is the system of equations (15). These equations have been obtained by assigning a Lagrangian density function and applying Hamilton’s variational procedure (see Appendix A), which yields a Cauchy stress consisting of a purely mechanical term, i.e. the classical Cauchy stress, and an electric one, identified with a generalised Maxwell stress tensor. The latter, in particular, is characterised by a constitutive expression that is strongly related to the way in which the mechanical and the electric descriptors (i.e. the displacement and the electric field, respectively) are coupled through the prescribed Lagrangian density function.

To study the effective behaviour of the composite under study, we have employed the asymptotic homogenisation technique, under the enforcement of a sharp separation between the length scales of the micro-structure and that of the composite as a whole. Under the hypothesis of micro-scale periodicity, we have obtained the boundary value problems (37) and (40), the auxiliary cell problems (44) and (45), each of which is defined on the elementary periodic cell, and the homogenised problem (65). Moreover, all the information concerning the micro-structure, i.e. the material coefficients of the two phases, the geometry of the periodic cell and the interaction of the two phases at the interface, is encapsulated in (62) through the effective stress tensors defined in (63) and (64).

One of the main results of our work, reported in Remark 4, is due to the fact that the formulation of the local equilibrium problem (15) does not call for the solvability condition that the theory of asymptotic homogenisation would invoke for equilibrium equations of the type \(\nabla \cdot {\varvec{\sigma }}+{\varvec{f}}={\varvec{0}}\). Indeed, in such cases, \({\varvec{f}}\) must fulfil the restriction \(\int \limits _{\Omega }{\varvec{f}}_{y}={\varvec{0}}\). However, since Eq. (15) features no external force of the same type as \({\varvec{f}}\), it is not necessary to look for any type of solvability. Rather, in the present framework, the only compatibility condition that has to be met is naturally supplied by Eq. (41). This is a direct consequence of the constitutive laws, of the fact that Eqs. (15a) and (15b) predict the vanishing of the divergence of the total stress tensor \({\varvec{\sigma }}^{\mathrm {(t)}}_{\alpha }\), and that the equilibrium condition at the interface (15c) involves \({\varvec{\sigma }}^{\mathrm {(t)}}_{\alpha }\) as a whole, rather than its mechanical part only, as in [64]. As such, our framework fully captures the role of the jumps of both the electric properties across the interfaces between different subphases in the composite, which could not be considered when specialising the work [64] to electrostriction.

Another result worth mentioning is that, with the introduction of an ansatz for the vector fields \({\varvec{\omega }}_{\mathrm {M}}\) and \({\varvec{\omega }}_{\mathrm {S}}\), we have converted the auxiliary cell problem (45) into a boundary value problem for the third-order tensor fields \({\varvec{\beta }}_{\mathrm {M}}\) and \({\varvec{\beta }}_{\mathrm {S}}\). Such differential problem, in fact, has the same structure as the auxiliary cell problem (44), with \({\varvec{\beta }}_{\alpha }\) sharing the same symmetry properties as \(\varvec{\chi }_{\alpha }\), with \(\alpha \in \{\mathrm {M}, \mathrm {S}\}\).

An advancement of our understanding of the theory of asymptotic homogenisation has been given by the introduction, for each boundary-value problem (see Eq. (45)), of its corresponding weak form. This has permitted to highlight the role of the elasticity tensors \({\mathbb {C}}_{\mathrm {M}}\) and \({\mathbb {C}}_{\mathrm {S}}\) in the determination of the auxiliary fields \(\varvec{\chi }_{\mathrm {M}}\) and \(\varvec{\chi }_{\mathrm {S}}\) (see Eqs. (46) and (47)), and the role of the contact force \([{\varvec{\sigma }}^{\mathrm (e)}_{\mathrm {M}}-{\varvec{\sigma }}^{\mathrm (e)}_{\mathrm {S}}]{\varvec{n}}\) on \({\varvec{\omega }}_{\mathrm {M}}\) and \({\varvec{\omega }}_{\mathrm {S}}\) (see Eqs. (48), (49) and (50)). Moreover, the use of the weak forms has improved our comprehension of the homogenised problem (65), especially for what concerns the nature of the “external forces” (see (66)).

The modelling hypotheses done in our work present some limitations. This fact makes this contribution open to generalisations from different points of view. Indeed, throughout this work, we have developed our model in the regime of small displacements and moderate electric fields. This choice permits to use, for the purposes of this work, a constitutive representation of the electric displacement that is linear in the electric field and independent of \(\nabla {\varvec{u}}_{\alpha }\). In other words, even though the “true” constitutive expression of \({\varvec{D}}_{\alpha }\) would couple \({\varvec{D}}_{\alpha }\) with \(\nabla {\varvec{u}}_{\alpha }\) (see e.g. Eq. (135a) in Appendix A), within our approximation we are allowed to use \({\varvec{D}}_{\alpha }={\varvec{\mathcal {E}}}_{\alpha }^{(0)}{\varvec{E}}\), with \({\varvec{\mathcal {E}}}_{\alpha }^{(0)}\) independent of \(\nabla {\varvec{u}}_{\alpha }\). Clearly, this approximation is an a priori estimate that should be confirmed a posteriori. Furthermore, in order to be compatible with Maxwell’s equations, the electric field has to be irrotational and solenoidal. In the same regime, if the electric field is not known from the outset, one has to solve the balance of linear momentum together with Maxwell’s equations. Within our approximation, the latter are completely decoupled from the former, and the electric field, computed by solving a Laplace equation for the scalar potential, plays the role of “an input” for the equations determining the mechanical displacement.

In a more general regime, it may be impossible to justify the approximation discussed so far, and the balance of linear momentum (i.e. the equation for \({\varvec{u}}_{\alpha }\)) should be solved together with Maxwell’s equation \(\nabla \cdot {\varvec{D}}_{\alpha }=0\).

In general, other modelling extensions of this work are possible. For example, one step ahead may consist in the possibility to admit also the existence of an imposed magnetic field, \({\varvec{H}}\), acting on the composite of interest. Such a situation would imply, on the one hand, the revision of the constitutive framework adopted so far and, on the other hand, the study of new cell problems to investigate the role played by each field in the determination of the effective quantities of the composite. Indeed, by activating the magnetic field \({\varvec{H}}\), these quantities would require to resolve the interactions among the three fields \({\varvec{u}}\), \({\varvec{E}}\) and \({\varvec{H}}\), thereby giving rise to a “game among three players”, in the jargon of [19] (see also [33]).

A step further in the research line of this work could be to assume the existence of a third scale, well separated from the other ones, and to switch to a three-scale-analysis of electro-sensitive composite materials [67,68,69]. In addition, another possible research direction could be to relax the hypothesis that the electric field is imposed from the outset, and to solve Maxwell’s equations in conjunction with the equilibrium problem (15). Doing this could pave the way towards the modelling of more realistic physical scenarios of biological and/or industrial interest. Finally, it is part of our research plans to formulate the presented framework in the context of finite deformations (see e.g. [12] and [11]) by following, for example, recent advances on homogenisation for nonlinear media, see also [18, 66]. Finally, the general results obtained in this work could be also exploited in contexts which are not necessarily related to electrostrictive composites but more generally to active composites, provided that suitable physical applications represented by the a priori given stress tensor are identified.