Homogenization of a coupled electrical and mechanical bidomain model for the myocardium

Abstract

We propose a coupled electrical and mechanical bidomain model for the myocardium tissue. The structure that we investigate possesses an elastic matrix with embedded cardiac myocytes. We are able to apply the asymptotic homogenization technique by exploiting the length scale separation that exists between the microscale where we see the individual myocytes and the overall size of the heart muscle. We derive the macroscale model which describes the electrical conductivity and elastic deformation of the myocardium driven by the existence of a Lorentz body force. The model comprises balance equations for the current densities and for the stresses, with the novel coefficients accounting for the difference in the electric potentials and elastic properties at different points in the microstructure. The novel coefficients of the model are to be computed by solving the periodic cell differential problems arising from application of the asymptotic homogenization technique. By combining both the mechanical and electrical behaviors, we obtain a macroscale model that highlights how the elastic deformation of the heart tissue is influenced and driven by the difference in the electric potentials at various points in the material.

1. Introduction

The heart is a muscular organ that pumps blood around the body. The muscular walls (myocardium) of the heart are formed from individual cardiac muscle cells which are connected together to form what is known as a functional electrical syncytium. This means there are strong electrical and mechanical connections between adjacent cardiac muscle cells in every direction allowing the myocardium to behave as a single contractile unit [1]. The cardiac muscle cells are known as myocytes and are connected at each end by the gap junctions at the intercalated disks. The heart contracts to pump blood around the body because of the electrical activation of the myocytes. In order to understand the physiology and electrical activity of the heart more thoroughly, we direct the reader to Katz [2], Opie [3], and Weidmann [4].

Many approaches have been made to model the electrophysiology and mechanical behavior of the heart [5–7]. One of the most important and prominent approaches is the constitutive nonlinear elastic modeling using the Holzapfel–Ogden law [8]. Holzapfel and Ogden [8] described the myocardium as a non-homogeneous, anisotropic, nonlinear elastic, and incompressible material. They then propose a strain energy function to describe the material. The parameters of the law are determined by biological measurements. This work has been subject to a variety of extensions that aim at understanding the behavior of the heart, such as in Wang et al. [9], and numerically in Pezzuto et al. [10].

Another very important approach is to consider an electrical bidomain model. This approach considers the propagation of cardiac action potentials and was first proposed in Barr and Jakobsson [11], Eisenberg [12], and Miller and Geselowitz [13]. The model is used to describe the development of the electrical potential through the heart muscle. The model accounts for the fact that at the microscopic scale, the action potentials arise due to the flow of specific types of ions through ion channels covering the cell membrane. We should note that the bidomain modeling approach is not exclusively considered for cardiac electrophysiology and can indeed be applied to the electrophysiology of neurons, e.g., Ellingsrud et al. [14] and Schreiner and Mardal [15].

Traditionally, the bidomain model was proposed for the electrical conductivity; however, there has been the development of the mechanical bidomain model to describe the elastic behavior of the heart, see Roth [16–19]. The mechanical bidomain model is able to account for the forces acting across the cell membrane which arise from the differences in the elastic displacement of the intracellular (inside the myocyte) and extracellular domains.

Each of the different bidomain models focuses on only either the mechanical or the electrical behavior and each has different advantages. The electrical bidomain model is used to simulate the electrical stimulation of the heart to understand how a voltage develops across the membrane between intracellular and extracellular domains. This is particularly useful to investigate cardiac pacing or defibrillation [20]. On the contrary, the mechanical bidomain model is used when considering mechanical forces are acting on and across the membrane, such as in the case of heart diseases and aging.

The heart has many structural features which are generally found at different microstructural levels (or scales). We are considering the heart muscle where we wish to look at the interactions between the myocytes and the extracellular matrix. Due to this, it is appropriate to consider a scale where we can see the myocytes and matrix clearly resolved. We call this scale the microscale. This is a fine microstructural level and has an associated length which is much smaller than the one characterising the whole heart. The complete heart muscle, however, can be described by a scale which we denote the macroscale.

The heart comprises multiple scales so we wish to create a model for it that can correctly characterize the effective properties but is also computationally feasible. To meet this goal, we require to relate the macroscale governing equations to the properties and interactions of the microscale. We approach this by creating a problem that describes the behavior and interactions of all constituents on the microscale. Once you have a problem of this type, it can be used in an upscaling process to obtain macroscale governing equations. There are a variety of techniques proposed in the literature that can be used for upscaling, and these are collectively known as homogenization techniques. Techniques such as mixture theory, effective medium theory, volume averaging, and asymptotic homogenization come under this category. Each technique has a variety of benefits and choice of method should be made depending on the application of the model and the information you wish to be encoded or available from the macroscale model. Each of the above-mentioned techniques has been reviewed and discussed in Hori and Nemat-Nasser [21, 22], while, e.g., Taffetani [23] provides an overview in context of mechanobiology.

Here, we make use of the asymptotic homogenization technique. Previously, this has been popularly used in modeling poroelasticity [24, 25], elastic composites [26, 27], and electroactive materials [28–30]. The technique has also proved useful in extending the theory such as to include growth and remodeling and vascularization of poroelastic materials [31, 32]. Recently, the technique was used in the derivation of the models of poroelastic composites and double poroelastic materials [33, 34]. To illustrate that the technique produces computationally feasible results, a micromechanical analysis of the effective stiffness of poroelastic composites has been investigated in Miller and Penta [35]. The technique has also previously been used in the context of heart modeling. In Miller and Penta [36], it has been used to investigate the structural changes involved in myocardial infarction and has even been used in the context of the electrical bidomain model (see Bader et al. [1] and Richardson and Chapman [37]).

This work involves applying the asymptotic homogenization technique to the problem that we have setup to describe the electrical and mechanical interactions between the cardiac myocytes and the surrounding matrix. We consider the heart muscle at a scale where the myocytes are distinct from the matrix. This scale is our microscale. The microscale has an associated length which is much smaller than the one of the entire heart muscle (where we cannot see individual myocytes) and so the scale of the heart is known as the macroscale. We then apply the asymptotic homogenization technique to upscale the problem, accounting for the continuity of current densities, stresses and elastic displacements as well as the difference in the electric potentials across the interface between the myocyte and the matrix. The novel macroscale model that is derived is a first attempt at an asymptotic homogenization model for a coupled electrical and mechanical bidomain model. The novel model comprises balance equations for the current densities and for the stresses as well as additional terms accounting for the difference in the electric potentials at different points in the microstructure. The coefficients of the model encode the properties of the microstructure and are computed by solving the microscale differential problems arising as a result of applying the asymptotic homogenization technique.

We have organized this work as follows. Section 2 focuses on introducing the mechanical and electrical problem that describes the microscale of the heart that we are considering. This problem describes the interactions between the intracellular domain (myocyte) and the extracellular domain (matrix). Then, in section 3, we perform the multiscale analysis of the problem. This procedure allows us to derive the new macroscale model governing the homogenized mechanical and electrical behavior of the heart tissue. In section 4, we present the new macroscale model and discuss the novel terms arising from considering both mechanical and electrical bidomain models coupled together. Then, in section 5, we consider the limit case where there is no change in electric potential at the interfaces between phases and reduce the model. We also provide a scheme that could be followed to solve the model numerically. In section 6, we summarize and draw conclusions to our work and provide further perspectives.

2. Electrical and mechanical problem

We begin by considering the microstructure of the myocardium which we call a set $\mathrm{\Omega }\in {\mathbb{R}}^3$ where $\mathrm{\Omega }$ is the union of the extracellular matrix ${\mathrm{\Omega }}_e$ and the intracellular myocyte ${\mathrm{\Omega }}_i$ , with $\overline{\mathrm{\Omega }}={\overline{\mathrm{\Omega }}}_i\cup {\overline{\mathrm{\Omega }}}_e$ . A sketch of a cross-section of the myocardium $\mathrm{\Omega }$ is shown in Figure 1.

Figure 1.

A 2D sketch representing a cross-section of the three-dimensional myocardium which we call our domain $\mathrm{\Omega }$ . The myocytes ${\mathrm{\Omega }}_i$ are shown in pink and the extracellular matrix ${\mathrm{\Omega }}_e$ is in blue.

We can then zoom in around one individual myocyte, and we see the clear microstructure which we have shown in Figure 2.

Figure 2.

A 2D sketch representing a cross-section of the three-dimensional domain $\mathrm{\Omega }$ . The intracellular domain (myocyte) ${\mathrm{\Omega }}_i$ is shown in pink and the extracellular domain ${\mathrm{\Omega }}_e$ is in blue. There is an interface $\mathrm{\Gamma }$ between the two domains.

We now wish to describe the equations that govern each domain and the chosen interface conditions to close the problem. We begin with the steady-state electrical bidomain equations:

$$\nabla ·({\mathsf{G}}_i\nabla {\varphi }_i)=\beta G({\varphi }_i-{\varphi }_e)\mathrm{in}{\mathrm{\Omega }}_i,$$

(1a)

$$\nabla ·({\mathsf{G}}_e\nabla {\varphi }_e)=-\beta G({\varphi }_i-{\varphi }_e)\mathrm{in}{\mathrm{\Omega }}_e,$$

(1b)

where ${\mathsf{G}}_i$ and ${\mathsf{G}}_e$ are the second rank conductivity tensors in the intracellular domain and extracellular domain, respectively, ${\varphi }_i$ and ${\varphi }_e$ are the scalar electric potentials of each phase, $\beta$ is the ratio of membrane area to tissue volume, and $G$ is the membrane conductance. Equations (1a) and (1b) are the balance equations for the electric current density for each phase given by:

$${\mathbf{j}}_i=-{\mathsf{G}}_i\nabla {\varphi }_i\mathrm{in}{\mathrm{\Omega }}_i,$$

(2a)

$${\mathbf{j}}_e=-{\mathsf{G}}_e\nabla {\varphi }_e\mathrm{in}{\mathrm{\Omega }}_e,$$

(2b)

which are Ohm’s law with the conductivity tensors ${\mathsf{G}}_i$ and ${\mathsf{G}}_e$ and the applied electric fields $\nabla {\varphi }_i$ and $\nabla {\varphi }_e$ .

In order to close the electrical part of the problem, we prescribe the following conditions on the interface $\mathrm{\Gamma }$ between the phases:

$${\mathsf{G}}_i\nabla {\varphi }_i·\mathbf{n}={\mathsf{G}}_e\nabla {\varphi }_e·\mathbf{n}\mathrm{on}\mathrm{\Gamma },$$

(3a)

$${\varphi }_i-{\varphi }_e=V\mathrm{on}\mathrm{\Gamma },$$

(3b)

where $V$ is a given and is the potential drop across the membrane [37], and n is the normal to the interface $\mathrm{\Gamma }$ pointing into the intracellular domain.

There are also the mechanical equations for each phase. We have the balance equations given by:

$$\nabla ·{\mathsf{T}}_i=K(\mathbf{u}-\mathbf{w})-{\mathsf{G}}_i\nabla {\varphi }_i\times \mathbf{B}\mathrm{in}{\mathrm{\Omega }}_i,$$

(4a)

$$\nabla ·{\mathsf{T}}_e=-K(\mathbf{u}-\mathbf{w})-{\mathsf{G}}_e\nabla {\varphi }_e\times \mathbf{B}\mathrm{in}{\mathrm{\Omega }}_e,$$

(4b)

where ${\mathsf{T}}_i$ and ${\mathsf{T}}_e$ are the stress tensors in each compartment, $\mathbf{u}$ and $\mathbf{w}$ are the elastic displacements of the intracellular and extracellular domains, respectively. Each domain is subject to the body forces given by the magnetic part of the Lorentz forces on the action potential currents, ${\mathbf{j}}_i$ and ${\mathbf{j}}_e$ , in a magnetic field $\mathbf{B}$ [17]. The body forces are therefore:

$${\mathbf{j}}_i\times \mathbf{B}=-{\mathsf{G}}_i\nabla {\varphi }_i\times \mathbf{B}\mathrm{and}{\mathbf{j}}_e\times \mathbf{B}=-{\mathsf{G}}_e\nabla {\varphi }_e\times \mathbf{B}.$$

