Fundamentals of electro-mechanically coupled cohesive zone formulations for electrical conductors

Motivated by the influence of (micro-)cracks on the effective electrical properties of material systems and components, this contribution deals with fundamental developments on electro-mechanically coupled cohesive zone formulations for electrical conductors. For the quasi-stationary problems considered, Maxwell’s equations of electromagnetism reduce to the continuity equation for the electric current and to Faraday’s law of induction, for which non-standard jump conditions at the interface are derived. In addition, electrical interface contributions to the balance equation of energy are discussed and the restrictions posed by the dissipation inequality are studied. Together with well-established cohesive zone formulations for purely mechanical problems, the present developments provide the basis to study the influence of mechanically-induced interface damage processes on effective electrical properties of conductors. This is further illustrated by a study of representative boundary value problems based on a multi-field finite element implementation.


Introduction
Material interfaces can occur at different material length scales and can significantly influence the effective constitutive response of the material (system) under consideration. Typical examples for such interfaces are grain boundaries in polycrystalline materials and transition zones between different phases in multi-phase materials. In this regard, the unique properties of the interfaces may significantly differ from those of the surrounding continuum and can be accounted for in simulations by the introduction of interface models. In these approaches, the physical interface (of finite thickness) is not geometrically resolved but rather approximated as a lower-dimensional object (e.g. a surface in a three-dimensional setting). Based on the assumed continuity of field quantities across the interface, different types of interface formulations may moreover be distinguished, see [15] and references cited therein. Elastic-interface models originate from the pioneering works of Gurtin and Murdoch on interface elasticity [11,23], and assume the displacement-type fields to be continuous, whereas traction-type quantities may exhibit jump-discontinuities across the interface. Classic cohesive interface models that date back to the seminal works by Barenblatt on quasi-brittle materials [2] and by Dugdale on ductile materials [6] on the other hand, assume traction continuity across the interface and allow for the modelling of jump-discontinuities in the displacement-type fields. Although seemingly different, it has recently been shown that classic cohesive zone formulations and interface elasticity formulations can be regarded as two extremes of the unifying theory of generalised imperfect interfaces elaborated in [15,16,27,31].
Since the works of Barenblatt and Dugdale [2,6], classic cohesive zone models have been subject of intense research with many contributions focusing on the consistency of cohesive zone formulations with fundamental requirements of continuum mechanics (particularly for finite deformations), e.g. [22,26,35,38], and on the elaboration of specific traction separation laws that account for irreversible processes such as damage and plasticity, e.g. [7,12,25,30,33,39]. In this regard, a salient feature of cohesive zone damage formulations as compared to continuum damage formulations is that the fracture energy does not depend on the finite element size such that no sophisticated regularisation schemes are required [3,24]. Moreover, cohesive zone formulations can conveniently be incorporated into finite element formulations by means of mixed-type boundary conditions or cohesive zone elements, see [24] and references cited therein.
In addition to the developments on cohesive zone formulations for purely mechanical problems, cohesive zone formulations for coupled multi-physics problems have been in the focus of intense research. In particular, thermomechanical coupling has been addressed in [8,9,28] and electro-mechanically coupled cohesive zone formulations for electro-active solids have been studied in [1,18,19,34,36,37]. Focussing on piezo-and ferroelectric effects in dielectric solids, the latter formulations differ significantly from the present developments on electrical conductors. In this regard, the principal (electrical) material parameters in the established formulations for electro-active solids are the piezoand dielectric moduli, whereas the present contribution is motivated by the influence of mechanically-induced interface damage processes on the electrical conductivity, as studied experimentally in, e.g., [4,5,10].
The contribution is organised as follows: Section 2 deals with the derivation of the fundamental set of balance equations and jump-conditions at material interfaces, with particular focus being on the electrical sub-problem. Based on these developments, a finite element-based implementation of the theory is discussed in Section 3. These fundamentals serve as the basis for the study of representative boundary value problems in Section 4 that show the applicability of the proposed formulation. In particular, the finite element implementation is validated by means of analytical solutions, and the influence of intergranular fracture processes on the electrical conductivity is studied numerically.

