A computational multiscale approach towards the modelling of microstructures with material interfaces in electrical conductors

Abstract

Motivated by the change of effective electrical properties grain or phase boundaries, a computational multiscale framework for continua with interfaces at the microscale is proposed. Cohesive-type interfaces are considered at the microscale, such that displacement and electrical potential jumps are accounted for. The governing equations for materials with interfaces under mechanical and electrical loads are provided. Based on these, a computational multiscale formulation is proposed. The coupling between the electrical and mechanical subproblem is established by the constitutive equations at the material interface. In order to investigate deformation-induced property changes at the microscale, the evolution of interface damage is elaborated. The proposed multiscale framework is further examined through various representative boundary value problems so as to identify its key properties.

1. Introduction

Every material in nature exhibits heterogeneous behaviour at a certain scale. In a material system, defects such as pores, grain boundaries, phase boundaries, secondary phases and particles can be the reasons for heterogeneity. The effective material response as observed on a macroscopic scale is a manifestation of these microscale heterogeneities. In particular, interfaces can affect the overall response of the material under consideration as illustrated in Figure 1. Experimental findings show, for instance, that grain boundaries may significantly affect the electrical properties [1, 2] such that interfaces at the microscale need to be accounted for in simulations in order to model the macroscopic behaviour realistically.

Figure 1.

(a) Backscattered electron (BSE) image of a pre-BM NbCo(Pt)Sn sample. (b) Corresponding energy-dispersive X-ray spectroscopy (EDX) map of Pt. Position layout of the four needles used for the local electrical measurements — (c) within the grain interior and (d) crossing a high-angle grain boundary. Source: Reprinted from “Dopant-segregation to grain boundaries controls electrical conductivity of n-type NbCo(Pt)Sn half-Heusler alloy mediating thermoelectric performance” by Luo et al. [2].

Understanding the properties and the behaviour of interfaces is of utmost importance because interfaces can significantly influence the effective constitutive response of the material under consideration. Different interface formulations can be introduced based on the assumed continuity condition of primary and secondary fields. Classic cohesive interface models date back to the seminal works of Barenblatt [3] and Dugdale [4]. In these classic models, traction continuity is assumed, whereas the displacement field may exhibit a jump discontinuity across the interfaces. Moreover, the cohesive interface has been proven to be an essential tool when coupled with the finite element method to model debonding and delamination processes of different materials since established in Hillerborg et al. [5]. An extensive amount of literature has been devoted to the development of this particular type of interface with various contributions focusing on purely mechanical problems [6–10] as well as on coupled multiphysics problems. In particular, thermomechanical coupling has been addressed in Willam et al. [11], Hattiangadi and Siegmund [12], Steinmann and Häsner [13], Fagerström and Larsson [14], and Özdemir et al. [15] and electro-mechanical coupling has been studied in Arias et al. [16], Utzinger et al. [17], Verhoosel et al. [18], Kozinov and Kuna [19], and Kaiser and Menzel [20]. Another type of interface is the elastic interface which dates back to the developments on surface elasticity by Gurtin and Murdoch [21]. Unlike the cohesive interface model, the elastic interface model allows traction jumps across the interface while a continuous displacement field is assumed [22, 23]. More recently, it was shown in [24] that cohesive and elastic interfaces are two extremes of generalised imperfect interfaces that account for discontinuities in both the displacement and traction field [25, 26].

Heterogeneous materials show more complex behaviour compared to their associated constituents. Homogenisation procedures can be used to determine the effective macroscopic properties of a heterogeneous material by analysing the response of its underlying microstructure [27–29]. These methods can be categorised into analytical and computational approaches. Analytical homogenisation typically requires certain simplifications and, in general, cannot resolve complex microstructures which motivates the usage of computational methods [30, 31]. Computational multiscale formulations, on the contrary, allow complex microscale processes to be accurately resolved in multiscale simulations and have successfully been applied over the past years to study multiphysics–multiscale problems, e.g., Miehe et al. [32], Kouznetsova et al. [33], Geers et al. [34], Fritzen et al. [35], Keip et al. [36], and Kochmann et al. [37]. For a comprehensive review of the subject, the reader is referred to Zohdi and Wriggers [38], Geers et al. [39], Schneider et al. [40], and Gierden et al. [41].

Classic mean field analysis techniques have been extended to cohesive interfaces in order to predict the macroscopic response in the case of debonding/dewetting of particles [42–46]. Other than analytical approaches, computational homogenisation methods have been employed in the presence of cohesive interfaces [18, 47–49]. In these methods, multiscale formulations provide the traction-separation relations that are governed by the evolution of microscopic cracks. More recently, the work [50] established a computational multiscale formulation for continua with generalised imperfect interfaces at the microscale and different aspects, such as damage evolution, are discussed [51–54].

The present contribution focuses on establishing an electro-mechanically coupled computational multiscale framework for continua with material interfaces at the microscale. Even though multiscale formulations for electro-active solids have already been elaborated in Schröder et al. [55] and Keip et al. [36], a framework for electrical conductors has been established only recently [56]. In this regard, and based on the cohesive zone formulation for electrical conductors [20], a small-strain multiscale framework that accounts for interfaces at the microscale is developed in the present contribution. The research question thus reads: Can a computational multiscale formulation be developed to capture the effect of cohesive-type (lowly-conducting) interfaces at the microscale on the macroscopic (effective) response of electrical conductors?

This paper is organised as follows: In Section 2, the thermodynamic fundamentals of continua with interfaces are summarised and the constitutive restrictions resulting from the dissipation inequality are investigated. Section 3 deals with the theoretical aspects of the electro-mechanical computational multiscale approach in the presence of cohesive interfaces. Specifically speaking, averaging theorems for the mechanical and electrical subproblems are discussed in Section 3.1, a generalised Hill–Mandel condition is proposed in Section 3.2 and suitable boundary conditions are presented. With the fundamental set of balance equations at hand, a weak form of the problem is derived and a finite element implementation is proposed in Section 4. Section 5 deals with the specific constitutive relations for the bulk and the interface, before representative boundary value problems are studied to demonstrate the capabilities of the proposed framework. The findings are summarised and concluding remarks are drawn in Section 6.

1.1. Notation

Throughout this paper, tensor notation is employed. Let $\mathit{\alpha }$ , $\mathit{\beta }$ , $\mathit{\gamma }$ , $\mathit{\delta }$ denote the first-order tensor-valued quantities. The standard dyadic product is represented as ⊗ and single and double tensor contractions are expressed as follows: $[\mathit{\alpha }\otimes \mathit{\beta }]·[\mathit{\gamma }\otimes \mathit{\delta }]=[\mathit{\beta }·\mathit{\gamma }][\mathit{\alpha }\otimes \mathit{\delta }]$ and $[\mathit{\alpha }\otimes \mathit{\beta }]:[\mathit{\gamma }\otimes \mathit{\delta }]=[\mathit{\alpha }·\mathit{\gamma }][\mathit{\beta }·\mathit{\delta }]$ . In addition, the non-standard dyadic product $[\mathit{\alpha }\otimes \mathit{\beta }]\overline{\otimes }[\mathit{\gamma }\otimes \mathit{\delta }]=[\mathit{\alpha }\otimes \mathit{\gamma }]\otimes [\mathit{\beta }\otimes \mathit{\delta }]$ is introduced. Tensor $\mathit{I}$ refers to the second-order identity tensor, $\mathit{ϵ}$ denotes the third-order permutation tensor, and $\nabla •$ , $\nabla ·•$ , denote the (right-)gradient and (right-)divergence operations, respectively.

Superscripts ${•}^{-}$ and ${•}^{+}$ indicate the quantities at opposing sides of an interface, where the jump of a quantity across the interface is given by $\left[[•]\right]={•}^{+}-{•}^{-}$ . The interfacial mean value is introduced as $\{\{•\}\}=\frac{1}{2}[{•}^{+}+{•}^{-}]$ . With these definitions at hand, the relation

$$\left[[•*]\right]=\left[[•]\right]\{\{*\}\}+\{\{•\}\}\left[[*]\right]$$

(1)

between the jump and the mean value of quantities at the interface is established.

1.2. Subdomains and interface

Control volume $\mathcal{V}$ consists of two disjoint subdomains ${\mathcal{V}}^{+}$ and ${\mathcal{V}}^{-}$ with material interface $\mathcal{I}$ in between (see Figure 2). The interface is a surface connecting the two sides of the bulk. The opposing sides of the interface can be defined as the intersection of the interface with the boundary of each subdomain, namely, ${\mathcal{I}}^{+}=\partial {\mathcal{V}}^{+}\cap \mathcal{I}$ and ${\mathcal{I}}^{-}=\partial {\mathcal{V}}^{-}\cap \mathcal{I}$ . The unit normal vector to the surface of the body $\partial \mathcal{V}$ is $\mathit{n}$ . The unit normal vectors to the plus side ${\mathcal{I}}^{+}$ and minus side ${\mathcal{I}}^{-}$ of interface $\mathcal{I}$ are ${\mathit{n}}^{+}$ and ${\mathit{n}}^{-}$ . Furthermore, the interface unit normal is defined as $\tilde{\mathit{n}}$ which points from ${\mathcal{I}}^{-}$ to ${\mathcal{I}}^{+}$ , i.e.,

