Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

Surface electrostatics: theory and computations G. Chatzigeorgiou, A. Javili and P. Steinmann Proc. R. Soc. A 2014 470, 20130628, published 5 February 2014

References

This article cites 35 articles

Subject collections

Articles on similar topics can be found in the following collections

http://rspa.royalsocietypublishing.org/content/470/2164/2013062 8.full.html#ref-list-1

materials science (143 articles) mechanical engineering (203 articles)

Email alerting service

Receive free email alerts when new articles cite this article - sign up in the box at the top right-hand corner of the article or click here

To subscribe to Proc. R. Soc. A go to: http://rspa.royalsocietypublishing.org/subscriptions

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

Surface electrostatics: theory and computations G. Chatzigeorgiou, A. Javili and P. Steinmann rspa.royalsocietypublishing.org

Research Cite this article: Chatzigeorgiou G, Javili A, Steinmann P. 2014 Surface electrostatics: theory and computations. Proc. R. Soc. A 470: 20130628. http://dx.doi.org/10.1098/rspa.2013.0628 Received: 23 September 2013 Accepted: 10 January 2014

Subject Areas: materials science, mechanical engineering Keywords: surface electricity, electrostatics, Maxwell stress, polarization stress, ponderomotive stress Author for correspondence: P. Steinmann e-mail: [email protected]

Chair of Applied Mechanics, University of Erlangen–Nuremberg, Egerlandstrasse 5, 91058 Erlangen, Germany The objective of this work is to study the electrostatic response of materials accounting for boundary surfaces with their own (electrostatic) constitutive behaviour. The electric response of materials with (electrostatic) energetic boundary surfaces (surfaces that possess material properties and constitutive structures different from those of the bulk) is formulated in a consistent manner using a variational framework. The forces and moments that appear due to bulk and surface electric fields are also expressed in a consistent manner. The theory is accompanied by numerical examples on porous materials using the finite-element method, where the influence of the surface electric permittivity on the electric displacement, the polarization stress and the Maxwell stress is examined.

1. Introduction Boundary surfaces typically present different properties than those of the bulk materials because of to various reasons (surface oxidation, ageing, coating, atoming rearrangement, termination of atomic bonds). Surface triggered phenomena are dominant mechanisms in micro and nanomaterials where the surface-to-volume ratio is large. Phenomenological models that consider surface effects on the mechanical behaviour of materials and introducing surfaces with their own structure have been studied firstly by Gibbs. In more recent contributions, Gurtin & Murdoch [1] described surface effects with the help of tensorial surface stresses. Moeckel [2] developed an alternative approach at the same period, based on the concept of a moving thermomechanical two-sided surface. Later, Daher & Maugin [3] invoked the method of virtual power to endow a surface with its own thermodynamic constituents. For applications of surface elasticity theory

2014 The Author(s) Published by the Royal Society. All rights reserved.

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

Direct notation is adopted throughout. Occasional use is made of index notation, the summation convention for repeated indices being implied. The scalar product of two vectors a and b is denoted a · b = [a]i [b]i . The tensor product of two vectors a and b is a second-order tensor D = a ⊗ b with [D]ij = [a]i [b]j . The dot product of a vector a with a second-order tensor A is given by [a · A]j = [A]ij [a]i . In the following, ε denotes the third-order permutation tensor with components εijk . Kronecker delta is denoted as δij while i denote the second-order identity tensor, i.e. [i]ij = δij . The latin letter indices take values 1, 2 or 3, while the greek letter indices take values 1 or 2. The preliminary mathematical tools from the differential geometry essential for this work are contained in table 3 of appendix A. The free space electric permittivity constant is denoted 0 . Unless stated otherwise quantities defined on the surface are differentiated from those in the bulk by a hat placed above the symbol. That is, {ˆ•} refers to a variable corresponding to the surface, its counterpart in the bulk being {•}. The basic variables used in this work are introduced in table 1.

...................................................

(a) Notation and definitions

2

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

in nanomaterials, see e.g. Cammarata [4], Dingreville & Qu [5], He & Lilley [6], Duan et al. [7] among many others. An important extension of the surface elasticity model to account for the flexural resistance of the surface was developed by Steigmann & Ogden [8]. Park et al. [9], Park & Klein [10] developed an alternative continuum framework based on the surface Cauchy–Born model, an extension of the classical Cauchy–Born model, to include surface stresses. Our own contributions in this field include [11] whereby the balance laws for a mechanical problem accounting for energetic boundaries were obtained in a variational setting based on potential energies. Later Javili & Steinmann [12] studied the thermomechanical response of materials with energetic surfaces and specified restrictions on the surface constitutive laws arising from the second law of thermodynamics. Javili & Steinmann [13–15] proposed a novel finite-element framework to model the surface elasticity and surface thermoelasticity theories. A unifying review of various approaches accounting for mechanical and/or thermal mechanisms of surface, interface and curve energies was presented in Javili et al. [16]. The interfaces between particles and matrix material play a significant role in the relative permittivity [17,18]. It has been shown experimentally that the electric performance of a particle can be improved with surface modification [19–22]. Also it has been observed that the increase in the area-to-volume ratio (i.e. decrease in the particle and/or grain size) can influence substantially the dielectric constant [23,24]. With regard to nanoporous materials, studies have shown that the dielectric constant is overpredicted by the classical micromechanics approaches [25,26]. A possible explanation for such behaviour could be the simplified approximations of the local electric field at the void surface [26]. Moreover, the influence of the independent surface behaviour on piezoelectric materials is studied in Dai et al. [27]. The aim of this work is to study, from a phenomenological point of view, the electric response of materials with boundary surfaces that possess their own electrostatic constitutive behaviour. The introduction of an energetic boundary allows to account for the influence of the surfaceto-volume ratio. Conceptually speaking, the energetic boundary surface resembles a thin layer of dielectric material with large dielectric constant. To this end, surface electric phenomena are described in the spirit of surface mechanics. The structure of the manuscript is as follows. First, notation and definitions are briefly introduced. Then in §2 the theoretical aspects of electric behaviour of continua with surface structure is elaborated. The presented theory in §2 is based on a variational approach. This theory is underpinned by the Maxwell equations given in appendix A. That is, it is possible to alternatively establish the same theory using pillbox arguments. This procedure is carefully explained in appendix A. The theory is then elucidated via a series of numerical examples in §3. Section 4 concludes the paper and discusses further extensions of the current work.

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