Notation
With α, β, γ , δ denoting tensor-valued quantities of first order and with the standard dyadic product denoted by ⊗, single and double tensor contractions take the form In addition, the non-standard dyadic prod- are introduced to shorten notation. Moreover, I denotes the second order identity tensor, the third order permutation tensor and ∇•, ∇ · •, ∇ × •, indicate (right-)gradient, (right-)divergence and (right-)curl operations, respectively. Furthermore, quantities at the two opposing sides of an interface will be indicated by superscripts • − and • + such that the jump of a quantity across the Fig. 1 Specification of quantities in the control volume V and on the interface I. The interface is characterised by its unit surface normal vector n and quantities at opposing sides of the interface are indicated by superscripts • − and • + interface follows as the interfacial mean value reads and the identity holds.

Continuum thermodynamics
This section deals with the thermodynamic fundamentals of a continuum with material interfaces that is subjected to mechanical and electrical loads, see Figure 1. For the convenience of the reader, and to show similarities and differences in the derivations of the local field equations, the mechanical sub-problem is briefly recapitulated in Section 2.1 before the electrical sub-problem is studied in detail in Section 2.2. Based on these fundamentals, an extended form of the balance equation of energy is proposed in Section 2.3, and constitutive restrictions that are posed by the extended dissipation inequality are discussed in Section 2.4. The ensuing derivations are based on the assumption of quasi-static deformation processes and of a quasi-stationary electrical problem.

Mechanical sub-problem
The control volume V is assumed to be loaded by volume distributed forces f that scale in the mass density (per unit volume) ρ, by surface distributed forces f that act on the interface I and scale in the mass density (per unit surface) ρ, and by tractions that are related to the small strain stress tensor σ via the outward unit surface normal vector n. Based on the latter assumptions, and by neglecting line-distributed (tangential) forces on ∂I, the extended integral form of the (quasi-static) balance equation of linear momentum reads The third summand on the right-hand side of (5) can further be rewritten as where the classic divergence theorem was used along with the definitions and with the geometric constraint on the unit surface normal vectors at opposing sides of the interface The localisation of the ensuing equation that results from the insertion of (6) into (5), i.e., for control volumes that contain and for control volumes that do not contain material interfaces, eventually gives rise to the local form of the balance equation of linear momentum in the bulk and at the interface, namely t + ρ f = 0 (at the interface).
Under the same assumptions that gave rise to (5), the extended integral form of the balance equation of angular momentum follows as with x denoting the position vector of a particle and with x Δ = x − x ref denoting the difference vector to a fixed but otherwise arbitrary reference point x ref . In a small deformation setting, the balance equations are evaluated with regard to the undeformed configuration such that x = x − = x + . Accordingly (11), takes the form By invoking (10), the localisation of (12) further simplifies to the classic symmetry condition of the bulk stress tensor σ = σ t (in the bulk).
In a finite deformation setting, the balance equation of angular momentum additionally stipulates a coaxiality constraint between {{ t}} and the jump in the displacement field across the interface u , see for instance [16] for a detailed derivation and the discussion in [26,38].