(5)

The choice of these body forces is for a variety of reasons. The first being there has been previous use of this force in Puwal and Roth [17] where they examine a mechanical bidomain model and the effect of magnetic forces on action currents associated with a propagating action potential wavefront. By considering the same force, this will allow for future comparison with the numerical results obtained within Puwal and Roth [17]. There are also important applications such as the imaging of the elastic displacement due to Lorentz force which has been recently proposed as a potential use of MRI (see conclusions in section 6 for further discussion on this application). For further details of electroelastic materials and applied electric body forces, see, e.g., Dorfmann and Ogden [38, 39] and Maugin [40]. We should also note that the intracellular and extracellular domains are not separate from each other and interact through a variety of transmembrane proteins. This means that any displacement in the intracellular domain will cause displacements in the extracellular domain and vice versa. To cope with this, we follow the approach taken by Roth [16–18] where these interactions are described using Hooke’s law for the displacement difference and so $K$ can be thought of as the spring constant of the law. We assume that both the phases are anisotropic linear elastic materials with the stress tensors:

$${\mathsf{T}}_i={\mathbb{C}}_i\nabla \mathbf{u}\mathrm{in}{\mathrm{\Omega }}_i,$$

(6a)

$${\mathsf{T}}_e={\mathbb{C}}_e\nabla \mathbf{w}\mathrm{in}{\mathrm{\Omega }}_e,$$

(6b)

where ${\mathbb{C}}_i$ and ${\mathbb{C}}_e$ are the fourth rank elasticity tensors with corresponding components $(C_i{)}_{\alpha \delta \gamma \nu }$ and $(C_e{)}_{\alpha \delta \gamma \nu }$ , for $\alpha ,\delta ,\gamma ,\nu =1,2,3$ . We have that ${\mathbb{C}}_i$ and ${\mathbb{C}}_e$ are equipped with right minor and major symmetries, namely:

$$(C_i{)}_{\alpha \delta \gamma \nu }=(C_i{)}_{\alpha \delta \nu \gamma };(C_e{)}_{\alpha \delta \gamma \nu }=(C_e{)}_{\alpha \delta \nu \gamma },$$

(7a)

$$(C_i{)}_{\alpha \delta \gamma \nu }=(C_i{)}_{\gamma \nu \alpha \delta };(C_e{)}_{\alpha \delta \gamma \nu }=(C_e{)}_{\gamma \nu \alpha \delta },$$

(7b)

and therefore also left minor symmetries follow by combining equations (7a) and (7b). In particular, by applying right minor symmetries, we can equivalently rewrite the stress equations (6a) and (6b) as:

$${\mathsf{T}}_i={\mathbb{C}}_i\xi (\mathbf{u})\mathrm{in}{\mathrm{\Omega }}_i,$$

(8a)

$${\mathsf{T}}_e={\mathbb{C}}_e\xi (\mathbf{w})\mathrm{in}{\mathrm{\Omega }}_e,$$

(8b)

where

$$\xi (•)=\frac{\nabla (•)+{(\nabla (•))}^{\mathrm{T}}}{2},$$

(9)

is the symmetric part of the gradient operator.

Finally, we need to close the problem by prescribing conditions on the interface $\mathrm{\Gamma }$ , these are continuity of stresses and continuity of elastic/mechanical displacements:

$${\mathsf{T}}_i·\mathbf{n}={\mathsf{T}}_e·\mathbf{n}\mathrm{on}\mathrm{\Gamma },$$

(10a)

$$\mathbf{u}=\mathbf{w}\mathrm{on}\mathrm{\Gamma },$$

(10b)

where again $\mathbf{n}$ is the normal to the interface $\mathrm{\Gamma }$ pointing into the intracellular domain. In the next section, we perform a multiscale analysis. This involves (1) non-dimensionalizing the partial differential equations (PDEs) that we have described in this section, (2) introducing two well-separated length scales, (3) applying the asymptotic homogenization technique to the resulting non-dimensional systems of PDEs, and (4) deriving the effective governing equations for the material as a whole.

3. Multiscale analysis

We now perform a multiscale analysis of the problem introduced in the previous section, which is summarized below:

$$\nabla ·({\mathsf{G}}_i\nabla {\varphi }_i)=\beta G({\varphi }_i-{\varphi }_e)\mathrm{in}{\mathrm{\Omega }}_i,$$

(11a)

$$\nabla ·({\mathsf{G}}_e\nabla {\varphi }_e)=-\beta G({\varphi }_i-{\varphi }_e)\mathrm{in}{\mathrm{\Omega }}_e,$$

(11b)

$${\mathbf{j}}_i=-{\mathsf{G}}_i\nabla {\varphi }_i\mathrm{in}{\mathrm{\Omega }}_i,$$

(11c)

$${\mathbf{j}}_e=-{\mathsf{G}}_e\nabla {\varphi }_e\mathrm{in}{\mathrm{\Omega }}_e,$$

(11d)

$${\mathsf{G}}_i\nabla {\varphi }_i·\mathbf{n}={\mathsf{G}}_e\nabla {\varphi }_e·\mathbf{n}\mathrm{on}\mathrm{\Gamma },$$

(11e)

$${\varphi }_i-{\varphi }_e=V\mathrm{on}\mathrm{\Gamma },$$

(11f)

$$\nabla ·{\mathsf{T}}_i=K(\mathbf{u}-\mathbf{w})-{\mathsf{G}}_i\nabla {\varphi }_i\times \mathbf{B}\mathrm{in}{\mathrm{\Omega }}_i,$$

(11g)

$$\nabla ·{\mathsf{T}}_e=-K(\mathbf{u}-\mathbf{w})-{\mathsf{G}}_e\nabla {\varphi }_e\times \mathbf{B}\mathrm{in}{\mathrm{\Omega }}_e,$$

(11h)

$${\mathsf{T}}_i={\mathbb{C}}_i\xi (\mathbf{u})\mathrm{in}{\mathrm{\Omega }}_i,$$

(11i)

$${\mathsf{T}}_e={\mathbb{C}}_e\xi (\mathbf{w})\mathrm{in}{\mathrm{\Omega }}_e,$$

(11j)

$${\mathsf{T}}_i·\mathbf{n}={\mathsf{T}}_e·\mathbf{n}\mathrm{on}\mathrm{\Gamma },$$

(11k)

$$\mathbf{u}=\mathbf{w}\mathrm{on}\mathrm{\Gamma }.$$

(11l)

The problem (11a)–(11l) is closed by prescribing suitable external boundary conditions on $\partial \mathrm{\Omega }$ . We make the assumption that there exist two different length scales in the system. We denote the average size of the whole domain by $L$ (the macroscale), while $d$ refers to the microscale, which here is assumed to be comparable with the distance between adjacent myocytes. In order to emphasize the difference between such scales, it is helpful to perform a non-dimensional analysis of the system of PDEs (11a)–(11l).

3.1. Non-dimensionalization of the problem

We wish to formulate our model in non-dimensional form so as to clarify the mutual weight of each of the relevant fields in the problem. The intention is that the derived model will be general in regard to the heart microstructure so as to allow specification at a later stage to certain conditions or diseases. For this reason, we do not motivate the non-dimensionalization via specific parameters but instead perform a formal non-dimensionalization that will highlight the proper asymptotic behavior of each of the relevant fields. We make the choice to scale the spatial variable and the elastic displacement, and therefore the stresses and elasticity tensors, by the characteristic length scale L of the domain and an elastic scaling parameter E, ${\mathrm{\Phi }}_0$ is the typical potential drop, and ${\mathsf{G}}_0$ is the typical conductance. We therefore carry out the non-dimensional analysis by assuming the following:

$$\begin{array}{c}\mathbf{x}=L{\mathbf{x}}^{\prime },{\mathbb{C}}_i=\mathrm{EL}{{\mathbb{C}}^{\prime }}_i,{\mathbb{C}}_e=\mathrm{EL}{{\mathbb{C}}^{\prime }}_e,{\mathsf{T}}_i=\mathrm{EL}{\mathsf{T}}_i^{\prime },\hfill \\ \multicolumn{1}{c}{{\mathsf{T}}_e=\mathrm{EL}{{\mathsf{T}}^{\prime }}_e,\mathbf{u}=L{\mathbf{u}}^{\prime },\mathbf{w}=L{\mathbf{w}}^{\prime },{\varphi }_i={\mathrm{\Phi }}_0{{\varphi }^{\prime }}_i,}\\ \multicolumn{1}{c}{{\varphi }_e={\mathrm{\Phi }}_0{{\varphi }^{\prime }}_e,V={\mathrm{\Phi }}_0{V}^{\prime },{\mathsf{G}}_i={\mathsf{G}}_0{{\mathsf{G}}^{\prime }}_i,}\\ \multicolumn{1}{c}{{\mathsf{G}}_e={\mathsf{G}}_0{{\mathsf{G}}^{\prime }}_e,\mathbf{B}=\frac{L}{{\mathsf{G}}_0{\mathrm{\Phi }}_0}{\mathbf{B}}^{\prime }.}\end{array}$$

(12)

We are then able to use equation (12) and observe that:

$$\nabla =\frac{1}{L}{\nabla }^{\prime }.$$

(13)

These are substituted into equations (11a)–(11l) to obtain the following non-dimensional form of the system of PDEs:

$$\nabla ·({\mathsf{G}}_i\nabla {\varphi }_i)=\hat{\beta }({\varphi }_i-{\varphi }_e)\mathrm{in}{\mathrm{\Omega }}_i,$$

(14a)

$$\nabla ·({\mathsf{G}}_e\nabla {\varphi }_e)=-\hat{\beta }({\varphi }_i-{\varphi }_e)\mathrm{in}{\mathrm{\Omega }}_e,$$

(14b)

$${\mathbf{j}}_i=-{\mathsf{G}}_i\nabla {\varphi }_i\mathrm{in}{\mathrm{\Omega }}_i,$$

(14c)

$${\mathbf{j}}_e=-{\mathsf{G}}_e\nabla {\varphi }_e\mathrm{in}{\mathrm{\Omega }}_e,$$

(14d)

$${\mathsf{G}}_i\nabla {\varphi }_i·\mathbf{n}={\mathsf{G}}_e\nabla {\varphi }_e·\mathbf{n}\mathrm{on}\mathrm{\Gamma },$$

(14e)

$${\varphi }_i-{\varphi }_e=V\mathrm{on}\mathrm{\Gamma },$$

(14f)

$$\nabla ·{\mathsf{T}}_i=\hat{K}(\mathbf{u}-\mathbf{w})-{\mathsf{G}}_i\nabla {\varphi }_i\times \mathbf{B}\mathrm{in}{\mathrm{\Omega }}_i,$$

(14g)

$$\nabla ·{\mathsf{T}}_e=-\hat{K}(\mathbf{u}-\mathbf{w})-{\mathsf{G}}_e\nabla {\varphi }_e\times \mathbf{B}\mathrm{in}{\mathrm{\Omega }}_e,$$

(14h)

$${\mathsf{T}}_i={\mathbb{C}}_i\xi (\mathbf{u})\mathrm{in}{\mathrm{\Omega }}_i,$$

(14i)

$${\mathsf{T}}_e={\mathbb{C}}_e\xi (\mathbf{w})\mathrm{in}{\mathrm{\Omega }}_e,$$

(14j)

$${\mathsf{T}}_i·\mathbf{n}={\mathsf{T}}_e·\mathbf{n}\mathrm{on}\mathrm{\Gamma },$$

(14k)

$$\mathbf{u}=\mathbf{w}\mathrm{on}\mathrm{\Gamma },$$

(14l)

where

$$\hat{\beta }=\frac{\beta G{L}^2}{G_0}\mathrm{and}\hat{K}=\frac{\mathrm{LK}}{E},$$

(15)

are the dimensionless parameters.

In the next section, we introduce the asymptotic homogenization technique which is used to upscale the non-dimensional system of PDEs (14a)–(14l) by formally assuming that the microscale and the macroscale are well separated.

3.2. The asymptotic homogenization technique