3

n

Figure 1. Material with energetic boundary. (Online version in colour.)

Table 1. Notation used in this manuscript. ..........................................................................................................................................................................................................

position vector



n ϕ

surface normal vector electric scalar potential

 f

bound charge density free charge density

p ppon

polarization stress ponderomotive stress

e d

electric field electric displacement

p 3

polarization vector three-dimensional Euclidean space

bpon cpon

ponderomotive body forces ponderomotive body moments

x

pmax

total charge density b

pol

Maxwell stress

..........................................................................................................................................................................................................

Consider a continuum body that occupies the configuration B ⊂ 3 within the free space occupying the open set O. The boundary surface of the body is denoted by ∂B, while the boundary unit vector is denoted n as illustrated in figure 1. The jump of a bulk quantity {•} across the surface and average of {•} over the surface are identified as [[{•}]] = {•}+ − {•}− and {{{•}}} = 12 [{•}+ + {•}− ], respectively. The terms {•}+ and {•}− are quantities in the free space O and in the body B adjacent to the boundary surface, respectively. The following identities hold for two bulk quantities a and b: [[a + b]] = [[a]] + [[b]], {{a + b}} = {{a}} + {{b}}, [[a · b]] = {{a}} · [[b]] + [[a]] · {{b}}, [[a ⊗ b]] = {{a}} ⊗ [[b]] + [[a]] ⊗ {{b}}.

2. Electric behaviour of bodies with energetic surface The electrostatic behaviour of a body under the action of an electric field is well identified in the electromagnetic theory [28–31]. With regard to an energetic surface attached to a body though, the development of conservation and constitutive laws requires to specify key assumptions on the electric variables that lie on the surface. Thus, we proceed with the following main hypotheses: (1) The boundary surface electric variables (electric field, electric displacement and polarization) have only tangential components. (2) We consider that the energetic nature of the boundary does not contribute to the free space electric energy. The first hypothesis states that eˆ ≡ eˆ · iˆ , and furthermore, it implies that

dˆ ≡ dˆ · iˆ

and pˆ ≡ pˆ · iˆ ,

eˆ := e · iˆ .

(2.1)

(2.2)

...................................................

configuration

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

x

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

In appendix A, we show that the non-existence of a normal component to the surface electric field causes the tangential part of the (bulk) electric field to be continuous on the boundary surface, i.e.

and therefore equation (2.2) follows from (2.1). This hypothesis is also supported theoretically by the work of Barham et al. [32] in the case of magnetoelasticity, when the magnetic field is connected with a scalar potential and the general three-dimensional problem is reduced to twodimensions through a membrane approximation. The second hypothesis (no boundary surface free space electric energy) states that the purely electric part of the free energy is only due to the free space and is not affected by surface electricity, which is exclusively due to the material nature of the boundary surface. Thus, the total internal energy of the energetic surface does not include a free space contribution, as opposed to the bulk material.

(a) Variational approach The variational approach that we develop in this section allows us to compute the field equations by identifying stationary points of a total potential energy functional. For better insight on energy formulations of electrostatics and energy formulations of continuum magneto-electroelasticity for identifying the relationship between the field equations, variational principles and thermodynamics, the interested reader is referred to Liu [33,34]. Based on the electric scalar potential ϕ, one can identify the variational format of the electric problem with (electrostatic) energetic surface. Recalling the usual relation for bulk materials that connects electric field and scalar potential via e = −grad ϕ,

(2.4)

we can similarly define the analogous relation for the boundary surface  ϕˆ eˆ = −grad ϕˆ · iˆ = −grad

with ϕˆ = ϕ|∂B .

(2.5)

Using thermodynamic considerations, we define the total potential energy as the sum of internal and external potential energies including the free space contribution in the body. The total potential energy for the bulk material is expressed as V tot (ϕ, e) = V int (e) + V ext (ϕ)

in B,

(2.6)

on ∂B.

(2.7)