Electrical sub-problem
The quasi-stationary processes of electrodynamics that are in the focus of this contribution are governed by the continuity equation for the electric current density and by Faraday's law of induction with e denoting the electric field vector, see e.g. [14,17] for a derivation based on Maxwell's equations of electromagnetism.
Following the same procedure that led to (10), the surface term in (14) is rewritten by making use of the classic divergence theorem The localisation of (16) for control volumes that contain and for control volumes that do not contain material interfaces stipulates the local form of the continuity equation for the electric current in the bulk and at the interface, namely Equation (18b) states that the normal jump of the electric current density vector across a material interface vanishes.
In the derivation of this continuity condition it has tacitly been assumed that the tangential component of the interfacial electric current density vector is negligible. The relaxation of this assumption would lead to a generalised formulation that includes contributions akin to those known from interface elasticity formulations in the mechanical case. In order to localise Faraday's law of induction (15) under the assumption of strong discontinuities in the electric field quantities, it is instructive to resolve the electrical processes at the interface as depicted in Figure 2. To this end, assume for now that it is possible to define the electric field as the (negative) gradient of an electric potential field φ, as is customary for quasi-stationary electrical problems. By integrating the electric field along a path T connecting two opposing points of the interface one arrives at (19) The mathematic analogue to the latter physically motivated derivation is the integration of a δ-distribution since the electric field becomes infinite in the limit case of vanishing interface thickness. With (19) at hand, the evaluation of (15) for an arbitrary control surface A that intersects the material interface, see Figure 3, yields where the Kelvin-Stokes theorem in its classic form has been used. The localisation of (20) for a control surface that does not intersect a material interface yields the local version of Faraday's law of induction in the bulk which can naturally be fulfilled by the introduction of an electric potential for the electric field vector. The localisation of (20) for a control surface that intersects the interface yields no additional condition for the studied cohesive zone formulation since, by invoking (21), the last three summands in (20) cancel out. It is noted that the latter statement tacitly assumes negligible tangential electric currents in the interface and hence neglects the tangential electric field in the interface. In an extended formulation that accounts for the latter quantities in the spirit of interface elasticity, an additional condition would occur that can be fulfilled by a proper definition of the interfacial electric potential. This discussion lies, however, beyond the scope of the present contribution.
Remark 1 (Continuity condition in the case of weak discontinuities). In the derivation of (20), strong discontinuities, i.e. jumps in the electric potential across material interfaces, were considered. Due to this assumption, the classic continuity condition for the tangential component of the electric field vector with h denoting tangential vectors to the interface, was not recovered. If only weak discontinuities, i.e. jumps in the electric field vector, were considered, (20) would take the form which would imply (22) in addition to (21).

Balance equation of energy
This section summarises the balance equation of energy for a continuum that is subjected to mechanical, thermal and electrical loads. The thermal problem is not in the focus of the present contribution but included in the presentation to highlight the fundamental differences between electrical and thermal cohesive zone formulations with regard to the evaluation of the dissipation inequality, cf. Section 2.4. By introducing the mass-specific internal energy densities of the bulk e and of the interface e, the (quasi-static) balance equation of energy reads The powers exerted by mechanical, thermal and electrical loads on the control volume are denoted by P u , P θ and P φ , respectively. In accordance with the assumptions made in Section 2.1, P u takes the form withu denoting the velocity field. By additionally assuming that the interface remains positioned in the middle of the two opposing surfaces, i.e.˙ u = {{u}}, and by invoking (10), (25) simplifies to The thermal power is given by with the classic heat flux vector q and with r denoting massspecific heat sources in the bulk. It is noted that heat sources and the tangential heat flux vector in the interface have been neglected in (27) for the sake of simplicity.
The electrical contribution to (24) is of purely dissipative type if polarisation and magnetisation effects are neglected. In particular, the dissipation caused by the electric current in the bulk is assumed to take the classic form In order to derive the energy contribution related to the electric current across the interface, a simple electrical resistor connecting the two opposing sides of the interface is analysed, as schematically depicted in Figure 2. According to (18b), the electric current density does not change across the interface, i.e. i + = i − = {{ i}}, such that is proposed, see also Remark 2. It is noted that tangential electric currents in the interface are neglected in (29), in accordance with the derivations presented in Section 2.2. By inserting (26)-(29) into (24) and by localising the ensuing equation in the bulk and at the interface one eventually arrives at where the definition of the small strain deformation tensor ε = ∇ sym u was used.
Remark 2 (Alternative derivation of (29)). The specific form of the dissipation related to the electric current across the interface, (29), was motivated in Section 2.3 based on physical arguments. In this remark, a more elaborated way to derive (29) is shown. First when focusing on a control volume that does not contain interfaces, it is observed that by virtue of (18a), (28) can alternatively be expressed in terms of the electric current density and of the electric potential at the boundaries, namely Analogously to the mechanical sub-problem (25), the evaluation of the surface integral-based representation of the electric dissipation for a control volume that contains interfaces yields where tangential electric currents in the interface have been neglected in accordance with Section 2.2. From a physics point of view, (32) states that two contributions are associated with the dissipation caused by the movement of electric charges from regions of high electric potential to regions of low electric potential. Whereas the first contribution is related to the classic electrical resistance of the bulk, the second contribution accounts for the electrical resistance of the interface. By additionally taking (3) and (18b) into account, the interface contribution in (32) can further be rewritten as