$${\mathit{n}}^{+}=-\tilde{\mathit{n}},{\mathit{n}}^{-}=\tilde{\mathit{n}}.$$

(2)

Figure 2.

Control volume $\mathcal{V}={\mathcal{V}}^{-}\cup {\mathcal{V}}^{+}$ under consideration. The interface $\mathcal{I}$ is characterised by its unit surface normal vector $\tilde{\mathit{n}}$ .

2. Continuum thermodynamics

The thermodynamic fundamentals of continua featuring cohesive-type material interfaces in small strain setting are briefly recapitulated in this section. In particular, focus is laid on the mechanical subproblem in Section 2.1, on the electrical subproblem in Section 2.2, and on the balance equation of energy and the dissipation inequality in Section 2.3.

2.1. Mechanical subproblem

The governing equations of control volume $\mathcal{V}$ with material interface $\mathcal{I}$ for electro-mechanical problems are derived in Kaiser and Menzel [20]. Based on the former contribution, tangential forces on $\partial \mathcal{I}$ are neglected. The integral form of the balance of linear momentum for a quasi-static case reads

$$0={\int }_{\mathcal{V}}\rho \mathit{f}\mathrm{d}v+{\int }_{\mathcal{I}}\tilde{\rho }\tilde{\mathit{f}}\mathrm{d}a+{\int }_{\partial \mathcal{V}}\mathit{t}\mathrm{d}a,$$

(3)

where $\rho$ and $\tilde{\rho }$ are the mass density of the bulk and interface. The volume under consideration is subjected to surface tractions

$$\mathit{t}=\mathit{\sigma }·\mathit{n},$$

(4)

and loaded by volume distributed forces $\rho \mathit{f}$ in the bulk as well as surface distributed forces $\tilde{\rho }\tilde{\mathit{f}}$ at the interface. In equation (4), $\mathit{\sigma }$ and $\mathit{n}$ correspond to the stress tensor and the outward unit normal vector, respectively. Using the divergence theorem and carrying out localisation of equation (3) for continua featuring material interfaces, with the definitions

$${\tilde{\mathit{t}}}^{+}={\mathit{\sigma }}^{+}·\tilde{\mathit{n}}\mathrm{and}{\tilde{\mathit{t}}}^{-}={\mathit{\sigma }}^{-}·\tilde{\mathit{n}},$$

(5)

the local form of the balance of linear momentum is obtained as

$$\nabla ·\mathit{\sigma }+\rho \mathit{f}=0(\mathrm{in}\mathrm{the}\mathrm{bulk}),$$

(6a)

$$\left[\left[\tilde{\mathit{t}}\right]\right]+\tilde{\rho }\tilde{\mathit{f}}=0(\mathrm{at}\mathrm{the}\mathrm{interface}).$$

(6b)

Under the same assumptions, localisation of conservation of angular momentum gives rise to the classic symmetry condition of the stress tensor

$$\mathit{ϵ}:\mathit{\sigma }=0\Rightarrow \mathit{\sigma }={\mathit{\sigma }}^{\mathrm{t}}(\mathrm{in}\mathrm{the}\mathrm{bulk}).$$

(7)

Cohesive zone formulations for small strain settings yield no additional condition regarding angular momentum at the interface.

2.2. Electrical subproblem

The integral form of the continuity equation for the electrical charge for quasi-stationary processes reads

$$0={\int }_{\partial \mathcal{V}}\mathit{j}·\mathit{n}\mathrm{d}a$$

(8)

with

$$i=\mathit{j}·\mathit{n}$$

(9)

where $\mathit{j}$ corresponds to the electric current density vector and $i$ denotes the projected current density vector. Following the same procedure as in Kaiser and Menzel [20] and introducing the definitions

$${\tilde{i}}^{+}={\mathit{j}}^{+}·\tilde{\mathit{n}}\mathrm{and}{\tilde{i}}^{-}={\mathit{j}}^{-}·\tilde{\mathit{n}},$$

(10)

the local version of the continuity equation for the electric charge is obtained as

$$\nabla ·\mathit{j}=0(\mathrm{in}\mathrm{the}\mathrm{bulk}),$$

(11a)

$$\left[\left[\tilde{i}\right]\right]=0(\mathrm{at}\mathrm{the}\mathrm{interface}).$$

(11b)

The introduction of a scalar-valued electric potential field $\varphi$ can fulfil the requirements of Faraday’s law of induction under the assumption of a (quasi-)stationary state, namely

$$\mathit{e}=-\nabla \varphi (\mathrm{in}\mathrm{the}\mathrm{bulk}),$$

(12)

where $\mathit{e}$ corresponds to the electric field vector. The localisation of Faraday’s law of induction at the interface yields no additional condition for cohesive zone formulations.

2.3. Conservation of energy and dissipation inequality

Taking into account mechanical, thermal and electrical contributions, the (quasi-static) balance equation of energy for a continuum with interfaces is given by

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\mathcal{V}}\rho e\mathrm{d}v+\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\mathcal{I}}\tilde{\rho }\tilde{e}\mathrm{d}a={\mathcal{P}}^{\mathit{u}}+{\mathcal{P}}^{\theta }+{\mathcal{P}}^{\varphi }$$

(13)

where $e$ and $\tilde{e}$ correspond to the mass-specific internal energy densities of the bulk and the interface, respectively, and $t$ denotes the time. The respective contributions to equation (13) are specified as

$${\mathcal{P}}^{\mathit{u}}={\int }_{\mathcal{V}\backslash \mathcal{I}}\nabla \stackrel{·}{\mathit{u}}:\mathit{\sigma }\mathrm{d}v+{\int }_{\mathcal{I}}\left[\left[\stackrel{·}{\mathit{u}}\right]\right]·\{\{\tilde{\mathit{t}}\}\}\mathrm{d}a,$$

(14a)

$${\mathcal{P}}^{\varphi }={\int }_{\mathcal{V}\backslash \mathcal{I}}\mathit{j}·\mathit{e}\mathrm{d}v+{\int }_{\mathcal{I}}-\{\{\tilde{i}\}\}\left[\left[\varphi \right]\right]\mathrm{d}a,$$

(14b)

$${\mathcal{P}}^{\theta }={\int }_{\mathcal{V}\backslash \mathcal{I}}\rho r-\nabla ·\mathit{q}\mathrm{d}v+{\int }_{\mathcal{I}}-\left[[\mathit{q}·\tilde{\mathit{n}}]\right]\mathrm{d}a.$$

(14c)

In the set of equations (14), $\stackrel{·}{\mathit{u}}$ denotes the velocity field whereas $\mathit{q}$ and $\mathit{r}$ correspond to the heat flux vector and mass-specific heat source in the bulk. By inserting equations (14a)–(14c) into equation (13), and by localising the ensuing equation, the local form of the balance equation of energy for a continuum with interfaces takes the form

$$\rho \stackrel{·}{e}=\mathit{\sigma }:\stackrel{·}{\epsilon }+\rho r-\nabla ·\mathit{q}+\mathit{j}·\mathit{e}(\mathrm{in}\mathrm{the}\mathrm{bulk}),$$

(15a)

$$\tilde{\rho }\tilde{e}=\left[\left[\stackrel{·}{\mathit{u}}\right]\right]·\{\{\tilde{\mathit{t}}\}\}-\left[[\mathit{q}·\tilde{\mathit{n}}]\right]-\{\{\tilde{i}\}\}\varphi (\mathrm{at}\mathrm{the}\mathrm{interface}),$$

(15b)

where the definitions of the small strain deformation tensor $\mathit{\epsilon }={\nabla }^{\mathrm{sym}}\mathit{u}$ with the material time derivative $\stackrel{·}{•}$ are used. Furthermore, localisation of dissipation inequality

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\mathcal{V}}\rho s\mathrm{d}v+\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\mathcal{I}}\tilde{\rho }\tilde{s}\mathrm{d}a\ge {\int }_{\mathcal{V}}\frac{\rho r}{\theta }\mathrm{d}v-{\int }_{\partial \mathcal{V}}\frac{\mathit{q}·\mathit{n}}{\theta }\mathrm{d}a,$$

(16)

with $\theta$ , $s$ and $\tilde{s}$ denoting the absolute temperature, the mass-specific entropy density of the bulk, and the mass-specific entropy density of the interface, respectively, yields

$$\mathit{\sigma }:\stackrel{·}{\mathit{\epsilon }}-\rho [\stackrel{·}{\psi }+s\stackrel{·}{\theta }]-\frac{1}{\theta }\mathit{q}·\nabla \theta +\mathit{j}·\mathit{e}\ge 0(\mathrm{in}\mathrm{the}\mathrm{bulk})$$

(17a)

and