and for the boundary surface as ˆ e) ˆ = Vˆ int (e) ˆ + Vˆ ext (ϕ) ˆ Vˆ tot (ϕ, We also consider the free space electric energy MO (e) = − 12 0 e · e in O.

(2.8)

According to the second hypothesis, the free space electric energy in the boundary surface is considered zero, i.e.1 ˆ = 0 on ∂B. M (2.9) The reason for omitting this term is because it is already included in the free space energy (2.8) of the bulk (recall that eˆ is the projection of e on the boundary according to (2.2)). 1 If one would like to account for an independent free space electric energy in the boundary, then it could be expressed in a similar manner to (2.8) as ˆ = M( ˆ e) M ˆ = − 12 ˆ0 eˆ · eˆ on ∂B,

where ˆ0 should be a surface parameter different than 0 . In such formalism, our second hypothesis implies that the parameter ˆ0 is zero.

...................................................

(2.3)

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

[[e]] × n = 0 on ∂B,

4

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

B

∂B

This last relation is expressed as    δI = [−δe · dO ] dv + [−δe · d + δϕf ] dv + O

B

∂B

[−δ eˆ · dˆ + δϕ ˆ f ] ds,

(2.12)

where we considered the following constitutive equations: dO = − d=−

∂MO ∂e

∂V int ∂e

∂ Vˆ int and dˆ = − ∂ eˆ

in O, in B on ∂B

and f =

∂V ext ∂ϕ

and ˆ f =

in B

∂ Vˆ ext ∂ϕ

⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ on ∂B.⎪

(2.13)

From the relations (2.4) and (2.5), the variation can be written only with respect to the electric potential, and renders     δϕ · dˆ + δϕ ˆ f ] ds. δI = grad δϕ · dO dv + [grad δϕ · d + δϕf ] dv + [grad (2.14) O

B

∂B

Using the bulk and surface divergence theorems (table 3) and taking into account that (i) dˆ is a tangential vector and (ii) the boundary surface is closed, we eventually obtain     dˆ − ˆ f ] ds. δϕ div dO dv − δϕ[div d − f ] dv − δϕ[[dO − d] · n + div δI = − O

B

∂B

Since δI vanishes for any arbitrary δϕ, we obtain the governing equations in local form div dO = 0

⎫ ⎪ ⎪ ⎪ ⎬

in O,

div d − f = 0 in B and

⎪ ⎪ ⎪  dˆ − ˆ = 0 on ∂B.⎭ [dO − d] · n + div

(2.15)

f

Relations (2.15) constitute the general system of electric equations that describe the free space, the bulk material and the boundary surface. By neglecting the surface electric field, equations (2.15)3 reduce to the classical boundary conditions in terms of the normal jump of the electric displacement. The same conclusions can be established if alternatively we examine the Maxwell equations using pillbox arguments (for details see appendix A). As a final remark we note that, even if we consider that there is a normal component of the electric displacement, this component will not participate in equation (2.15) (see appendix A, or in the case of thin films, eqn (223) of [36]). 2 3

Mathematically speaking, one could think of ϕ as the only purely independent variable and refine the derivations here, too.

In a similar manner as in Kankanala & Triantafyllidis [35] approach for magnetostatics, we split the electric displacement in the free space (and consequently the electric field and the Maxwell stress) in two parts, the applied electric displacement and the perturbed electric displacement owing to the presence of the solid. The perturbed electric quantities in the free space are assumed to decay and eventually vanish far away from the solid. In the sequel, we will omit the terms that are related to the applied electric field at ||x|| → ∞.

...................................................

O

for which we want to identify a stationary point for all admissible variations of the independent variables, i.e. δI = 0 ∀ δϕ, ∀ δe, ∀ δ e. ˆ (2.11)

5

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

Next, considering as independent variables2 the potential ϕ, the bulk electric field e and the surface electric field e, ˆ we identify the total potential energy functional3    tot I(ϕ, e, e) ˆ = MO (e) dv + V (ϕ, e) dv + ˆ ds, (2.10) Vˆ tot (ϕ, e)

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

(b) Polarization

6

in B.

(2.16)

For the boundary surface, we consider that the polarization, the electric field and the electric displacement are connected through the relation4 dˆ − pˆ = 0 on ∂B,

(2.17)

which means that the surface electric displacement is equal to the surface polarization. The motivation for equation (2.17) arises from the tangential continuity of the electric field across the boundary surface. Using equation (A 12) of appendix A in conjunction with (2.16), we have that n × 0 [[e]] = n × [[d]] − n × [[h]] = 0 on ∂B, which states that the tangential parts of the bulk electric displacement and the bulk polarization across the surface should be equal.

(c) Resultant bound charge and electric moment Following classical approaches [31], the total charge density inside a body is considered as the sum of the bound charge density, which appears due to the material, and the free charge density which is due to moving charges. Thus, using equations (2.15)–(2.17), we have

b

b := −div p

in B,

 pˆ ˆ b := ˆ bd + ˆ bs = −[[p]] · n − div

f := div d

in B,

f f  dˆ ˆ f := ˆ d + ˆ s = [[d]] · n + div

f

and  :=  +  = 0 div e

in B,

b

ˆ := ˆ + ˆ

f

= 0 [[e]] · n = ˆ bd

on ∂B, on ∂B

f + ˆ d

= ˆ d

⎫ ⎪ ⎪ ⎪ ⎬

⎪ ⎪ ⎪ ⎭ on ∂B, (2.18)

where the subscript d denotes the jump of a bulk quantity and the subscript s denotes the contribution of the energetic surface. Using the divergence theorems, the equations (2.18) and noticing that (i) the polarization has a meaning only inside the material, i.e. [[p]] = −p, (ii) the boundary surface polarization pˆ is a tangential vector and (iii) the boundary is a closed surface, we can show that: 1. The presence of a boundary surface with different constitutive behaviour than the material does not influence the general requirement that the total bound charge in an isolated body vanishes, i.e. it holds      p] b dv + ˆ b ds = − div p dv − [[[p]] · n + div ˆ ds B

∂B

 =

4

B

∂B

∂B

[−p · n + p · n] ds +

 ∂B

kˆpˆ · n ds = 0.

(2.19)

Strictly speaking, the analogous relation of (2.16) for surfaces is written as ˆ0 eˆ = dˆ − pˆ

on ∂B,

where ˆ0 is a surface parameter different than 0 . As we already mentioned, the second hypothesis implies that ˆ0 is zero, thus the last relation leads to (2.17).

...................................................

0 e = d − p

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

The relationship between the polarization, the electric field and the electric displacement inside the bulk is written as

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

∂B

B

 =

B

∂B



+  =  =

B B



[p − div(x ⊗ p)] dv +

∂B

∂B

[x ⊗ p] · n ds

 ⊗ p)] [pˆ · iˆ − div(x ˆ ds 

p dv +  p dv +

 ∂B ∂B

pˆ ds +

∂B

kˆ [x ⊗ p] ˆ · n ds

pˆ ds.

(2.20)

(d) Ponderomotive force and moment densities The ponderomotive body forces bpon in the bulk material and the tractions (or alternatively surface forces) bˆ pon at the boundary surface satisfy     e dv + {{e}} ˆ ds = bpon dv + (2.21) bˆ pon ds. B

∂B

∂B

B

The term in the second integral of the left-hand side represents the contribution of the boundary as a surface with discontinuity in the tangential part of the electric field [37].5 After proper calculations, the resultant ponderomotive force densities are expressible as f bpon = grad e · p + f e and bˆ pon = d {{e}} + [[[e]] ⊗ {{p}}] · n.

(2.22)

Note that in expression (2.22)2 the boundary surface ponderomotive force depends not on the f

total boundary surface free charge density ˆ f , but only on its part ˆ d that arises from the normal discontinuity of d across the boundary. The ponderomotive body force is connected with a ponderomotive stress tensor, which consists of two parts, the Maxwell stress and the stress due to the polarization. Thus, we can express the bulk ponderomotive stress tensor in the form ppon = pmax + ppol

and pmax = 0 e ⊗ e − 12 0 [e · e]i

and ppol = e ⊗ p.

(2.23)

Using the expressions (2.23), we obtain div ppon = bpon

and

[[ppon ]] · n = bˆ pon .

Following the same approach, the ponderomotive moment densities satisfy     x × [e] dv + x × [{{e}}] ˆ ds = cpon dv + cˆpon ds, B

∂B

B

∂B

(2.24)

where cpon and cˆpon denote the ponderomotive bulk and boundary surface moments, respectively. Thus,     [x × bpon + p × e] dv + x × bˆ pon ds = cpon dv + cˆpon ds, B

∂B

leading to cpon = x × bpon + p × e,

B

∂B

cˆpon = x × bˆ pon .

(2.25)

A summary of the equations and expressions that are related to both the bulk and the boundary surface is given in table 2. 5

Since {{e}} ˆ = {{e ˆ · iˆ⊥ }} + ˆ eˆ

on ∂B,

the contribution of the boundary surface to the ponderomotive force is already included in equation (2.21).

...................................................

B

7

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

2. The material nature of the boundary influences the electric moment, since it includes the contribution of the bulk polarization p and the surface polarization p, ˆ i.e.      p]x b x dv + ˆ b x ds = − x div p dv − [[[p]] · n + div ˆ ds

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

Table 2. Bulk material and energetic boundary electrostatic equations. on boundary (∂B) ϕ eˆ = −grad ..........................................................................................................................................................................................................  dˆ first Maxwell equation f = div d ˆ f = [ d]] · n + div ..........................................................................................................................................................................................................

second Maxwell equation

curl e = 0

constitutive law

d = 0 e + p = −∂e V int

n × [ e]] = 0 dˆ = pˆ = −∂eˆ Vˆint