In this section, we introduce the two-scale asymptotic homogenization technique which is used to derive a macroscale model for equations (14a)–(14l). We first assume that the microscale (where individual myocytes are clearly resolved from the matrix), denoted by $d$ , is small compared to average size of the domain $L$ , i.e.:

$$ϵ<<1.$$

(16)

We then introduce a local spatial variable to capture microscale variations of the field via setting:

$$\mathbf{y}=\frac{\mathbf{x}}{ϵ}.$$

(17)

The spatial variables $\mathbf{x}$ and $\mathbf{y}$ are to be considered formally independent and represent the macroscale and the microscale, respectively. The gradient operator then transforms as:

$$\nabla \to {\nabla }_{\mathbf{x}}+\frac{1}{ϵ}{\nabla }_{\mathbf{y}}.$$

(18)

We further assume that all the fields $\mathbf{u}$ , $\mathbf{w}$ , ${\mathsf{T}}_i$ , ${\mathsf{T}}_e$ , ${\varphi }_i$ , ${\varphi }_e$ , $V$ , as well as the elasticity tensors ${\mathbb{C}}_i$ and ${\mathbb{C}}_e$ , are the functions of both $\mathbf{x}$ and $\mathbf{y}$ . We also assume that the fields $\mathbf{u}$ , $\mathbf{w}$ , ${\mathsf{T}}_i$ , ${\mathsf{T}}_e$ , ${\varphi }_i$ , ${\varphi }_e$ , and $V$ can be represented in terms of a series expansion in powers of $ϵ$ , i.e.:

$${\phi }^{ϵ}(\mathbf{x},\mathbf{y},t)=\sum _{l=0}^{\infty }{\phi }^{(l)}(\mathbf{x},\mathbf{y},t){ϵ}^l,$$

(19)

where $\phi$ collectively denotes each field involved in the present analysis. We do, however, make the assumption that the difference in the electric potentials $V$ is given and has the following expansion:

$$V=V^{(0)}(\mathbf{x},t)+\cdots ,$$

(20)

i.e., it depends only on the macroscale at order zero. This assumption has also been made in Richardson and Chapman [37].

We make the following two assumptions in the next sections of this work.

Remark 1. Microscale periodicity. To simplify the analysis in this work, we restrict our attention to a single subset of the domain which we call the periodic cell. For this to be possible, we assume that every field ${\phi }^{(l)}$ in our problem (14a)–(14l) is $\mathbf{y}$ -periodic. By making this assumption, we can solve the microscale differential problems arising from using the asymptotic homogenization technique on just a finite bounded subset of the material. This assumption need not be made and the analysis can be carried out by assuming local boundedness of fields only. This approach, however, only allows us to determine the functional form of the macroscale model and the model coefficients are related to microscale problems that need to be solved on the whole microstructure of the material. This makes solving the model very computationally expensive when using the local boundedness of fields approach in comparison with microscale periodicity. Some examples of this are found in Burridge and Keller [24] and Penta and Gerisch [24].

Remark 2. Macroscopic uniformity. At individual macroscale points, it is clear that the microstructure can vary. This variation has been investigated by Burridge and Keller [24], Penta and Gerisch [26], Penta et al. [31], Holmes [42], and Dalwadi et al. [43]. The dependence of the microscale on the macroscale adds additional terms to the final model by proper application of the Reynolds transport theorem. In this work, we will neglect this dependence to simplify the derivation of the model. This means we assume that at every macroscale point, the microstructure will be the same. That is, the microscale geometry does not depend on $\mathbf{x}$ . This property is known as macroscopic uniformity. We will make this assumption in this work. We therefore have the simple differentiation under the integral sign:

$${\int }_{\mathrm{\Omega }}{\nabla }_{\mathbf{x}}·(•)\mathrm{dy}={\nabla }_{\mathbf{x}}·{\int }_{\mathrm{\Omega }}(•)\mathrm{dy},$$

(21)

where $(•)$ is a tensor or a vector quantity.

3.3. Multiple scales expansion

We apply the asymptotic homogenization assumptions (18) and (19) to equations (14a)–(14l) to obtain, accounting for periodicity, the following multiscale system of PDEs:

$${ϵ}^2{\nabla }_{\mathbf{x}}·({\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{ϵ})+ϵ{\nabla }_{\mathbf{x}}·({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{ϵ})+ϵ{\nabla }_{\mathbf{y}}·({\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{ϵ})+{\nabla }_{\mathbf{y}}·({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{ϵ})={ϵ}^2\hat{\beta }({\varphi }_i^{ϵ}-{\varphi }_e^{ϵ}),$$

(22a)

$${ϵ}^2{\nabla }_{\mathbf{x}}·({\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{ϵ})+ϵ{\nabla }_{\mathbf{x}}·({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{ϵ})+ϵ{\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{ϵ})+{\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{ϵ})=-{ϵ}^2\hat{\beta }({\varphi }_i^{ϵ}-{\varphi }_e^{ϵ}),$$

(22b)

$${\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{ϵ}·\mathbf{n}+ϵ{\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{ϵ}·\mathbf{n}={\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{ϵ}·\mathbf{n}+ϵ{\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{ϵ}·\mathbf{n},$$

(22c)

$$ϵ{\mathbf{j}}_i^{ϵ}=-ϵ{\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{ϵ}-{\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{ϵ},$$

(22d)

$$ϵ{\mathbf{j}}_e^{ϵ}=-ϵ{\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{ϵ}-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{ϵ},$$

(22e)

$${\varphi }_i^{ϵ}-{\varphi }_e^{ϵ}=V^{ϵ},$$

(22f)

$$ϵ{\nabla }_{\mathbf{x}}·{\mathsf{T}}_i^{ϵ}+{\nabla }_{\mathbf{y}}·{\mathsf{T}}_i^{ϵ}=ϵ\hat{K}({\mathbf{u}}^{ϵ}-{\mathbf{w}}^{ϵ})-ϵ{\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{ϵ}\times {\mathbf{B}}^{ϵ}-{\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{ϵ}\times {\mathbf{B}}^{ϵ},$$

(22g)

$$ϵ{\nabla }_{\mathbf{x}}·{\mathsf{T}}_e^{ϵ}+{\nabla }_{\mathbf{y}}·{\mathsf{T}}_e^{ϵ}=-ϵ\hat{K}({\mathbf{u}}^{ϵ}-{\mathbf{w}}^{ϵ})-ϵ{\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{ϵ}\times {\mathbf{B}}^{ϵ}-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{ϵ}\times {\mathbf{B}}^{ϵ},$$

(22h)

$$ϵ{\mathsf{T}}_i^{ϵ}=ϵ{\mathbb{C}}_i{\xi }_{\mathbf{x}}({\mathbf{u}}^{ϵ})+{\mathbb{C}}_i{\xi }_{\mathbf{y}}({\mathbf{u}}^{ϵ}),$$

(22i)

$$ϵ{\mathsf{T}}_e^{ϵ}=ϵ{\mathbb{C}}_e{\xi }_{\mathbf{x}}({\mathbf{w}}^{ϵ})+{\mathbb{C}}_e{\xi }_{\mathbf{y}}({\mathbf{w}}^{ϵ}),$$

(22j)

$${\mathsf{T}}_i^{ϵ}·\mathbf{n}={\mathsf{T}}_e^{ϵ}·\mathbf{n},$$

(22k)

$${\mathbf{u}}^{ϵ}={\mathbf{w}}^{ϵ}.$$

(22l)

We can now substitute power series of the type $(19)$ into the relevant fields in equations (22a)–(22l). Then, by equating the coefficients of ${ϵ}^l$ for $l=0,1,...$ , we derive the macroscale model for the material in terms of the relevant leading (zeroth) order fields. Whenever a component in the asymptotic expansion retains a dependence on the microscale, we can take the integral average, which we define as:

$$〈\phi {〉}_k=\frac{1}{|\mathrm{\Omega }|}{\int }_{{\mathrm{\Omega }}_i}\phi (\mathbf{x},\mathbf{y},t)d\mathbf{y}k=i,e,$$

(23)

where the integral average can be performed over one representative cell due to $\mathbf{y}$ -periodicity and $|\mathrm{\Omega }|$ is the volume of the domain and the integration is performed over the microscale. We note that $|\mathrm{\Omega }|=|{\mathrm{\Omega }}_i|+|{\mathrm{\Omega }}_e|$ . Due to the assumption of $\mathbf{y}$ -periodicity, the integral average can be performed over one representative cell. Therefore, equation (23) represents a cell average.

Equating coefficients of ${ϵ}^0$ :

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(0)})=0,$$

(24a)

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(0)})=0,$$

(24b)

$${\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(0)}·\mathbf{n}={\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(0)}·\mathbf{n},$$

(24c)

$${\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(0)}=0,$$

(24d)

$${\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(0)}=0,$$

(24e)

$${\varphi }_i^{(0)}-{\varphi }_e^{(0)}=V^{(0)},$$

(24f)

$${\nabla }_{\mathbf{y}}·{\mathsf{T}}_i^{(0)}=-{\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(0)}\times {\mathbf{B}}^{(0)},$$

(24g)

$${\nabla }_{\mathbf{y}}·{\mathsf{T}}_e^{(0)}=-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(0)}\times {\mathbf{B}}^{(0)},$$

(24h)

$${\mathbb{C}}_i{\xi }_{\mathbf{y}}({\mathbf{u}}^{(0)})=0,$$

(24i)

$${\mathbb{C}}_e{\xi }_{\mathbf{y}}({\mathbf{w}}^{(0)})=0,$$

(24j)

$${\mathsf{T}}_i^{(0)}·\mathbf{n}={\mathsf{T}}_e^{(0)}·\mathbf{n},$$

(24k)

$${\mathbf{u}}^{(0)}={\mathbf{w}}^{(0)}.$$

(24l)

From equations (24i) and (24j), we see that ${\mathbf{u}}^{(0)}$ and ${\mathbf{w}}^{(0)}$ are rigid body motions and therefore, by $\mathbf{y}$ -periodicity, do not depend on the microscale variable $\mathbf{y}$ . That is:

$${\mathbf{u}}^{(0)}={\mathbf{u}}^{(0)}(\mathbf{x},t),$$

(25a)

$${\mathbf{w}}^{(0)}={\mathbf{w}}^{(0)}(\mathbf{x},t).$$

(25b)

Since we have the interface condition (24l) ${\mathbf{u}}^{(0)}={\mathbf{w}}^{(0)}$ on $\mathrm{\Gamma }$ , we can use that:

$${\mathbf{u}}^{(0)}={\mathbf{w}}^{(0)},$$

(26)

throughout the following sections.

Due to the assumption (20) that $V^{(0)}$ is given and does not depend on the microscale variable $\mathbf{y}$ , we can write the ${ϵ}^0$ problem for ${{\varphi }_i}^{(0)}$ and ${{\varphi }_e}^{(0)}$ . We first define a new variable:

$${\overline{\varphi }}_e^{(0)}={\varphi }_e^{(0)}+V^{(0)}.$$

(27)

We then write the ${ϵ}^0$ problem in terms of ${{\varphi }_i}^{(0)}$ and our new variable ${\overline{\varphi }}_e^{(0)}$ :

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(0)})=0,$$

(28a)

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\overline{\varphi }}_e^{(0)})=0,$$

(28b)

$${\varphi }_i^{(0)}={\overline{\varphi }}_e^{(0)},$$

(28c)

$${\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(0)}·\mathbf{n}={\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\overline{\varphi }}_e^{(0)}·\mathbf{n}.$$

(28d)

The boundary value problem (28a)–(28d) is a linear-elastic type problem with no source terms in equations (28a) and (28b). It is equipped with the jump condition (28d) between the current densities and with the continuity condition between the zero-order electric potentials (28c). For problems of the type (28a)–(28d), it has been proved that the only solutions are constant with respect to the macroscale variable $\mathbf{y}$ (see Bakhvalov and Panasenko [44] and Cioranescu and Donato [45]). Therefore, the problem admits the solution that ${{\varphi }_i}^{(0)}$ and ${\overline{\varphi }}_e^{(0)}$ do not depend on the microscale variable $\mathbf{y}$ . It then follows that both ${{\varphi }_i}^{(0)}$ and ${{\varphi }_e}^{(0)}$ do not depend on that microscale, so we write:

$${{\varphi }_i}^{(0)}={{\varphi }_i}^{(0)}(\mathbf{x},t)$$

(29a)

$${{\varphi }_e}^{(0)}={{\varphi }_e}^{(0)}(\mathbf{x},t).$$

(29b)

This means that balance equations (24g) and (24h) can be rewritten as:

$${\nabla }_{\mathbf{y}}·{\mathsf{T}}_i^{(0)}=0,$$

(30a)

$${\nabla }_{\mathbf{y}}·{\mathsf{T}}_e^{(0)}=0.$$

(30b)

Similarly, we can now equate the coefficients of ${ϵ}^1$ :

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)})+{\nabla }_{\mathbf{y}}·({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(1)})=0,$$

(31a)

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)})+{\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(1)})=0,$$