Dissipation inequality
With θ and s denoting the absolute temperature and the massspecific entropy density of the bulk, and with s denoting the mass-specific entropy density of the interface, the dissipation inequality reads In order to localise (34), the volume integral on the right-hand side is rewritten in a first step. More specifically speaking holds, where (3) was used. To proceed, the convex-concave Legendre(-Fenchel) transformations with energetic duals are invoked. In the latter set of equations, ψ and ψ denote the mass-specific Helmholtz free energy density functions of the bulk and of the interface, and θ denotes the interface temperature. By inserting (30) and (35) into (34), and by making use of (36) one eventually arrives at the local forms of the dissipation inequality The application of a Coleman-Noll-type procedure to evaluate (38) naturally gives rise to the constitutive restrictions for the electrical sub-problem j · e ≥ 0 (in the bulk), From a physics point of view, (39) implies that the flow of electric charges is opposed to the electric potential gradient. In contrast to the electrical sub-problem, the thermal subproblem at the interface cannot further be evaluated without additional assumptions on the relation between the bulk and the interface temperature.

Weak form of equilibrium
In order to derive the weak form of equilibrium for the electro-mechanically coupled continuum under consideration, field equations (10a) and (18a) are multiplied by test functions η u , respectively η φ (41) and integrated over (the subdomains of) body B. Likewise, (10b) and (18b) are multiplied by test functions and integrated over the interface(s) I. By additionally applying the classic divergence theorem and by making use of (3), the weak form of the balance equation of linear momentum (10), i.e.
takes the form Following the same procedure, the weak form of the continuity equation for the electric current density can be rewritten as

Finite element implementation
The finite element implementation of the proposed cohesive zone formulation is based on the discretisation of primary fields and test functions by means of Lagrange polynomials, namely, with n en denoting the number of element nodes, with N • denoting shape functions and with nodal values being indicated by subscripts, e.g. u C . By inserting (47) into the weak forms of the field equations, (44) and (46), one arrives at the residual-type format that is based on the (generalised) internal, surface and volume forces vectors defined in Appendix A. The linearisation of (48) at some iteration step q in an iterative, gradient-based solution method moreover results in the update relation for the global lists of nodal degrees of freedom, u and φ. In (49), K • * and K • * denote bulk and interface contributions to the generalised global stiffness matrix. The specific forms of these quantities are provided in Appendix B.

Representative simulation results
This section focuses on the study of representative boundary value problems. To this end, the electro-mechanical material behaviour is specified in Section 4.1 and Section 4.2.
In particular, a coupling of the mechanical and electrical field equations due to the influence of deformation-induced interface damage processes on the electrical conductivity is considered. With the specific material model at hand, analytical solutions for a one-dimensional tensile problem are derived in Section 4.3 and compared with finite elementbased simulation results in Section 4.4. Finally, tension and compression tests of a test specimen with geometrically resolved microstructure are studied in Section 4.5. It is emphasised that the boundary value problems studied in the proof of concept-type simulations are of academic nature and that, in particular, neither the material models nor the material parameters are adapted to a specific material (system). Moreover, it is noted that the derivations presented in Section 2 and Section 3 hold independently of the spatial dimensions such that the application of the proposed formulation to general three-dimensional boundary value problems is possible.