bound charge density

b = −div p

 pˆ ˆ b = −[[p]] · n − div

total charge density

 = 0 div e

ˆ = 0 [ e]] · n

polarization stress

p =e⊗p

pˆpol = 0

Maxwell stress

pmax = 0 e ⊗ e − 1 0 [e · e]i

pˆmax = 0

ponderomotive stress

ppon = pmax + ppol

pˆpon = 0

.......................................................................................................................................................................................................... .......................................................................................................................................................................................................... .......................................................................................................................................................................................................... ..........................................................................................................................................................................................................

pol

.......................................................................................................................................................................................................... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..........................................................................................................................................................................................................

Before we close the theoretical section we would like to point out that, from a thermodynamic point of view, our theory assumes that the energetic surface supplies electrical power without any mechanical effect, as opposed to the bulk material, in which we have also mechanical power due to the polarization. A proper thermodynamic treatment though of the current problem would require to write our equations taking into account large deformation formulations, i.e. identifying the deformation gradient and its Jacobian for connecting spatial and material description. This would exceed the scope of the current contribution.

3. Numerical examples In order to carry out a finite-element implementation of the proposed theory, one would need to employ the weak form of balance equations (2.15). The procedure of obtaining the weak form from a strong form is standard in the context of the finite-element method [38]. Here, the weak form is obtained by multiplying the strong form (2.15) by appropriate test functions of proper Sobolev spaces in the bulk and on the surface and applying the (bulk) divergence theorem in the bulk and surface divergence theorem on the surface. This procedure would lead to a weak form resembling the variational form (2.14) and therefore, we do not repeat it. The numerical implementation of the problem accounting for the surface conservation laws and constitutive equations is similar to the one developed in the case of thermoelastic solids [39]. Thus in this work, we do not describe the finite-element formulation. In the sequel, we present numerical examples for a two-dimensional toy example. Clearly, the theory is not limited to two dimensions and holds for arbitrary geometries in three dimensions. This simple example is chosen to illustrate key features of the proposed theory. Consider a (linear) electric body as shown in figure 2 with square shape of length  that includes a circular void at its centre. The surface of the void is considered to present linear electric behaviour. The (thermodynamically consistent) constitutive relations that are used for both the bulk and the surface of the void are expressed as ˆ d = 0 e and dˆ = 0 ˆe.

(3.1)

In the above expressions, 0 is the electric permittivity of free space, while  and ˆ denote the relative dielectric constants of the body and the energetic surface, or rather interface, respectively. In the numerical analyses, we consider a solid material with  = 10, i.e. 10 times the free space. The dielectric constant of the void is clearly the same as the free space, i.e. equal to one. The surface relative dielectric constant varies from 0 to 100. In terms of boundary conditions, we choose zero normal electric displacement on the left and right side, while on the top and the

...................................................

in bulk (B) e = −grad ϕ

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

equation/expression electric field

8

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

non-homogeneous Dirichlet-type boundary conditions homogeneous Neumann-type boundary conditions

homogeneous Neumann-type boundary conditions

=1 y