$$\left[\left[\stackrel{·}{\mathit{u}}\right]\right]·\{\{\tilde{\mathit{t}}\}\}-\tilde{\rho }[\tilde{\psi }+\tilde{s}\tilde{\theta }]-\{\{\tilde{i}\}\}\varphi +\tilde{\theta }[\{\{\mathit{q}·\tilde{\mathit{n}}\}\}\left[\left[{\theta }^{-1}\right]\right]+[\{\{{\theta }^{-1}\}\}-{\tilde{\theta }}^{-1}]\left[[\mathit{q}·\tilde{\mathit{n}}]\right]]\ge 0(\mathrm{at}\mathrm{the}\mathrm{interface}).$$

(17b)

In order to derive the dissipation inequality in the form (17), the Legendre(–Fenchel) transformation, relating the internal and free energy densities of the bulk and the interface as

$$\psi (\mathit{\epsilon },\theta ,•)=\underset{s}{inf}\{e(\mathit{\epsilon },s,•)-\theta s\},$$

(18a)

$$\tilde{\psi }(\left[\left[\mathit{u}\right]\right],\tilde{\theta },•)=\underset{\tilde{s}}{inf}\{\tilde{e}(\left[\left[\mathit{u}\right]\right],\tilde{s},•)-\tilde{\theta }\tilde{s}\},$$

(18b)

was invoked, where $\tilde{\theta }$ denotes the interface temperature. Evaluation of the Coleman–Noll procedure yields the constitutive relations for the mechanical subproblem

$$\mathit{\sigma }=\rho \frac{\partial \psi }{\partial \mathit{\epsilon }}(\mathrm{in}\mathrm{the}\mathrm{bulk}),$$

(19a)

$$\{\{\tilde{\mathit{t}}\}\}=\tilde{\rho }\frac{\partial \tilde{\psi }}{\partial \left[\left[\mathit{u}\right]\right]}(\mathrm{at}\mathrm{the}\mathrm{interface}).$$

(19b)

For the electrical subproblem, the constitutive restrictions for the electric current density read

$$\mathit{j}·\mathit{e}\ge 0(\mathrm{in}\mathrm{the}\mathrm{bulk}),$$

(20a)

$$-\{\{\tilde{i}\}\}\left[\left[\varphi \right]\right]\ge 0(\mathrm{at}\mathrm{the}\mathrm{interface}).$$

(20b)

3. Computational homogenisation

This section focuses on computational multiscale formulations for electro-mechanical problems in the presence of material interfaces at the microscale. The derivations in this section are based on, and a direct extension of, the multiscale framework for materials with interfaces documented in Javili et al. [50] and the work on multiscale methods for electrical conductors [56]. Whereas the focus in Kaiser and Menzel [50] was on purely mechanical problems, the focus in the present contribution is laid on the extension to electro-mechanical coupling. Likewise, material interfaces have not been considered in [56].

Constitutive relations at the microscale are assumed to be known and the objective is to compute the macroscopic response through homogenising the response of the underlying microstructure. In particular, interfaces between different constituents at the microscale are considered and assumed to significantly affect the overall behaviour of the material (see Figure 3). The governing equations at the microscale are solved subject to the assumptions of quasi-statics, quasi-stationary, and negligible body forces.

Figure 3.

The effective macroscale material response is determined by the underlying heterogeneous microstructure. In each material point on the macroscale, a representative volume element (RVE) that characterises the material microstructure is considered.

3.1. Averaging theorems

Macroscopic quantities are related to their microscopic counterparts through volume averages over the representative volume element (RVE) and additional jump contributions originating from the material interfaces. In particular, it is observed that the macroscopic strain tensor ${\mathit{\epsilon }}_{\mathrm{M}}$ , the electric field vector ${\mathit{e}}_{\mathrm{M}}$ , the stress tensor ${\mathit{\sigma }}_{\mathrm{M}}$ , and the electric current density vector ${\mathit{j}}_{\mathrm{M}}$ can be expressed as surface integrals of microscopic quantities. The definition $v=\mathrm{vol}\left(\mathcal{V}\right)$ holds. The effective macroscopic strain tensor may be specified in terms of boundary displacements, namely

$$\begin{array}{c}{\mathit{\epsilon }}_{\mathrm{M}}=\frac{1}{v}{\int }_{\mathcal{V}}\mathit{\epsilon }\mathrm{d}v+\frac{1}{v}{[{\int }_{\mathcal{I}}\left[\left[\mathit{u}\right]\right]\otimes \tilde{\mathit{n}}\mathrm{d}a]}^{\mathrm{sym}}\hfill \\ \multicolumn{1}{c}{=\frac{1}{v}{[{\int }_{\partial \mathcal{V}}\mathit{u}\otimes \mathit{n}\mathrm{d}a]}^{\mathrm{sym}}+\frac{1}{v}{[{\int }_{\mathcal{I}}[{\mathit{u}}^{+}\otimes {\mathit{n}}^{+}+{\mathit{u}}^{-}\otimes {\mathit{n}}^{-}]+\left[\left[\mathit{u}\right]\right]\otimes \tilde{\mathit{n}}\mathrm{d}a]}^{\mathrm{sym}}}\\ \multicolumn{1}{c}{=\frac{1}{v}{[{\int }_{\partial \mathcal{V}}\mathit{u}\otimes \mathit{n}\mathrm{d}a]}^{\mathrm{sym}}.}\end{array}$$

(21)

Analogous to equation (21), the effective macroscopic electric field vector

$$\begin{array}{c}{\mathit{e}}_{\mathrm{M}}=\frac{1}{v}{\int }_{\mathcal{V}}\mathit{e}\mathrm{d}v+\frac{1}{v}{\int }_{\mathcal{I}}-\left[\left[\varphi \right]\right]\tilde{\mathit{n}}\mathrm{d}a\hfill \\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\partial \mathcal{V}}-\varphi \mathit{n}\mathrm{d}a+\frac{1}{v}{\int }_{\mathcal{I}}-[{\varphi }^{+}{\mathit{n}}^{+}+{\varphi }^{+}{\mathit{n}}^{-}]-\left[\left[\varphi \right]\right]\tilde{\mathit{n}}\mathrm{d}a}\\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\partial \mathcal{V}}-\varphi \mathit{n}\mathrm{d}a}\end{array}$$

(22)

is obtained. In a small strain setting $(\tilde{\mathit{x}}={\mathit{x}}^{+}={\mathit{x}}^{-})$ and in the absence of interfacial stresses (see Javili et al. [50]), the effective macroscopic stress tensor can be expressed as

$$\begin{array}{c}{\mathit{\sigma }}_{\mathrm{M}}=\frac{1}{v}{\int }_{\mathcal{V}}\mathit{\sigma }\mathrm{d}v\hfill \\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\partial \mathcal{V}}\mathit{t}\otimes \mathit{x}\mathrm{d}a+\frac{1}{v}{\int }_{\mathcal{I}}{\mathit{t}}^{+}\otimes {\mathit{x}}^{+}+{\mathit{t}}^{-}\otimes {\mathit{x}}^{-}\mathrm{d}a}\\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\partial \mathcal{V}}\mathit{t}\otimes \mathit{x}\mathrm{d}a-\frac{1}{v}{\int }_{\mathcal{I}}\underset{=0}{\underbrace{\tilde{\mathit{t}}}}\otimes \tilde{\mathit{x}}\mathrm{d}a=\frac{1}{v}{\int }_{\partial \mathcal{V}}\mathit{t}\otimes \mathit{x}\mathrm{d}a}\end{array}$$

(23)

where use was made of equation (6b). Analogously and in view of equation (11b), the effective electric current density results in

$$\begin{array}{c}{\mathit{j}}_{\mathrm{M}}=\frac{1}{v}{\int }_{\mathcal{V}}\mathit{j}\mathrm{d}v\hfill \\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\partial \mathcal{V}}\mathit{x}i\mathrm{d}a+\frac{1}{v}{\int }_{\mathcal{I}}{\mathit{x}}^{+}{i}^{+}+{\mathit{x}}^{-}{i}^{-}\mathrm{d}a}\\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\partial \mathcal{V}}\mathit{x}i\mathrm{d}a-\frac{1}{v}{\int }_{\mathcal{I}}\tilde{\mathit{x}}\underset{=0}{\underbrace{i}}\mathrm{d}a=\frac{1}{v}{\int }_{\partial \mathcal{V}}\mathit{x}i\mathrm{d}a.}\end{array}$$

(24)

3.2. Hill–Mandel conditions

The Hill–Mandel condition, also known as the micro–macro energy equivalence condition, is a criterion used in homogenisation theory to ensure that the effective material properties of a heterogeneous solid are physically meaningful and well-defined. This criterion postulates that the virtual work on the macroscale is equivalent to the average virtual work on the microscale. In view of equation (13), the classic Hill–Mandel condition for the mechanical problem (6) is extended by interface contributions and, by using equation (1), takes the form