(31b)

$${\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(1)}·\mathbf{n}+{\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}·\mathbf{n}={\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(1)}·\mathbf{n}+{\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}·\mathbf{n},$$

(31c)

$${\mathbf{j}}_i^{(0)}=-{\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}-{\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(1)},$$

(31d)

$${\mathbf{j}}_e^{(0)}=-{\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(1)},$$

(31e)

$${\varphi }_i^{(1)}-{\varphi }_e^{(1)}=V^{(1)},$$

(31f)

$${\nabla }_{\mathbf{x}}·{\mathsf{T}}_i^{(0)}+{\nabla }_{\mathbf{y}}·{\mathsf{T}}_i^{(1)}=\hat{K}({\mathbf{u}}^{(0)}-{\mathbf{w}}^{(0)})-{\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}\times {\mathbf{B}}^{(0)}-{\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(1)}\times {\mathbf{B}}^{(0)},$$

(31g)

$${\nabla }_{\mathbf{x}}·{\mathsf{T}}_e^{(0)}+{\nabla }_{\mathbf{y}}·{\mathsf{T}}_e^{(1)}=-\hat{K}({\mathbf{u}}^{(0)}-{\mathbf{w}}^{(0)})-{\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}\times {\mathbf{B}}^{(0)}-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(1)}\times {\mathbf{B}}^{(0)},$$

(31h)

$${\mathsf{T}}_i^{(0)}={\mathbb{C}}_i{\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)})+{\mathbb{C}}_i{\xi }_{\mathbf{y}}({\mathbf{u}}^{(1)}),$$

(31i)

$${\mathsf{T}}_e^{(0)}={\mathbb{C}}_e{\xi }_{\mathbf{x}}({\mathbf{w}}^{(0)})+{\mathbb{C}}_e{\xi }_{\mathbf{y}}({\mathbf{w}}^{(1)}),$$

(31j)

$${\mathsf{T}}_i^{(1)}·\mathbf{n}={\mathsf{T}}_e^{(1)}·\mathbf{n},$$

(31k)

$${\mathbf{u}}^{(1)}={\mathbf{w}}^{(1)}.$$

(31l)

Now that we have equated the coefficients of powers 0 and 1 of $ϵ$ , we can use these equations to form problems for the order 1 electric potentials and elastic displacements.

3.4. Problem for electric potentials ${\varphi }_i^{(1)}$ and ${\varphi }_e^{(1)}$

Using equations (31a)–(31c) and (31f), we can form the following problem for the electric potentials ${\varphi }_i^{(1)}$ and ${\varphi }_e^{(1)}$ . That is:

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)})+{\nabla }_{\mathbf{y}}·({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(1)})=0\mathrm{in}{\mathrm{\Omega }}_i,$$

(32a)

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)})+{\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(1)})=0\mathrm{in}{\mathrm{\Omega }}_e,$$

(32b)

$${\varphi }_i^{(1)}-{\varphi }_e^{(1)}=V^{(1)}\mathrm{on}\mathrm{\Gamma },$$

(32c)

$$({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(1)}-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(1)})·\mathbf{n}=({\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}-{\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)})·\mathbf{n}\mathrm{on}\mathrm{\Gamma }.$$

(32d)

Exploiting the linearity of problem (32a)–(32d), we can propose the following ansatz:

$${\varphi }_i^{(1)}={\mathrm{\Phi }}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}+{\hat{\mathrm{\Phi }}}_i{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}+{\tilde{\varphi }}_i,$$

(33a)

$${\varphi }_e^{(1)}={\mathrm{\Phi }}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}+{\hat{\mathrm{\Phi }}}_e{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}+{\tilde{\varphi }}_e,$$

(33b)

where ${\mathrm{\Phi }}_i$ , ${\mathrm{\Phi }}_e$ , ${\hat{\mathrm{\Phi }}}_i$ , and ${\hat{\mathrm{\Phi }}}_e$ are the vectors, and ${\tilde{\varphi }}_i$ and ${\tilde{\varphi }}_e$ are scalars. The auxiliary fields ${\mathrm{\Phi }}_i$ , ${\mathrm{\Phi }}_e$ , ${\hat{\mathrm{\Phi }}}_i$ , ${\hat{\mathrm{\Phi }}}_e$ , ${\tilde{\varphi }}_i$ , and ${\tilde{\varphi }}_e$ satisfy the following cell problems:

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_i)+{\nabla }_{\mathbf{y}}·{\mathsf{G}}_i=0\mathrm{in}{\mathrm{\Omega }}_i,$$

(34a)

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\hat{\mathrm{\Phi }}}_e)=0\mathrm{in}{\mathrm{\Omega }}_e,$$

(34b)

$${\mathrm{\Phi }}_i={\hat{\mathrm{\Phi }}}_e\mathrm{on}\mathrm{\Gamma },$$

(34c)

$$({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_i-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\hat{\mathrm{\Phi }}}_e)·\mathbf{n}=-{\mathsf{G}}_i·\mathbf{n}\mathrm{on}\mathrm{\Gamma },$$

(34d)

and

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\hat{\mathrm{\Phi }}}_i)=0\mathrm{in}{\mathrm{\Omega }}_i,$$

(35a)

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_e)+{\nabla }_{\mathbf{y}}·{\mathsf{G}}_e=0\mathrm{in}{\mathrm{\Omega }}_e,$$

(35b)

$${\hat{\mathrm{\Phi }}}_i={\mathrm{\Phi }}_e\mathrm{on}\mathrm{\Gamma },$$

(35c)

$$({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\hat{\mathrm{\Phi }}}_i-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_e)·\mathbf{n}={\mathsf{G}}_e·\mathbf{n}\mathrm{on}\mathrm{\Gamma },$$

(35d)

and

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\tilde{\varphi }}_i)=0\mathrm{in}{\mathrm{\Omega }}_i,$$

(36a)

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\tilde{\varphi }}_e)=0\mathrm{in}{\mathrm{\Omega }}_e,$$

(36b)

$${\tilde{\varphi }}_i-{\tilde{\varphi }}_e=V^{(1)}\mathrm{on}\mathrm{\Gamma },$$

(36c)

$$({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\tilde{\varphi }}_i)·\mathbf{n}=({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\tilde{\varphi }}_e)·\mathbf{n}\mathrm{on}\mathrm{\Gamma }.$$

(36d)

where periodic conditions also apply on the boundary $\partial \mathrm{\Omega }\backslash \mathrm{\Gamma }$ , and we require a further condition on the auxiliary fields ${\mathrm{\Phi }}_i$ , ${\mathrm{\Phi }}_e$ , ${\hat{\mathrm{\Phi }}}_i$ , ${\hat{\mathrm{\Phi }}}_e$ , ${\tilde{\varphi }}_i$ , and ${\tilde{\varphi }}_e$ for the solution to be unique. We propose zero average on their own domains, i.e.:

$$\begin{array}{c}{〈{\mathrm{\Phi }}_i〉}_i=0,{〈{\mathrm{\Phi }}_e〉}_e=0,{〈{\hat{\mathrm{\Phi }}}_i〉}_i=0,{〈{\hat{\mathrm{\Phi }}}_e〉}_e=0,\hfill \end{array}$$

(37)

Since we now have expressions for ${\varphi }_i^{(1)}$ and ${\varphi }_e^{(1)}$ from the ansatz (33a) and (33b), then it is possible to write Ohm’s law for the current densities ${\mathbf{j}}_i^{(0)}$ and ${\mathbf{j}}_e^{(0)}$ (equations (31d) and (31e)) as:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{j}}_i^{(0)}=-{\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}-{\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(1)}}\\ \multicolumn{1}{c}{=-{\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}-{\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}-{\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\hat{\mathrm{\Phi }}}_i{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}}\\ \multicolumn{1}{c}{-{\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\tilde{\varphi }}_i}\\ \multicolumn{1}{c}{=-({\mathsf{G}}_i+{\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_i){\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}-{\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\hat{\mathrm{\Phi }}}_i{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}}\\ \multicolumn{1}{c}{-{\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\tilde{\varphi }}_i}\\ \multicolumn{1}{c}{=-({\mathsf{G}}_i+{\mathsf{G}}_i{\mathsf{R}}_i){\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}-({\mathsf{G}}_i{\mathsf{Q}}_i){\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}}\\ \multicolumn{1}{c}{-{\mathsf{G}}_i{\mathbf{s}}_i,}\end{array}$$

(38)

and

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{j}}_e^{(0)}=-{\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(1)}}\\ \multicolumn{1}{c}{=-{\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\hat{\mathrm{\Phi }}}_e{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}}\\ \multicolumn{1}{c}{-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\tilde{\varphi }}_e}\\ \multicolumn{1}{c}{=-({\mathsf{G}}_e+{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_e){\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\hat{\mathrm{\Phi }}}_e{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}}\\ \multicolumn{1}{c}{-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\tilde{\varphi }}_e}\\ \multicolumn{1}{c}{=-({\mathsf{G}}_e+{\mathsf{G}}_e{\mathsf{R}}_e){\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}-({\mathsf{G}}_e{\mathsf{Q}}_e){\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}}\\ \multicolumn{1}{c}{-{\mathsf{G}}_e{\mathbf{s}}_e,}\end{array}$$

(39)

where we have used the notation:

$$\begin{array}{c}{\mathsf{R}}_i={\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_i,{\mathsf{R}}_e={\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_e,{\mathsf{Q}}_i={\nabla }_{\mathbf{y}}{\hat{\mathrm{\Phi }}}_i,\hfill \\ \multicolumn{1}{c}{{\mathsf{Q}}_e={\nabla }_{\mathbf{y}}{\hat{\mathrm{\Phi }}}_e,{\mathbf{s}}_i={\nabla }_{\mathbf{y}}{\tilde{\varphi }}_i,{\mathbf{s}}_e={\nabla }_{\mathbf{y}}{\tilde{\varphi }}_e.}\end{array}$$

(40)

We now require a balance equation for the current densities. This requires us to equate the powers of ${ϵ}^2$ in equations (22a) and (22b) and use the ${ϵ}^2$ definition of equations (22d) and (22e). That is, we need:

$$\begin{array}{c}\multicolumn{1}{c}{{\nabla }_{\mathbf{x}}·({\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)})+{\nabla }_{\mathbf{x}}·({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(1)})+{\nabla }_{\mathbf{y}}·({\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(1)})}\\ \multicolumn{1}{c}{+{\nabla }_{\mathbf{y}}·({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(2)})=\hat{\beta }({\varphi }_i^{(0)}-{\varphi }_e^{(0)}),}\end{array}$$

(41a)

$$\begin{array}{c}\multicolumn{1}{c}{{\nabla }_{\mathbf{x}}·({\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)})+{\nabla }_{\mathbf{x}}·({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(1)})+{\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(1)})}\\ \multicolumn{1}{c}{+{\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(2)})=-\hat{\beta }({\varphi }_i^{(0)}-{\varphi }_e^{(0)}),}\end{array}$$

(41b)

with the ${ϵ}^2$ expansions of Ohm’s law:

$${\mathbf{j}}_i^{(1)}=-{\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(1)}-{\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(2)},$$

(42a)

$${\mathbf{j}}_e^{(1)}=-{\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(1)}-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(2)}.$$