Mechanical material models
For the sake of simplicity, the mechanical response of the bulk is assumed to be isotropic, linear-elastic and to be governed by the Helmholtz free energy density function with the elasticity tensor being parametrised in terms of Young's modulus E and Poisson's ratio ν. For the interface, a linear-elastic mechanical material response combined with a brittle damage formulation is used, which represents a direct extension of the model proposed in [29]. More specifically speaking, the area specific interface Helmholtz free energy density function Ψ = ρ ψ is formulated as a function of the displacement jump u and of the interface damage variable d. Furthermore, it is assumed that interface damage processes manifest themselves in a reduction of the mechanical stiffness E in the tensile region whereas the interface is assumed to regain its normal stiffness under compressive loadings With (50) and (52) at hand, the evaluation of (38) yields the specific form of the bulk stress tensor and of the traction separation law under tensile loadings respectively under compressive loadings By virtue of (38)-(40) and (52a), the reduced form of the interface dissipation inequality in the tensile region takes the form which stipulates that the damage process is intrinsically energy driven and that no self-healing occurs. By additionally assuming that damage evolution only takes place under tensile loadings (i.e.˙ d = 0 if u · n < 0), the failure function with admissibility conditions is proposed. In accordance with the developments presented in [12,29], d is moreover assumed to be related to the displacement jump-type internal variable at time t * , In the light of (57), χ 0 takes the interpretation of the critical interface opening that is related to the initial fracture strength Q 0 according to Furthermore, G F can be interpreted as the fracture toughness, see [29] and Remark 3.

Remark 3 (Fracture toughness)
. The fracture toughness of the brittle material under consideration is defined by the dissipation that is associated with the crack formation. In particular, time integration of (56) yields

Electrical material models
The electrical material models are chosen in accordance with the restrictions posed by the dissipation inequalities (39). In particular, a linear relation between the electric field vector and the electric current density vector, namely is assumed in the bulk, with denoting the (positive definite) electrical conductivity tensor that is defined in terms of the scalar-valued electrical conductivity parameter κ. Moreover, the deformation-induced damage processes at the interface are assumed to influence the effective electrical conductivity. The specific form of the constitutive equation for the electric current density with κ denoting the idealised conductivity of the interface, thus establishes a coupling between the electrical and mechanical field equations. Equation (65) is based on the assumption that the effective conducting interface area in the tensile region is reduced due to the damage processes, whereas the interface recovers its initial conductivity under compressive loads. Conceptually speaking, the electrical problem is thus treated similarly to the mechanical problem with the damage evolution being associated with a reduction in the number of springs, respectively a reduction in the number of wires, connecting the opposing sides of the interface, cf. Figure 2.

Analytical solution
This section deals with the derivation of an analytical solution for the one-dimensional sample boundary value problem sketched in Figure 4, subject to the assumption of monotonic tensile loadings. The total elongation of the bar results from the elongation of the regions B − and B + , i.e. 2 Δu, and from the displacement jump across the interface u . Moreover, the latter kinematic quantities are related to the axial stress in the bulk via and determine the stresses at the interface For monotonic loadings, d (χ ) can equivalently been expressed as a function of the displacement jump u , such that (67) establishes an analytical relation between the displacement jump across the interface and the total elongation of the bar. In this regard, it is noted that, by virtue of (59) and (60), Expressions for the effective mechanical and electrical properties of the bar may further be derived by interpreting the bar in Figure 4 as a series connection of springs and resistors. More specifically speaking, the axial stress σ as a function of the displacement jump takes the form and the axial electric current density j as a function of the displacement jump and of the prescribed electric potential difference Δφ reads With parametrisations of the total elongation, of the axial stress and of the electric current density in terms of the displacement jump at hand, (implicit) relations between the total elongation and the axial stress, respectively between the total elongation and the electric current density, are established.
Remark 4 (Size effects). According to (69), the electrical conductance of the bar for a given deformation state reads Likewise, by taking (69) into account and by neglecting snapback-type effects the mechanical (unloading) tensile stiffness The specific forms of (70) and (71) indicate a size-dependent theory and motivate the introduction of internal length scales which can be used to derive a unified form of (70) and (71). More specifically speaking, by indicating the reference (unloading) tensile stiffness and the reference conductance of the virgin material by subscripts • 0 one arrives at Relation (73) suggests a similar evolution of the mechanical (unloading) stiffness and of the electrical conductance with damage evolution, and is plotted in Figure 5 for various values of the length scale ratio • /l. In particular, it is observed that the extent to which interface damage processes manifest themselves in a reduction of the mechanical stiffness and of the electrical conductance significantly depends on the length scale ratio • /l which may take different values for the mechanical and the electrical problem.