$$\begin{array}{c}{\mathit{\sigma }}_{\mathrm{M}}:\delta {\mathit{\epsilon }}_{\mathrm{M}}{=}^{!}\frac{1}{v}{\int }_{\mathcal{V}}\mathit{\sigma }:\delta \mathit{\epsilon }\mathrm{d}v+\frac{1}{v}{\int }_{\mathcal{I}}\{\{\tilde{\mathit{t}}\}\}·\delta \left[\left[\mathit{u}\right]\right]\mathrm{d}a\hfill \\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\mathcal{V}}\nabla ·[\delta \mathit{u}·\mathit{\sigma }]-\delta \mathit{u}·\underset{=0}{\underbrace{[\nabla ·\mathit{\sigma }]}}\mathrm{d}v+\frac{1}{v}{\int }_{\mathcal{I}}\{\{\tilde{\mathit{t}}\}\}·\delta \left[\left[\mathit{u}\right]\right]\mathrm{d}a}\\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\partial \mathcal{V}}\mathit{t}·\delta \mathit{u}\mathrm{d}a-\frac{1}{v}{\int }_{\mathcal{I}}\left[[\tilde{\mathit{t}}·\delta \mathit{u}]\right]\mathrm{d}a+\frac{1}{v}{\int }_{\mathcal{I}}\{\{\tilde{\mathit{t}}\}\}·\delta \left[\left[\mathit{u}\right]\right]\mathrm{d}a}\\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\partial \mathcal{V}}\mathit{t}·\delta \mathit{u}\mathrm{d}a.}\end{array}$$

(25)

Analogous to equation (25), the extended Hill–Mandel condition for the electrical subproblem can be derived as

$$\begin{array}{c}{\mathit{j}}_{\mathrm{M}}·\delta {\mathit{e}}_{\mathrm{M}}{=}^{!}\frac{1}{v}{\int }_{\mathcal{V}}\mathit{j}·\delta \mathit{e}\mathrm{d}v+\frac{1}{v}{\int }_{\mathcal{I}}-\{\{\tilde{i}\}\}\delta \left[\left[\varphi \right]\right]\mathrm{d}a\hfill \\ \multicolumn{1}{c}{=-\frac{1}{v}{\int }_{\mathcal{V}}\nabla ·\left[\delta \varphi \mathit{j}\right]\mathrm{d}v+\frac{1}{v}{\int }_{\mathcal{V}}-\delta \varphi \underset{=0}{\underbrace{[\nabla ·\mathit{j}]}}\mathrm{d}v+\frac{1}{v}{\int }_{\mathcal{I}}-\{\{\tilde{i}\}\}\delta \left[\left[\varphi \right]\right]\mathrm{d}a}\\ \multicolumn{1}{c}{=-\frac{1}{v}{\int }_{\partial \mathcal{V}}i\delta \varphi \mathrm{d}a+\frac{1}{v}{\int }_{\mathcal{I}}\left[\left[\tilde{i}\delta \varphi \right]\right]\mathrm{d}a+\frac{1}{v}{\int }_{\mathcal{I}}-\{\{\tilde{i}\}\}\delta \left[\left[\varphi \right]\right]\mathrm{d}a}\\ \multicolumn{1}{c}{=-\frac{1}{v}{\int }_{\partial \mathcal{V}}i\delta \varphi \mathrm{d}a.}\end{array}$$

(26)

In accordance with Javili et al. [50], the representation of the Hill–Mandel condition are expressed as a surface integral, which remains equivalent to its classic form even in the presence of an interface. This has advantages in computational implementations because existing frameworks and subroutines developed for classic homogenisation, i.e., homogenisation without interfacial contributions, can be applied.

3.3. Affine boundary conditions

Affine displacement boundary conditions are prescribed based on the macroscopic deformation state, namely

$$\mathit{u}={\mathit{\epsilon }}_{\mathrm{M}}·\mathit{x}\mathrm{on}\partial \mathcal{V}.$$

(27)

The affine displacement boundary condition fulfils the averaging theorem for the macroscopic strain tensor (21), i.e., inserting equation (27) into equation (21) results in

$$\begin{array}{c}\frac{1}{v}{[{\int }_{\partial \mathcal{V}}\mathit{u}\otimes \mathit{n}\mathrm{d}a]}^{\mathrm{sym}}=\frac{1}{v}{\int }_{\partial \mathcal{V}}[{\mathit{\epsilon }}_{\mathrm{M}}·\mathit{x}]\otimes \mathit{n}\mathrm{d}a\hfill \\ \multicolumn{1}{c}{=\frac{1}{v}[{\mathit{\epsilon }}_{\mathrm{M}}·{\int }_{\mathcal{V}}\nabla ·[\mathit{x}\otimes \mathit{I}]\mathrm{d}v]={\mathit{\epsilon }}_{\mathrm{M}}.}\end{array}$$

(28)

For the electrical subproblem, the particular form of microscopic electric potential field

$$\varphi =-{\mathit{e}}_{\mathrm{M}}·\mathit{x}\mathrm{on}\partial \mathcal{V}$$

(29)

is chosen. Analogous to the mechanical problem (28), consistency of the boundary condition (29) with the scale-bridging relation (22) is shown by inserting equation (29) in equation (22), namely

$$\begin{array}{c}\frac{1}{v}{\int }_{\partial \mathcal{V}}-\varphi \mathit{n}\mathrm{d}a=\frac{1}{v}{\int }_{\partial \mathcal{V}}-[-{\mathit{e}}_{\mathrm{M}}·\mathit{x}]\mathit{n}\mathrm{d}a\hfill \\ \multicolumn{1}{c}{=\frac{1}{v}[{\mathit{e}}_{\mathrm{M}}·{\int }_{\mathcal{V}}\nabla ·[\mathit{x}\otimes \mathit{I}]\mathrm{d}v]={\mathit{e}}_{\mathrm{M}}.}\end{array}$$

(30)

By inserting equation (27) into the Hill–Mandel condition for the mechanical problem (25)

$$\begin{array}{c}\frac{1}{v}{\int }_{\mathcal{V}}\mathit{\sigma }:\delta \mathit{\epsilon }\mathrm{d}v+\frac{1}{v}{\int }_{\mathcal{I}}\{\{\tilde{\mathit{t}}\}\}·\delta \left[\left[\mathit{u}\right]\right]\mathrm{d}a=\frac{1}{v}{\int }_{\partial \mathcal{V}}\mathit{t}·\delta \mathit{u}\mathrm{d}a\hfill \\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\partial \mathcal{V}}\mathit{t}·[\delta {\mathit{\epsilon }}_{\mathrm{M}}·\mathit{x}]\mathrm{d}a}\\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\partial \mathcal{V}}\mathit{t}\otimes \mathit{x}\mathrm{d}a:\delta {\mathit{\epsilon }}_{\mathrm{M}}={\mathit{\sigma }}_{\mathrm{M}}:\delta {\mathit{\epsilon }}_{\mathrm{M}}}\end{array}$$

(31)

is obtained. Furthermore, inserting equation (29) into equation (26) shows the consistency of the boundary condition (29) with the Hill–Mandel condition for the electrical problem, i.e.,

$$\begin{array}{c}\frac{1}{v}{\int }_{\mathcal{V}}\mathit{j}·\delta \mathit{e}\mathrm{d}v-\frac{1}{v}{\int }_{\mathcal{I}}\{\{\tilde{i}\}\}\delta \left[\left[\varphi \right]\right]\mathrm{d}a=-\frac{1}{v}{\int }_{\partial \mathcal{V}}i\delta \varphi \mathrm{d}a\hfill \\ \multicolumn{1}{c}{=-\frac{1}{v}{\int }_{\partial \mathcal{V}}i[-\delta {\mathit{e}}_{\mathrm{M}}·\mathit{x}]\mathrm{d}a}\\ \multicolumn{1}{c}{=+\frac{1}{v}{\int }_{\partial \mathcal{V}}i\mathit{x}\mathrm{d}a·\delta {\mathit{e}}_{\mathrm{M}}={\mathit{j}}_{\mathrm{M}}·\delta {\mathit{e}}_{\mathrm{M}}.}\end{array}$$

(32)

3.4. Uniform flux boundary conditions

Uniform (constant) traction conditions are prescribed at the boundary of the RVE as

$$\mathit{t}={\mathit{\sigma }}_{\mathrm{M}}·\mathit{n}\mathrm{on}\partial \mathcal{V}.$$

(33)

By inserting equation (33) into equation (23), consistency of the uniform flux boundary condition with scale-bridging relation (23) can be shown, namely

$$\begin{array}{c}\frac{1}{v}{\int }_{\mathcal{V}}\mathit{\sigma }\mathrm{d}v=\frac{1}{v}{\int }_{\partial \mathcal{V}}\mathit{t}\otimes \mathit{x}\mathrm{d}a\hfill \\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\partial \mathcal{V}}{\mathit{\sigma }}_{\mathrm{M}}·\mathit{n}\otimes \mathit{x}\mathrm{d}a-\frac{1}{v}{\int }_{\mathcal{V}}\underset{=0}{\underbrace{\nabla ·{\mathit{\sigma }}_{\mathrm{M}}}}\otimes \mathit{x}\mathrm{d}v-\frac{1}{v}{\int }_{\mathcal{I}}\underset{=0}{\underbrace{\left[\left[\tilde{\mathit{t}}\right]\right]}}\otimes \tilde{\mathit{x}}\mathrm{d}a}\\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\mathcal{V}}{\mathit{\sigma }}_{\mathrm{M}}\mathrm{d}v={\mathit{\sigma }}_{\mathrm{M}}.}\end{array}$$

(34)