(42b)

We can then use equations (42a) and (42b), along with equations (31d) and (31e) in equations (41a) and (41b) to obtain:

$${\nabla }_{\mathbf{x}}·{\mathbf{j}}_i^{(0)}+{\nabla }_{\mathbf{y}}·{\mathbf{j}}_i^{(1)}=\hat{\beta }({\varphi }_i^{(0)}-{\varphi }_e^{(0)}),$$

(43a)

$${\nabla }_{\mathbf{x}}·{\mathbf{j}}_e^{(0)}+{\nabla }_{\mathbf{y}}·{\mathbf{j}}_e^{(1)}=-\hat{\beta }({\varphi }_i^{(0)}-{\varphi }_e^{(0)}),$$

(43b)

Now, we can take the integral average of the sum of equations (43a) and (43b) to obtain:

$$\begin{array}{c}\multicolumn{1}{c}{{\int }_{{\mathrm{\Omega }}_i}{\nabla }_{\mathbf{x}}·{\mathbf{j}}_i^{(0)}d\mathbf{y}+{\int }_{{\mathrm{\Omega }}_i}{\nabla }_{\mathbf{y}}·{\mathbf{j}}_i^{(1)}d\mathbf{y}-{\int }_{{\mathrm{\Omega }}_i}\hat{\beta }({\varphi }_i^{(0)}-{\varphi }_e^{(0)})d\mathbf{y}}\\ \multicolumn{1}{c}{+{\int }_{{\mathrm{\Omega }}_e}{\nabla }_{\mathbf{x}}·j_e^{(0)}d\mathbf{y}+{\int }_{{\mathrm{\Omega }}_e}{\nabla }_{\mathbf{y}}·{\mathbf{j}}_e^{(1)}d\mathbf{y}+{\int }_{{\mathrm{\Omega }}_e}\hat{\beta }({\varphi }_i^{(0)}-{\varphi }_e^{(0)})d\mathbf{y}=0.}\end{array}$$

(44)

This can be rewritten as the following when applying macroscopic uniformity to the first and fourth integrals and the divergence theorem to the second and fifth integrals:

$$\begin{array}{c}\multicolumn{1}{c}{{\nabla }_{\mathbf{x}}·〈{\mathbf{j}}_i^{(0)}{〉}_i+{\nabla }_{\mathbf{x}}·〈{\mathbf{j}}_e^{(0)}{〉}_e+{\int }_{\partial {\mathrm{\Omega }}_i\backslash \mathrm{\Gamma }}{\mathbf{j}}_i^{(1)}{\mathbf{n}}_{\partial {\mathrm{\Omega }}_i}\mathrm{dS}+}\\ \multicolumn{1}{c}{{\int }_{\mathrm{\Gamma }}{\mathbf{j}}_i^{(1)}\mathbf{n}\mathrm{dS}+{\int }_{\partial {\mathrm{\Omega }}_e\backslash \mathrm{\Gamma }}{\mathbf{j}}_e^{(1)}{\mathbf{n}}_{\partial {\mathrm{\Omega }}_e}\mathrm{dS}-{\int }_{\mathrm{\Gamma }}{\mathbf{j}}_e^{(1)}\mathbf{n}\mathrm{dS}-}\\ \multicolumn{1}{c}{{\int }_{{\mathrm{\Omega }}_i}\hat{\beta }({\varphi }_i^{(0)}-{\varphi }_e^{(0)})\mathrm{d}\mathbf{y}+{\int }_{{\mathrm{\Omega }}_e}\hat{\beta }({\varphi }_i^{(0)}-{\varphi }_e^{(0)})\mathrm{d}\mathbf{y}=0.}\end{array}$$

(45)

Due to periodicity, the terms on the external boundaries will cancel out, and using the ${ϵ}^2$ coefficient of equation (22c) with equations (42a) and (42b), we see that the terms on $\mathrm{\Gamma }$ will also disappear leaving:

$${\nabla }_{\mathbf{x}}·〈{\mathbf{j}}_i^{(0)}{〉}_i+{\nabla }_{\mathbf{x}}·〈{\mathbf{j}}_e^{(0)}{〉}_e-{\int }_{{\mathrm{\Omega }}_i}\hat{\beta }({\varphi }_i^{(0)}-{\varphi }_e^{(0)})\mathrm{d}\mathbf{y}+{\int }_{{\mathrm{\Omega }}_e}\hat{\beta }({\varphi }_i^{(0)}-{\varphi }_e^{(0)})\mathrm{d}\mathbf{y}=0.$$

(46)

Since we have equations (24f) and (29a) and (29b), we can rewrite equation (46) as:

$${\nabla }_{\mathbf{x}}·〈{\mathbf{j}}_i^{(0)}{〉}_i+{\nabla }_{\mathbf{x}}·〈{\mathbf{j}}_e^{(0)}{〉}_e-\hat{\beta }{V}^{(0)}(|{\mathrm{\Omega }}_e|-|{\mathrm{\Omega }}_i|)=0,$$

(47)

where

$$\begin{array}{c}\multicolumn{1}{c}{〈{\mathbf{j}}_i^{(0)}{〉}_i=-〈{\mathsf{G}}_i+{\mathsf{G}}_i{\mathsf{R}}_i{〉}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}-〈{\mathsf{G}}_i{\mathsf{Q}}_i{〉}_i{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}}\\ \multicolumn{1}{c}{-〈{\mathsf{G}}_i{\mathbf{s}}_i{〉}_i}\end{array}$$

(48a)

$$\begin{array}{c}\multicolumn{1}{c}{〈{\mathbf{j}}_e^{(0)}{〉}_e=-〈{\mathsf{G}}_e+{\mathsf{G}}_e{\mathsf{R}}_e{〉}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}-〈{\mathsf{G}}_e{\mathsf{Q}}_e{〉}_e{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}}\\ \multicolumn{1}{c}{-〈{\mathsf{G}}_e{\mathbf{s}}_e{〉}_e.}\end{array}$$

(48b)

3.5. Problem for the elastic displacements ${\mathbf{u}}^{(1)}$ and ${\mathbf{w}}^{(1)}$

We now require a problem for ${\mathbf{u}}^{(1)}$ and ${\mathbf{w}}^{(1)}$ , so by taking equations (30a), (30b), and (31k) with equations (31i) and (31j) along with the interface condition (31l), we can write:

$${\nabla }_{\mathbf{y}}·({\mathbb{C}}_i{\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)}))+{\nabla }_{\mathbf{y}}·({\mathbb{C}}_i{\xi }_{\mathbf{y}}({\mathbf{u}}^{(1)}))=0,$$

(49a)

$${\nabla }_{\mathbf{y}}·({\mathbb{C}}_e{\xi }_{\mathbf{x}}({\mathbf{w}}^{(0)}))+{\nabla }_{\mathbf{y}}·({\mathbb{C}}_e{\xi }_{\mathbf{y}}({\mathbf{w}}^{(1)}))=0,$$

(49b)

$${\mathbf{u}}^{(1)}={\mathbf{w}}^{(1)},$$

(49c)

$$({\mathbb{C}}_i{\xi }_{\mathbf{y}}({\mathbf{u}}^{(1)})-{\mathbb{C}}_e{\xi }_{\mathbf{y}}({\mathbf{w}}^{(1)}))·\mathbf{n}=({\mathbb{C}}_e{\xi }_{\mathbf{x}}({\mathbf{w}}^{(0)})-{\mathbb{C}}_i{\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)}))·\mathbf{n}.$$

(49d)

Since we have condition (26), ${\mathbf{u}}^{(0)}={\mathbf{w}}^{(0)}$ , then we can replace the ${\mathbf{w}}^{(0)}$ s in equations (49a)–(49d) to obtain:

$${\nabla }_{\mathbf{y}}·({\mathbb{C}}_i{\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)}))+{\nabla }_{\mathbf{y}}·({\mathbb{C}}_i{\xi }_{\mathbf{y}}({\mathbf{u}}^{(1)}))=0,$$

(50a)

$${\nabla }_{\mathbf{y}}·({\mathbb{C}}_e{\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)}))+{\nabla }_{\mathbf{y}}·({\mathbb{C}}_e{\xi }_{\mathbf{y}}({\mathbf{w}}^{(1)}))=0,$$

(50b)

$${\mathbf{u}}^{(1)}={\mathbf{w}}^{(1)},$$

(50c)

$$({\mathbb{C}}_i{\xi }_{\mathbf{y}}({\mathbf{u}}^{(1)})-{\mathbb{C}}_e{\xi }_{\mathbf{y}}({\mathbf{w}}^{(1)}))·\mathbf{n}=({\mathbb{C}}_e-{\mathbb{C}}_i){\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)})·\mathbf{n}.$$

(50d)

We can exploit the linearity of the problem to propose the ansatz:

$${\mathbf{u}}^{(1)}={\mathcal{A}}_i{\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)}),$$

(51a)

$${\mathbf{w}}^{(1)}={\mathcal{A}}_e{\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)}),$$

(51b)

where ${\mathcal{A}}_i$ and ${\mathcal{A}}_e$ are the third rank tensors which solve the following cell problem:

$${\nabla }_{\mathbf{y}}·({\mathbb{C}}_i{\xi }_{\mathbf{y}}({\mathcal{A}}_i))+{\nabla }_{\mathbf{y}}·{\mathbb{C}}_i=0,$$

(52a)

$${\nabla }_{\mathbf{y}}·({\mathbb{C}}_e{\xi }_{\mathbf{y}}({\mathcal{A}}_e))+{\nabla }_{\mathbf{y}}·{\mathbb{C}}_e=0,$$

(52b)

$${\mathcal{A}}_i={\mathcal{A}}_e,$$

(52c)

$$({\mathbb{C}}_i{\xi }_{\mathbf{y}}({\mathcal{A}}_i)-{\mathbb{C}}_e{\xi }_{\mathbf{y}}({\mathcal{A}}_e))·\mathbf{n}=({\mathbb{C}}_e-{\mathbb{C}}_i)·\mathbf{n}.$$

(52d)

For uniqueness of solution, we require an additional condition on the auxiliary tensors ${\mathcal{A}}_i$ and ${\mathcal{A}}_e$ , so we propose:

$$〈{\mathcal{A}}_i{〉}_i=0,\mathrm{and}〈{\mathcal{A}}_e{〉}_e=0.$$

(53)

Using the expressions, we have just obtained for ${\mathbf{u}}^{(1)}$ and ${\mathbf{w}}^{(1)}$ , we can write the leading order stresses (31i) and (31j) as:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathsf{T}}_i^{(0)}={\mathbb{C}}_i{\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)})+{\mathbb{C}}_i{\xi }_{\mathbf{y}}({\mathcal{A}}_i){\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)})}\\ \multicolumn{1}{c}{=({\mathbb{C}}_i+{\mathbb{C}}_i{\mathbb{M}}_i){\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)}),}\end{array}$$

(54)

and

$$\begin{array}{c}\multicolumn{1}{c}{{\mathsf{T}}_e^{(0)}={\mathbb{C}}_e{\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)})+{\mathbb{C}}_e{\xi }_{\mathbf{y}}({\mathcal{A}}_e){\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)})}\\ \multicolumn{1}{c}{=({\mathbb{C}}_e+{\mathbb{C}}_e{\mathbb{M}}_e){\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)}),}\end{array}$$

(55)

where we have used the notation:

$${\mathbb{M}}_i={\xi }_{\mathbf{y}}({\mathcal{A}}_i),{\mathbb{M}}_e={\xi }_{\mathbf{y}}({\mathcal{A}}_e).$$

(56)

We now want to consider a balance equation that takes into consideration each compartment. Taking the integral average of equations (31g) and (31h) gives:

$$\begin{array}{c}\multicolumn{1}{c}{{\int }_{{\mathrm{\Omega }}_i}({\nabla }_{\mathbf{x}}·{\mathsf{T}}_i^{(0)})\mathrm{d}\mathbf{y}+{\int }_{{\mathrm{\Omega }}_i}({\nabla }_{\mathbf{y}}·{\mathsf{T}}_i^{(1)})\mathrm{d}\mathbf{y}-{\int }_{{\mathrm{\Omega }}_i}\hat{K}({\mathbf{u}}^{(0)}-{\mathbf{w}}^{(0)})\mathrm{d}\mathbf{y}}\\ \multicolumn{1}{c}{+{\int }_{{\mathrm{\Omega }}_i}({\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}+{\int }_{{\mathrm{\Omega }}_i}({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(1)}\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}}\\ \multicolumn{1}{c}{+{\int }_{{\mathrm{\Omega }}_e}({\nabla }_{\mathbf{x}}·{\mathsf{T}}_e^{(0)})\mathrm{d}\mathbf{y}+{\int }_{{\mathrm{\Omega }}_e}({\nabla }_{\mathbf{y}}·{\mathsf{T}}_e^{(1)})\mathrm{d}\mathbf{y}}\\ \multicolumn{1}{c}{+{\int }_{{\mathrm{\Omega }}_e}\hat{K}({\mathbf{u}}^{(0)}-{\mathbf{w}}^{(0)})\mathrm{d}\mathbf{y}+{\int }_{{\mathrm{\Omega }}_e}({\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}}\\ \multicolumn{1}{c}{+{\int }_{{\mathrm{\Omega }}_e}({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(1)}\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}=0.}\end{array}$$

(57)

Since we have that ${\mathbf{u}}^{(0)}={\mathbf{w}}^{(0)}$ (26), we can cancel the third and eighth integral, we can then apply the divergence theorem to the second and seventh integral and use the expressions (33a) and (33b) that we have for ${\varphi }_i^{(1)}$ and ${\varphi }_e^{(1)}$ to write:

$$\begin{array}{c}\multicolumn{1}{c}{{\nabla }_{\mathbf{x}}·〈{\mathsf{T}}_i^{(0)}{〉}_i+{\nabla }_{\mathbf{x}}·〈{\mathsf{T}}_e^{(0)}{〉}_e+{\int }_{\partial {\mathrm{\Omega }}_i\backslash \mathrm{\Gamma }}({\mathsf{T}}_i^{(1)}·{\mathbf{n}}_{\partial {\mathrm{\Omega }}_i})\mathrm{dS}}\\ \multicolumn{1}{c}{+{\int }_{\mathrm{\Gamma }}({\mathsf{T}}_i^{(1)}·\mathbf{n})\mathrm{dS}+{\int }_{\partial {\mathrm{\Omega }}_e\backslash \mathrm{\Gamma }}({\mathsf{T}}_e^{(1)}{\mathbf{n}}_{\partial {\mathrm{\Omega }}_e})\mathrm{dS}}\\ \multicolumn{1}{c}{+{\int }_{\mathrm{\Gamma }}({\mathsf{T}}_e^{(1)}·\mathbf{n})\mathrm{dS}+{\int }_{{\mathrm{\Omega }}_i}({\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}}\\ \multicolumn{1}{c}{+{\int }_{{\mathrm{\Omega }}_e}({\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}}\\ \multicolumn{1}{c}{+{\int }_{{\mathrm{\Omega }}_i}({\mathsf{G}}_i{\nabla }_{\mathbf{y}}({\mathrm{\Phi }}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}+{\hat{\mathrm{\Phi }}}_i{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}+{\tilde{\varphi }}_i)\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}}\\ \multicolumn{1}{c}{+{\int }_{{\mathrm{\Omega }}_e}({\mathsf{G}}_e{\nabla }_{\mathbf{y}}({\mathrm{\Phi }}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}+{\hat{\mathrm{\Phi }}}_e{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}+{\tilde{\varphi }}_e)\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}=0.}\end{array}$$

(58)

The terms of the external boundaries cancel due to periodicity and the terms on $\mathrm{\Gamma }$ cancel due to equation (31k) so that we have:

$$\begin{array}{c}\multicolumn{1}{c}{{\nabla }_{\mathbf{x}}·〈{\mathsf{T}}_i^{(0)}{〉}_i+{\nabla }_{\mathbf{x}}·〈{\mathsf{T}}_e^{(0)}{〉}_e+{\int }_{{\mathrm{\Omega }}_i}({\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}}\\ \multicolumn{1}{c}{+{\int }_{{\mathrm{\Omega }}_e}({\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}+{\int }_{{\mathrm{\Omega }}_i}({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}}\\ \multicolumn{1}{c}{+{\int }_{{\mathrm{\Omega }}_i}({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\hat{\mathrm{\Phi }}}_i{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}+{\int }_{{\mathrm{\Omega }}_i}({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\tilde{\varphi }}_i\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}}\\ \multicolumn{1}{c}{+{\int }_{{\mathrm{\Omega }}_e}({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}}\\ \multicolumn{1}{c}{+{\int }_{{\mathrm{\Omega }}_e}({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\hat{\mathrm{\Phi }}}_e{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}}\\ \multicolumn{1}{c}{+{\int }_{{\mathrm{\Omega }}_e}({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\tilde{\varphi }}_e\times {\mathbf{B}}^{(0)})\mathrm{d}\mathbf{y}=0.}\end{array}$$

(59)

Collecting integrals together and using the notation (40), we have:

$$\begin{array}{c}\multicolumn{1}{c}{{\nabla }_{\mathbf{x}}·〈{\mathsf{T}}_i^{(0)}{〉}_i+{\nabla }_{\mathbf{x}}·〈{\mathsf{T}}_e^{(0)}{〉}_e+〈{\mathsf{G}}_i{\mathsf{R}}_i+{\mathsf{G}}_i{〉}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}\times 〈{\mathbf{B}}^{(0)}{〉}_i}\\ \multicolumn{1}{c}{+〈{\mathsf{G}}_e{\mathsf{R}}_e+{\mathsf{G}}_e{〉}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}\times 〈{\mathbf{B}}^{(0)}{〉}_e}\\ \multicolumn{1}{c}{+〈{\mathsf{G}}_i{\mathsf{Q}}_i{〉}_i{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}\times 〈{\mathbf{B}}^{(0)}{〉}_i+〈{\mathsf{G}}_i{\mathbf{s}}_i\times {\mathbf{B}}^{(0)}{〉}_i}\\ \multicolumn{1}{c}{+〈{\mathsf{G}}_e{\mathsf{Q}}_{\mathsf{e}}{〉}_e{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}\times 〈{\mathbf{B}}^{(0)}{〉}_e+〈{\mathsf{G}}_e{\mathbf{s}}_e\times {\mathbf{B}}^{(0)}{〉}_e=0.}\end{array}$$

(60)

This can then be written as:

$$\begin{array}{c}\multicolumn{1}{c}{{\nabla }_{\mathbf{x}}·〈{\mathsf{T}}_i^{(0)}{〉}_i+{\nabla }_{\mathbf{x}}·〈{\mathsf{T}}_e^{(0)}{〉}_e}\\ \multicolumn{1}{c}{=-〈{\mathsf{G}}_i{\mathsf{R}}_i+{\mathsf{G}}_i{〉}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}\times 〈{\mathbf{B}}^{(0)}{〉}_i}\\ \multicolumn{1}{c}{-〈{\mathsf{G}}_e{\mathsf{R}}_e+{\mathsf{G}}_e{〉}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}\times 〈{\mathbf{B}}^{(0)}{〉}_e}\\ \multicolumn{1}{c}{-〈{\mathsf{G}}_i{\mathsf{Q}}_i{〉}_i{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}\times 〈{\mathbf{B}}^{(0)}{〉}_i-〈{\mathsf{G}}_e{\mathsf{Q}}_{\mathsf{e}}{〉}_e{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}\times 〈{\mathbf{B}}^{(0)}{〉}_e}\\ \multicolumn{1}{c}{-〈{\mathsf{G}}_i{\mathbf{s}}_i\times {\mathbf{B}}^{(0)}{〉}_i-〈{\mathsf{G}}_e{\mathbf{s}}_e\times {\mathbf{B}}^{(0)}{〉}_e.}\end{array}$$

(61)

We have now derived all the equations required to be able to state our macroscale model.

4. Macroscale model

The macroscale equations describe the effective behavior of the heart in terms of the leading order elastic displacement ${\mathbf{u}}^{(0)}$ and the leading order electric potentials ${\varphi }_i^{(0)}$ and ${\varphi }_e^{(0)}$ . The model is given by:

$${\nabla }_{\mathbf{x}}·〈{\mathbf{j}}_i^{(0)}{〉}_i+{\nabla }_{\mathbf{x}}·〈{\mathbf{j}}_e^{(0)}{〉}_e=\hat{\beta }{V}^{(0)}(|{\mathrm{\Omega }}_e|-|{\mathrm{\Omega }}_i|),$$

(62a)

$$\begin{array}{c}\multicolumn{1}{c}{{\nabla }_{\mathbf{x}}·〈{\mathsf{T}}_i^{(0)}{〉}_i+{\nabla }_{\mathbf{x}}·〈{\mathsf{T}}_e^{(0)}{〉}_e=-〈{\mathsf{G}}_i{\mathsf{R}}_i+{\mathsf{G}}_i{〉}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}\times 〈{\mathbf{B}}^{(0)}{〉}_i}\\ \multicolumn{1}{c}{-〈{\mathsf{G}}_e{\mathsf{R}}_e+{\mathsf{G}}_e{〉}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}\times 〈{\mathbf{B}}^{(0)}{〉}_e-〈{\mathsf{G}}_i{\mathsf{Q}}_i{〉}_i{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}\times 〈{\mathbf{B}}^{(0)}{〉}_i}\\ \multicolumn{1}{c}{-〈{\mathsf{G}}_e{\mathsf{Q}}_{\mathsf{e}}{〉}_e{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}\times 〈{\mathbf{B}}^{(0)}{〉}_e-〈{\mathsf{G}}_i{\mathbf{s}}_i\times {\mathbf{B}}^{(0)}{〉}_i}\\ \multicolumn{1}{c}{-〈{\mathsf{G}}_e{\mathbf{s}}_e\times {\mathbf{B}}^{(0)}{〉}_e,}\end{array}$$

(62b)

$${\varphi }_i^{(0)}-{\varphi }_e^{(0)}=V^{(0)},$$

(62c)

where we have the averaged leading order current densities:

$$\begin{array}{c}\multicolumn{1}{c}{〈{\mathbf{j}}_i^{(0)}{〉}_i=-〈{\mathsf{G}}_i+{\mathsf{G}}_i{\mathsf{R}}_i{〉}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}-〈{\mathsf{G}}_i{\mathsf{Q}}_i{〉}_i{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}}\\ \multicolumn{1}{c}{-〈{\mathsf{G}}_i{\mathbf{s}}_i{〉}_i,}\end{array}$$

(63a)

$$\begin{array}{c}\multicolumn{1}{c}{〈{\mathbf{j}}_e^{(0)}{〉}_e=-〈{\mathsf{G}}_e+{\mathsf{G}}_e{\mathsf{R}}_e{〉}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}-〈{\mathsf{G}}_e{\mathsf{Q}}_e{〉}_e{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)}}\\ \multicolumn{1}{c}{-〈{\mathsf{G}}_e{\mathbf{s}}_e{〉}_e,}\end{array}$$

(63b)

and the averaged leading order solid stresses:

$$〈{\mathsf{T}}_i^{(0)}{〉}_i=〈{\mathbb{C}}_i+{\mathbb{C}}_i{\mathbb{M}}_i{〉}_i{\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)}),$$

(64a)

$$〈{\mathsf{T}}_e^{(0)}{〉}_i=〈{\mathbb{C}}_e+{\mathbb{C}}_e{\mathbb{M}}_e{〉}_e{\xi }_{\mathbf{x}}({\mathbf{u}}^{(0)}).$$

(64b)