Comparisson of simulation results and analytical solutions
This section focuses on the validation of the finite element implementation proposed in Section 3 based on the analytical solution derived in Section 4.3. To this end, the elementary boundary value problem depicted in Figure 4 is discretised by means of four-node quadrilateral bulk and four-node (two nodes on each side) linear interface elements, e.g. [32,33], with the same discretisation being used for the geometry and the primary electrical and mechanical field quantities. Occurring integrals are evaluated numerically by means of Gaussian quadrature schemes featuring four quadrature points, respectively one quadrature point.
Since the two-dimensional finite element calculations are carried out subject to the simplifying assumption of a plane strain deformation setting, the material parameter E in the analytical solutions derived in Section 4.3 needs to be substituted by its plane strain equivalent E ps , i.e.
The finite element results and the analytical solutions for the axial stress σ and for the axial electric current density j are provided in Figure 6 for a bar of length 2 l = 100mm. The finite element simulations are based on the material parameters provided in Table 1, whereas, according to (74), E ps ≈ 230769Nmm −2 and E = 210000Nmm −2 are used in the evaluation of the analytical solutions.

Intergranular crack propagation
Motivated by cracks that propagate along grain boundaries, the two-dimensional sample boundary value problem depicted in Figure 7 is studied in this section. In particular, the individual grains of the microstructure are resolved in the ensuing finite element calculations by making use of the proposed electro-mechanically coupled cohesive zone formulation. By doing so, transgranular fracture processes where cracks grow through the material grains are neglected for now and the set of admissible crack paths is significantly reduced. Furthermore, the material models are chosen according to Section 4.1 and Section 4.2, a plane strain setting is assumed and the set of material parameters provided in Table 1 is adopted. In accordance with Section 4.4, the finite element discretisation is based on four-node quadrilateral bulk and four-node (two nodes on each side) linear interface elements, and Gaussian quadrature schemes with four quadrature points, respectively one quadrature point, are used. Whereas the vertical displacement at the two ends of the test sample is fixed in the simulations, a relative displacement   Table 1 and assume a plane strain deformation setting. By virtue of (74), the analytical solution is evaluated for E ps ≈ 230769Nmm −2 and E = 210000Nmm −2  in horizontal direction is prescribed that results in a total elongation Δu. Moreover, an electric potential difference Δφ between the two ends is prescribed and the electric potential at the left boundary is set to 0.0mV. The potential difference Δφ = 0.1mV is applied in the first load step (A→B) and kept constant during tensile loading (B→C), unloading (C→D) and compressive loading (D→E). For the latter load path, the horizontal reaction force F · e 1 and the electric current I are depicted in Figure 8 as functions of deformation. Due to the evolution of interface damage, significant reductions in both the mechanical reaction force and the electric current are observed under tensile loading (B→C). In the (purely elastic) unloading step (C→D), no damage evolution occurs so that the overall mechanical stiffness and electrical conductance remain constant. Under compressive loading (D→E), the damage evolution is suppressed and the interfaces are assumed to regain their normal stiffness and conductivity. Hence, jumps in the horizontal reaction force-elongation curve and in the electric current-elongation curve are observed at the end of the unloading step (C→D). In addition to the generalised reaction force-displacement curves, the electric potential field and the electric current density field are depicted in Figs. 9 and 10 for various load steps. Since the material interfaces pose additional obstacles for the electric current, a complex, inhomogeneous electric current density distribution is observed. Moreover, due to the evolution of interface damage and its influence on the electrical conductivity, the electric current density distribution and the electric potential field significantly change with deformation. Electric potential field φ at the centre region of the boundary value problem depicted in Figure 7 for deformation states according to Figure 8. In particular, state B resembles a sample in an undeformed state, state C a sample subjected to tensile loading and state E a sample subjected to compressive loading. It is noted that the material interface in state E is damaged due to the deformation history. Moreover, (normalised) electric current density vectors are indicated by black arrows Finally, an additional study of the influence of the finite element discretisation on the simulation results is provided in Appendix C.