homogeneous Dirichlet-type boundary conditions

x

Figure 2. Model and boundary conditions. Mesh quality, dimensions and material properties are shown. (Online version in colour.)

= 0.0

= 0.1

= 0.2

= 0.3

= 0.4

electric potential j

= 0.5

= 1.0

= 2.0

= 5.0

= 10

10.0 9.00 8.00 7.00 6.00 5.00 4.00 3.00 2.00 1.00 0

Figure 3. Distribution of the scalar electric potential. (Online version in colour.)

bottom side we apply 0 and 10 V electric potential, respectively. The domain is discretized using 2000 bilinear quadrilateral elements. The mesh quality of the finite-element method and details about the boundary conditions are presented in figure 2. In order to better understand the influence of the surface electric behaviour on the overall response of the body, in a series of numerical examples, distributions of the quantities are illustrated in figures 3–11. Note that / ˆ has the dimension of length, i.e. reducing the size  by a factor f is equivalent to increase the / ˆ by the same factor f . More importantly, for all stresses the results depend on the length squared 2 which indicates a strong size effect. The results here are given for  = 1 and normalized by multiplying with 2 for stresses and with  for electric fields due to similarity reasons. Also, all stress values are divided by 0 . The results in figures 3–11 are carefully produced in a self-explanatory manner. Here, for more clarity some of the key features are explained. In figure 3, we observe that the distribution of the scalar potential in the area close to the void depends strongly on the surface dielectric constant. When the surface of the void does not have independent material behaviour (ˆ = 0) the potential lines tend to enter the void. In all the following results, this (top-left) illustration shall be understood as the classical response of the example of interest with no additional surface constitutive behaviour. Increasing the surface dielectric constant ˆ varies the distributions of electric potential lines. For / ˆ = 10, one observes a very strong surface effect such that the

...................................................

0.4

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

matrix = 10

9

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

= 0.0

= 0.1

= 0.2

= 0.3

= 0.4

= 1.0

= 2.0

= 5.0

= 10

Figure 4. Distribution of the electric field in x-direction. (Online version in colour.)

= 0.0

= 0.1

= 0.2

= 0.3

= 0.4

electric field y×

= 0.5

= 1.0

= 2.0

= 5.0

= 10

21.16 19.04 16.93 14.81 12.70 10.58 8.46 6.35 4.23 2.12 0

Figure 5. Distribution of the electric field in y-direction. (Online version in colour.)

= 0.0

= 0.1

= 0.2

= 0.3

= 0.4

Maxwell stress 2 pmax xx ×

= 0.5

= 1.0

= 2.0

= 5.0

= 10

231.06 207.86 184.66 161.46 138.26 115.06 91.85 68.65 45.45 22.25 –0.95

Figure 6. Distribution of the Maxwell stresses in xx-direction. (Online version in colour.)

potential lines tend to move away from the void.6 We also see that for / ˆ = 0.2 the void becomes (almost) invisible to the scalar potential, leading to an (almost) homogeneous behaviour (see also figures 4 and 5). Thus, electrostatic energetic surfaces allow to create ‘invisible’ metamaterials. Figure 4 demonstrates the distribution of the electric field in the direction normal to the applied electric potential difference. The presence of a material void surface alters significantly the ex close 6

The surface dielectric constant is given in terms of the bulk matrix dielectric constant, i.e. / ˆ = 10 corresponds to ˆ = 100.

...................................................

= 0.5

9.68 7.74 5.81 3.87 1.94 0.00 –1.94 –3.87 –5.81 –7.74 –9.68

10

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

electric field x×

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

= 0.0

= 0.1

= 0.2

= 0.3

= 0.4

11

= 1.0

= 2.0

= 5.0

= 10

Figure 7. Distribution of the Maxwell stresses in yy-direction. (Online version in colour.)

= 0.0

= 0.1

= 0.2

= 0.3

= 0.4

polarization stress 2 ppol xx ×

= 0.5

= 1.0

= 2.0

= 5.0

= 10

837.56 753.49 669.41 585.34 501.26 417.19 333.12 249.04 164.97 80.89 –3.18

Figure 8. Distribution of the polarization stresses in xx-direction. (Online version in colour.)

= 0.0

= 0.1

= 0.2

= 0.3

= 0.4

polarization stress 2 ppol yy ×

= 0.5

= 1.0

= 2.0

= 5.0

= 10

4027.80 3624.27 3220.73 2817.20 2413.66 2010.13 1606.59 1203.06 799.52 395.99 –7.55

Figure 9. Distribution of the polarization stresses in yy-direction. (Online version in colour.)

to the void, since for / ˆ = 0 to / ˆ = 0.2 the electric field tends to become homogeneous and for / ˆ = 0.2 to / ˆ = 10 it becomes again inhomogeneous close to the void with opposite signs compared to / ˆ = 0. Similar behaviour is observed for the electric field in the direction parallel to the applied electric potential difference in figure 5. Here, we also observe that for / ˆ = 0.2 the electric field ey inside the void almost vanishes. Regarding the Maxwell stress distributions, a strong void surface in terms of electric properties causes almost disappearance of stress inside the void and a very high stress concentration close to

...................................................

= 0.5

671.62 604.32 537.01 469.71 402.40 335.10 267.79 200.49 133.18 65.88 –1.43

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

Maxwell stress 2 pmax yy ×

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

= 0.0

= 0.1

= 0.2

= 0.3

= 0.4

= 1.0

= 2.0

= 5.0

= 10

Figure 10. Distribution of the ponderomotive stresses in xx-direction. (Online version in colour.)

= 0.0

= 0.1

= 0.2

= 0.3

= 0.4

ponderomotive stress 2 ppon yy ×

= 0.5

= 1.0

= 2.0

= 5.0

= 10