The novel model comprises the balance equation for the leading order current densities (62a). The current densities (63a) and (63b) comprise both the electric fields of each compartment premultiplied by second rank tensors that are to be obtained by solving the cell problems (34a)–(34d) and (35a)–(35d) and a vector term that arises from premultiplying the solutions to the cell problem (36a)–(36d) by the conductivity tensors. These coefficients arising from the cell problem solutions account for the differences in the electric potentials at each point in the microstructure and encode these in the model. We also have the balance equation for the solid stresses (62b), where the stresses are (64a) and (64b) with the standard effective elasticity tensors ${\mathbb{C}}_i+{\mathbb{C}}_i{\mathbb{M}}_i$ and ${\mathbb{C}}_e+{\mathbb{C}}_e{\mathbb{M}}_e$ found for elastic composite problems such as in Penta and Gerisch [26, 27]. These coefficients encode the difference in elastic properties at different points in the microstructure. Equation (62b) also has terms that show how the electric potentials and Lorentz forces influence and drive the elastic deformations of the material. These terms again arise by solving the cell problems (34a)–(34d) and (35a)–(35d) and (36a)–(36d) which encode the difference in electric potentials in the microstructure.

5. Limit case and scheme to solve the model

In this section, we will consider some modifications to the model setup and the outcome they will have on the cell problems and overall model. We will also propose a scheme that can be used to solve the macroscale model.

5.1. Limit case here $V=0$

Throughout this work, we have assumed that there is a potential drop $V$ across the interface between the intracellular and extracellular domains. We can here, however, consider what will happen to the model (62a)–(62c) and the corresponding cell problems in the case in which $V=0$ .

Under this assumption, the macroscale equation (62c) will become:

$${\varphi }_i^{(0)}={\varphi }_e^{(0)}.$$

(65)

We can then consider the problem (32a)–(32d) assuming $V=0$ , i.e.:

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)})+{\nabla }_{\mathbf{y}}·({\mathrm{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(1)})=0\mathrm{in}{\mathrm{\Omega }}_i,$$

(66a)

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)})+{\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(1)})=0\mathrm{in}{\mathrm{\Omega }}_e,$$

(66b)

$${\varphi }_i^{(1)}={\varphi }_e^{(1)}\mathrm{on}\mathrm{\Gamma },$$

(66c)

$$({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(1)}-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(1)})·\mathbf{n}=({\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }_e^{(0)}-{\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }_i^{(0)})·\mathbf{n}\mathrm{on}\mathrm{\Gamma }.$$

(66d)

We can use the fact that ${\varphi }_i^{(0)}$ and ${\varphi }_e^{(0)}$ do not depend on the microscale variable $y$ and using equation (24f) under the assumption $V=0$ , we have:

$${\varphi }_i^{(0)}(x,t)={\varphi }_e^{(0)}(x,t)={\varphi }^{(0)}(x,t).$$

(67)

We can then propose the following ansatz:

$${\varphi }_i^{(1)}={\mathrm{\Phi }}_i{\nabla }_{\mathbf{x}}{\varphi }^{(0)},$$

(68a)

$${\varphi }_e^{(1)}={\mathrm{\Phi }}_e{\nabla }_{\mathbf{x}}{\varphi }^{(0)},$$

(68b)

where ${\mathrm{\Phi }}_i$ and ${\mathrm{\Phi }}_e$ are the vectors. The auxiliary fields ${\mathrm{\Phi }}_i$ , ${\mathrm{\Phi }}_e$ , satisfy the following cell problem:

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_i)+{\nabla }_{\mathbf{y}}·{\mathsf{G}}_i=0\mathrm{in}{\mathrm{\Omega }}_i,$$

(69a)

$${\nabla }_{\mathbf{y}}·({\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_e)+{\nabla }_{\mathbf{y}}·{\mathsf{G}}_e=0\mathrm{in}{\mathrm{\Omega }}_e,$$

(69b)

$${\mathrm{\Phi }}_i={\mathrm{\Phi }}_e\mathrm{on}\mathrm{\Gamma },$$

(69c)

$$({\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_i-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\mathrm{\Phi }}_e)·\mathbf{n}=({\mathsf{G}}_e-{\mathsf{G}}_i)·\mathbf{n}\mathrm{on}\mathrm{\Gamma }.$$

(69d)

This is a vector problem which is driven by the difference in the conductivity tensors on the normal to the interface between the phases.

We can also use this ansatz to determine the current densities ${\mathbf{j}}_i^{(0)}$ and ${\mathbf{j}}_e^{(0)}$ in the case where $V=0$ . We can write equations (31d) and (31e):

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{j}}_i^{(0)}=-{\mathsf{G}}_i{\nabla }_{\mathbf{x}}{\varphi }^{(0)}-{\mathsf{G}}_i{\nabla }_{\mathbf{y}}{\varphi }_i^{(1)}}\\ \multicolumn{1}{c}{=-({\mathsf{G}}_i+{\mathsf{G}}_i{\mathsf{R}}_i){\nabla }_{\mathbf{x}}{\varphi }^{(0)},}\end{array}$$

(70)

and

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{j}}_e^{(0)}=-{\mathsf{G}}_e{\nabla }_{\mathbf{x}}{\varphi }^{(0)}-{\mathsf{G}}_e{\nabla }_{\mathbf{y}}{\varphi }_e^{(1)}}\\ \multicolumn{1}{c}{=-({\mathsf{G}}_e+{\mathsf{G}}_e{\mathsf{R}}_e){\nabla }_{\mathbf{x}}{\varphi }^{(0)}.}\end{array}$$

(71)

The macroscale balance equation for the current densities (62a) can be simplified to:

$${\nabla }_{\mathbf{x}}·〈{\mathbf{j}}_i^{(0)}{〉}_i+{\nabla }_{\mathbf{x}}·〈{\mathbf{j}}_e^{(0)}{〉}_e=0.$$

(72)

The leading order solid stresses will remain unchanged as in equations (64a) and (64b) while the balance of the solid stresses (62b) will simplify to:

$$\begin{array}{c}\multicolumn{1}{c}{{\nabla }_{\mathbf{x}}·〈{\mathsf{T}}_i^{(0)}{〉}_i+{\nabla }_{\mathbf{x}}·〈{\mathsf{T}}_e^{(0)}{〉}_e}\\ \multicolumn{1}{c}{=-〈{\mathsf{G}}_i{\mathsf{R}}_i+{\mathsf{G}}_i{〉}_i{\nabla }_{\mathbf{x}}{\varphi }^{(0)}\times 〈{\mathbf{B}}^{(0)}{〉}_i}\\ \multicolumn{1}{c}{-〈{\mathsf{G}}_e{\mathsf{R}}_e+{\mathsf{G}}_e{〉}_e{\nabla }_{\mathbf{x}}{\varphi }^{(0)}\times 〈{\mathbf{B}}^{(0)}{〉}_e.}\end{array}$$

(73)

This means that we can summarize the simplified macroscale model using equations (72), (73), and (65) with amended leading order current densities (70) and (71) and the unchanged leading order solid stresses (64a) and (64b).

5.2. Scheme for solving the macroscale model

We now propose a scheme that can be implemented to solve the macroscale model (62a)–(62c). We progress in a step-by-step manner explaining how to find the effective model coefficients. We also explain how these coefficients will be used in later steps when solving the macroscale model (62a)–(62c). The model coefficients encode the structural details such as geometry, elastic, and electrical properties. Since we have assumed macroscopic uniformity of the material (see Remark 2) and assumed that the two scales are fully decoupled, then we can propose the following steps to solve the model. The process is as follows:

• 1. The first step is to set the original material properties of the intracellular and the extracellular domains at the microscale. We choose to make the assumption that each of the domains is isotropic. This means we require two parameters for each domain. These are two independent elastic constants that can be Poisson’s ratio and Young’s modulus (or alternatively the Lamé constants). We could, however, not make this assumption and just provide the elasticity tensor with up to 81 components depending on the symmetry that exists in each phase. We also must fix the original electrical properties such as the conductivity tensors in each phase with up to nine components and the potential drop $V$ across the interface. The last property we need to fix is the magnetic field $\mathbf{B}$ .

• 2. We must determine the microscale geometry; this includes defining the specific geometry of a single periodic cell and the volume of each of the phases.

• 3. We now begin the process that allows us to determine the macroscale model coefficients. We will begin with the elastic coefficients. We are able to solve the elastic-type cell problem (52a)–(52d) which has the solution auxiliary tensors ${\mathbb{M}}_e$ and ${\mathbb{M}}_i$ . These tensors then appear in the macroscale model coefficients, such as in the leading order stress. The cell problem to be solved is, in components:

$$\frac{\partial }{\partial y_{\delta }}({(C_i)}_{\alpha \delta \tau \kappa }{\xi }_{\tau \kappa }^{\gamma \nu }({\mathcal{A}}_i))+\frac{\partial {(C_i)}_{\alpha \delta \gamma \nu }}{\partial y_{\delta }}=0\mathrm{in}{\mathrm{\Omega }}_i,$$

(74a)

$$\frac{\partial }{\partial y_{\delta }}({(C_e)}_{\alpha \delta \tau \kappa }{\xi }_{\tau \kappa }^{\gamma \nu }({\mathcal{A}}_e))+\frac{\partial {(C_e)}_{\alpha \delta \gamma \nu }}{\partial y_{\delta }}=0\mathrm{in}{\mathrm{\Omega }}_e,$$

(74b)

$$({\mathcal{A}}_i{)}_{\alpha \gamma \nu }=({\mathcal{A}}_e{)}_{\alpha \gamma \nu }\mathrm{on}\mathrm{\Gamma },$$

(74c)

$$(C_i{)}_{\alpha \delta \tau \kappa }{\xi }_{\tau \kappa }^{\gamma \nu }({\mathcal{A}}_i)n_{\delta }-(C_e{)}_{\alpha \delta \tau \kappa }{\xi }_{\tau \kappa }^{\gamma \nu }({\mathcal{A}}_e)n_{\delta }=(C_e-C_i{)}_{\alpha \delta \gamma \nu }{n}_{\delta }\mathrm{on}\mathrm{\Gamma }.$$

(74d)

The solution to problem (74a)–(74d) is found by fixing the couple of indices $\tau ,\kappa =1,2,3$ such that we are solving six elastic-type cell problems. When we do this, we obtain the strains ${\xi }_{\tau \kappa }^{\gamma \nu }({\mathcal{A}}_i)$ and ${\xi }_{\tau \kappa }^{\gamma \nu }({\mathcal{A}}_e)$ so that for each fixed couple of indices $\tau$ , $\kappa$ , we have a linear elastic problem. For other examples of where this process has been carried out and utilized (see also Miller and Penta [35] and Dehghani [46, 47]). We have used the notation:

$$\begin{array}{c}{\xi }_{\tau \kappa }^{\gamma \nu }({\mathcal{A}}_i)=\frac{1}{2}(\frac{\partial {({\mathcal{A}}_i)}_{\tau \gamma \nu }}{\partial y_{\kappa }}+\frac{\partial {({\mathcal{A}}_i)}_{\kappa \gamma \nu }}{\partial y_{\tau }});\hfill \\ \multicolumn{1}{c}{{\xi }_{\tau \kappa }^{\gamma \nu }({\mathcal{A}}_e)=\frac{1}{2}(\frac{\partial {({\mathcal{A}}_e)}_{\tau \gamma \nu }}{\partial y_{\kappa }}+\frac{\partial {({\mathcal{A}}_e)}_{\kappa \gamma \nu }}{\partial y_{\tau }}).}\end{array}$$

(75)

• 4. We now will solve the problems that will allow us to determine the tensors and vectors that contribute to the coefficients in the electrical current densities. We have the problems in components:

$$\frac{\partial }{\partial y_{\delta }}({({\mathsf{G}}_i)}_{\alpha \delta }{\xi }_{\alpha \delta }({\mathrm{\Phi }}_i))+\frac{\partial {({\mathsf{G}}_i)}_{\alpha \delta }}{\partial y_{\delta }}=0\mathrm{in}{\mathrm{\Omega }}_i,$$

(76a)

$$\frac{\partial }{\partial y_{\delta }}({({\mathsf{G}}_e)}_{\alpha \delta }{\xi }_{\alpha \delta }({\hat{\mathrm{\Phi }}}_e))=\mathrm{in}{\mathrm{\Omega }}_e,$$

(76b)

$$({\mathrm{\Phi }}_i{)}_{\alpha }=({\hat{\mathrm{\Phi }}}_e{)}_{\alpha }\mathrm{on}\mathrm{\Gamma },$$