Uniform electric current density boundary conditions are given by

$$i={\mathit{j}}_{\mathrm{M}}·\mathit{n}\mathrm{on}\partial \mathcal{V}.$$

(35)

By inserting equation (35) into equation (24) and by analogy with the mechanical subproblem, one arrives at

$$\begin{array}{c}\frac{1}{v}{\int }_{\mathcal{V}}\mathit{j}\mathrm{d}v=\frac{1}{v}{\int }_{\partial \mathcal{V}}i\mathit{x}\mathrm{d}a\hfill \\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\partial \mathcal{V}}\mathit{x}\otimes {\mathit{j}}_{\mathrm{M}}·\mathit{n}\mathrm{d}a-\frac{1}{v}{\int }_{\mathcal{V}}\underset{=0}{\underbrace{\nabla ·{\mathit{j}}_{\mathrm{M}}}}\mathit{x}\mathrm{d}v-\frac{1}{v}{\int }_{\mathcal{I}}\tilde{\mathit{x}}\underset{=0}{\underbrace{\left[\left[i\right]\right]}}\mathrm{d}a}\\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\mathcal{V}}{\mathit{j}}_{\mathrm{M}}\mathrm{d}v={\mathit{j}}_{\mathrm{M}}}\end{array}$$

(36)

which shows the consistency of the uniform electric current density boundary condition with scale-bridging relation (24). In view of equations (33) and (35), the evaluation of the Hill–Mandel condition for the mechanical subproblem (25) results in

$$\begin{array}{c}\frac{1}{v}{\int }_{\mathcal{V}}\mathit{\sigma }:\delta \mathit{\epsilon }\mathrm{d}v+\frac{1}{v}{\int }_{\mathcal{I}}\{\{\tilde{\mathit{t}}\}\}·\delta \left[\left[\mathit{u}\right]\right]\mathrm{d}a=\frac{1}{v}{\int }_{\partial \mathcal{V}}\mathit{t}·\delta \mathit{u}\mathrm{d}a\hfill \\ \multicolumn{1}{c}{=\frac{1}{v}{\int }_{\partial \mathcal{V}}\delta \mathit{u}·[{\mathit{\sigma }}_{\mathrm{M}}·\mathit{n}]\mathrm{d}a}\\ \multicolumn{1}{c}{={\mathit{\sigma }}_{\mathrm{M}}:[\frac{1}{v}{\int }_{\partial \mathcal{V}}\delta \mathit{u}\otimes \mathit{n}\mathrm{d}a]={\mathit{\sigma }}_{\mathrm{M}}:\delta {\mathit{\epsilon }}_{\mathrm{M}}}\end{array}$$

(37)

and for the electrical subproblem (26)

$$\begin{array}{c}\frac{1}{v}{\int }_{\mathcal{V}}\mathit{j}·\delta \mathit{e}\mathrm{d}v-\frac{1}{v}{\int }_{\mathcal{I}}\{\{\tilde{i}\}\}\delta \left[\left[\varphi \right]\right]\mathrm{d}a=-\frac{1}{v}{\int }_{\partial \mathcal{V}}i\delta \varphi \mathrm{d}a\hfill \\ \multicolumn{1}{c}{=-\frac{1}{v}{\int }_{\partial \mathcal{V}}\delta \varphi [{\mathit{j}}_{\mathrm{M}}·\mathit{n}]\mathrm{d}a}\\ \multicolumn{1}{c}{={\mathit{j}}_{\mathrm{M}}:[\frac{1}{v}{\int }_{\partial \mathcal{V}}-\delta \varphi \mathit{n}\mathrm{d}a]={\mathit{j}}_{\mathrm{M}}:\delta {\mathit{e}}_{\mathrm{M}}}\end{array}$$

(38)

is obtained.

4. Finite element implementation

This section addresses the finite element implementation of the proposed multiscale formulation for electrical conductors in the presence of material interfaces. In particular, the focus is on the weak form of the microscale boundary value problem in Section 4.1 before discrete closed-form solutions for the homogenised macroscale fields, and the associated tangent stiffness contributions are derived in Sections 4.2 and 4.3, respectively.

4.1. Weak form of the coupled problem

To derive a weak form of the mechanical problem, equation (6a) is multiplied with test function ${\mathit{\eta }}^{\mathit{u}}$ , integrated over domain $\mathcal{V}$ (which consists of subdomains ${\mathcal{V}}^{-}$ and ${\mathcal{V}}^{+}$ )

$$0={\int }_{{\mathcal{V}}^{+}}{\mathit{\eta }}^{\mathit{u}}·[\nabla ·\mathit{\sigma }]\mathrm{d}v+{\int }_{{\mathcal{V}}^{-}}{\mathit{\eta }}^{\mathit{u}}·[\nabla ·\mathit{\sigma }]\mathrm{d}v,$$

(39)

and use is made of the divergence theorem and equation (1) to arrive at

$$\begin{array}{c}0=-{\int }_{{\mathcal{V}}^{+}\cup {\mathcal{V}}^{-}}\nabla {\mathit{\eta }}^{\mathit{u}}:\mathit{\sigma }\mathrm{d}v+{\int }_{\partial \mathcal{V}}{\mathit{\eta }}^{\mathit{u}}·\mathit{t}\mathrm{d}a-{\int }_{\mathcal{I}}\left[\left[{\mathit{\eta }}^{\mathit{u}}\right]\right]·\{\{\tilde{\mathit{t}}\}\}+\{\{{\mathit{\eta }}^{\mathit{u}}\}\}·\left[\left[\tilde{\mathit{t}}\right]\right]\mathrm{d}a.\hfill \end{array}$$

(40)

By inserting equation (6b) into equation (40), the governing equation for the mechanical subproblem in weak form reads

$$0=-{\int }_{{\mathcal{V}}^{+}\cup {\mathcal{V}}^{-}}\nabla {\mathit{\eta }}^{\mathit{u}}:\mathit{\sigma }\mathrm{d}v+{\int }_{\partial \mathcal{V}}{\mathit{\eta }}^{\mathit{u}}·\mathit{t}\mathrm{d}a-{\int }_{\mathcal{I}}\left[\left[{\mathit{\eta }}^{\mathit{u}}\right]\right]·\{\{\tilde{\mathit{t}}\}\}\mathrm{d}a.$$

(41)

For the electrical problem, the continuity equation of electric charge (11a) is multiplied with test function ${\eta }^{\varphi }$

$$0={\int }_{{\mathcal{V}}^{+}}{\eta }^{\varphi }\nabla ·\mathit{j}\mathrm{d}v+{\int }_{{\mathcal{V}}^{-}}{\eta }^{\varphi }\nabla ·\mathit{j}\mathrm{d}v,$$

(42)

and the application of the divergence theorem together with equation (1)

$$0=-{\int }_{{\mathcal{V}}^{+}\cup {\mathcal{V}}^{-}}\nabla {\eta }^{\varphi }·\mathit{j}\mathrm{d}v+{\int }_{\partial \mathcal{V}}{\eta }^{\varphi }\mathit{j}·\mathit{n}\mathrm{d}a-{\int }_{\mathcal{I}}\left[\left[{\eta }^{\varphi }\right]\right]\{\{\tilde{i}\}\}+\{\{{\eta }^{\varphi }\}\}\tilde{i}\mathrm{d}a$$

(43)

specifies, in virtue of equation (11b), the weak form

$$0=-{\int }_{{\mathcal{V}}^{+}\cup {\mathcal{V}}^{-}}\nabla {\eta }^{\varphi }·\mathit{j}\mathrm{d}v+{\int }_{\partial \mathcal{V}}{\eta }^{\varphi }\mathit{j}·\mathit{n}\mathrm{d}a-{\int }_{\mathcal{I}}\left[\left[{\eta }^{\varphi }\right]\right]\{\{\tilde{i}\}\}\mathrm{d}a.$$

(44)

4.2. Homogenisation

The finite element implementation of continua with interfaces is adopted from [20]. The primary fields and test functions are discretised as follows

$${\mathit{u}}^{\mathrm{h}}=\sum _{C=1}^{{\mathrm{n}}_{\mathrm{en}}}{\mathit{u}}_C{N}_C,{\mathit{\eta }}^{\mathit{u}\mathrm{h}}=\sum _{A=1}^{{\mathrm{n}}_{\mathrm{en}}}{\mathit{\eta }}^{\mathit{u}}{N}_A,$$

(45a)

$${\varphi }^{\mathrm{h}}=\sum _{D=1}^{{\mathrm{n}}_{\mathrm{en}}}{\varphi }_D{N}_D,{\eta }^{\varphi \mathrm{h}}=\sum _{B=1}^{{\mathrm{n}}_{\mathrm{en}}}{\eta }_B^{\varphi }{N}_B,$$

(45b)

where $n_{\mathrm{en}}$ and $N_{•}$ are the number of element nodes and the shape functions, respectively. The linearised system of equations for the coupled problem with affine boundary conditions takes the form

