## Abstract

An original boundary integral formulation is proposed for the problem of a semi-infinite crack at the interface between two dissimilar elastic materials in the presence of heat flows and mass diffusion. Symmetric and skew-symmetric weight function matrices are used together with a generalized Betti's reciprocity theorem in order to derive a system of integral equations that relate the applied loading, the temperature and mass concentration fields, the heat and mass fluxes on the fracture surfaces and the resulting crack opening. The obtained integral identities can have many relevant applications, such as for the modelling of crack and damage processes at the interface between different components in electrochemical energy devices characterized by multi-layered structures (solid oxide fuel cells and lithium ions batteries).

## 1. Introduction

Modelling of interfacial crack problems in elastic bimaterials in presence of thermodiffusion represents an important issue for many engineering applications. In particular, it is crucial for studying fracture initiation and propagation at the interface between different components of electrochemical energy devices which are characterized by multi-layered structures, such as solid oxide fuel cells and lithium ions batteries. Indeed, because of the high operational temperature that can be reached and to the intense particle fluxes that are required for maintaining the electrical current, the components of these devices are subject to severe thermomechanical stresses as well as stresses induced by the particle diffusion [1,2]. This can cause damage and crack formation compromising their performances in terms of power generation and energy conversion efficiency [3–5]. For this reason, the modelling of fracture processes at the interface between the components of such battery devices is fundamental to the prediction of these phenomena, and subsequently enables successful manufacture and reliable performances of the systems.

Owing to their distinctive feature of reducing by one the dimension of the considered problems, so that only the boundary or surface of the domain needs to be modelled, boundary integral formulations are particularly suitable for studying multiphysics phenomena such as dynamic and static fracture processes in thermodiffusive elastic materials. Several boundary integral approaches have been proposed for solving crack problems in linear elastic, thermoelastic and thermodiffusive materials, such as analytical techniques based on the method of singular integral equations [6–8], and numerical techniques based on the boundary element method [9–13]. In all these approaches, the displacements and stress fields are defined by integral relations involving the Green's functions, which need to be derived analytically in explicit form [14] or computed numerically [15]. Although Green's functions have been derived for several crack problems in linear thermoelastic and thermodiffusive elastic materials [16–19], their utilization for calculating physical displacements and stress fields on the crack faces requires challenging numerical estimation of integrals whose convergence should be assessed carefully. Moreover, the approach based on the Green's function method works when the tractions and the thermal and diffusive stresses acting on the discontinuity surface are symmetric, but not in the case where asymmetric mechanical and thermodiffusive loading distributions are applied on the crack faces.

In this paper, the problem of a semi-infinite quasi-static crack at the interface between two dissimilar thermodiffusive elastic materials is addressed by means of an original boundary integral formulation which avoids the use of the Green's functions and the challenging computations connected, without any assumptions regarding the symmetry of the loading and of the temperature and mass concentration profile at the interface. The general approach recently proposed in Piccolroaz & Mishuris [20], Morini *et al.* [21], Vellender *et al.*[22] and Mishuris *et al.* [23] for interfacial crack problems in isotropic and anisotropic elastic bimaterials, based on Betti's reciprocal theorem and weight functions theory, is extended in order to study fracture processes in presence of thermodiffusion. The volume integral terms present in the reciprocity identity, associated with the temperature and mass concentration effects [24], are converted into surface integrals through an exact transformation based on the notion of Lamé elastic potentials [25] while assuming that the temperature and mass concentration fields are harmonic in the domain. The derived original form of Betti's identity is used together with symmetric and skew-symmetric weight function matrices derived by [20] for formulating the considered crack problem in thermodiffusive bimaterials in terms of singular integral equations.

The article is organized as follows: in §2, the static governing equations for a linear elastic thermodiffusive media are formulated. Reciprocity identities are introduced and the volume integral terms associated with the temperature and mass concentration fields are transformed into boundary surface integrals by introducing Lamé elastic potentials and applying second Green's theorem. In §3, the problem of a quasi-static crack at the interface between two dissimilar elastic thermodiffusive media is introduced. The weight functions, defined as a special singular solution of the homogeneous traction-free problem are used together with the obtained Betti's identity for formulating the problem in terms of boundary integral equations. The case of a plane strain crack is analysed in §4. Lamé potentials in the weight functions space are derived in closed form, and explicit integral identities relating the applied mechanical loading, the profiles of temperature, mass concentration, heat and mass fluxes on the interface and the resulting crack opening are obtained. Finally, in §5, the integral identities are used to study some illustrative examples of plane crack problems in the presence of thermodiffusion. Exact expressions for the crack opening and tractions ahead of the crack tip associated with the introduced temperature and heat flux distributions are derived, and the corresponding stress intensity factors are calculated in closed form.

## 2. Preliminary results: governing equations and reciprocity theorem