(76c)

$$(({\mathsf{G}}_i{)}_{\alpha \delta }{\xi }_{\alpha \delta }({\mathrm{\Phi }}_i)-({\mathsf{G}}_e{)}_{\alpha \delta }{\xi }_{\alpha \delta }({\hat{\mathrm{\Phi }}}_e))n_{\delta }=-({\mathsf{G}}_i{)}_{\alpha \delta }{n}_{\delta }\mathrm{on}\mathrm{\Gamma },$$

(76d)

and

$$\frac{\partial }{\partial y_{\delta }}(({\mathsf{G}}_i{)}_{\alpha \delta }{\xi }_{\alpha \delta }({\hat{\mathrm{\Phi }}}_i))=0\mathrm{in}{\mathrm{\Omega }}_i,$$

(77a)

$$\frac{\partial }{\partial y_{\delta }}(({\mathsf{G}}_e{)}_{\alpha \delta }{\xi }_{\alpha \delta }({\mathrm{\Phi }}_e))+\frac{\partial {({\mathsf{G}}_e)}_{\alpha \delta }}{\partial y_{\delta }}=0\mathrm{in}{\mathrm{\Omega }}_e,$$

(77b)

$$({\hat{\mathrm{\Phi }}}_i{)}_{\alpha }=({\mathrm{\Phi }}_e{)}_{\alpha }\mathrm{on}\mathrm{\Gamma },$$

(77c)

$$(({\mathsf{G}}_i{)}_{\alpha \delta }{\xi }_{\alpha \delta }({\hat{\mathrm{\Phi }}}_i)-({\mathsf{G}}_e{)}_{\alpha \delta }{\xi }_{\alpha \delta }({\hat{\mathrm{\Phi }}}_e))n_{\delta }=({\mathsf{G}}_e{)}_{\alpha \delta }{n}_{\delta }\mathrm{on}\mathrm{\Gamma },$$

(77d)

and

$$\frac{\partial }{\partial y_{\delta }}\left({({\mathsf{G}}_i)}_{\alpha \delta }\frac{\partial ({\tilde{\varphi }}_i)}{\partial y_{\delta }}\right)=0\mathrm{in}{\mathrm{\Omega }}_i,$$

(78a)

$$\frac{\partial }{\partial y_{\delta }}\left({({\mathsf{G}}_e)}_{\alpha \delta }\frac{\partial ({\tilde{\varphi }}_e)}{\partial y_{\delta }}\right)=0\mathrm{in}{\mathrm{\Omega }}_e,$$

(78b)

$$({\tilde{\varphi }}_i)-({\tilde{\varphi }}_e)=V^{(1)}\mathrm{on}\mathrm{\Gamma },$$

(78c)

$$\left({({\mathsf{G}}_i)}_{\alpha \delta }\frac{\partial ({\tilde{\varphi }}_i)}{\partial y_{\delta }}\right)n_{\alpha }=\left({({\mathsf{G}}_e)}_{\alpha \delta }\frac{\partial ({\tilde{\varphi }}_e)}{\partial y_{\delta }}\right)n_{\alpha }\mathrm{on}\mathrm{\Gamma },$$

(78d)

where we have used the notation:

$${\xi }_{\alpha \delta }({\mathrm{\Phi }}_i)=\frac{1}{2}(\frac{\partial {({\mathrm{\Phi }}_i)}_{\alpha }}{\partial y_{\delta }}+\frac{\partial {({\mathrm{\Phi }}_i)}_{\delta }}{\partial y_{\alpha }});$$

(79a)

$${\xi }_{\alpha \delta }({\hat{\mathrm{\Phi }}}_i)=\frac{1}{2}(\frac{\partial {({\hat{\mathrm{\Phi }}}_i)}_{\alpha }}{\partial y_{\delta }}+\frac{\partial {({\hat{\mathrm{\Phi }}}_i)}_{\delta }}{\partial y_{\alpha }}),$$

(79b)

$${\xi }_{\alpha \delta }({\mathrm{\Phi }}_e)=\frac{1}{2}(\frac{\partial {({\mathrm{\Phi }}_e)}_{\alpha }}{\partial y_{\delta }}+\frac{\partial {({\mathrm{\Phi }}_e)}_{\delta }}{\partial y_{\alpha }});$$

(79c)

$${\xi }_{\alpha \delta }({\hat{\mathrm{\Phi }}}_e)=\frac{1}{2}(\frac{\partial {({\hat{\mathrm{\Phi }}}_e)}_{\alpha }}{\partial y_{\delta }}+\frac{\partial {({\hat{\mathrm{\Phi }}}_e)}_{\delta }}{\partial y_{\alpha }}).$$

(79d)

The problems (76a)–(76d) and (77a)–(77d) are the vector problems. These have driving forces of conductivity tensor of each phase applied to the normal of the interface. The problem (78a)–(78d) is a scalar problem.

• 5. When solving the cell problems, we need the solution to be unique. We therefore require one additional condition. We make the choice that the cell averages of the auxiliary variables are zero. That is:

$$\begin{array}{c}{〈{\mathrm{\Phi }}_i〉}_i=0,{〈{\mathrm{\Phi }}_e〉}_e=0,{〈{\hat{\mathrm{\Phi }}}_i〉}_i=0,\hfill \\ \multicolumn{1}{c}{{〈{\hat{\mathrm{\Phi }}}_e〉}_e=0,{〈{\mathcal{A}}_i〉}_i=0,{〈{\mathcal{A}}_e〉}_e=0.}\end{array}$$

(80)

• 6. The auxiliary second rank tensors and vectors arising from the cell problems ( ${\mathsf{Q}}_e$ , ${\mathsf{Q}}_i$ , ${\mathsf{R}}_e$ , ${\mathsf{R}}_i$ , ${\mathbf{s}}_e$ , and ${\mathbf{s}}_i$ ) are then substituted in where they appear in the macroscale equations and this leads to a macroscale model with coefficients that encode the microstructural details.

• 7. The structure and geometry of the macroscale must be set. This includes giving boundary conditions for the homogenized cell. We also must provide initial conditions for the macroscale solid elastic displacement and the electric potential drop $V$ .

• 8. Finally, after carrying out all the above steps, our macroscale model (62a)–(62c) can then be solved.

6. Conclusion

In this work, we have derived a novel system of PDEs describing the effective electrical and mechanical behavior of an elastic composite designed to represent the myocardium tissue. Our structure comprises an elastic matrix, known as the extracellular domain, and embedded elastic subphases, known as the intracellular domain. This structure has been designed to account for the cardiac myocytes that are surrounded by a matrix which creates the heart muscle.

In order to derive the new model, we set up a microscale problem that describes the electrical and mechanical interactions between the cardiac myocytes and the surrounding matrix. We close the problem by accounting for the continuity of current densities, stresses and elastic displacements as well as the difference in the electric potentials across the interface between the myocyte and the matrix. We consider the heart muscle at a scale where the myocytes are distinct from the matrix. The microscale has an associated length much smaller than the one of the entire heart muscle (where we cannot see individual myocytes) and so the scale of the heart is known as the macroscale. Since we have this sharp scale separation, we are then able to apply the asymptotic homogenization technique to upscale the problem. The novel macroscale model that is derived is a first attempt at an asymptotic homogenization model for a coupled electrical and mechanical bidomain model. The novel model comprises balance equations for the current densities and for the stresses as well as additional terms accounting for the difference in the electric potentials at different points in the microstructure. The coefficients of the model encode the properties of the microstructure and are computed by solving the microscale differential problems arising as a result of applying the asymptotic homogenization technique.

The novel model obtained in this work is a first example of combining both the electrical and mechanical bidomain models of the heart via asymptotic homogenization. This can be considered as the next natural step in deriving computationally feasible electro-mechanical models for the heart based on the underlying microstructure. The key novelty of this work resides in considering how the mechanical deformations of the heart are influenced by the different electric fields arising from the electric potentials in each of the domains. The coupled cell problems (34a)–(34d), (35a)–(35d), as well as (36a)–(36d) encode the details of the difference in electric potentials at different points in the microstructure. The problem (52a)–(52d) is to be solved in order to obtain the tensors ${\mathbb{M}}_i$ and ${\mathbb{M}}_e$ which encode the details of the geometry and stiffness of the microstructure in the macroscale equations.

This model is designed to describe the behavior of the heart and can be used in both healthy and diseased scenarios including growth and remodeling of cardiac tissue. We could use the model to investigate the displacement of the heart tissue due to the body forces that can be applied, such as the Lorentz force that we used in this work. This could provide some important results as the imaging of the elastic displacement due to Lorentz force has been recently proposed as a potential use of MRI [48]. This would allow traditional MRI to be used to measure bio-currents and differs from functional MRI that need to use blood oxygenation level-dependent signaling to detect currents [48].

The differences in intracellular and extracellular displacements that the new model includes can be considered to determine the role they play in mechanotransduction. This investigation would suggest whether a bidomain formulation, of the type in this work, is necessary when considering growth and remodeling. The mechanical element of the model can be used to predict tissue changes in the border zone surrounding a region of ischemia following myocardial infarction. It may also be considered when trying to explain why the complete heart muscle thickens during hypertrophy. In Puwal [49], the mechanical bidomain model has been used to predict how the heart responds to elevated blood pressure and other structural abnormalities.

Having a model that can describe the electrical activity of the heart is a very useful tool to aid the understanding of how the electrical function is impaired or changed by various diseases or disorders of the heart conduction system. With the novel model, we have derived here we can investigate how the structural changes caused by myocardial ischemia induce differences in the heart electrophysiology. The ischemia causes a change in the action potentials and the membrane potential increases with a larger uptake of potassium ions [50].

The current model assumes that both the matrix and myocytes are anisotropic linear elastic materials. We could, however, allow for these materials to be hyperelastic. This would involve a method similar to that carried out in Miller and Penta [51], Collis et al. [52], and Ramírez-Torres et al. [53]. Using nonlinear elasticity for the phases, we complicate the numerical simulations involved in computing the model coefficients and the macroscale model solution. The complication is due to the fact that the two length scales in the system remain coupled and massively increase the computational load. Recently, there has been the emergence of techniques to try to overcome this coupling of the scales (see Dehghani and Zilian [54, 55]).

By allowing for the matrix and myocytes to be hyperelastic, we could approach using nonlinear elasticity theory [56] and the Holzapfel–Ogden law [8]. This would allow us to treat the myocardium as a non-homogeneous, anisotropic, nonlinear elastic, incompressible material, with a strain energy function to describe the material. The parameters of the law would be to be determined by biological measurements using techniques such as those that have been investigated in Ogden et al. [57] and Gao [58]. Using such an approach, we could even consider the influence of residual stresses as investigated by Wang et al. [9] or different filling phases of the individual heart chambers [59].

In this work, we have assumed that both the myocyte and the matrix are anisotropic elastic solids. This could, however, be modified to provide slightly more realistic behavior of the phases. In Puwal and Roth [17], the assumption that the matrix is isotropic is made. This would be an easy enough modification to be made in this work. Puwal and Roth [17] assumed that the myocyte is not isotropic and therefore used a fluid-fiber continuum to represent the mechanical properties. They make the assumption that the myocardial fibers are directed along unit vectors which are in general a function of position. We can therefore modify the stress equation to comprise the intracellular pressure and the fiber tension. The fiber tension has two parts [60]: an active isometric tension at zero strain and a passive linear stress–strain relationship characterized by Young’s modulus. This allows for curving fibers in the myocyte to be accounted for as in Chadwick [60] and Holzapfel et al. [61]. These changes to the properties of the elastic phases could be incorporated in our model to improve the applicability to the heart and to provide a framework that can be used with experimental data.

The current work could be further developed in many ways, however, potentially the most important and useful of these would be to carry out the numerical simulations. This would mean obtaining solutions to the model using a specific microstructure with the parameters chosen by real-world data for the heart, myocytes, and electrical conductivity. Comparison between the numerical results and experimental data will allow for model validation and will give an insight into the predictive capabilities of the model as a potential diagnostic tool.