$${\left[\begin{array}{cc}{\mathbf{K}}^{\mathit{u}\mathit{u}}+{\tilde{\mathbf{K}}}^{\mathit{u}\mathit{u}}\hfill & {\mathbf{K}}^{\mathit{u}\varphi }+{\tilde{\mathbf{K}}}^{\mathit{u}\varphi }\hfill \\ \multicolumn{1}{c}{{\mathbf{K}}^{\varphi \mathit{u}}+{\tilde{\mathbf{K}}}^{\varphi \mathit{u}}}& {\mathbf{K}}^{\varphi \varphi }+{\tilde{\mathbf{K}}}^{\varphi \varphi }\hfill \end{array}\right]}_q·{\left[\begin{array}{c}\mathrm{\Delta }\hat{\mathit{u}}\hfill \\ \multicolumn{1}{c}{\mathrm{\Delta }\hat{\mathit{\varphi }}}\end{array}\right]}_q={\left[\begin{array}{c}\mathrm{\Delta }{\mathit{f}}^{\mathit{u}}+\mathrm{\Delta }{\tilde{f}}^{\mathit{u}}\hfill \\ \multicolumn{1}{c}{\mathrm{\Delta }{\mathit{f}}^{\varphi }+\mathrm{\Delta }{\tilde{f}}^{\varphi }}\end{array}\right]}_q,$$

(46)

for the Newton iteration step $q$ , where $\hat{\mathit{u}}$ and $\hat{\mathit{\varphi }}$ denote the global lists of nodal degrees of freedom. In equation (46), ${\mathbf{K}}^{•*}$ and ${\tilde{\mathbf{K}}}^{•*}$ are the generalised stiffness matrices of bulk and interface. Similarly, ${\mathit{f}}^{•}$ and ${\tilde{f}}^{•}$ are the contributions of bulk and interface to the generalised reaction force vector. By collecting the contributions of the individual fields in the matrix $\mathbf{K}$ and the vectors $\mathcal{X}$ and $\mathit{f}$ , system (46) can be written as

$$\mathbf{K}·\mathrm{\Delta }\mathcal{X}=\mathrm{\Delta }\mathit{f}$$

(47)

In the case of periodic boundary conditions, the generalised stiffness matrices and generalised reaction force vectors are substituted by reduced versions that result from enforcing linear constraints between opposing RVE boundaries (see, e.g., Kaiser and Menzel [56]). The proposed finite element implementation uses scale-bridging relations that rely on equations (23) and (24). Accordingly, the discrete versions of the macroscopic stress tensor

$${\mathit{\sigma }}_{\mathrm{M}}=\frac{1}{v}{\int }_{\partial \mathcal{V}}\mathit{t}\otimes \mathit{x}\mathrm{d}a\approx \frac{1}{v}{\sum _{i=1}^{{\mathrm{n}}_{\mathrm{pn}}}}^{(i)}{\mathit{f}}^{\mathit{u}}{\otimes }^{(i)}\mathit{x}$$

(48)

and the electric current density vector

$${\mathit{j}}_{\mathrm{M}}=\frac{1}{v}{\int }_{\partial \mathcal{V}}i\mathit{x}\mathrm{d}a\approx \frac{1}{v}{\sum _{i=1}^{{\mathrm{n}}_{\mathrm{pn}}}}^{(i)}{f^{\varphi }}^{(i)}\mathit{x}$$

(49)

are obtained. In equation (48), ${\mathrm{n}}_{\mathrm{pn}}$ denotes the number of Dirichlet nodes, and ${\mathit{x}}^{(i)}$ and ${}^{(i)}{\mathit{f}}^{\mathit{u}}$ are the position and the reaction force vector of node $i$ , corresponding to the mechanical problem. Similarly, ${}^{(i)}{f}^{\varphi }$ represents the generalised reaction force of node $i$ for the electrical problem.

4.3. Consistent tangent stiffness tensors

In order to derive consistent algorithmic tangent stiffness operators, the linear system (47) is partitioned into Dirichlet ${\mathcal{X}}_{\mathrm{d}}$ and free ${\mathcal{X}}_{\mathrm{f}}$ degrees of freedom

$$\left[\begin{array}{cc}{\mathbf{K}}_{\mathrm{dd}}\hfill & {\mathbf{K}}_{\mathrm{df}}\hfill \\ \multicolumn{1}{c}{{\mathbf{K}}_{\mathrm{fd}}}& {\mathbf{K}}_{\mathrm{ff}}\hfill \end{array}\right]·\left[\begin{array}{c}\mathrm{\Delta }{\mathcal{X}}_{\mathrm{d}}\hfill \\ \multicolumn{1}{c}{\mathrm{\Delta }{\mathcal{X}}_{\mathrm{f}}}\end{array}\right]=\left[\begin{array}{c}\mathrm{\Delta }{\mathit{f}}_{\mathrm{d}}\hfill \\ \multicolumn{1}{c}{0}\end{array}\right].$$

(50)

Solving the second line of the system of equations for ${\mathcal{X}}_{\mathrm{f}}$ and inserting the ensuing relation in the first line, the closed-form expression

$$[{\mathbf{K}}_{\mathrm{dd}}-{\mathbf{K}}_{\mathrm{df}}·{\mathbf{K}}_{\mathrm{ff}}^{-1}·{\mathbf{K}}_{\mathrm{fd}}]·\mathrm{\Delta }{\mathcal{X}}_{\mathrm{d}}=\mathrm{\Delta }{\mathit{f}}_{\mathrm{d}}$$

(51)

is obtained. For affine and periodic boundary conditions, the kinematic relations are given as

$${\mathrm{\Delta }}^{(j)}{\mathit{u}}_{\mathrm{d}}=\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{M}}·{}^{(j)}\mathit{x}$$

(52a)

$${\mathrm{\Delta }}^{(j)}{\varphi }_{\mathrm{d}}=-\mathrm{\Delta }{\mathit{e}}_{\mathrm{M}}·{}^{(j)}\mathit{x}.$$

(52b)

By partitioning equation (50) into mechanical and electrical contributions and after inserting equation (52), the relations

$${\mathrm{\Delta }}^{(i)}{\mathit{f}}^{\mathit{u}}={\sum _{j=1}^{{\mathrm{n}}_{\mathrm{pn}}}}^{(\mathrm{ij})}{\hat{\mathbf{K}}}^{\mathit{u}\mathit{u}}·\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{M}}·{}^{(j)}\mathit{x}-{\sum _{j=1}^{{\mathrm{n}}_{\mathrm{pn}}}}^{(\mathrm{ij})}{\hat{\mathbf{K}}}^{\mathit{u}\varphi }\otimes \mathrm{\Delta }{\mathit{e}}_{\mathrm{M}}·{}^{(j)}\mathit{x},$$

(53)

and

$${\mathrm{\Delta }}^{(i)}{f}^{\varphi }={\sum _{j=1}^{{\mathrm{n}}_{\mathrm{pn}}}}^{(\mathrm{ij})}{\hat{\mathbf{K}}}^{\varphi \mathit{u}}·\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{M}}·{}^{(j)}\mathit{x}-{\sum _{j=1}^{{\mathrm{n}}_{\mathrm{pn}}}}^{(\mathrm{ij})}{\hat{\mathbf{K}}}^{\varphi \varphi }\mathrm{\Delta }{\mathit{e}}_{\mathrm{M}}·{}^{(j)}\mathit{x},$$

(54)

between increments in the generalised reaction force vectors at the microscale and changes in ${\mathit{\epsilon }}_{\mathrm{M}}$ and ${\mathit{e}}_{\mathrm{M}}$ are obtained.

Finally, by inserting equation (53) into the discrete representation of the stress tensor (48), one arrives at

$$\begin{array}{c}\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{M}}\approx \stackrel{=\frac{\mathrm{d}{\mathit{\sigma }}_{\mathrm{M}}}{\mathrm{d}{\mathit{\epsilon }}_{\mathrm{M}}}}{\overbrace{\frac{1}{v}\sum _{i=1}^{{\mathrm{n}}_{\mathrm{pn}}}{\sum _{j=1}^{{\mathrm{n}}_{\mathrm{pn}}}}^{(\mathrm{ij})}{\hat{\mathbf{K}}}^{\mathit{u}\mathit{u}}\overline{\otimes }[{}^{(i)}\mathit{x}\otimes {}^{(j)}\mathit{x}]}}:\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{M}}\hfill \\ \multicolumn{1}{c}{\underset{=\frac{\mathrm{d}{\mathit{\sigma }}_{\mathrm{M}}}{\mathrm{d}{\mathit{e}}_{\mathrm{M}}}}{\underbrace{-\frac{1}{v}\sum _{i=1}^{{\mathrm{n}}_{\mathrm{pn}}}\sum _{j=1}^{{\mathrm{n}}_{\mathrm{pn}}}{}^{(\mathrm{ij})}{\hat{\mathbf{K}}}^{\mathit{u}\varphi }\otimes [{}^{(i)}\mathit{x}\otimes {}^{(j)}\mathit{x}]}}·\mathrm{\Delta }{\mathit{e}}_{\mathrm{M}}.}\end{array}$$

(55)

Analogously, inserting equation (54) into the discrete representation of the electric current density vector (49) yields