In this section, constitutive relations and static balance equation for linear thermodiffusive elastic media in infinitesimal deformations are introduced. Betti's integral identities are derived by means of the theorem of the reciprocity of work, and volume integral terms associated with heat flows and mass diffusion are converted into surface integrals using an exact transformation based on the introduction of Lamé elastic potentials.

### (a) Governing equations

For a linear isotropic elastic body where temperature changes and mass diffusion are considered, the constitutive relationship between the stress ** σ** and the strain

**is given by [24]:**

*ε**μ*are Lamé's constants,

*θ*=

*T*−

*T*

_{0}is the temperature of the medium with respect to the reference state (

**=0,**

*ε**T*=

*T*

_{0}

*C*=

*C*

_{0}),

*χ*=

*C*−

*C*

_{0}is the concentration of diffusing particles with respect to the natural state,

**is the identity matrix,**

*I**γ*

_{t}=(3λ+2

*μ*)

*α*

_{t}and

*γ*

_{c}=(3λ+2

*μ*)

*α*

_{c}with

*α*

_{t}and

*α*

_{c}coefficients of thermal and diffusive expansions, respectively. Then in the natural state of the system we have

**=0,**

*ε**θ*=0,

*χ*=0 and in order to describe the material using linear thermoelastic diffusion theory [24,26], the values of

*θ*and

*χ*are assumed such that |

*θ*/

*T*

_{0}|≪1 and |

*χ*/

*C*

_{0}|≪1. In the static case, the equations of equilibrium are given by

**represents the body forces. Substituting the constitutive relation (2.1) in (2.2), the Navier's equations in terms of displacements**

*b***are obtained**

*u*The conservation of energy and mass in steady-state conditions leads to
*W* is the heat generated for unit of time and volume by internal sources. For the considered linear isotropic thermodiffusive media, the heat flux ** q** and particles current

**can be expressed as**

*j**k*

_{t}is the thermal conductivity coefficient and

*D*

_{c}is the mass diffusivity of the material. Using relations (2.5) in (2.4), the following Poisson's and Laplace's equations for the steady-state temperature and mass concentration fields are derived:

### (b) Reciprocity theorem: transformation of the volume integrals

Let us consider a thermodiffusive elastic body of volume *V* subject to the action of body forces *b*^{(1)}, tractions *t*^{(1)}, heat sources *W*^{(1)}, surface heating to the temperature *u*^{(1)}(** x**), temperature

*θ*

^{(1)}(

**) and mass concentration**

*x**χ*

^{(1)}(

**). The set of causes and effects characterizing the state of the medium can be written symbolically in the compact form**

*x*

*σ*^{(1)}and the strains

*ε*^{(1)}are connected by the relationship (2.1), where the dilatation is defined as tr

**=∇⋅**

*ε***, and they are assumed to be continuous together with their first derivatives. The displacements, temperature and mass concentration are also continuous and moreover have continuous derivatives up to the second order for**

*u***∈**

*x**V*+∂

*V*. These functions satisfy the fields equations (2.2), (2.6)

_{(1)}and (2.6)

_{(2)}with boundary conditions:

**denotes the unit outward normal to the surface ∂**

*n**V*.

Introducing another set of causes *u*^{(1)},*σ*^{(1)} and *θ*^{(1)}, respectively, are derived as functions of the test quantities *u*^{(2)},*σ*^{(2)} and *θ*^{(2)}, which in general are assumed to be fundamental solutions of the evolution equations or weight functions [9,10]. Several analytical and numerical transformations have been proposed for reducing the volume integrals associated with body forces and coupling between mechanical strains and temperature to surface integrals [29,30].

Here, we assume zero body forces and heat sources acting on the system, and we introduce an exact procedure for transforming coupling volume terms involving temperature and mass concentration in equation (2.10). From Helmholtz's decomposition theorem, as both the displacements *u*^{(1)} and *u*^{(2)} are solutions of the equilibrium equation (2.3), they can be represented by *scalar displacement potentials* *φ*^{(1)} and *φ*^{(2)} and *vector displacement potentials* *ψ*^{(1)} and *ψ*^{(2)} [25]:

Applying the Green's second identity to the volume integral terms in (2.14) [30], using equations (2.6) and remembering that zero heating sources are assumed, we get