Closure
Motivated by the influence of mechanically-induced interface damage processes on the electrical conductivity, this contribution dealt with the fundamentals of electro-mechanically coupled cohesive zone formulations for electrical conductors. At the outset of the theory, local forms and jump conditions of the governing set of balance equations for the electrical sub-problem, that reduced to Faraday's law of induction and to the continuity equation for the electric current, were discussed. In addition, an extended form of the balance equation of energy that accounted for additional electrical interface contributions was proposed and the associated restrictions stipulated by the dissipation inequality were studied. Next, the electro-mechanically coupled cohesive zone formulation was implemented into a multifield finite element framework. The coupling between the  Table 1 and assume a plane strain deformation setting electrical and mechanical sub-problem was established by the constitutive equations at the material interface -i.e. mechanically-induced damage processes were assumed to reduce the interface conductivity in the sample boundary value problems studied. More specifically speaking, the finite element formulation was validated by means of analytical solutions in a first step, before the influence of intergranular crack propagation on the effective electrical properties was studied.
The sample boundary value problems have been of academic nature to show the principle applicability of the proposed formulation. However, to truly understand and predict the electro-mechanical interaction at material interfaces, tailored material models (including damage, plasticity and their influence on the electrical conductivity) need to be furnished for a specific material system in close relation to experiments. In this regard, sophisticated computational multiscale approaches for material interfaces, as elaborated for instance in [13,20,21], may be invoked to take detailed (lower scale) information of the material interfaces into account. Put into perspective, the general modelling framework for electro-mechanically coupled cohesive zone formulations Comparison of electric potential field φ and electric current density field j at the centre region of the boundary value problem depicted in Figure 7 for different finite element discretisations and defor-mation states according to Figure 8, respectively Figure 11. In particular, state B resembles a sample in an undeformed state and state C a sample subjected to tensile loading that has been established in the present work provides a basis for these developments.
In conclusion, it is remarked that tangential electric currents in the interface have been neglected in the proposed formulation. From a modelling point of view, these processes share similarities with the theory of interface elasticity known from mechanical problems. An extension of the present formulation in the spirit of generalised imperfect interfaces is thus aspired in future works. in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Generalised nodal force vectors
The bulk and interface contributions to the (generalised) internal, surface and volume forces vectors that result from the discretisation of (44) and (46) are briefly summarised in this appendix. In the ensuing equations, assembly operator A will frequently be used to establish the relation between local quantities and their global analogues. In this regard, n el and n el denote the total number of bulk and interface elements, respectively.
For the mechanical subproblem, the bulk and interface contribution to the internal force vector are given by

B Tangent stiffness contributions
This appendix focuses on the specific form of the tangent stiffness contributions that occur in the linearisation of (48). For the sake of simplicity, functional dependencies of the volume distributed forces and of the (generalised) surface forces on the primary field variables are neglected in the following. The contributions to the global tangent stiffness matrix that results from the linearisation of (75) and (76) with respect to the displacement field are given by The respective sensitivities with respect to the electric potential field read Likewise, the linearisation of (80) and (81) with respect to the displacement field give rise to whereas the linearisation with respect to the electric potential field results in

C Influence of the finite element discretisation
This appendix refers to the influence of the finite element discretisation on the solution of the boundary value problem depicted in Figure 7. For the material models, material parameters and load path introduced in Section 4.5 and for two different finite element discretisations, the resulting (generalised) reaction force-elongation curves are depicted in Figure 11. In addition, the electric potential field and the electric current density field are provided in Figure 12. The comparison of the simulation results eventually reveals a comparably small influence of the finite element discretisation on the simulation results.