$$\begin{array}{c}\mathrm{\Delta }{\mathit{j}}_{\mathrm{M}}\approx \stackrel{=\frac{\mathrm{d}{\mathit{j}}_{\mathrm{M}}}{\mathrm{d}{\mathit{\epsilon }}_{\mathrm{M}}}}{\overbrace{\frac{1}{v}\sum _{i=1}^{{\mathrm{n}}_{\mathrm{pn}}}{\sum _{j=1}^{{\mathrm{n}}_{\mathrm{pn}}}}^{(i)}\mathit{x}\otimes {}^{(\mathrm{ij})}{\hat{\mathbf{K}}}^{\varphi \mathit{u}}\otimes {}^{(j)}\mathit{x}}}:\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{M}}\hfill \\ \multicolumn{1}{c}{\underset{=\frac{\mathrm{d}{\mathit{j}}_{\mathrm{M}}}{\mathrm{d}{\mathit{e}}_{\mathrm{M}}}}{\underbrace{-\frac{1}{v}\sum _{i=1}^{{\mathrm{n}}_{\mathrm{pn}}}{\sum _{j=1}^{{\mathrm{n}}_{\mathrm{pn}}}}^{(\mathrm{ij})}{\hat{\mathbf{K}}}^{\varphi \varphi }{}^{(i)}\mathit{x}\otimes {}^{(j)}\mathit{x}}}·\mathrm{\Delta }{\mathit{e}}_{\mathrm{M}}.}\end{array}$$

(56)

5. Numerical examples

This section deals with a study of representative boundary value problems. Constitutive relations that are employed in the current work for the bulk and interface are discussed in Section 5.1. In particular, the electrical and mechanical material models are adopted from Kaiser and Menzel [20] for demonstration purposes. An analytical solution for the effective conductivity in a (quasi-)one-dimensional setting is provided in Section 5.2. Various microstructures are studied, and the effective macroscopic conductivity tensors are obtained and compared in Section 5.3. Moreover, the internal length scale is discussed in Section 5.4, and the coupled electro-mechanical response due to mechanically induced microscale damage processes is exemplarily studied in Section 5.5. The specific material parameters are summarised in Table 1, a plane strain setting is assumed, and quadrilateral bulk and linear interface elements are used.

Table 1.

Material parameters used in the analytical solution and finite element–based simulations.

$E$ , $\tilde{E}$	210,000 N/mm 2	210,000 N/mm 3
$\nu$	0.3	–
$\kappa$ , $\tilde{\kappa }$	1450 A/[Vmm]	1450 A/[Vmm 2]
$Q_0$	–	500 N/mm 2
$G_{\mathrm{F}}$	–	200 N/mm 2

5.1. Constitutive relations

The mechanical response of the bulk material is assumed to be governed by the volume specific Helmholtz free energy density function

$$\rho \psi =\mathit{\psi }\left(\mathit{\epsilon }\right)=\frac{K}{2}\mathrm{t}{\mathrm{r}}^2\left(\mathit{\epsilon }\right)+G{\mathit{\epsilon }}_{\mathrm{dev}}:{\mathit{\epsilon }}_{\mathrm{dev}}.$$

(57)

In equation (57), ${\mathit{\epsilon }}_{\mathrm{dev}}$ is the deviatoric part of the strain tensor, $K=\frac{E}{3[1-2\nu ]}$ is the bulk modulus, and $G=\frac{E}{2[1+\nu ]}$ is the shear modulus. Evaluation of the constitutive equation (19) for the volume-specific Helmholtz free energy density (57) results in

$$\mathit{\sigma }=\frac{\partial \mathit{\psi }}{\partial \mathit{\epsilon }}=K\mathrm{tr}\left(\mathit{\epsilon }\right)\mathit{I}+2G{\mathit{\epsilon }}_{\mathrm{dev}}.$$

(58)

The interface incorporates a model that accounts for brittle damage, along with a linear-elastic material response. More specifically speaking, the area-specific interface Helmholtz free energy density is taken as a function of displacement jump $\left[\left[\mathit{u}\right]\right]$ and interface damage variable $\tilde{d}$ . The mechanical response at the interface is assumed to show different behaviour patterns under tensile and compressive loadings. The interface elastic stiffness $\tilde{E}$ is assumed to reduce due to damage evolution in the tensile region $(\left[\left[\mathit{u}\right]\right]·\tilde{\mathit{n}}>0)$

$$\begin{array}{c}\tilde{\rho }\tilde{\psi }=\tilde{\mathit{\psi }}(\left[\left[\mathit{u}\right]\right],\tilde{d})=\frac{1}{2}[1-\tilde{d}]\tilde{E}\left[\left[\mathit{u}\right]\right]·\mathit{I}·\left[\left[\mathit{u}\right]\right],\hfill \end{array}$$

(59a)

whereas the interface is assumed to regain its normal stiffness under compressive loadings $(\left[\left[\mathit{u}\right]\right]·\tilde{\mathit{n}}\le 0)$

$$\begin{array}{c}\tilde{\rho }\tilde{\psi }=\tilde{\mathit{\psi }}(\left[\left[\mathit{u}\right]\right],\tilde{d})=\frac{1}{2}\tilde{E}\left[\left[\mathit{u}\right]\right]·[1-\tilde{d}][\mathit{I}-\tilde{\mathit{n}}\otimes \tilde{\mathit{n}}]·\left[\left[\mathit{u}\right]\right]\hfill \\ \multicolumn{1}{c}{+\frac{1}{2}\tilde{E}\left[\left[\mathit{u}\right]\right]·\tilde{\mathit{n}}\otimes \tilde{\mathit{n}}·\left[\left[\mathit{u}\right]\right].}\end{array}$$

(59b)

The corresponding traction–separation law at the interface