4699.40 4228.56 3757.72 3286.89 2816.05 2345.21 1874.37 1403.53 932.70 461.86 –8.98

Figure 11. Distribution of the ponderomotive stresses in yy-direction. (Online version in colour.)

the surface. This phenomenon holds for both normal and parallel directions to the applied electric potential difference (figures 6 and 7). The distribution of the polarization stress follows a similar pattern with the distribution of the electric fields according to figures 8 and 9. The distribution of the ponderomotive stress, i.e. the sum of the Maxwell and the polarization stress, is shown in figures 10 and 11, where we observe that for / ˆ > 2 there is a significant stress concentration around the void.

4. Conclusion In this paper, we have identified the governing equations that describe the behaviour of a body with (electrostatic) energetic boundary surface. The assumption of surface electric quantities that are tangential to the surface has the following advantages: (i) it allows the surface electric field to be described in terms of the usual electric scalar potential and thus to connect it with the bulk electric field, (ii) it provides a correct representation of the surface bound charges satisfying the classical principle in electrostatics that the total bound charge in an isolated body vanishes and (iii) the surface polarization pˆ represents electric moment per unit area. In the developed framework, we also observe that the surface electric variables do not provide additional Maxwell and polarization surface forces. An important application of the proposed framework is to study the overall response of microand nanoporous materials whereby the surface plays a significant role due to a more dominant surface-to-volume ratio. In order to do so, one needs to set up a micro-to-macro transition

...................................................

= 0.5

1028.60 925.33 822.05 718.78 615.50 512.23 408.96 305.68 202.41 99.13 –4.14

12

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

ponderomotive stress 2 ppon xx ×

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

Table 3. Tensor properties in bulk material and on energetic surface. surface

divergence theorem   ∂L {•} · m ds = L div {•} dv

divergence theorem    ˆ  ˆ ∂I {•} · ms dl = I div {•} ds + I k{•} · n ds

Stoke’s theorem   ∂S u · mt dl = S curl u · mn ds

Stoke’s theorem    ∂S u · mt dl = S curl u · mn ds

..........................................................................................................................................................................................................

..........................................................................................................................................................................................................

framework to calculate the overall response from the microscopic behaviour. This extension will be the subject of our further contributions to this field. Funding statement. The support of this work by the ERC Advanced Grant MOCOPOLY is gratefully acknowledged. The first author would like to acknowledge the support by the King Abdullah University of Science and Technology (KAUST) project NumPor.

Appendix A. Maxwell laws It is enlightening to show how the governing equations (2.15) can be alternatively obtained by examining the Maxwell equations using pillbox arguments. Table 3 summarizes certain key definitions and properties of differential geometry necessary in the following derivations (for further details, see [16] and references therein). Let xq denote the Cartesian coordinates, while θˆα denotes the curvilinear surface coordinates. Moreover, gq and gq denote Cartesian covariant and contravariant base vectors in the bulk, while gˆα and gˆα denote curvilinear covariant and contravariant base vectors on the surface, respectively. Some of the vectors used in the table are shown in figure 12.

(a) First Maxwell equation We examine the following cases for the first Maxwell equation. 1. We consider the control region of figure 12a. The canonical control volume region is defined as one which includes a part of the material, i.e. L ⊂ B. L has a smooth boundary ∂L with unit normal vector m. In the case of the electrostatic problem, the second Maxwell equation reads   d · m ds = f dv. (A 1) ∂L

L

From the volume divergence theorem, we obtain  [div d − f ] dv = 0, L

(A 2)

which upon localization leads to the local equation div d = f

in B.

(A 3)

...................................................

q ˆα gˆα := ∂x , gˆα := ∂ θ gq := ∂xq , gq := ∂x∂x α . . . . . . . . . . . . .∂x . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .∂ . . θ.ˆ. . . . . . . . . . . . . . . . . . . . . ∂x ..................................................................... q δs = gq · gs δ α = gˆα · gˆβ , gˆ3 = gˆ1 × gˆ2 , n = gˆ1 × gˆ2 /|gˆ1 × gˆ2 | . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .β. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . q i := δs gq ⊗ gs = gq ⊗ gq iˆ := δ α gˆα ⊗ gˆβ = gˆα ⊗ gˆα = i − n ⊗ n, iˆ⊥ := n ⊗ n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .β ............................................................................................. q  n grad{•} := ∂{•} ⊗ g = grad{•} · i grad{•} := ∂{•}α ⊗ gˆα = grad{•} · iˆ , kˆ := −grad q . . . . . . . . . . . . . . . . . . . . . . . ∂x . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .∂. .θˆ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ∂{•} ∂{•}  {•} :=  n div{•} := q · gq = grad{•} : i div · gˆα = grad{•} : iˆ , kˆ := −div α . . . . . . . . . . . . . . . . . . . . .∂x . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .∂. .θ.ˆ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ∂{•} ∂{•}  := −  {•} : ε, ˆ εˆ = ε curl{•} := − q × gq = −grad {•} : ε curl{•} × gˆα = −grad α . . . . . . . . . . . . . . . . . . . . . . . . . . ∂x . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .∂ . .θ .ˆ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

bulk

13

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

(a)

(b)

14 ...................................................

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

I

I

m

n

ˆt m

L mˆ s

L

(c)

(d) n mn

n mt

nˆ

n

Figure 12. (a) Control volume region in the bulk, (b) control surface region on the boundary, (c) control surface region in the bulk, (d) control line region on the boundary. (Online version in colour.)

2. In order to identify the continuity condition between the bulk and the boundary surface, we follow the procedure described by Kovetz [31]. We consider a small cylinder passing through the boundary surface and including both bulk material and free space (figure 12b). The cylinder crosses the boundary on the surface I, bounded by the curve ˆ s (tangent to the surface) and tangent unit ∂I. The curve ∂I has normal unit vector m ˆ t . We apply the second Maxwell equation on this cylinder and we consider the vector m limit as it collapses onto I. In this case, we have    ˆ s dl = ˆ f ds. [[d]] · n ds + (A 4) dˆ · m I