By means of the proposed general procedure, the volume integral terms associated with thermal and diffusive stresses in the reciprocity identity (2.14) has been reduced to surface integrals. Remembering the fluxes definitions (2.5), expression (2.19) can be written as follows:
*β*_{t}=*γ*_{t}/*k*_{t} and *β*_{c}=*γ*_{t}/*D*_{c}. Assuming *W*=0, reciprocity identities (2.11) and (2.12) can also be expressed in terms of fluxes:

In the cases where it is possible to determine explicitly the elastic potentials (2.13), using expression (2.20) together with (2.21), boundary integral formulation of static crack problems in thermodiffusive solids can be obtained avoiding numerical estimation of volume integrals [31]. In the next sections, the obtained relation (2.20) will be extensively applied in order to derive integral identities for the modelling of fracture phenomena at the interface between dissimilar elastic materials in the presence of thermodiffusion.

## 3. Interfacial cracks in thermodiffusive media

A static semi-infinite crack between two dissimilar elastic materials in the presence of thermodiffusion is considered, the geometry of the system is shown in figure 1. No body forces, heat and mass sources are assumed. The crack is situated in the half-plane *f*〉, and the jump, [[*f*]], of a function *f* across the plane containing the crack, *x*_{2}=0.

In the proposed model, thermal conduction and diffusion are both assumed to be isotropic, and then only normal stresses are induced by temperature changes and mass diffusion (see constitutive relation (2.1)).

Further in the paper, using the Betti identities (2.20) and (2.21), the considered crack problem will be formulated in terms of boundary integral equations. In order to correctly apply expressions (2.20) and (2.21), a preliminary discussion regarding the behaviour of the temperature and mass concentration functions in the upper and lower elastic thermodiffusive half-spaces is reported in next §3a, and boundary conditions for these fields and the fluxes (2.5) are introduced.

### (a) Temperature and mass concentration in thermodiffusive half-spaces

As zero heat and mass sources are assumed, the temperature *θ* and mass concentration fields are harmonic functions determined by the solution of a Laplace's equation assuming the form of (2.6)_{(2)}. Considering the upper elastic thermodiffusive half-space these equations become

The values of *θ* and *χ* for *x*_{2}=0^{+}, corresponding to the plane containing both the bounded interface and the crack, are given by
*θ*^{+} with *θ*^{−} and *χ*^{+} with *χ*^{−} a couple of Laplace's equations identical to the (3.2) and defined for *x*_{2}<0 are obtained for the lower thermodiffusive half-space, and then for *x*_{2}=0^{−} boundary conditions analogous to (3.3) can be written.

The average and the jump of the temperature and concentration are defined on the whole plane *x*_{2}=0:

Further in the text, we will refer to (3.4) and (3.5) as to imperfect thermodiffusive interface conditions. The values of the components of the fluxes normal to the interface plane, *q*_{2}=** q**⋅

*e*_{2}and

*j*

_{2}=

**⋅**

*j*

*e*_{2}, for

*x*

_{2}=0

^{+}are

*x*

_{2}=0

^{−}.

The average and the jump of *q*_{2} and *j*_{2} are also defined on the whole plane *x*_{2}=0:

Both in the upper and lower half-spaces, the temperatures and mass concentration fields are assumed to vanish at infinity:
*q*_{2} and *j*_{2} must satisfy the self-balance conditions on the boundary. For the considered upper thermodiffusive half-plane this means that
*x*=0 have been utilized. Note that the self-balance conditions for the fluxes on *x*_{2}=0^{−} are given by relations similar to the (3.10). Remembering the definitions (3.7) and (3.8), integral balance conditions analogues to the (3.10) can be derived for the average and jump of the fluxes across the plane *x*_{2}=0.

### (b) Boundary integral equations

Considering the geometry of the model shown in figure 1, the Betti identities (2.20) and (2.21) are applied to a semi-spherical domain of radius *r* in the upper and in the lower half spaces *u*^{(1)} and *u*^{(2)} decay suitably fast at infinity and remembering conditions (3.9) for the temperature and mass concentration fields, the reciprocity relation (2.20) in the upper half-space becomes
** u**=[

*u*

_{1},

*u*

_{2},

*u*

_{3}]

^{T}and

*σ*_{2}denotes the traction vector acting on the plane

*x*

_{2}=0:

*σ*_{2}=

*σ*

*e*_{2}. In the same limit

Considering the boundary *x*=0^{−} instead of *x*=0^{+}, three similar expressions are derived for the lower half-space. In the cases where elastic potentials can be derived explicitly, integral equations (3.11)–(3.13) can be used together with their analogues in the lower half-plane in order to obtain expressions for the physical fields associated with the crack problem

### (c) Betti formula and weight functions

Following the approach proposed in several boundary element formulations of thermoelasticity [30,31], for the auxiliary solutions system we assume *u*^{(2)} is represented by the non-trivial singular solution of the homogeneous traction-free crack problem, commonly known as the weight function [33,34], and defined as follows:
** R** is the rotation matrix:

The transformation (3.14) corresponds to a change of coordinates consisting in a rotation of an angle *π* around the *x*_{2}-axis, and it is connected to the fact that the weight function ** U** is defined in a different domain with respect to physical displacement, where the crack is placed along the positive semi-plane

*x*

_{1}=0,

*x*

_{2}>0 [35]. The stress tension

**associated with the displacements**

*Σ***is introduced:**

*U***(**

*U**x*

_{1},

*x*

_{2},

*x*

_{3}) with

**(**

*U**x*′

_{1}−

*x*

_{1},

*x*

_{2},

*x*′

_{3}−

*x*

_{3}) which corresponds to a shift within the plane (

*x*

_{1},

*x*

_{3}), for the upper half-plane we obtain:

^{+}is replaced with 0

^{−}. Subtracting this expression from (3.17), we obtain

*Φ*

_{γt},

*Φ*

_{γc},

*Φ*

_{βt}and

*Φ*

_{βc}indicates the following normalized potentials:

**,**

*u*

*σ*_{2},

*θ*,

*q*

_{2},

*χ*and

*j*

_{2}involved in the (3.18) can be represented as follows:

*H*denotes the Heaviside function. The reciprocity identity (3.18) then becomes

*x*

_{1}and

*x*

_{3}, which is defined as follows [36]:

On the basis of representation (3.20), 〈*σ*_{2}〉^{(+)} is the traction along the interface, ahead of the crack tip, and [[** u**]]

^{(−)}is the crack opening (displacement discontinuity across the crack faces). [[

**]] and 〈**

*U***〉 are the symmetrical and skew- symmetrical weight functions matrices defined and derived in closed form in [32], whereas the term 〈**

*U*

*Σ*_{2}〉 stands for the traction along the

*x*

_{1}-axis corresponding to the singular auxiliary displacements

**.**

*U*The integral equation (3.22) is the generalization to thermoelastic diffusive media of the Betti identity derived in [20,21], and it relates the physical solution ** u**,

*σ*_{2},

*θ*,

*χ*to the auxiliary singular solution

**,**

*U*

*Σ*_{2}. Note that (3.22) is valid for the most general case of static interfacial crack problems in presence of thermal and diffusive effects, and includes also the case of imperfect mechanical interface conditions, associated with the discontinuity of tractions and displacements at the interface ahead of the crack tip [22]. In the sequel, perfect contact conditions at the interface for

*x*

_{1}>0 will be assumed for mechanical fields, then: [[

**]]**

*u*^{(+)}=[[

*σ*_{2}]]

^{(+)}=0, whereas according to definitions reported in §3a discontinuity of the thermodiffusive quantities along the whole plane

*x*

_{2}=0 will be considered (imperfect thermodiffusive interface conditions). Moreover, the notation (3.1) is used to indicate the symmetrical and skew-symmetrical mechanical loads applied at the crack faces, then 〈

*σ*_{2}〉

^{(−)}=〈

**〉 and [[**

*p*

*σ*_{2}]]

^{(−)}=[[

**]], respectively. The integral identity (3.22) becomes**

*p**θ*〉,〈

*χ*〉,〈

*q*

_{2}〉 and 〈

*j*

_{2}〉 stand, respectively, for the mean values of the temperature, mass concentration, heat flux and particle current on the plane

*x*

_{2}=0, defined, respectively, by expressions (3.4) and (3.7). Similarly, [[

*θ*]],[[

*χ*]],[[

*p*

_{2}]] and [[

*j*

_{2}]] denotes the jumps of these functions across the same plane

*x*

_{2}=0, and defined by (3.5) and (3.8).

The Betti's identity (3.24) will be used further in the text in order to formulate the illustrated crack problem at the interface between dissimilar thermodiffusive materials in terms of singular integral equations. Explicit integral identities relating the loading applied at the crack faces, the temperature, the mass concentration, the heat and mass fluxes to the resulting crack opening and traction ahead of the crack tip are derived for the two-dimensional case.

## 4. Integral identities

In this section, starting from the general identity (3.24), the problem of a two-dimensional interfacial crack in presence of thermodiffusion is reduced to a system of explicit integral equations. The solution of these equations provides exact expressions for the crack opening and for the tractions ahead of the tip corresponding to an arbitrary loading configuration and to the temperature and concentration profiles of the system, obtained by solving the Poisson's equations (2.6). As thermal conduction and diffusion are both assumed to be isotropic, only normal stresses are induced by temperature changes and mass diffusion (see constitutive relation (2.1)). For this reason, in the considered two-dimensional problem anti-plane stress and deformations are not affected by thermodiffusion, so that only the plane strain case is studied in order to estimate the contribution of thermal and diffusive stresses on interface fracture phenomena.

For plane strain deformations, the Betti identity (3.24) relating the physical solution ** u**=[

*u*

_{1},

*u*

_{2}]

^{T}and

*σ*_{2}=[

*σ*

_{21},

*σ*

_{22}]

^{T}with the weight function

**and**

*U*

*Σ*_{2}becomes

**〉=[〈**

*p**p*

_{1}〉,〈

*p*

_{2}〉]

^{T}, [[

**]]=[[[**

*p**p*

_{1}]],[[

*p*

_{2}]]]

^{T}and the symbol * denotes the convolution with respect to the variable

*x*

_{1}. Here and in the sequel of the article we will used the following 2×2 matrices:

For plane elastic bimaterials, two linearly independent weight functions,

In order to express the weight function matrix (4.3)_{(1)} in terms of elastic Lamé potentials (4.17), it is important to note that in the considered plane strain case the vector potential ** Ψ** possesses only one non-zero component directed along

*x*

_{3}-axis, and then:

**=**

*Ψ**Ψ*

_{3}

*e*_{3}=

*Ψ*

*e*_{3}. Two couples of elastic potentials (

*Φ*

^{1},

*Ψ*

^{1}) and (

*Φ*

^{2},

*Ψ*

^{2}) are introduced, and the components of

**are expressed as follows:**

*U*Let us introduce the Fourier transform of a generic function *f* with respect to the variable *x*_{1}:
*Φ*_{βc}=[*β*_{c}*Φ*^{1}, *β*_{c}*Φ*^{2}]^{T}, whereas the superscripts + and − denote ‘+’ and ‘−’ functions, respectively. In order to obtain from expression (4.6) a system of integral equations relating the crack opening and the traction ahead of the tip with applied loading and thermodiffusive quantities, explicit expressions for the elastic potentials *Φ*^{1} and *Φ*^{2} are needed.

### (a) Elastic potentials

Applying the Fourier transform to equations (4.4), the following system of non-homogeneous ODEs is derived:
*x*_{2}. As both the pairs of transformed elastic potentials

The Fourier transform of the singular displacements *et al.* [32]. Substituting the expressions obtained by Piccolroaz *et al.* [32] into (4.7), the solution of the linear system yields expressions for the transformed elastic potentials in both the upper and the lower half-plane. The explicit derivation of these expressions for elastic potentials is reported in detail in the electronic supplementary material. The traces of the elastic potentials on the plane *x*_{2}=0 containing the crack are used together with the definition (3.19) for deriving the jump and the average of the Fourier transform of the normalized potentials *Φ*_{βt}:
*et al.* [37]) and
*h*_{t},ℓ_{t},*h*_{c} and ℓ_{c} are reported in the supplementary material.

In order to derive explicit integral identities by equation (4.6), also expressions for the jump and the average across of the plane containing the crack of the derivatives of the normalized elastic potentials with respect to *x*_{2} are required. As *x*_{2}, *x*_{2} can be finally written as
*m*_{t},*n*_{t},*p*_{t},*q*_{t},*m*_{c},*n*_{c},*p*_{c} and *q*_{c} are defined in the supplementary material.

The obtained expressions for the jump and the average of the Fourier transform of elastic potentials and fluxes across of the plane *x*_{2}=0 are now used together with the Betti identity (4.6) for reducing the interface crack problem to a system of explicit singular integral equations.

### (b) Explicit integral identities

Substituting expressions (4.8) and (4.11) in equation (4.6), we obtain:
** A** and

**are the following matrices, identical to those introduced in Piccolroaz & Mishuris [20] and Morini**

*B**et al.*[21] for the elastic problem without thermodiffusion:

*g*_{t}〉,[[

*g*_{t}]],〈

*d*_{t}〉 and [[

*d*_{t}]], associate to thermal effects, are given:

*g*_{c}〉,[[

*g*_{c}]],〈

*d*_{c}〉 and [[

*d*_{c}]], corresponding to the contribution of the mass diffusion, are defined as follows:

**is given by**

*M*The matrices ** A**,

**and**

*B***can be easily computed using the symmetric and skew-symmetric weight functions derived by [32]:**

*M**α*and

*d*are the so-called Dundurs parameters,

*b*and

*γ*are bimaterial constants defined in the Appendix A in supplementary material. Using expressions (4.18) together with the definitions (4.14) and (4.17), we get

Substituting (4.10) and (4.17) in equation (4.15), the explicit expressions for the vectors [[*d*_{t}]],〈*d*_{t}〉,[[*g*_{t}]] and 〈*g*_{t}〉 and for 〈*d*_{c}〉, 〈*d*_{c}〉, 〈*g*_{c}〉 and 〈*g*_{c}〉 are obtained:

As explicit expressions for all the terms of the identity (4.13) have been derived, we can now apply the inverse Fourier transform to this equation, obtaining two distinct relationships corresponding to the two cases *x*_{1}<0 and *x*_{1}>0:
*f* can be both the jump and the average of *θ*,*χ*,*q*_{2} and *j*_{2}, and *γ*_{eul} is the Euler's gamma constant. Using the inverse transforms (4.25) and (4.26) in equations (4.23) and (4.24), we finally obtain the explicit integral identities for plane cracks problems at the interface between dissimilar thermodiffusive elastic materials:
*et al.* [21] by introducing the singular operator of the Cauchy type

The vector operators *x*_{2}=0, depend by

Equations (4.27) and (4.28) form a system of integral equations relating the crack opening and the traction ahead of the tip to the mechanical loading applied at the crack faces and to the values of temperature, mass concentration, heat and mass fluxes on the plane containing the fracture. The solution of (4.27) by the inversion of the matrix operator ** u**]]

^{(−)}can then be used in (4.28) for evaluating the traction ahead of the crack tip.

Note that, in the case where the Dundurs parameter *d* vanishes, the operators (4.29) and (4.30) can be written as
*α*=0 for the operators (4.38), *q*_{t}=*m*_{t}=*q*_{c}=*m*_{c}=0 for the (4.33) and (4.34) and *h*_{t}=*h*_{c}=0 for the (4.35) and (4.36). Considering an homogeneous material and assuming symmetrical mechanical load ([[*p*_{1}]]=[[*p*_{2}]]=0) together with negligible thermodiffusive effects, it can be easily verified that equations (4.39) become identical to expressions derived by Rice [43].

Different from the general case (4.27), where the inversion of the matrix operator

## 5. Illustrative examples

In this section, the derived general integral identities are applied in order to study the effects of two particular temperature and heat flux profiles on crack opening and traction ahead of the tip. The temperature and heat flux functions used in the illustrative examples here proposed are derived in the electronic supplementary material by the solution of the Laplace's equation in a semi-plane considering two different sets of boundary conditions. Note that, because of the same form of the associate vector integral operators (see expressions (4.33) and (4.34) and (4.35) and (4.36)), the results here reported for the temperature field can be immediately generalized for studying the effects of mass concentration and mass flux variations at the interface.

### (a) Punctually localized heat flux profile at the interface

The following punctually localized profiles for the average and the jump of normal heat flux across the plane *x*_{2}=0, containing both the crack and the interface, are assumed:
*a*_{1},*a*_{2},*L*>0 and *a*_{1}>*a*_{2}. It can be easily verified that the functions (5.1) and (5.2) satisfy the integral balance conditions introduced in §3. The average and the jump of the temperature at *x*_{2}=0 corresponding to the heat flux profiles (5.1) and (5.2) can be evaluated solving the steady state heat conduction equation for both the upper and the lower thermodiffusive half-space. This procedure is reported in details in the electronic supplementary material, and the resulting average and jump of the temperature are given by
*x*_{1}/*L* assuming *a*_{1}/*L*=4 and *a*_{2}/*L*=2.

For simplicity, the Dundurs parameter *d* involved in the weight functions (4.28) is assumed to be zero, and then the integral identities (4.19) and (4.28) decouple and equations (4.27) take the form (4.39). Considering only the contributions due to the fluxes and temperature profiles, the uncoupled integral identities (4.39) become

Substituting the (5.5) and (5.6) into equations (5.4), we get
*et al.* [21]. This leads to

The tractions ahead of the crack tip are given by the solution of equations (4.28), which in this case assume the form:
*H*(*a*_{1}−*x*_{1}) and *H*(*a*_{2}−*x*_{1}) are Heaviside step functions (see the electronic supplementary material for details regarding the derivation of (5.13)). Substituting expressions (5.13) and (5.14) into equations (5.12), the tractions ahead of the crack tip finally become

In order to study the influence of the different thermal expansion coefficients of the materials on the crack opening and on the tractions ahead of the tip, the upper and lower half-space media are assumed to have identical values for the elastic moduli λ^{+}=λ^{−}=λ, *μ*^{+}=*μ*^{−}=*μ*, and for the thermal conductivity coefficients *k*^{+}_{t}=*k*^{−}_{t}=*k*_{t}, and dissimilar values for *γ*_{t}: *σ*_{21}〉^{(+)} vanishes independently of the results derived by applying the integral operators to the temperature profiles. This means that for the case of a plane crack in an homogeneous elastic thermodiffusive material, no shear stresses are present because of the temperature and heat flux profiles on the crack surface.

The variation of the normalized crack opening (5.18) is reported in figure 3 as a function of the spatial coordinate *x*_{1}/*L* assuming a single value of the ratio *a*_{1}/*L* and *a*_{2}/*L*. The values for these distances have been chosen such that *a*_{1}>*a*_{2}, as imposed in the definition of the heat flux profiles (5.1) and (5.2). The crack opening profiles shown in figure 3*a* have been computed considering the same value of *a*_{1}/*L*=0.8 and four different values of *a*_{2}/*L*={0.001,0.1,1,2}. It can be observed that, as *a*_{2}/*L* decreases, and then the contribution to the heat flux functions depending by *δ*((*x*_{1}−*a*_{2})/*L*) approaches the crack tip, the crack opening increases. In the authors' opinion, this is because of the fact that, for small values of *a*_{2}/*L*, the peak of the heat flux profiles (5.1) and (5.2) and the maximum of the temperature (5.3) are localized near to the crack tip. Consequently, the maximum of the thermal load provided by these distributions of heat flux and temperature approaches the crack tip, the crack process is favoured and an increase of the crack opening is produced, in analogy to what is observed for mechanical loading functions localized near to the tip [45].

In figure 3*b*, the crack opening [[*u*_{2}]]^{(−)} is reported for the same value of *a*_{2}/*L*=2 and four different values of *a*_{1}/*L*={2.1,2.5,3,4}. Differently from what is detected for the variation of *a*_{2}/*L*, it can be noted that the crack opening becomes larger as the normalized distance *a*_{1}/*L* increases. As we expect, [[*u*_{2}]]^{(−)} is almost zero for *a*_{1}/*L*=2.1, as in this case *a*_{1}≈*a*_{2}, and then observing expressions (5.1)–(5.3) it can be easily deduced that 〈*q*_{1}〉≈0, [[*q*_{1}]]≈0 and 〈*θ*〉≈0. Looking at the structure of profiles (5.1)–(5.3), it can be easily deduced that if the value of *a*_{2}/*L* is fixed and *a*_{1}/*L* increase, the effect of the temperature and heat fluxes becomes more relevant. Consequently, the loading action provided by the temperature and the flux at the interface increases and gives rise to higher values of the crack opening.

The variation of the normalized traction 〈*σ*_{22}〉^{(+)} is plotted in figure 4 as a function of the spatial coordinate *x*_{1}/*L* for the same single value of the ratio *a*_{1}/*L* and *a*_{2}/*L* assumed for the crack opening. In order to compute profiles of the traction 〈*σ*_{22}〉^{(+)}, the symbolic computation program Mathematica was used. As it can be expected, the tractions exhibit a singular behaviour at the crack tip, for *x*_{1}=0, and for *x*_{1}=*a*_{1} and *x*_{1}=*a*_{2}, where the temperature profile (5.3) diverges.

In figure 5, the crack opening [[*u*_{2}]]^{(−)} and the traction ahead of the tip are reported as functions of *x*_{1}/*L* for single values of *a*_{1}/*L*=4 and of *a*_{2}/*L*=2 and several values of the ratio

The traction expressions (5.20) can now be used for calculating the stress intensity factors. Applying the standard definition reported in Freund [45], we obtain
*a*, the variation of the normalized stress intensity factor *K*_{I} with the ratio *a*_{2}/*L* is reported for *a*_{1}/*L*=4 and several values the ratio *K*_{I}=0 for *a*_{2}/*L*=*a*_{1}/*L*=4. This is because of the fact that for *a*_{2}/*L*=*a*_{1}/*L* both the heat flux and temperature profiles vanish (see expressions (5.1)–(5.3)), and consequently no thermal load is present at the interface and no crack opening is generated. Moreover, it can be noted that as *a*_{2}/*L* decreases, *K*_{I} increases monotonically. This behaviour is in agreement with what is detected analysing the crack opening, and also in this case it can be explained observing that for small values of *a*_{2}/*L* the heat flux profiles (5.1) and (5.2) and the temperature (5.3) possess a peak localized at the crack tip. As a consequence, the maximum of the provided thermal load approaches the crack tip and the crack process is favoured. In figure 6*b*, *K*_{I} is plotted as a function of *a*_{1}/*L* assuming *a*_{2}/*L*=2 and the same values *a*, *K*_{I}=0 for *a*_{2}/*L*=*a*_{1}/*L*=2. The reported curves show that the stress intensity factor becomes larger as *a*_{1}/*L* increases. This behaviour is because of the characteristics of the heat flux and temperature expressions (5.1)–(5.3): looking at the structure of these profiles, it can be easily verified that if the value of *a*_{2}/*L* is fixed and *a*_{1}/*L*, the effects of temperature and heat fluxes becomes more relevant, then the loading action provided by thermal effects increases. Consequently, the crack process is favoured and the value of *K*_{I} becomes larger as does the crack opening (figure 3*b*). Observing both figure 6*a*,*b*, it can be noted that the stress intensity factor increases as the ratio

Considering the case of a homogeneous material, where we have *u*_{1}]]^{(−)}=〈*σ*_{21}〉^{(+)}=0, the crack opening (5.18) assumes the form

### (b) Symmetrically distributed temperature profile at the interface

The following distribution profile for the average temperature across the plane *x*_{2}=0, containing both the crack and the interface, is considered:
*θ*]](*x*_{1})=0 for *q*_{2} at the interface corresponding to the average temperature profile (5.26) and to the assumption [[*θ*]](*x*_{1})=0 has been evaluated solving the steady-state heat conduction equation for both the upper and the lower thermodiffusive half-space. This procedure is reported details in the electronic supplementary material, and the resulting average and jump of the of the normal heat flux are given by
*x*_{1}/*L* in figure 7.

Also for this example, the Dundurs parameter *d* involved in the weight functions (4.18) is assumed to be zero, and then considering the contributions due to the fluxes and temperature profiles, the uncoupled integral identities have the form (5.4). Applying the integral operator

where *Υ*_{t} is the same bimaterial parameter defined in §5a. Using the procedure reported in Piccolroaz & Mishuris [20] and Morini *et al.* [21], the singular integral operator

The results provided by the application of the operators *u*_{1}]]^{(−)} and [[*u*_{2}]]^{(−)} are derived by integrating expressions (5.31) and (5.32), respectively.

Applying the operators

Also for the example described in this section, the upper and lower half-space media are assumed to have identical values for the elastic moduli λ^{+}=λ^{−}=λ, *μ*^{+}=*μ*^{−}=*μ*, and for the thermal conductivity coefficients *k*^{+}_{t}=*k*^{−}_{t}=*k*_{t}, and dissimilar values for *γ*_{t}:

The normalized profiles of the crack opening components [[*u*_{1}]]^{(−)} and [[*u*_{2}]]^{(−)}, obtained by integrating expressions (5.36) and (5.37), are reported in figure 8 as function of *x*_{1}/*L*. Four different values of the ratio *u*_{1}]]^{(−)}=0. This behaviour can be directly deduced observing expression (5.36). Observing figure 8*a*, it can be noted that the crack opening component [[*u*_{1}]]^{(−)} changes sign at a determinate distance *a*/*L* and *b*/*L*. For *u*_{1}]]^{(−)} is positive near to the crack tip, and then becomes negative for *u*_{2}]]^{(−)} illustrated in figure 8*b* is positive everywhere for any value of *x*_{2}-direction, corresponding to the propagation Mode I we have crack opening for any value of the thermal expansion coefficients of the materials.

The variation of the normalized tractions 〈*σ*_{21}〉^{(+)} and 〈*σ*_{22}〉^{(+)} is plotted in figure 9 as a function of the spatial coordinate *x*_{1}/*L* for the same values of *σ*_{21}〉^{(+)} and 〈*σ*_{22}〉^{(+)} exhibit a singular behaviour at the crack tip, for *x*_{1}=0. This behaviour is associated with the leading term of the asymptotic expansion of expressions (5.38) and (5.39) for

Also in this case, exact expressions for the stress intensity factors can be evaluated by means of the tractions (5.38) and (5.39). Applying the standard definition reported in Freund [45], we obtain

## 6. Conclusion

The problem of a quasi-static semi-infinite interfacial crack between two dissimilar thermodiffusive elastic materials has been formulated in terms of boundary integral equations by means of Betti's reciprocity identity and weight functions. For the case of a plane strain crack, explicit integral identities relating the applied mechanical loading, the temperature, mass density, heat and mass flux distributions at the interface and the resulting crack opening have been obtained. The solution of this system of singular integral equations provides explicit expressions for the crack opening and the stress fields associate to any arbitrary harmonic temperature and mass density profile and self-balanced heat or mass flux distribution at the interface.

Illustrative examples of the application of the integral identities to crack problems characterized by punctually localized and distributed heat flux profiles on the interface were performed. The crack opening and tractions ahead of the tip corresponding to the assumed temperature and heat flux distributions were derived by analytical inversion of the singular integral operators. In addition, exact expressions for the stress intensity factors were obtained.

The proposed original method provides a general boundary integral formulation for interfacial crack problems in thermodiffusive bimaterials avoiding the use of the Green's function and the related challenging numerical calculations. The integral identities can be applied for studying the effects of any harmonic temperature and mass concentration profile as well as of any self-balanced heat and mass flux distributions on crack opening and traction fields ahead of the tip. This approach can have various relevant applications, especially in the modelling of fracture processes induced by thermal and diffusive stresses at the interface between different components of multi-layered electrochemical devices such as solid oxide fuel cells and lithium ions batteries. Moreover, the derived integral identities have their own value from the mathematical point of view, because, to the authors best knowledge, a similar explicit formulation in terms of singular integral equations seems to be unknown in the literature.

## Data accessibility

This work does not have any experimental data. Calculations were performed using Wolfram's MATHEMATICA v. 9.

## Authors' contributions

L.M. and A.P. developed the problem's formulation, carried out the analytical solution and the illustrative examples and co-wrote the paper.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by the Italian Ministry of Education, University and Research in the framework of the FIRB project 2010 ‘Structural mechanics models for renewable energy applications’.

- Received May 14, 2015.
- Accepted July 20, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.