$$\{\{\tilde{\mathit{t}}\}\}=\frac{\partial \tilde{\mathit{\psi }}}{\partial \left[\left[\mathit{u}\right]\right]}=\{\begin{array}{cc}[1-\tilde{d}]\tilde{E}\left[\left[\mathit{u}\right]\right]\hfill & \mathrm{if}\left[\left[\mathit{u}\right]\right]·\tilde{\mathit{n}}>0\hfill \\ \multicolumn{1}{c}{\tilde{E}[[1-\tilde{d}]\mathit{I}+\tilde{d}\tilde{\mathit{n}}\otimes \tilde{\mathit{n}}]·\left[\left[\mathit{u}\right]\right]}& \mathrm{if}\left[\left[\mathit{u}\right]\right]·\tilde{\mathit{n}}\le 0\hfill \end{array}$$

(60)

is derived by evaluating equation (19) for equation (59). For the damage evolution, an exponential form of the interface damage variable

$$\tilde{d}=1-\frac{{\chi }_0}{\chi }\exp (-[\chi -{\chi }_0]\frac{\tilde{E}{\chi }_0}{G_{\mathrm{F}}-\frac{1}{2}\tilde{E}{\chi }_0^2})$$

(61)

is assumed as proposed in Radulovic et al. [57] and Heitbreder et al. [58], where $G_{\mathrm{F}}$ , $\chi$ , and ${\chi }_0$ are the fracture toughness, a displacement jump-type internal variable, and the critical interface opening, respectively. Moreover, variable $\chi$ at time $t^{*}$ is defined as

$$\chi =max\{\underset{0<t\le t^{*}}{max}\{\mathrm{sgn}(\left[\left[\mathit{u}\right]\right]·\tilde{\mathit{n}})\parallel \left[\left[\mathit{u}\right]\right]\parallel \},{\chi }_0\}.$$

(62)

For the electrical subproblem, a linear relation between the electric field vector and the electric current density vector in the bulk is adopted in accordance with the restrictions posed by the dissipation inequality (20), i.e.,

$$\mathit{j}=\mathcal{S}·\mathit{e}.$$

(63)

For the sake of simplicity, an isotropic electrical conductivity tensor is assumed

$$\mathcal{S}=\kappa \mathit{I},$$

(64)

where $\kappa$ denotes a scalar-valued electrical conductivity parameter. Similar to the bulk, a linear relation between the jump of the electric potential and the mean value of the electric current is considered at the interface. The electrical conductivity is assumed to be influenced by damage evolution at the interface. To this end, a $[1-\tilde{d}]$ formulation is employed as

$$\{\{\tilde{i}\}\}=\{\begin{array}{cc}-[1-\tilde{d}]\tilde{\kappa }\varphi \hfill & \mathrm{if}\left[\left[\mathit{u}\right]\right]·\tilde{\mathit{n}}>0\hfill \\ \multicolumn{1}{c}{-\tilde{\kappa }\varphi }& \mathrm{if}\left[\left[\mathit{u}\right]\right]·\tilde{\mathit{n}}\le 0\hfill \end{array}$$

(65)

where $\tilde{\kappa }$ denotes the idealised conductivity of the interface.

5.2. Analytical solution

This section deals with the derivation of a (quasi-)one-dimensional analytical solution for the calculation of the effective conductivity of an RVE with a material interface, see Figure 4. For the electrical problem under periodic boundary conditions with a prescribed electric potential difference $\overline{\mathrm{\Delta }\varphi }$ , electric current density $j=\parallel \mathit{j}\parallel$ can be calculated as

$$j={[\frac{l/2}{\kappa }+\frac{1}{\tilde{\kappa }}+\frac{l/2}{\kappa }]}^{-1}\overline{\mathrm{\Delta }\varphi }$$

(66)

using elementary calculation rules for serially connected resistors in electrical circuits. Accordingly, the effective conductivity of the (quasi-)one-dimensional RVE in the direction of ${\mathit{e}}_1$ , i.e., ${\kappa }_{\mathrm{M}}={\mathit{e}}_1·{\mathcal{S}}_{\mathrm{M}}·{\mathit{e}}_1$ , reads

$${\kappa }_{\mathrm{M}}={[\frac{l}{\kappa }+\frac{1}{\tilde{\kappa }}]}^{-1}l.$$

(67)

Figure 4.

RVE of length $l$ and cross-sectional area $A=l\times 1.0\mathrm{mm}$ , subjected to periodic electrical boundary conditions. A material interface crossing the boundary is placed at the middle of the unit cell. The boundary conditions are applied in macroscopic problem such that $\mathrm{\Delta }\varphi$ is parallel to ${\mathit{e}}_1$ for RVE. Dimensions are given in $\mathrm{mm}$ .

Comparing the analytical (67) and finite element simulation (56) results verifies the proposed finite element framework. Moreover, the simulation results shown in Figure 4 clearly indicate a size-dependent macroscale material response due to the presence of material interfaces at the microscale. This size effect will be studied in detail in Section 5.4.

5.3. Effective macroscopic conductivity tensor

The objective of this section is to demonstrate the use of the proposed electro-mechanical multiscale formulation within a two-dimensional setting. Therefore, effective macroscopic conductivity tensors are calculated for the different periodic microstructures depicted in Figure 5. The results are normalised with respect to the idealised material parameters given in Table 1, i.e.

$${\left[{\mathcal{S}}_{\mathrm{M}}\right]}_{{\mathit{e}}_{1,2}}^{\mathrm{mat}}=\left[\begin{array}{cc}1.0000\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1.0000\hfill \end{array}\right]\kappa$$

(68a)

$${\left[{\mathcal{S}}_{\mathrm{M}}\right]}_{{\mathit{e}}_{1,2}}^{\mathrm{cir}}=\left[\begin{array}{cc}0.7000\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0.7000\hfill \end{array}\right]\kappa$$

(68b)

$${\left[{\mathcal{S}}_{\mathrm{M}}\right]}_{e_{1,2}}^{\mathrm{gro}}=\left[\begin{array}{cc}0.9096\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0.7322\hfill \end{array}\right]\kappa$$

(68c)

Figure 5.

Sketch of different two-dimensional RVEs that are analysed with the electro-mechanical multiscale finite element formulation. The microstructures were chosen in accordance with Kaiser and Menzel [56]. Dimensions are given in $\mathrm{mm}$ .

The effective macroscopic conductivity tensors corresponding to microstructures with circular (68b) and groove-shaped interfaces (68c) are compared with the conductivity tensor of an idealised material (68a). The presence of the circular cohesive interface causes a reduction in conductivity of approximately 30%. In accordance with theoretical predictions, the constitutive response remains isotropic. In comparison, the reduction in the 22-coefficient of the conductivity tensor caused by the presence of a groove-shaped interface is significantly higher than the reduction in the 11-coefficent such that the effective macroscopic material response is anisotropic.

5.4. Internal length scale

Classic first-order multiscale formulations do not account for size-dependent material behaviour. In other words, the macroscopic response does not change with respect to the RVE size. To solve this problem, second-order homogenisation methods that incorporate higher-order gradients are developed (see, for instance, Kouznetsova et al. [33]). On the contrary, it was observed in Section 5.2 that by accounting for material interfaces at the microscale, the continuum under consideration is endowed with an internal length scale such that size effects can naturally be accounted for. In accordance with Kaiser and Menzel [20], the internal length scale of the electrical subproblem is introduced as

$${\ell }_{\varphi }=\frac{\kappa }{\tilde{\kappa }}.$$

(69)

The ratio of the internal length scale to the RVE size $l$ is the length scale ratio ${\ell }_{\varphi }/l$ . By changing the length scale ratio (e.g., fixing the RVE size, $l=1.0\mathrm{mm}$ , while varying the internal length scale), the size effect induced by the presence of material interfaces can be studied. The matrix–inclusion problem with a cohesive interface shown in Figure 6 is exemplarily analysed, and the material response for various length scale ratios is depicted in Figure 7. The upper limit corresponds to the perfectly conducting RVE, with the effect of the interface being negligible. The lower limit, on the contrary, corresponds to the case of a pore, i.e., the material interface electrically insulates the inclusion.

Figure 6.

Sketch of two-dimensional multiphase microstructure analysed with the electro-mechanical multiscale finite element formulation. The conductivity of the matrix (indicated by grey colour) and the interface between inclusion and matrix are given in Table 1. The conductivity of the inclusion (indicated by green colour) is chosen as ${\kappa }_{\mathrm{inc}}=0.1\kappa$ . Dimensions are given in $\mathrm{mm}$ .

Figure 7.

Effective macroscopic conductivity of the multiphase RVE depicted in Figure 6. A size-dependent behaviour is observed with limit cases corresponding to a perfect conduction interface and a perfect insulation.

5.5. Electro-mechanical coupling

This section focuses on the change in the effective macroscopic conductivity tensor due to mechanically induced damage processes at the microscale. An inclusion with low conduction embedded in high conductor matrix is selected for this example. The material parameters $E_{\mathrm{inc}}=210,000\mathrm{N}/\mathrm{m}{\mathrm{m}}^2$ , ${\kappa }_{\mathrm{inc}}=0.3$ , and ${\kappa }_{\mathrm{inc}}=145\mathrm{A}/$ (Vmm) govern the constitutive response of the inclusion and the material parameters of the matrix and the interface parameters being chosen according to Table 1. In view of equation (65), the electrical and mechanical subproblems are coupled at the interface by means of damage variable $\tilde{d}$ . The RVE is loaded mechanically by

$${\mathit{\epsilon }}_{\mathrm{M}}=\alpha (t){\mathit{e}}_1\otimes {\mathit{e}}_1$$

(70a)

and electrically via

$${\mathit{e}}_{\mathrm{M}}=\beta \left(t\right){\mathit{e}}_1$$

(70b)

with parameters $\alpha \left(t\right)$ and $\beta \left(t\right)$ controlling the monotonic load history as a function of time $t$ .

The decrease in the $[{\mathcal{S}}_{\mathrm{M}}{]}_{11}$ and $[{\mathcal{S}}_{\mathrm{M}}{]}_{22}$ coefficients of the effective macroscopic conductivity tensor as a function of deformation in terms of the prescribed macroscale strain ${\mathit{\epsilon }}_{\mathrm{M}}$ is shown in Figure 8. The obtained values are normalised with respect to the conductivity of the matrix. A constant conductivity in ${\mathit{e}}_1$ - and ${\mathit{e}}_2$ -direction is observed up to the onset of damage evolution, i.e., for $\alpha \lesssim 0.01$ . The initial conductivity is reduced due to the inclusion with comparatively low conduction as well as the resistance of the cohesive interface. After damage evolution starts, a significantly different decrease of the electrical conductivity in the two spatial directions is observed. The norm of electric current density field $\parallel \mathit{j}\parallel =\sqrt{\mathit{j}·\mathit{j}}$ for different load cases is additionally illustrated in Figure 9 on the deformed configuration. With the displacement jump at the interface being closely related to the damage evolution at the interface and, accordingly, to the degradation of electrical properties, Figure 9 clearly reveals the cause of the deformation-induced anisotropy.

Figure 8.

Effective conductivity as a function of prescribed mechanical boundary conditions.

Figure 9.

Electric current density field $\parallel \mathit{j}\parallel$ of the microstructure depicted in Figure 6 under given loading conditions for $\beta =0.1\mathrm{mV}$ . Left: damage evolution is not initiated $(\alpha =0.002)$ . Right: deformation-induced damage evolution at the interface leads to a decrease in the effective conductivity $(\alpha =0.05)$ . Electric current density vectors are indicated by black arrows.

6. Conclusion

In this work, a computational multiscale formulation for electrical conductors that feature material interfaces at the microscale is proposed and a finite element implementation is discussed. Several numerical examples are studied in two-dimensional settings to show the applicability of the proposed formulation: (1) The finite element results are compared against analytical solutions in order to validate the proposed formulation. (2) Different microstructures are analysed to show the influence of the material interfaces on the effective macroscale conductivity tensor. (3) The characteristic size effect that is associated with the presence of material interfaces is revealed, and it is observed that (as opposed to classic first-order homogenisation approaches) the macroscale response is not independent of the geometric dimensions of the RVE. (4) The influence of mechanically induced interfacial degradation processes on the effective electrical conductivity is exemplarily studied to demonstrate the capabilities of the proposed formulation for coupled problems. In this regard, it is important to note that the constitutive relations used are purely academic at this stage and not based on experimental data. They were selected to show the fully coupled behaviour and the applicability of the proposed framework. The application to experiments will, accordingly, be the focus of future work.