∂I

I

ˆ even if it exists, does not contribute to the last relation, We note that a normal part of d, since the term dˆ · n at the top face of I cancels out with the dˆ · [−n] term at the bottom face of I. Using the surface divergence theorem, we can rewrite the last equation as    dˆ − ˆ f ] ds + kˆ dˆ · n ds = 0. [[[d]] · n + div (A 5) I

I

The second integral vanishes since dˆ is tangential, thus the arbitrariness of I allows us to identify the local format as  dˆ = ˆ f [[d]] · n + div

on ∂B.

(A 6)

The final result dictates that the boundary surface free charge density can be split in f

two parts: the first part ˆ d = [[d]] · n is due to the normal discontinuity of d across the f  dˆ is due to the independent behaviour of the boundary, while the second part ˆ s = div boundary surface.

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

(b) Second Maxwell equation

15

∂S

which, by using the Stokes’ theorem, rewrites as  curl e · mn ds = 0. S

(A 8)

Since S is an arbitrary surface in the bulk material and mn is an arbitrary vector, the local form of (A 8) leads to curl e = 0 in B. (A 9) 2. In order to identify the continuity condition between the bulk and the boundary surface, we follow the procedure described by Kovetz [31]. We consider a plane perpendicular to ˆ which is tangential ∂B (figure 12d). The orientation of the plane is defined by its normal n, to the surface. We apply the (first) Maxwell equation on this plane and we consider the limit as it collapses onto ∂B in a single curve C with tangent unit vector n˜ = nˆ × n and boundary points ∂C. In this case, we have  [[e]] · n˜ dl − eˆ · n = 0. (A 10) C

∂C

In the above relation, we note that in a closed loop around C the term [eˆ · iˆ ] · n˜ above the ˜ below the curve. Since eˆ is tangential, the curve is cancelled with the term [eˆ · iˆ ] · [−n] second term of equation (A 10) vanishes and the equation reads  [n × [[e]]] · nˆ dl = 0. (A 11) C

Since C is arbitrary, we derive the jump condition in local form as n × [[e]] = 0 on ∂B.

(A 12)

The last equation states that the tangential part of the electric field is not affected by the boundary surface.

References 1. Gurtin ME, Murdoch AI. 1975 A continuum theory of elastic material surfaces. Arch. Rational Mech. Anal. 57, 291–323. (doi:10.1007/BF00261375) 2. Moeckel GP. 1975 Thermodynamics of an interface. Arch. Rational Mech. Anal. 57, 255–280. (doi:10.1007/BF00280158) 3. Daher N, Maugin GA. 1986 The method of virtual power in continuum mechanics application to media presenting singular surfaces and interfaces. Acta Mechanica 60, 217–240. (doi:10.1007/BF01176354) 4. Cammarata RC. 1994 Surface and interface stress effects in thin films. Progr. Surf. Sci. 46, 1–38. (doi:10.1016/0079-6816(94)90005-1) 5. Dingreville R, Qu J. 2005 Surface free energy and its effect on the elastic behavior of nano-sized particles, wires and films. J. Mech. Phys. Solids 53, 1827–1854. (doi:10.1016/j.jmps.2005.02.012) 6. He J, Lilley CM. 2008 Surface effect on the elastic behavior of static bending nanowires. Nano Lett. 8, 1798–1802. (doi:10.1021/nl0733233) 7. Duan HL, Wang J, Karihaloo BL. 2009 Theory of elasticity at the nanoscale. Adv. Appl. Mech. 42, 1–68. (doi:10.1016/S0065-2156(08)00001-X)

...................................................

1. We assume a surface S with boundary curve ∂S and normal unit vector mn , exclusively embedded in the bulk material. The boundary curve of the bulk has tangent unit vector mt (figure 12c). Then the second Maxwell equation in the case of electrostatics reads  e · mt dl = 0, (A 7)

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

We examine the following cases for the second Maxwell equation.

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

16 ...................................................

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

8. Steigmann DJ, Ogden RW. 1999 Elastic surface–substrate interactions. Proc. R. Soc. Lond. A 455, 437–474. (doi:10.1098/rspa.1999.0320) 9. Park HS, Klein PA, Wagner GJ. 2006 A surface Cauchy–Born model for nanoscale materials. Int. J. Numer. Methods Eng. 68, 1072–1095. (doi:10.1002/nme.1754) 10. Park HS, Klein PA. 2007 Surface Cauchy–Born analysis of surface stress effects on metallic nanowires. Phys. Rev. B 75, 1–9. (doi:10.1103/PhysRevB.75.085408) 11. Steinmann P. 2008 On boundary potential energies in deformational and configurational mechanics. J. Mech. Phys. Solids 56, 772–800. (doi:10.1016/j.jmps.2007.07.001) 12. Javili A, Steinmann P. 2010 On thermomechanical solids with boundary structures. Int. J. Solids Struct. 47, 3245–3253. (doi:10.1016/j.ijsolstr.2010.08.009) 13. Javili A, Steinmann P. 2009 A finite element framework for continua with boundary energies. Part I: the two-dimensional case. Comp. Methods Appl. Mech. Eng. 198, 2198–2208. (doi:10.1016/j.cma.2009.02.008) 14. Javili A, Steinmann P. 2010 A finite element framework for continua with boundary energies. Part II: The three-dimensional case. Comp. Methods Appl. Mech. Eng. 199, 755–765. (doi:10.1016/j.cma.2009.11.003) 15. Javili A, Steinmann P. 2011 A finite element framework for continua with boundary energies. Part III: the thermomechanical case. Comp. Methods Appl. Mech. Eng. 200, 1963–1977. (doi:10.1016/j.cma.2010.12.013) 16. Javili A, McBride A, Steinmann P. 2013 Thermomechanics of solids with lower-dimensional energetics: on the importance of surface, interface and curve structures at the nanoscale. A unifying review. Appl. Mech. Rev. 65, 010802. (doi:10.1115/1.4023012) 17. Todd MG, Shi FG. 2003 Characterizing the interphase dielectric constant of polymer composite materials: effect of chemical coupling agents. J. Appl. Phys. 94, 4551–4557. (doi:10.1063/1.1604961) 18. Roy M, Nelson JK, MacCrone RK, Schadler LS, Reed CW, Keefe R, Zenger W. 2005 Polymer nanocomposite dielectrics—the role of the interface. IEEE Trans. Dielectrics Electrical Insulation 12, 629–643. (doi:10.1109/TDEI.2005.1511089) 19. Park KS, Son JT, Chung HT, Kim SJ, Lee CH, Kang KT, Kim HG. 2004 Surface modification by silver coating for improving electrochemical properties of LiFePO4 . Solid State Commun. 129, 311–314. (doi:10.1016/j.ssc.2003.10.015) 20. Ma D, Hugener TA, Siegel RW, Christerson A, Mårtensson E, Önneby C, Schadler LS. 2005. Influence of nanoparticle surface modification on the electrical behaviour of polyethylene nanocomposites. Nanotechnology 16, 724–731. (doi:10.1088/0957-4484/16/6/016) 21. Wang T, Li W, Luo L, Zhu Y. 2013 Ultrahigh dielectric constant composites based on the oleic acid modified ferroferric oxide nanoparticles and polyvinylidene fluoride. Appl. Phys. Lett. 102, 092904. (doi:10.1063/1.4795128) 22. Yu K, Niu Y, Zhou Y, Bai Y, Wang H. 2013 Nanocomposites of surface-modified BaTiO3 nanoparticles filled ferroelectric polymer with enhanced energy density. J. Am. Ceramic Soc. 96, 2519–2524. (doi:10.1111/jace.12338) 23. Lee H, Kim H. 1990 Ceramic particle size dependence of dielectric and piezoelectric properties of piezoelectric ceramicpolymer composites. J. Appl. Phys. 67, 2024–2028. (doi:10.1063/1.345584) 24. Papet P, Dougherty JP, Shrout TR. 1990 Particle and grain size effects on the dielectric behavior of the relaxor ferroelectric Pb(Mg1/3 Nb2/3 )O3 . J. Mater. Res. 5, 2902–2909. (doi:10.1557/JMR.1990.2902) 25. Hedrick JL, Miller RD, Hawker CJ, Carter KR, Volksen W, Yoon DY, Trollsås M. 1998 Templating nanoporosity in thin-film dielectric insulators. Adv. Mater. 10, 1049–1053. (doi:10.1002/(SICI)1521-4095(199809)10:133.0.CO;2-F) 26. Min CK, Wu TB, Yang WT, Chen CL. 2008 Functionalized mesoporous silica/polyimide nanocomposite thin films with improved mechanical properties and low dielectric constant. Comp. Sci. Technol. 68, 1570–1578. (doi:10.1016/j.compscitech.2007.09.021) 27. Dai S, Gharbi M, Sharma P, Park HS. 2011 Surface piezoelectricity: Size effects in nanostructures and the emergence of piezoelectricity in non-piezoelectric materials. J. Appl. Phys. 110, 104305. (doi:10.1063/1.3660431) 28. Toupin RA. 1956 The elastic dielectric. J. Rational Mech. Anal. 5, 849–914. 29. Pao YH. 1978 Electromagnetic forces in deformable continua. In Mechanics today, vol. 4, pp. 209–305. Oxford, UK: Pergamon Press.

Downloaded from rspa.royalsocietypublishing.org on May 16, 2014

17 ...................................................

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130628

30. Maugin GA. 1993 Material inhomogeneities in elasticity, vol. 3. London, UK: Chapman & Hall/CRC. 31. Kovetz A. 2000 Electromagnetic theory. New York, NY: Oxford. 32. Barham M, Steigmann DJ, White D. 2012 Magnetoelasticity of highly deformable thin films: theory and simulation. Int. J. Non-Linear Mech. 47, 185–196. (doi:10.1016/ j.ijnonlinmec.2011.05.004) 33. Liu L. 2013. On energy formulations of electrostatics for continuum media. J. Mech. Phys. Solids 61, 968–990. (doi:10.1016/j.jmps.2012.12.007) 34. Liu L. 2014 An energy formulation of continuum magneto-electro-elasticity with applications. J. Mech. Phys. Solids 63, 451–480. (doi:10.1016/j.jmps.2013.08.001) 35. Kankanala S, Triantafyllidis N. 2004 On finitely strained magnetorheological elastomers. J. Mech. Phys. Solids 52, 2869–2908. (doi:10.1016/j.jmps.2004.04.007) 36. Steigmann DJ. 2013 Surface waves in elastic half spaces coated with crystalline films. In CISM course on dynamic localization phenomena in elasticity, acoustics and electromagnetism, vol. 547 (eds R. Craster & J. Kaplunov), pp. 225–256. Vienna, Austria and New York, NY: Springer. 37. Griffiths DJ. 2013 Introduction to electrodynamics, 4th edn. Harlow, Essex, UK: Pearson. 38. Hughes TJR. 1987 The finite element method: linear static and dynamic finite element analysis, vol. 65. Englewood Cliffs, NJ: Prentice Hall, Inc. 39. Javili A, McBride A, Steinmann P. 2013 Numerical modelling of thermomechanical solids with highly conductive energetic interfaces. Int. J. Numer. Methods Eng. 93, 551–574. (doi:10.1002/nme.4402)