A thermo-elastoplastic self-consistent homogenization method for inter-granular plasticity with application to thermal ratcheting of TATB

Abstract

A novel thermo-elastoplastic self-consistent homogenization model for granular materials that exhibit inter-granular plasticity is presented. The model, TEPSCA, is made possible by identifying a new inter-granular plastic Eshelby-like tensor. A micromechanical model of interfacial yielding between grains of a Mohr–Coulomb type is provided, which is relatable to the description of imperfect interfaces within the paradigm of self-consistent homogenization. The local grain constitutive laws are consistent with the description of an interphase layer comprised of local pore volume between grains, such that inelastic inter-particle displacements are directly relatable to changes in bulk porosity, i.e., dilation. The model was developed for the purpose of modeling thermally induced plasticity—the phenomenon known as thermal ratcheting or “ratchet growth”—of composites made from the high explosive triaminotrinitrobenzene (TATB). Model simulations are compared to ratchet growth measurements during cyclic thermal loading of a TATB pellet under stress-free conditions.

Introduction

Many materials can be described as granular, being comprised of distinct solid particles (grains) and often possessing some degree of inter-granular porosity. Examples of granular materials include such ubiquitous materials as geologic materials (geomaterials), concretes, and manufactured materials formed from powders or larger grains (cf. [20]). The overall (effective) thermo-mechanical behavior of granular materials is generally attributable to mechanisms occurring at the grain scale (microscale), i.e., granular and inter-granular micromechanical mechanisms (cf. [30]). Determination of the effective thermo-mechanical response of an assembly of grains from micromechanical quantities is thus often of practical interest, especially when effective behaviors are sensitive to the microstructure, i.e., sensitive to the arrangement, morphology, and contact (interface) properties of grains (e.g., [5, 10, 62, 68, 70]).

The process of upscaling from the heterogeneous microscale to a homogeneous mesoscale of the material is generally called in the context of mechanics homogenization, and various specific methodologies have been proposed (cf. [17, 41]), falling broadly into two basic categories: direct numerical simulation of microstructure (e.g., [18, 35,36,37, 63]), and homogenization by analytical means, typically called micromechanical modeling or sometimes statistical mechanics (cf. [13, 32, 44, 48, 52, 60]). Some approaches consist of a mix of both numerical and analytical means (e.g., [56, 66]). Homogenized expressions of overall strength, i.e., yield criteria, derived from microscopic quantities are often also of interest (e.g., [8, 15, 25, 47]).

This work is concerned with analytical homogenization (micromechanical modeling) of bonded granular materials within the context of Self-Consistent Homogenization (SCH), which has its origin in statistical concepts of multi-scale mechanics originally elucidated by Hershey [26] and Kröner [33], generally credited to have been formalized by Hill [28, 29]. SCH lends itself well to assemblies of grains, especially when the grains are bonded [6, 39]. This is because the underlying idea behind SCH is that the local mechanical fields associated with each grain are solved (typically making use of Eshelby’s solution [21] for local fluctuations in the strain field) considering that the grain is embedded in a matrix that consists of all the other grains, including itself, thus the origin of the term “self-consistent” [11].

An incremental elastoplastic SCH method was originally proposed by Hill [28]; however, due to lack of accounting, roughly speaking, for relaxation of stresses caused by plastic stains, it was found to be in general too stiff, motivating the development of various modifications [7, 40, 53, 64]. Molinari et al. [42] (see also [43]) proposed a fundamental approach of treating elastic and (visco)plastic interactions between grains and their surroundings separately but contributing concurrently to the local fluctuations in the strain fields, which has since found much utility by others (e.g., [61]). An alternative approach to that of Molinari et al. [42] for elastic-plastic coupling makes use of certain projection operators and translation techniques [46, 54]. Some self-consistent homogenization methods inclusive of inter-granular plasticity can be found, though lacking generality in the sense of being restricted to certain materials/applications (e.g., [2, 23, 55]).

A new Thermo-ElastoPlastic Self-Consistent bonded-particle Assembly (TEPSCA) homogenization method and model is presented in this work. A general SCH solution for inter-granular yielding and plastic flow of granular assemblies is made possible by the identification of a novel inter-granular plasticity Eshelby-like tensor. Local yielding laws are viscoplastic; however, the provided framework is meant to be general in the sense of possible future alternative rate-independent yielding, e.g., in a limiting sense or with alternative plastic strain evolution laws. For this reason, we use the term plastic strains to mean in a general sense inelastic strains, although we note that the term viscoplastic (and elasto-viscoplastic) is commonly used in the literature for similar plastic strain evolution laws [34]. The model draws on basic concepts of relating inter-granular damage and porosity from a recently developed thermoelastic modified SCH scheme (TE-SCH) for bonded-particle assemblies [6], where herein extension to inter-particle (visco)plasticity and associated changes in porosity are developed. Inter-granular stress-strain relations from granular micromechanical theory are integrated with continuum matrix-inclusion approaches by defining elastic and inelastic parts of displacement (and velocity) discontinuities across grain boundaries. The TEPSCA model has been developed with the objective of modeling irreversible changes in porosity (bulk plasticity) under thermal loading attributable to heterogeneous and anisotropic thermal strains of the microstructure—what is sometimes called “ratchet growth”—of composites made from the high explosive triaminotrinitrobenzene (TATB), such as the Plastic Bonded eXplosive (PBX) 9502. The driving mechanism of ratchet growth is modeled from basic granular mechanics principles of inter-granular force-displacement relationships congruent with the classical notion of shear induced dilation [22] and non-associative flow [59]. Thus, a micromechanical description of inter-granular (visco)plasticity within the context of SCH is provided, where appropriately derived so-called interaction and concentration tensors are identified. The model is shown to be capable of predicting effective plastic strains induced under thermal loading even at a globally stress-free state, what we refer to herein also as “thermoplastic” strains synonymous with ratchet growth. Model predictions are compared to measurements by Woznick et al. [65] of the ratchet growth of a pressed TATB pellet subjected to thermal cycling, demonstrating the model’s effectiveness.

Notation and material scales

Tensors are indicated throughout in direct notation (boldface), e.g., $\varvec{A}$, although indicial notation (and summation convention), e.g., $A_{ij}$, is sometimes used for clarity. Contraction of indices is implied with dots between tensors, e.g., $a=\varvec{b}{\cdot }\varvec{c}$ is the scalar product of vectors $\varvec{b}$ and $\varvec{c}$ and $a = \varvec{1}{:}\varvec{A}$ is the trace of second order tensor $\varvec{A}$ with $\varvec{1}$ the second order identity. Symmetric dyadic products are specified by the $\overset{s}{\otimes }$ symbol, e.g., $\varvec{A}=\varvec{a}\overset{s}{\otimes }\varvec{b}$ implies that $\varvec{a}$ and $\varvec{b}$ are first order tensors (vectors) and $\varvec{A}$ is a symmetric second order tensor. Fourth order tensors are indicated with double-bold font, e.g., $\mathbb {I}$ is the fourth order symmetric identity, such that $\mathbb {I}{:}\varvec{A} = \varvec{A}$.

Average quantities are denoted with over-bars. For example, the ensemble average of $\varvec{A}$ over volume V is $\overline{\varvec{A}} :=1/V \int _V \varvec{A}\, dV \approx \langle \varvec{A} \rangle := \sum _{\alpha =1}^{\alpha =N} \phi _\alpha \varvec{A}_\alpha$, where the discretized form is exact if $\varvec{A}$ is uniform over each of the $\alpha ^{\mathrm{th}}$ (for $\alpha = 1 \text { to } N$) $\phi _\alpha$ partial volume fractions. The exception to this rule is global (homogenized) fields. The global stress and strain fields are average fields, but for the sake of clarity are demarked in uppercase. The global stress is defined ${\varvec{\Sigma }}:=\langle \overline{\varvec{\sigma }}\rangle$, i.e., the ensemble average of the mean fields over each fraction. Similarly, the global strain is defined $\varvec{E}:= \langle \overline{\varvec{\varepsilon }} \rangle$. The global tangents are distinguished from the local ones by an overbar even though they are not strictly averages of the local ones, e.g., $\varvec{\Sigma } = \overline{\mathbb {L}}^e {:}\varvec{E}^e$ is the global stress constitutive equation, where the e superscript denotes elastic strain and elastic stiffness. Rates are denoted with dots above them, e.g., $\dot{\varvec{E}}$ is the global strain rate, and fluctuations are denoted by tildes, e.g., $\tilde{\varvec{\varepsilon }}$ is the local fluctuation in the strain field and $\dot{\tilde{\overline{\varvec{\varepsilon }}}}$ is the local fluctuation in the average strain rate of a grain.

Fig. 1

Adapted from Bennett et al. [6]

Granular assembly: a slice through a representative volume element (RVE) showing local porosity and imperfect interface associated with a grain; b treatment of a grain within SCH scheme having imperfect interface represented as an interphase layer with thickness of grain’s local porosity, such that changes to it translate to changes in porosity.

Figure 1 shows schematically the division of material scales defined herein, where the mesoscale refers to the homogenized scale of the material at a material point, e.g., the scale of a representative volume element (RVE) having a distinct microstructure [9, 45]. The microscale is the heterogeneous sub-scale of the RVE. Distinction is made herein between the “local” description of the microscale and the “global” description of the mesoscale.

The local and global strain rates are decomposed into elastic, plastic, and thermal parts, respectively,

$$\begin{aligned} \dot{\varvec{\varepsilon }} =\dot{\varvec{\varepsilon }}^e + \dot{\varvec{\varepsilon }}^p +\dot{\varvec{\varepsilon }}^{th} \, , \qquad \dot{\varvec{E}} = \dot{\varvec{E}}^e + \dot{\varvec{E}}^p +\dot{\varvec{E}}^{th}. \end{aligned}$$

(1)

The local and global stress constitutive laws are described in terms of generally anisotropic linear elastic stiffness tensors in rate form by

$$\begin{aligned} \dot{\varvec{\sigma }} = \mathbb {L}^e{:} \dot{\varvec{\varepsilon }}^e \, , \qquad \dot{\varvec{\Sigma }} = \overline{\mathbb {L}}^{e}{:} \dot{\varvec{E}}^e , \end{aligned}$$

(2)

respectively. The plastic strain rate, alternatively, is related to the stress by

$$\begin{aligned} \dot{\varvec{\varepsilon }}^p = \mathbb {M}^p {:} \varvec{\sigma } \, , \qquad \dot{\varvec{E}}^p = \overline{\mathbb {M}}^p {:} \varvec{\Sigma }, \end{aligned}$$

(3)

where we use the general notation $\mathbb {M}^{\bullet } =\mathbb {L}^{\bullet ^{-1}}$. For example,

$$\begin{aligned} \dot{\varvec{\varepsilon }}^e = \mathbb {M}^e{:}\dot{\varvec{\sigma }} \, , \qquad \dot{\varvec{E}}^e = \overline{\mathbb {M}}^e{:} \dot{\varvec{\Sigma }}. \end{aligned}$$

(4)

The local and global plastic compliances, $\mathbb {M}^p$ and $\overline{\mathbb {M}}^p$, respectively, in general change in time for non-trivial cases (are unique over an increment of strain).

Micromechanical model with inter-granular plasticity

Motivated by the desire to treat classical inter-granular force-displacement relations of granular mechanics in a continuum sense, the imperfect interface between grains is treated as a smeared interface region according to the method described by Bennett and Luscher [5, 6] (see Fig. 1). At every point on the surface of each particle, the displacement discontinuity rate (velocity jump), surface traction, and surface traction rate can be computed. The velocity jump between the particle surface and its surroundings is defined,

$$\begin{aligned} \llbracket \varvec{\dot{u}} \rrbracket := \varvec{\dot{u}}^+ -\varvec{\dot{u}}^-, \end{aligned}$$

(5)

for $\varvec{x}$ on S (see Fig. 2). The surface traction rate $\dot{\varvec{t}}$ is assumed continuous across the interface, and resolved in the standard way such that $\dot{\varvec{t}}=\dot{\varvec{t}}_n + \dot{\varvec{t}}_r$, where $\dot{\varvec{t}}_n$ is the normal component and $\dot{\varvec{t}}_r$ is the tangential (see Fig. 3). The velocity jump vector is similarly resolved to its normal $\llbracket {\dot{\varvec{u}}_n}\rrbracket$ and tangential $\llbracket {\dot{\varvec{u}}_r}\rrbracket$ components. Moreover, the velocity jump is additively decomposed into elastic, $\llbracket {\dot{\varvec{u}}^e}\rrbracket$, and inelastic (plastic), $\llbracket {\dot{\varvec{u}}^p }\rrbracket$, parts, i.e.,

$$\begin{aligned} \llbracket \dot{\varvec{u}} \rrbracket =\llbracket {\dot{\varvec{u}}^e }\rrbracket +\llbracket {\dot{\varvec{u}}^p }\rrbracket , \end{aligned}$$

(6)

which can each also be decomposed into normal and tangential parts in the standard way. The magnitudes of the traction and velocity jump are given by $f={\Vert }{\varvec{t}}{\Vert }$ (force per unit area) and $\dot{\delta }={\Vert }{\llbracket {\dot{\varvec{u}}}\rrbracket }{\Vert }$ (speed), where the various components are distinguished with subscripts, e.g., $f = \sqrt{f_n^2 + f_r^2}$ and $\dot{\delta }^p=\sqrt{(\dot{\delta }^p_n)^2 + (\dot{\delta }^{p}_r)^2}$, etc.

Fig. 2

Imperfect interface between grains comprising interphase region with displacement discontinuity (and traction continuity) across it: a SCH treatment of grain interfacing with matrix, and b idealized grain/grain interface reconciled with SCH description as mobilized volume of interface associated with discretized particle surface (see also Fig. 3)

Fig. 3

Discretization of surface with tangent planes: a shown schematically as 2D unit circle with a standard normal orientation for discretization into 16 equal parts of angle $2\varphi$ (dotted lines); b surface patch of ellipsoid with associated tangent plane and surface traction vector, $\varvec{t}$, with normal and tangent parts also shown

The constitutive law relating the elastic part of the displacement discontinuity to the traction can be expressed in either total or incremental (rate) form. The rate form is given by

$$\begin{aligned} \llbracket {\dot{\varvec{u}}^e }\rrbracket :=\varvec{\eta } {{\cdot }} \dot{\varvec{\sigma }}{{\cdot }}\varvec{n} = \varvec{\eta } {{\cdot }} \dot{\varvec{t}}, \end{aligned}$$

(7)

where $\dot{{\varvec{\sigma }}}{{\cdot }}\varvec{n} = \dot{\varvec{t}}$ is the traction increment and the non-evolving in time second order tensor $\varvec{\eta }$ represents the elastic interface compliance [50, 51] such that $\varvec{\eta }\rightarrow {0}$ recovers an elastically perfect interface, and $\varvec{\eta }\rightarrow \infty$ represents complete debonding. The traction-separation law of Eq. (7) is analogous to the traction-separation laws used within cohesive finite elements [1], and (it will be subsequently shown in Sect. 3.2) Eq. (6) is relatable in an integrated sense to additive decomposition of the strain rate into elastic and plastic parts. The elastic interface compliance tensor is given by

$$\begin{aligned} \varvec{\eta } =\alpha \varvec{1} + (\beta -\alpha )\varvec{n} \otimes \varvec{n} , \end{aligned}$$

(8)

where strictly non-negative $\alpha$ and $\beta$ represent tangential and normal compliance, respectively. Bennett et al. [6] described the interphase effective shear modulus $G^{\mathrm{int}}$ and bulk modulus $K^{\mathrm{int}}$ conjugate to the interphase thickness, h, such that

$$\begin{aligned} \alpha = \frac{h}{G^{\mathrm{int}}}, \quad \beta =\frac{h}{K^{\mathrm{int}}}, \end{aligned}$$

(9)

where $G^{\mathrm{int}}$ and $K^{\mathrm{int}}$ are effective properties of an interphase region surrounding the particle that phenomenologically describe the heterogeneous bonding and void space distribution at the inter-particle interfaces. Note that the traction is then an effective traction, and we emphasize that it represents what may be in actuality a variable traction distribution among the grains represented by a single statistically representative grain (SRG) because of the heterogeneous nature of the interface (e.g., including regions of stress concentration). The interphase thickness, h, can be determined from bulk porosity measurements according to the method described by Bennett et al. [6] with the porosity assumed to be completely evenly dispersed. Although it is a subtle distinction, the identification of effective interphase properties notably represents the identification of another, smaller, interfacial scale of the material, which is dealt with in a purely phenomenological sense with respect to inter-granular constitutive laws in what follows.

Plastic yield

Inelastic inter-particle displacement occurs when the traction locally exceeds the strength criteria. A combined cohesive and frictional strength criterion between particles, i.e., a Mohr–Coulomb criterion, is defined

$$\begin{aligned} \mathcal {F}:= f_r -c_r - f_n \xi (q), \end{aligned}$$

(10)

where $c_r$ is the cohesive frictional strength between particles and $\xi (q)$ is a hardening function of the stress-like internal state variable (ISV) q. A hyperbolic hardening law is posed,

$$\begin{aligned} \xi := \frac{ q\tan [\phi ] }{f_n \tan [\phi ] + q}, \end{aligned}$$

(11)

where $\phi$ is the friction angle between particles (also called the internal friction angle). It can be seen from Eq. (11) that the initial slope of the hyperbolic curve is $q/f_n$, and the value of $\xi$ asymptotically approaches that of the apparent inter-particle friction angle (i.e., $\tan [\phi ]$).

The strain-like ISV energy conjugate to q is z, such that $q=Hz$, where H is the hardening modulus. A common option for the evolution of z is to take $z \equiv z_0 + \delta _r^p$, where $z_0$ describes an initial (reference) state (cf. [16]). An alternative option is taken here, which is to provide a separate evolution equation for z. The motivation for this choice is to provide a framework consistent with ISV thermodynamic theory (cf. [4]) that is advantageous for ongoing development of this research and model.

A local non-associative flow rule is posed,

$$\begin{aligned} \llbracket {\dot{\varvec{u}}^p}\rrbracket = \dot{\gamma } \frac{\partial \mathcal {F}}{\partial \varvec{t}_r} +\lambda \dot{\gamma } \varvec{n}, \end{aligned}$$

(12)

which in implicit integrated (incremental) form is,

$$\begin{aligned} \llbracket {\Delta \varvec{u}^p}\rrbracket = \Delta {\gamma } \left. \frac{\partial \mathcal {F}}{\partial \varvec{t}_r} \right| _{(n+1)}+ \lambda \Delta {\gamma } \varvec{n}_{(n+1)}, \end{aligned}$$

(13)

where $\Delta {\gamma }$ is the strictly positive incremental plastic multiplier over an increment of time, $\Delta t = t_{(n+1)} - t_{(n)}$, and $\lambda$ is a dilation coefficient that relates the amount of dilation/contraction (by normal displacement) to the amount of shear sliding. The partial derivative with respect to tangential traction gives the direction of sliding, which is simply the tangential traction’s direction, i.e., $\partial _{\varvec{t}_r} \mathcal {F}= \varvec{t}_r/f_r$. It is emphasized that Eq. (13) holds for any point on the surface; however, the yield function has a distinct value at each point on the surface. It is convenient to rearrange Eq. (13) in order to identify the plastic compliance of the surface interface,

$$\begin{aligned} \llbracket {\Delta \varvec{u}^p}\rrbracket= & {} \Delta {\gamma } \left. \frac{\varvec{t}_r}{f_r}\right| _{(n+1)} +\lambda \Delta {\gamma } \varvec{n}_{(n+1)} \nonumber \\= & {} \left. \Delta {\gamma }\frac{\varvec{t}_r}{f_r}\right| _{(n+1)} + \lambda \Delta {\gamma } \left. \frac{\varvec{t}_n}{f_n}\right| _{(n+1)} \nonumber \\= & {} \varvec{\zeta }_{(n+1)} {\cdot } \varvec{t}_{(n+1)}, \end{aligned}$$

(14)

where, explicitly noting the dependence on (surface) position and dropping the $n+1$ subscript for brevity,

$$\begin{aligned} \varvec{\zeta }(\varvec{x}):= \Delta {\gamma }(\varvec{x}) \Big ( f_r^{-1}(\varvec{x}) \varvec{1} + \left( f_n^{-1}(\varvec{x}) \lambda -f_r^{-1}(\varvec{x})\right) \varvec{n}(\varvec{x})\otimes \varvec{n}(\varvec{x}) \Big ), \end{aligned}$$

(15)

for $\varvec{x} \in S$. The point on the surface where inter-particle slip is predicted to initiate, say $\varvec{x}^*$, is the tangent plane inclined at the angle $\psi ^*$ where the so-called obliquity is greatest, i.e., at $\psi ^* = \text {max}[ \tan ^{-1}[f_r/f_n]]$. For example, the approximation is often made in two dimensions that $\psi ^* \approx 45^{\circ } + \phi /2$ measured from perpendicular to the axis of loading, which is exact for isotropic materials (cf. [16]). For the numerical examples subsequently provided, we find the maximum value of the yield function (Eq. 10) on the surface of each grain in 3D in order to identify $\varvec{x}^*$. In our experience, we have found that $\varvec{x}^*$ remains fixed (or nearly so) for each grain during monotonic loading (and cyclically monotonic loading). In other words, at least for the case of small strain, the location of the tangent planes (the so-called slip-planes) depend on the microstructure of the assembly and the thermo-mechanical loading.

The average plastic modulus of the interphase is identified,

$$\begin{aligned} \overline{\varvec{\zeta }}:= m\varvec{\zeta }, \end{aligned}$$

(16)

where the mobilization factor m is defined as the relative amount of the interphase volume mobilized by the slip, i.e., $m:=v^m/v^{\mathrm{int}}$, with $v^m$ the mobilized volume.

Fig. 4

Inter-granular contact between agglomerate grains: a SEM image of TATB microstructure showing interface between two agglomerate grains (from [6]), and b Idealized 2D agglomerate spherical grain comprised of perfectly aligned and bonded TATB crystal platelets. Slip-plane at ${\overline{n}}^{*}$ is described relative to the local grain axes aligned with crystals’ basal plane normal by angle $\theta$

The mobilized volume is obtained by identifying the mobilized area, $a^m$, of yielding such that $v^m=a^mh$, where h is the interphase thickness. A characteristic particle interaction angle $\varphi$ describes $a^m$ as a function of its mean normal $\overline{\varvec{n}}=\varvec{n}(\varvec{x}^*)$, i.e., the normal to the tangent plane at $\varvec{x}^*$ (also called the standard normal), depicted in Fig. 4b. The interaction angle $\varphi$ is the angle between the mean normal $\overline{\varvec{n}}$ and normals along the perimeter of $a^m$ (see Fig. 3). The mobilized area is related to the discretized tangent plane area $\delta _A$ by

$$\begin{aligned} a^m = a^m_{0} + \delta _A, \end{aligned}$$

(17)

where $a^m_{0}$ describes phenomenologically excess area greater than the discretized tangent plane being mobilized.¹ For a sphere, $\delta _A$ is constant, corresponding to half the angle uniformly discretizing the surface (see Fig. 3). However, for ellipsoidal shaped particles, $\delta _A$ varies with $\varvec{x}^*$ with constant $\varphi$. The relation between the mobilized area on the surface of the ellipsoid $\delta _A$ and the characteristic angle $\varphi$ is provided by the Gauss map between differential areas at $\overline{\varvec{n}}$ and that corresponding to a uniformly discretized sphere in the standard way (cf. [49]), written in terms of Gaussian curvature, K,

$$\begin{aligned} \delta A = \oint _{A_\sigma } {|}{K}{|}d\mathcal {A}_\sigma , \end{aligned}$$

(18)

where $A_\sigma$ is a surface patch (see Fig. 3b).

This allows a local elastoplastic compliance relating the overall displacement rate and the traction to be formed from Eqs. (6), (7), and (14) as

$$\begin{aligned} \llbracket {\varvec{u}}\rrbracket ^e_n + \llbracket {\Delta \varvec{u}}\rrbracket = \underbrace{(\varvec{\eta } +\overline{\varvec{\zeta }})}_{=:\varvec{m}}{\cdot }\varvec{t}_{(n + 1)} = \varvec{m} {\cdot } \varvec{t}_{(n+1)}, \end{aligned}$$

(19)

where $\varvec{m}:=\varvec{\eta } + \overline{\varvec{\zeta }}$ is the (nonlinear) elastoplastic compliance, and which implies

$$\begin{aligned} \llbracket {\Delta \varvec{u}}\rrbracket = \varvec{m} {\cdot } \varvec{t}_{(n+1)} - \varvec{\eta } \cdot \varvec{t}_{(n)} = \overline{\varvec{\zeta }} \cdot \varvec{t}_{(n)} +\varvec{m}\cdot \Delta {\varvec{t}}. \end{aligned}$$

(20)

The evolution of the hardening ISV is given by

$$\begin{aligned} \dot{z} = -\dot{\gamma } \frac{\partial \mathcal {F}}{\partial q}, \implies \Delta {z} =-\Delta \gamma \frac{\partial \mathcal {F}}{\partial q}, \end{aligned}$$

(21)

where

$$\begin{aligned} \frac{\partial \mathcal {F}}{\partial q} = \frac{\xi }{q} \left( 1-\frac{\xi }{\tan [\phi ]}\right) f_n. \end{aligned}$$

(22)

Resolution of local elastic and plastic fields

At each time step, the thermoelastic SCH (TE-SCH) is first solved with an initial guess of a perfectly elastic increment (assuming no particles slip). In this way, a trial stress state,

$$\begin{aligned} \varvec{\sigma }^\alpha _{(tr)} = \varvec{\sigma }^\alpha _{(n)} +\Delta \varvec{\sigma }^\alpha _{(tr)}, \end{aligned}$$

(23)

is obtained for each particle through their respective trial concentration tensors, $\{ \overline{\mathbb {B}}_\alpha ^{(tr)}, \overline{\varvec{b}}_\alpha ^{(tr)}\}$, described in “Localization tensors” section. The subscripts in parentheses denote the trial value, (tr), and the previous value, (n), for the current $(n+1)$ time step. The yield criteria is then evaluated at each particle on the surface tangent plane of maximum obliquity, and all the local rates of elastic and plastic displacement discontinuity are solved for along with all the local interaction and concentration tensors. The TEPSCA homogenization scheme, subsequently described, is then solved with the updated values to obtain new local and global stress values. The process is then repeated ad finitum, i.e., until there is negligible difference between subsequent effective predictions of plastic strain rate and the plastic tangent values.

For each local solution of the traction-jump relation over a time step, a viscoplastic update to the plastic multiplier is obtained. The rate of plastic shearing on a given particle’s critical slip plane is assumed to depend on the current stress state through the power law

$$\begin{aligned} \dot{\gamma } = \dot{\gamma _0} \left( \frac{\mathcal {F}}{\tau }\right) ^{n}, \end{aligned}$$

(24)

where $\dot{\gamma _0}$ is a normalization factor with dimensions of 1 over time, $\tau$ is the so-called fluidity with dimensions of stress, and n is the rate-sensitivity exponent. The (implicit) integrated form of Eq. (24) is

$$\begin{aligned} \Delta \gamma = \Delta \gamma _0 \left( \frac{\mathcal {F}|_{(n+1)}}{\tau }\right) ^{n}. \end{aligned}$$

(25)

This closes the microscale constitutive and evolution equations. The connection to local stress/strain fields is made possible by recognizing from Eq. (14) that $\varvec{m}$ of Eq. (19) can be written

$$\begin{aligned} \varvec{m}=\tilde{\alpha }\varvec{1} +(\tilde{\beta }-\tilde{\alpha })\varvec{n}\otimes \varvec{n}, \end{aligned}$$

(26)

where the coefficients of $\varvec{m}$ are given by

$$\begin{aligned} \tilde{\alpha }&:= \alpha + \Delta {\gamma } f_r^{-1},\nonumber \\ \tilde{\beta }&:= \beta + \lambda \Delta {\gamma }f_n^{-1}. \end{aligned}$$

(27)

The fourth order constitutive tensor relating the interphase strain to the particle stress is then derived in a similar way to the purely elastic one first described by Qu [50, 51],

$$\begin{aligned} \varvec{\varepsilon }_{\llbracket {u}\rrbracket }^{e,(n)} +\Delta \varvec{\varepsilon }_{\llbracket {u}\rrbracket }= & {} \oint _S \left( \llbracket {\varvec{u}}\rrbracket ^e_n + \llbracket {\Delta \varvec{u}}\rrbracket \right) \overset{s}{\otimes }\varvec{n} \, dS\nonumber \\= & {} \oint _S (\varvec{m}{\cdot }\varvec{\sigma }{\cdot }\varvec{n})\overset{s}{\otimes }\varvec{n} \, dS\nonumber \\= & {} \overline{\varvec{\sigma }}{:}\underbrace{\oint _S\varvec{n}\overset{s}{\otimes }\varvec{m}\overset{s}{\otimes }\varvec{n} \, dS}_{=: \mathbb {R}}\nonumber \\= & {} \mathbb {R}{:}\overline{\varvec{\sigma }}, \end{aligned}$$

(28)

where it is clear from Eq. (14) that $\mathbb {R}(\varvec{m}) = \mathbb {R}^e(\varvec{\eta }) +\mathbb {R}^p(\overline{\varvec{\zeta }})$. A constitutive update equation for the traction can then be derived from Eq. (14),

$$\begin{aligned} \varvec{t}= & {} \varvec{m}^{-1}{\cdot } \left( \llbracket {\varvec{u}}\rrbracket ^e_n +\llbracket {\Delta \varvec{u}}\rrbracket \right) \nonumber \\ \varvec{t}= & {} \varvec{m}^{-1}{\cdot }\mathbb {R}{:}\overline{\varvec{\sigma }}{\cdot }\varvec{x}^*. \end{aligned}$$

(29)

The rate form of Eq. (28) can further be split into elastic and plastic parts. The elastic part is the rate form of that used in Bennett et al. [6],

$$\begin{aligned} \dot{\varvec{\varepsilon }}_{\llbracket {u}\rrbracket }^e&= \oint _S \left( \llbracket {\dot{\varvec{u}}^e}\rrbracket \right) \overset{s}{\otimes }\varvec{n} \, dS\nonumber \\&=\oint _S (\varvec{\eta }{\cdot }\dot{\varvec{\sigma }}{\cdot }\varvec{n}) \overset{s}{\otimes }\varvec{n} \, dS\nonumber \\&\approx \dot{\overline{\varvec{\sigma }}}{:}\underbrace{\oint _S\varvec{n} \overset{s}{\otimes }\varvec{\eta }\overset{s}{\otimes }\varvec{n} \, dS}_{=: \mathbb {R}^e}\nonumber \\&=\mathbb {R}^e{:}\dot{\overline{\varvec{\sigma }}}. \end{aligned}$$

(30)

The plastic part is given by

$$\begin{aligned} \dot{\varvec{\varepsilon }}_{\llbracket {u}\rrbracket }^p&= \oint _S \left( \llbracket {\dot{\varvec{u}}^p}\rrbracket \right) \overset{s}{\otimes }\varvec{n} \, dS\nonumber \\&=\oint _S (\overline{\varvec{\zeta }}{\cdot }\varvec{\sigma }{\cdot } \varvec{n})\overset{s}{\otimes }\varvec{n} \, dS\nonumber \\&\approx \overline{\varvec{\sigma }}{:}\underbrace{\oint _S\varvec{n} \overset{s}{\otimes }\overline{\varvec{\zeta }}\overset{s}{\otimes }\varvec{n} \, dS}_{=: \mathbb {R}^p}\nonumber \\&=\mathbb {R}^p{:}\overline{\varvec{\sigma }}. \end{aligned}$$

(31)

Modified Eshelby tensors for elastoplastic imperfect interface

It is emphasized that the interior of the solid particles are modeled as thermoelastic, while the plastic strain occurs only within the surrounding interphase layer, i.e., the plasticity is only inter-granular (not intra-granular). The modification of Eshelby’s tensor due to elastic imperfection of the interface was described by Qu [50, 51],

$$\begin{aligned} \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{te} = \mathbb {S}{:}\dot{\varvec{e}} +\mathbb {S}{:}\mathbb {R}^e{:}\overline{\mathbb {L}}^e{:}(\mathbb {I} -\mathbb {S}){:}\dot{\varvec{e}} = \overline{\mathbb {S}}^M{:}\dot{\varvec{e}}, \end{aligned}$$

(32)

where $\overline{\mathbb {S}}^M:= \mathbb {S} +\mathbb {S}{:}\mathbb {R}^e{:}\overline{\mathbb {L}}^e{:}(\mathbb {I} -\mathbb {S})$. The modification due to the plastic part of the imperfection is novel and is provided in detail within Appendix A; however, the novel plastic interface part is summarized within this section to provide an overall understanding of the relevant equations that make the homogenization of the subsequently described micromechanical model possible.

The approach taken here is to add an additional irreversible contribution to the imperfection (plastic imperfection) quantified by the irreversible part of the velocity jump (Eq. 6). The fluctuation in the strain rate field due to sliding and separation at the interface, $\dot{\tilde{\varvec{\varepsilon }}}^{\delta p}$, is then given by an analogous surface integral as the one identified by Qu for elastic imperfection, but with here the additional plastic part of the velocity jump integrated over the surface,

$$\begin{aligned} \dot{\tilde{\varvec{\varepsilon }}}^{\delta p} = - \int _{S} \mathbb {G}{:}\overline{\mathbb {L}}^e{:}\llbracket \dot{\varvec{u}^p} \rrbracket \overset{s}{\otimes }\varvec{n}\, dS, \end{aligned}$$

(33)

where $\mathbb {G}$ is the second gradient of Green’s function (see Appendix A). Notably, there is an approximation made here that the matrix surrounding the grain is elastic ($\overline{\mathbb {L}}^e$), which allows the contribution from plastic displacements (and their rates) to be derived as an added part of the imperfection (see also further details in Appendix A). It could be argued that the matrix should be in fact elastoplastic, requiring the identification of a single elastoplastic stiffness; however, this would add considerable complexity for which a tractable solution (alternative to the asymptotic solution obtained herein) is not readily evident to us. Moreover, we advocate the method described herein because it obviates the need to identify a single effective elastoplastic tangent, it reduces to the elastic imperfection case for vanishing plasticity (and the classical perfect interface case for totally vanishing imperfection), and it seems to work well in comparison with measurements shown in the “Comparison with measurements” section.

Making use of the continuous-traction/discontinuous-velocity constitutive law of Eq. (14), $\llbracket {\dot{\varvec{u}}^p}\rrbracket = \overline{\varvec{\zeta }} {\cdot } \varvec{t} = \overline{\varvec{\zeta }} {\cdot } \varvec{\sigma }{\cdot } \varvec{n}$, allows the integral to be written in terms of traction, i.e.,

$$\begin{aligned} \dot{\tilde{\varvec{\varepsilon }}}^{\delta p} = - \int _S \mathbb {G}{:}\overline{\mathbb {L}}^e{:}(\overline{\varvec{\zeta }}{\cdot } {\varvec{\sigma }}{\cdot } \varvec{n})\overset{s}{\otimes }\varvec{n} \, dS. \end{aligned}$$

(34)

The traction term can be eliminated making use of Hooke’s law within the inclusion ${\varvec{\sigma }}=\mathbb {L}^e{:}( {\tilde{{{\,\mathrm{\varvec{\varepsilon }}\,}}}}^e -{\varvec{e}})$, to obtain

$$\begin{aligned} \dot{\tilde{\varvec{\varepsilon }}}^{\delta p}= & {} - \int _S \mathbb {G}{:}\overline{\mathbb {L}}^e{:}\mathbb {D}{:}\overline{\mathbb {L}}^e{:} {\tilde{{{\,\mathrm{\varvec{\varepsilon }}\,}}}}^e \, dS \nonumber \\&+ \int _S \mathbb {G}{:}\overline{\mathbb {L}}^e{:} \mathbb {D}{:}\overline{\mathbb {L}}^e{:} {\varvec{e}} \, dS, \end{aligned}$$

(35)

where $\mathbb {D}:=\text {sym}[\varvec{n}\otimes \overline{\varvec{\zeta }} \otimes \varvec{n}]$.

The asymptotic solution for small imperfection of Eq. (35) similar to that provided by Qu [50, 51] is obtained by assigning ${\tilde{{{\,\mathrm{\varvec{\varepsilon }}\,}}}}^{e(0)} =\mathbb {S}{:}{\varvec{e}}$ corresponding to the perfect interface state, in which case

$$\begin{aligned} \dot{\tilde{\varvec{\varepsilon }}}^{\delta p}= & {} - \int _S \mathbb {G}{:}\overline{\mathbb {L}}^e{:}\mathbb {D}{:}\overline{\mathbb {L}}^e{:} {\tilde{{{\,\mathrm{\varvec{\varepsilon }}\,}}}}^{e(0)} \, dS \nonumber \\&+\int _S \mathbb {G}{:}\overline{\mathbb {L}}^e{:} \mathbb {D}{:}\overline{\mathbb {L}}^e{:}\mathbb {S}^{-1}{:} {\tilde{{{\,\mathrm{\varvec{\varepsilon }}\,}}}}^{e(0)} \, dS\nonumber \\= & {} - \int _S \mathbb {G}{:}\overline{\mathbb {L}}^e{:} \mathbb {D}{:}\overline{\mathbb {L}}^e{:}(\mathbb {S}-\mathbb {I}){:} {\varvec{e}} \, dS\nonumber \\= & {} \mathbb {Y}(\varvec{x}){:}(\mathbb {I}-\mathbb {S}){:} {\varvec{e}}, \end{aligned}$$

(36)

where the definition is made, again making explicit note of the dependence on position,

$$\begin{aligned} \mathbb {Y}(\varvec{x}):= \int _S \mathbb {G}(\varvec{x},\varvec{\xi }){:}\overline{\mathbb {L}}^e{:} \mathbb {D}(\varvec{\xi }){:}\overline{\mathbb {L}}^e \, dS(\varvec{\xi }). \end{aligned}$$

(37)

Equation (36) gives rise to the notion of an Eshelby-like tensor, defining

$$\begin{aligned} \mathbb {T}(\varvec{x}):= \mathbb {Y}(\varvec{x}){:} \left( \mathbb {I} - \mathbb {S} \right) . \end{aligned}$$

(38)

Integration provides the relation of the average fields related by $\overline{\mathbb {T}}$, i.e.,

$$\begin{aligned} \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} = \overline{\mathbb {T}}{:}\varvec{e}, \end{aligned}$$

(39)

where

$$\begin{aligned} \overline{\mathbb {T}} :=\mathbb {S}{:}\mathbb {R}^p{:} \overline{\mathbb {L}}^e{:} \left( \mathbb {I} -\mathbb {S}\right) . \end{aligned}$$

(40)

Taking the time derivative of the (average) plastic Eshelby tensor provides,

$$\begin{aligned} \dot{\overline{\mathbb {T}}} :=\mathbb {S}{:}\dot{\mathbb {R}}^p{:} \overline{\mathbb {L}}^e{:} \left( \mathbb {I} -\mathbb {S}\right) , \end{aligned}$$

(41)

where $\dot{\mathbb {R}}^p$ can be evaluated in (implicit) integrated form, i.e., by solving $\mathbb {R}^p_{(n\,+\,1)}$ and storing $\mathbb {R}^p_{n}$.

Making use of the plastic interface Eshelby-like tensor of Eq. (39) and the thermoelastic modified Eshelby tensor of Eq. (32), the total average fluctuation strain rate is

$$\begin{aligned} \dot{\tilde{\overline{\varvec{\varepsilon }}}} = \overline{\mathbb {S}}^M {:} \dot{\bar{\varvec{e}}} +\overline{\mathbb {T}} {:}\bar{\varvec{e}}. \end{aligned}$$

(42)

The additive form of Eqs. (42) and (43) into elastic and plastic parts relating to eigenstrain rate and total eigenstrain, respectively, are consistent with thinking of the plasticity at the interface as an added imperfection, i.e., plastic imperfection. This motivates the additive decomposition of the relations between strain rate fluctuations with stress and stress rate fluctuations described in the following section, i.e., additive decomposition of the so-called interaction equations.

Self-consistent thermo-elastoplastic approach

Interaction and localization tensors

Key to elastoplastic SCH is the interaction equation that relates the local fluctuations in the strain rate to the stress rate (and stress) fluctuation fields. Since the relation between the fluctuation in strain rate is related to the eigenstrain and eigenstrain rate through modified Eshelby tensor $\mathbb {S}^M$ and Eshelby-like tensor $\mathbb {T}$, we assume existence of two corresponding interaction tensors in the interaction equation. The average local fields over each grain are thus assumed to be related in an additive form similar to that suggested by Molinari et al. [42], but different in the sense that here what is added to the thermoelastic part is the contribution from the interface plasticity (not plasticity of the solid grain itself), as described previously in “Modified Eshelby tensors for elastoplastic imperfect interface” section. The additive decomposition of the interaction equation is thus assumed to be of the form,

$$\begin{aligned} \dot{\tilde{\overline{\varvec{\varepsilon }}}} = -\widetilde{\mathbb {M}}^e {:}\dot{\tilde{\overline{\varvec{\sigma }}}} -\widetilde{\mathbb {M}}^{\delta p} {:} \tilde{\overline{\varvec{\sigma }}}, \end{aligned}$$

(43)

where the $\widetilde{\mathbb {M}}^\bullet$ are the interaction tensors (also called the overall constraint tensors by Hill [27]),

$$\begin{aligned} \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{te} = -\widetilde{\mathbb {M}}^{te} {:}\dot{\tilde{\overline{\varvec{\sigma }}}} \, , \qquad \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} =-\widetilde{\mathbb {M}}^{\delta p} {:} \tilde{\overline{\varvec{\sigma }}}. \end{aligned}$$

(44)

The te superscript is used to denote thermoelastic, e.g., $\dot{\tilde{\varvec{\varepsilon }}}^{te} = \dot{\tilde{\varvec{\varepsilon }}}^{e} +\dot{\tilde{\varvec{\varepsilon }}}^{th}$. The fluctuation fields are defined as

$$\begin{aligned} \dot{\tilde{\overline{\varvec{\varepsilon }}}} := \dot{\overline{\varvec{\varepsilon }}} - \dot{\varvec{E}}\, , \qquad \dot{\tilde{\overline{\varvec{\sigma }}}}:= \dot{\overline{\varvec{\sigma }}} - \dot{\varvec{\Sigma }} \, , \qquad \tilde{\overline{\varvec{\sigma }}} := \overline{\varvec{\sigma }} -\varvec{\Sigma }. \end{aligned}$$

(45)

The interaction tensors are solved making use of the Eshelby tensor [21] relations between fluctuation in strain and eigenstrain, $\varvec{e}$, e.g., by $\tilde{\varvec{\varepsilon }}^{te} =\mathbb {S}{:}\varvec{e}$. For example, the thermoelastic interaction tensor for perfect inter-particle interface (cf. [27]) is given by

$$\begin{aligned} \widetilde{\mathbb {M}}^{te}|_{\mathrm{perfect interface}} =(\mathbb {I}-{\mathbb {S}})^{-1}{:}{\mathbb {S}}{:}\overline{\mathbb {M}}^{e}, \end{aligned}$$

(46)

and the more general thermoelastic interaction tensor allowing also for imperfect interface (cf. [6]) is given by

$$\begin{aligned} \widetilde{\mathbb {M}}^{te}=(\mathbb {I}-\overline{\mathbb {S}}^M)^{-1}{:} \overline{\mathbb {S}}^M{:}\overline{\mathbb {M}}^{e}, \end{aligned}$$

(47)

where $\overline{\mathbb {S}}^M$ is the average value over the particle of the so-called modified Eshelby tensor inclusive of imperfect interface first derived by Qu [50, 51], such that $\tilde{\overline{\varvec{\varepsilon }}}^{te} =\overline{\mathbb {S}}^M{:}\varvec{e}$.

The interaction equation for the local fluctuation in the strain rate from plastic imperfection of the interface is solved for making use of the classical notion of an equivalent inclusion [21] considering here solely the additional part arising from the plastic imperfection of the interface, $\widetilde{\mathbb {M}}^{te}$ already being solved for in Eq. (47) (and provided in detail in Bennett et al. [6]). For such case, the sequence of derivation steps are as follows:

$$\begin{aligned} \overline{\varvec{\sigma }}&= \overline{\mathbb {L}}^p {:} \dot{\overline{\varvec{\varepsilon }}}' \qquad \text {Equivalent inclusion}\end{aligned}$$

(48a)

$$\begin{aligned} \tilde{\overline{\varvec{\sigma }}} + \varvec{\Sigma }&= \overline{\mathbb {L}}^p {:} \left( \dot{\overline{\varvec{\varepsilon }}}^{ p} -\dot{\varvec{e}} \right) \qquad \text {where } \dot{\overline{\varvec{\varepsilon }}}^{ p} = \dot{\overline{\varvec{\varepsilon }}}' + \dot{\varvec{e}} \end{aligned}$$

(48b)

$$\begin{aligned} \tilde{\overline{\varvec{\sigma }}} + \varvec{\Sigma }&= \overline{\mathbb {L}}^p {:} \left( \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} + \dot{\varvec{E}}^p -\dot{\varvec{e}} \right) \qquad \text {where }\dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} =\dot{\overline{\varvec{\varepsilon }}}^{ p} - \dot{\varvec{E}}^p\end{aligned}$$

(48c)

$$\begin{aligned} \tilde{\overline{\varvec{\sigma }}}&= \overline{\mathbb {L}}^p {:} \left( \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} -\dot{\varvec{e}} \right) \qquad \text {where } \varvec{\Sigma }=\overline{\mathbb {L}}^p{:}\dot{\varvec{E}}^p\end{aligned}$$

(48d)

$$\begin{aligned} \overline{\mathbb {M}}^p{:}\tilde{\overline{\varvec{\sigma }}}&= \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} -\dot{\varvec{e}} \qquad \text {where } \overline{\mathbb {M}}^p{:}\overline{\mathbb {L}}^p=\mathbb {I} \end{aligned}$$

(48e)

$$\begin{aligned} \overline{\mathbb {M}}^p{:}\tilde{\overline{\varvec{\sigma }}}&= \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} -\frac{\varvec{e}_{(n+1)} -\varvec{e}_{(n)}}{\Delta t} \qquad \text {where } \dot{\varvec{e}} =\frac{\varvec{e}_{(n+1)} -\varvec{e}_{(n)}}{\Delta t}\end{aligned}$$

(48f)

$$\begin{aligned} \overline{\mathbb {M}}^p{:}\tilde{\overline{\varvec{\sigma }}}&= \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} - \frac{1}{\Delta t} \overline{\mathbb {T}}_{(n+1)}^{-1} {:} \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} + \frac{1}{\Delta t} \overline{\mathbb {T}}_{(n)}^{-1} {:} \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p}_{(n)} \qquad \text {where } \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} = \overline{\mathbb {T}}{:}\varvec{e}\end{aligned}$$

(48g)

$$\begin{aligned} \overline{\mathbb {M}}^p{:}\tilde{\overline{\varvec{\sigma }}}&= \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} -\dot{\overline{\mathbb {T}}}^{-1}{:} \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} \qquad {\left\{ \begin{array}{ll} \text {where }&{}(\overline{\mathbb {T}}_{(n+1)}^{-1} {:} \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} - \overline{\mathbb {T}}_{(n)}^{-1} {:} \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p}_{(n)})/\Delta t \\ &{}\approx \dot{\overline{\mathbb {T}}}^{-1} {:} \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} \end{array}\right. }, \end{aligned}$$

(48h)

which is rearranged to provide

$$\begin{aligned} \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} = -\underbrace{\left( \mathbb {I} -\dot{\overline{\mathbb {T}}} \right) ^{-1}{:}\dot{\overline{\mathbb {T}}}{:} \overline{\mathbb {M}}^p}_{=:\widetilde{\mathbb {M}}^{\delta p}} {:} \tilde{\overline{\varvec{\sigma }}}, \end{aligned}$$

(49)

defining,

$$\begin{aligned} \widetilde{\mathbb {M}}^{\delta p} := \left( \mathbb {I} -\dot{\overline{\mathbb {T}}} \right) ^{-1} {:}\dot{\overline{\mathbb {T}}}{:}\overline{\mathbb {M}}^p, \end{aligned}$$

(50)

such that

$$\begin{aligned} \dot{\tilde{\overline{\varvec{\varepsilon }}}}^{\delta p} = -\widetilde{\mathbb {M}}^{\delta p}{:}\tilde{\overline{\varvec{\sigma }}}. \end{aligned}$$

(51)

Some discussion on Eq. (48) is warranted. Backward Euler integration was utilized to time-discretize the eigenstrain rate over a time-step from $t_{(n)}$ to $t_{(n+1)}$. The approximation made within the RHS of Eq. (48h) may not be appropriate for large strain increments, which suggests greater accuracy when limiting to relatively small plastic strain increments.

Localization tensors

From the interaction equation (43) and the constitutive laws (Eqs. 2–4) are derived the localization tensors (also called concentration tensors) that relate the local to global stress fields. The method proposed by Zecevic and Lebensohn [69] for perfect interfaces is followed (cf. details provided therein), where the total stress and stress rate localization tensors are generally different,

$$\begin{aligned} \overline{\varvec{\sigma }} = \mathbb {B}{:}\varvec{\Sigma } + \varvec{b} \, , \qquad \dot{\overline{\varvec{\sigma }}} = \hat{\mathbb {B}}{:}\dot{\varvec{\Sigma }} + \hat{\varvec{b}}. \end{aligned}$$

(52)

The definition of both total stress and stress rate localization tensors is necessary for derivation of effective properties [69]. The localization equation in terms of the total stress provides the forms of $\mathbb {B}$ and $\varvec{b}$ for the case of imperfect interfaces:

$$\begin{aligned} \mathbb {B} := \left[ \frac{1}{\Delta t} (\mathbb {M}^e +\widetilde{\mathbb {M}}^e) + (\mathbb {R}^p+ \widetilde{\mathbb {M}}^{\delta p}) \right] ^{-1}{:} \left[ \frac{1}{\Delta t} (\overline{\mathbb {M}}^e + \widetilde{\mathbb {M}}^e) + (\overline{\mathbb {M}}^p + \widetilde{\mathbb {M}}^{\delta p}) \right] , \end{aligned}$$

(53)

and

$$\begin{aligned} \varvec{b}:= & {} \left[ \frac{1}{\Delta t} (\mathbb {M}^e + \widetilde{\mathbb {M}}^e) + (\mathbb {R}^p + \widetilde{\mathbb {M}}^{\delta p}) \right] ^{-1} \nonumber \\&: \left[ \frac{1}{\Delta t} (\mathbb {M}^e + \widetilde{\mathbb {M}}^e){:}\overline{\varvec{\sigma }}_{(n)} - \frac{1}{\Delta t} (\overline{\mathbb {M}}^e + \widetilde{\mathbb {M}}^e) {:}\varvec{\Sigma }_{(n)} -(\dot{\overline{\varvec{\varepsilon }}}^{th} - \dot{\varvec{E}}^{th}) \right] . \end{aligned}$$

(54)

It is found that $\hat{\mathbb {B}}\equiv \mathbb {B}$, but $\hat{\varvec{b}} \ne \varvec{b}$ (cf. [69]). Rather,

$$\begin{aligned} \hat{\varvec{b}}:= & {} \frac{1}{\Delta t} \left[ \frac{1}{\Delta t}(\mathbb {M}^e + \widetilde{\mathbb {M}}^e) + (\mathbb {R}^p +\widetilde{\mathbb {M}}^{\delta p})\right] ^{-1}\nonumber \\&: \left[ (\overline{\mathbb {M}}^p + \widetilde{\mathbb {M}}^p){:} \varvec{\Sigma }_{(n)} - (\mathbb {R}^p + \widetilde{\mathbb {M}}^{\delta p}){:} \overline{\varvec{\sigma }}_{(n)} - ( \dot{\overline{\varvec{\varepsilon }}}^{th} - \dot{\varvec{E}}^{th})\right] . \end{aligned}$$

(55)

For the self-consistent solution allowing for different grain shapes and orientations, we obtain

$$\begin{aligned} \overline{\varvec{\sigma }} = \overline{\mathbb {B}}{:}\varvec{\Sigma } +\overline{\varvec{b}} \, , \qquad \dot{\overline{\varvec{\sigma }}} =\hat{\overline{\mathbb {B}}}{:}\dot{\varvec{\Sigma }} + \hat{\overline{\varvec{b}}}, \end{aligned}$$

(56)

where

$$\begin{aligned} \overline{\mathbb {B}}:=\mathbb {B}\langle \mathbb {B}\rangle ^{-1} \, , \qquad \hat{\overline{\mathbb {B}}}:=\hat{\mathbb {B}}\langle \hat{\mathbb {B}}\rangle ^{-1}, \end{aligned}$$

(57)

and

$$\begin{aligned} \overline{\varvec{b}} := \varvec{b} - \overline{\mathbb {B}}{:}\langle \varvec{b} \rangle \, , \qquad \hat{\overline{\varvec{b}}} := \hat{\varvec{b}} -\hat{\overline{\mathbb {B}}}{:}\langle \hat{\varvec{b}} \rangle . \end{aligned}$$

(58)

SCH equations

Plastic part

Multiplying the stress concentration Eq. (56) by $\mathbb {R}^p$, we obtain the expression for the global plastic strain rate,

$$\begin{aligned} \begin{aligned} \mathbb {R}^p{:}\overline{\varvec{\sigma }}&= \mathbb {R}^p{:}\overline{\mathbb {B}}{:}\varvec{\Sigma } +\mathbb {R}^p{:}\overline{\varvec{b}}\\ \dot{\overline{\varvec{\varepsilon }}}^{p}&=\mathbb {R}^p{:}\overline{\mathbb {B}}{:}\varvec{\Sigma } +\mathbb {R}^p{:}\overline{\varvec{b}}\\ \underbrace{\langle \dot{\overline{\varvec{\varepsilon }}}^{p} \rangle }_{=\dot{\varvec{E}}^p}&=\underbrace{\langle \mathbb {R}^p{:} \overline{\mathbb {B}}\rangle }_{=\overline{\mathbb {M}}^p} {:}\varvec{\Sigma } + \underbrace{\langle \mathbb {R}^p{:}\overline{\varvec{b}}\rangle }_{=:\dot{\varvec{E}}^p_0}\\ \dot{\varvec{E}}^p&= \overline{\mathbb {M}}^p {:}\varvec{\Sigma } + \dot{\varvec{E}}^p_0 \end{aligned} \end{aligned}$$

(59)

where the ensemble average was used, and the global plastic compliance is identified

$$\begin{aligned} \overline{\mathbb {M}}^p := \langle \mathbb {R}^p{:}\overline{\mathbb {B}} \rangle . \end{aligned}$$

(60)

Since plastic strain is only inter-granular under the stated assumptions, Eq. (59) is the total global plastic strain rate. Eq. (59) states that plasticity is possible under thermal loading even in the absence of a global stress field, quantified by what we refer to here as the “ratchet growth” or “thermoplastic” strain rate, $\dot{\varvec{E}}^p_0$.

Thermoelastic part

The proposed total stress concentration tensors in Zecevic and Lebensohn [69] were used for derivation of viscoplastic effective properties, while the stress rate concentration tensors were used for derivation of elastic effective properties in the context of elasto-viscoplastic self-consistent method with perfect interfaces. For the case of imperfect interfaces with thermal strains we adopt a similar approach. The global thermoelastic (and elastic) strain rates can be decomposed into the contribution from the solid grains (exclusive of the interphases), , and the contribution from the interphase compliance, $\dot{\varvec{E}}^{\llbracket {te}\rrbracket }$, i.e.,

(61)

respectively. Obviously, and $\dot{\varvec{E}}^{\llbracket {te}\rrbracket } = \dot{\varvec{E}}^{\llbracket {e}\rrbracket } + \dot{\varvec{E}}^{\llbracket {th}\rrbracket }$.

The expression can be obtained for the contribution to the global elastic strain from the interface compliance by multiplying the stress rate concentration Eq. (56) by $\mathbb {R}^e$,

$$\begin{aligned} \begin{aligned} \mathbb {R}^e{:}\dot{\overline{\varvec{\sigma }}}&= \mathbb {R}^e{:}\hat{\overline{\mathbb {B}}}{:}\dot{\varvec{\Sigma }} +\mathbb {R}^e{:}\hat{\overline{\varvec{b}}}\\ \dot{\overline{\varvec{\varepsilon }}}^{\llbracket {te}\rrbracket }&=\mathbb {R}^e{:}\hat{\overline{{\mathbb {B}}}}{:}\dot{\varvec{\Sigma }} +\mathbb {R}^e{:}\hat{\overline{\varvec{b}}}\\ \underbrace{\langle \dot{\overline{\varvec{\varepsilon }}}^{\llbracket {te}\rrbracket } \rangle }_{=\dot{\varvec{E}}^{\llbracket {te}\rrbracket }} =&\underbrace{\langle \mathbb {R}^e{:}\hat{\overline{\mathbb {B}}}\rangle }_{=\overline{\mathbb {M}}^{\llbracket {e}\rrbracket }} {:}\dot{\varvec{\Sigma }} + \underbrace{\langle \mathbb {R}^e{:}\hat{\overline{\varvec{b}}}\rangle }_{=:\dot{\varvec{E}}^{\llbracket {th}\rrbracket }}\\ \dot{\varvec{E}}^{\llbracket {te}\rrbracket }&= \overline{\mathbb {M}}^{\llbracket {e}\rrbracket }{:}\dot{\varvec{\Sigma }} + \dot{\varvec{E}}^{\llbracket {th}\rrbracket }\\ \dot{\varvec{E}}^{\llbracket {te}\rrbracket }&= \dot{\varvec{E}}^{\llbracket {e}\rrbracket } + \dot{\varvec{E}}^{\llbracket {th}\rrbracket } \end{aligned} \end{aligned}$$

(62)

The thermoelastic strain rate exclusive of the contribution from the interphase compliance is given by

(63)

We note that the expressions for and have essentially the same form as the ones proposed by Zecevic and Lebensohn [69] in the context of elasto-viscoplastic self-consistent method with perfect interfaces, but with thermal strains replacing the so-called back-extrapolated term described in that work.

The global elastic tangent is thus identified through the elastic strain expression,

(64)

where from Eqs. (62) and (63) it is clear that the elastic modulus is solved by

$$\begin{aligned} \overline{\mathbb {M}}^e = \langle ( \mathbb {M}^e + \mathbb {R}^e){:} \hat{\overline{\mathbb {B}}} \rangle . \end{aligned}$$

(65)

The thermal strain rate is identified from Eqs. (62) and (63),

(66)

completely specifying the self consistent equations.

Discussion of SCH equations

Table 1 summarizes the SCH equations. Significantly, the thermoplastic strain rate, $\dot{\varvec{E}}^p_0$, (also called the “ratchet growth rate”) has been identified, which, as can be seen by its dependence on $\mathbb {R}^p$ and $\overline{\varvec{b}}$, is non-zero only in the presence of heterogeneous thermal strains sufficient to cause irreversible inter-particle displacements, irrespective of the global stress-state. In the absence of plasticity, the provided framework reduces exactly to the Modified Self-Consistent (M-SCH) thermoelastic framework provided by Bennett et al. [6]. This is evident by seeing that if $\mathbb {R}^p \rightarrow 0$, then $\dot{\varvec{E}}^p \rightarrow 0$ and $\overline{\mathbb {M}}^p \rightarrow 0$. Similarly, if the interfaces are completely perfect also in the elastic sense, the classical thermoelastic SCH (TE-SCH) solution is recovered. This can be seen by the fact that thermoelastic modified Eshelby tensor of Eq. (32) reduces to the classical Eshelby tensor in the absence of interface imperfection and as $\mathbb {R}^e \rightarrow 0$ (and $\mathbb {R}^p \rightarrow 0$). For perfect interfaces, the localization tensors thus become the classical thermoelastic ones, and the effective compliance and thermal strains are in the limit of perfect interface the classical TE-SCH solutions.

Table 1 Summary of homogenization equations

The numerical implementation of TEPSCA is described in Algorithm Boxs 1 and 2. The stress and temperature driven thermo-mechanical problem is considered. At time $t=t_{(n+1)}$ an increment in stress $\Delta \varvec{\Sigma }$ and temperature $\Delta \theta$ are given. At $t=t_{(n)}$, the initial global state $\theta _{(n)}, \varvec{E}^p_{(n)}, \varvec{E}^e_{(n) }, \varvec{E}^{th}_{(n)}, \varvec{\Sigma }_{(n)}$, and the set of local states $\{\bar{\varvec{\sigma }}, \bar{\varvec{\varepsilon }}^p, \tilde{\bar{\varvec{\varepsilon }}}^{th}, z\}_{(n)}$ (for all grains) are known. The global thermo-elastoplastic constitutive response $\varvec{E}^e_{(n+1)}, \varvec{E}^{th}_{(n+1)}, \varvec{E}^p_{(n+1)}$ are required. In solving the constitutive response, we also determine for numerical analysis the global tangents $\overline{\mathbb {L}}^e$ and $\overline{\mathbb {L}}^p$ and the updated local states $\{\bar{\varvec{\sigma }}, \bar{\varvec{\varepsilon }}^p, \tilde{\bar{\varvec{\varepsilon }}}^{th}, z\}_{(n+1)}$. On first call, a purely thermoelastic increment is assumed with $\{\Delta \gamma =0\}$ to determine trial values.

Comparison with measurements

TEPSCA was used to simulate the thermal cycling of a pressed TATB pellet and compared with the measurements of Woznick et al. [65] in order to demonstrate the applicability of the model. Ratchet growth measurements of a neat pressed TATB specimen were chosen over those of specimens inclusive of polymer binder (e.g., PBX 9502) because PBX are known to exhibit complex nonlinear thermal expansion behavior due to the presence of the dispersed polymer binder (cf. [6]), which was desired to be avoided in exemplifying the novel developments of TEPSCA. Future work may incorporate inter-particle binder effects similar to those discussed by Bennett et al. [6], which we suspect especially influence observed creep behavior of PBX 9502 and are known to reduce somewhat the overall susceptibility to ratchet growth (cf. [58]).

TATB has a triclinic crystal structure, and its single particle morphology can be described as platelet-like (or graphite-like) corresponding to the (001) planes [12, 14]. TATB crystal thermal expansion is approximately 20 times higher along the c-axis (basal pole) than for the a- and b-lattice directions [31, 67] and the crystals possess strong elastic anisotropy as well [3, 67]. During uniaxial pressing, the platelet-like TATB crystals tend to become preferentially aligned perpendicular to the axis of pressing, resulting in effectively transverse isotropy of the pellet (cf. [6, 12, 38]). Previous texture measurements from identically pressed TATB pellets reported by Yeager et al. [67] were utilized according to the procedure described in Bennett et al. [6] for determining effective thermoelastic anisotropy. Figure 5a shows the texture measurements of Yeager et al. [67] expressed as the probability distribution of basal plane normals aligned with the material axis of transverse isotropy.

Fig. 5

Pole figures showing crystal orientation by basal plane normal distribution with sample axis of transverse isotropy aligned with the center of pole figure: a TATB neat-pressed pellet texture after measurements of Yeager et al. [67], and b mobilized grain orientations at onset of plasticity where approximately 1/4 of the grains (252 of them) are slipping. Units are “times random” probability

Figure 4a shows TATB micostructure, where the single platelet-like TATB crystals stack together to form agglomerate particles, not unlike graphite [57] and clay [24] particles are known to do. Figure 4b shows schematically how the agglomerate particles are treated within TEPSCA. The agglomerates are assumed to be composed of sequentially stacked TATB crystals with perfect bonding (perfect interface) between crystals, forming spherical shaped agglomerate grains. In principle, more complex crystal structure of the agglomerates could be included by performing TE-SCH for agglomerate scale properties; however, for simplicity, the agglomerates are assigned single crystal thermoelastic properties corresponding to the mentioned perfectly interfacing stacked TATB crystals. Inter-granular plastic slip is thus considered to occur between the agglomerate particles only, i.e., the solid particle agglomerates themselves are modeled as thermoelastic. The motivation for this is twofold, stemming from the highly anisotropic thermal expansion properties of TATB: (i) stacked TATB crystals will not develop strong interface forces (since they thermally expand essentially together), and (ii) shear traction forces sufficient for yield are predicted to develop on the thin edges of the TATB crystals where the displacement from thermal expansion is greatest (i.e., towards the left and right edges of the agglomerate in Fig. 4b), such that agglomerate interfaces like that shown in Fig. 4a are predicted to develop large shear forces over a relatively large area of the agglomerate grain surface, which we hypothesize is the root cause of ratchet growth.

The ratchet growth measurements consisted of cyclic thermal loading at a rate of $1\,^{\circ }\text {C}$/min of an unconfined and stress-free specimen from room temperature ($23\,^{\circ }\text {C}$) up to a maximum temperature of $153\,^{\circ }\text {C}$. Hold times at maximum and minimum temperatures for each thermal load cycle were 10 min each, with a total of 18 load cycles performed over approximately three and a half days. Thermally induced strain, consisting of recoverable (thermal) strain and irrecoverable (thermoplastic) strain were measured along the pressing axis of the pellet (coaxial with the axis of transverse isotropy). Table 2 lists the specimen material properties that are directly input into TEPSCA (without calibration). Complete details of the experimental methods and procedure are available in Woznick et al. [65].

The single crystal TATB elastic constants were taken from Bedrov et al. [3] to be

$$\begin{aligned} {[}\mathbb {L}^e]=\begin{bmatrix} 65.7 &{}\quad 18.5 &{}\quad 4.0 &{}\quad -0.2 &{}\quad -1.0 &{}\quad 1.0 \\ &{}\quad 62.0 &{}\quad 5.0 &{}\quad 0.6 &{}\quad -0.5 &{}\quad 1.0 \\ &{} &{}\quad 18.3 &{}\quad 0.2 &{}\quad -0.4 &{}\quad -0.4 \\ &{} \text {sym.} &{} &{}\quad 1.4 &{}\quad 0.1 &{}\quad 0.3 \\ &{} &{} &{} &{}\quad 0.7 &{}\quad 0.4 \\ &{} &{} &{} &{} &{} \quad 21.6 \end{bmatrix} \text {GPa}. \end{aligned}$$

Single crystal anisotropic coefficients of thermal expansion for TATB were taken from the fit by Luscher et al. [38] to the measurements of Kolb and Rizzo [31] reported in Table 3. The anisotropic elastic and thermal expansion constants ($\mathbb {L}^e$ and $\varvec{\varepsilon }^{th}(\theta )$) are direct inputs into the model without any need for calibration. Elastic interface-damage parameters $\alpha$ and $\beta$ were calibrated to best fit the unload curves in Fig. 6b. The plastic interface parameters were then calibrated to fit the nonlinear accumulation of plastic strain. All calibrated TEPSCA model parameters are provided in Table 4.

Table 2 Specimen material properties

Fig. 6

Comparison of TEPSCA simulations of the ratchet growth of a TATB pellet to the measurements of Woznick et al. [65]: a temperature load was cycled from 23 to 153 $^{\circ }\text {C}$ at a rate of $1\,^{\circ }\text {C}$/min a total of 18 times with 10 min hold times at maximum and minimum temperatures, and b TEPSCA simulation results are compared with measured strain response

Fig. 7

TEPSCA simulation results: a break up of global strain into thermal and thermoplastic (“ratchet growth”) parts, and b evolution of strain-like ISV z for all 1000 grains of RVE showing some grains slip (thus harden) more than others

Fig. 8

Location of mobilized planes relative to TATB basal plane within agglomerate ($\theta$) (compare with Fig. 4b). Notably, most particles are predicted to yield in-plane with the basal plane, where shear stresses and strains are greatest due to thermal expansion anisotropy of TATB crystals (near $\theta = 0$)

Figure 6 shows the thermal loading time history and compares TEPSCA simulation results with the measurements of Woznick et al. [65]. Figure 7a compares simulation thermoelastic (reversible) and plastic (accumulated) strains during the cyclic thermal loading, and Fig. 7b shows the evolution of the ISV z for all grains. The variability of evolution of z between grains is due to the heterogeneity, where some grains are predicted to slip before and more than others, evident in Fig. 7b, where the overall (effective) strain is the cumulative effect of the many slipping grains. Only approximately 1/4 of the grains become mobilized at the onset of plasticity. From Fig. 5b, it appears that a more or less uniform sampling of grains become mobilized, i.e., no certain grain orientation exhibit strong preferential slip (although some subtle preferential slip of particles aligned with the material axis arguably may be evident). Figure 8 shows the location of slip planes relative to each grain’s local coordinates (see also Fig. 4b). The grains clearly exhibit preferential slip in the direction perpendicular to the basal plane, supportive of our hypothesis that the anisotropic thermal expansion of the TATB crystals leads to the build up of strong shear forces along these interfaces, driving ratchet growth.

Table 3 Polynomial fitting coefficients of the form $c_0 + c_1T +c_2 T^2$, where T is temperature, for single crystal TATB lattice parameter variation fit to the measurements of [31]

Table 4 Calibrated micromechanical parameters

Discussion of comparison with measurements

As can be seen in Fig. 6, the simulation results match the last unload curve very closely, but do not match the first quite as well. Notably, Eq. (9) states that the accumulation of porosity (plastic volume strain) increases the interface compliance (where herein $\dot{h}\ne 0$); however, although this effect is included in the model, the changes are negligible and the change in slope of the unload curves between the first and last cycles evident in the measurements is not captured well by the model (also see uniform elastic strain cycles of Fig. 7), which would require either additional softening or reversal of plastic flow under reduction of thermal load (cooling). We believe what is missing from the model in order to match more precisely the unloading measurements are two things: (i) evolving damage due to breaking of cohesive bonds, and (ii) plastic slip between grains due to reversal of traction forces during cooling. Since these phenomena are not presently modeled explicitly, we suspect that the elastic interface compliance parameters ($G^\mathrm{{int}}$ and $K^\mathrm{{int}}$) were calibrated to be overly stiff in order to match the thermal strain on the unload. Future work may extend TEPSCA to model evolving damage and plasticity under cooling; however, it appears that the present model provides a reasonably good match to the ratchet growth measurements as is, which exemplifies its effectiveness for modeling inter-granular plasticity arising from heterogeneous and anisotropic thermal expansion of the microstructure.

The viscoplastic behavior evident in the measurements appears to be captured well by the viscoplastic constitutive equations described within Sect. 3. Accumulation of plastic strain exhibits roughly hyperbolic hardening with cycles—see Eq. (11) and compare with measurements in Fig. 6b and also model predictions of plastic strain in Fig. 7a. The measurement data in Fig. 6b does not exactly exhibit hyperbolic hardening, where the trend appears slightly linear out to 18 cycles. Notably, this is somewhat unusual behavior (cf. [58]), and we expect that the hyperbolic hardening law of Eq. (11) should overall perform well in future applications modeling ratchet growth. In our view, a hyperbolic law is appropriate because the accumulation of volumetric plastic strain must asymptote to some maximum value and is also consistent with classical notions of granular plasticity [19], despite this expected behavior not being clearly evident in these measurements.

Fitting the overall strain measurements of Fig. 6 required a large number of trial simulations varying the model parameters reported in Table 4 in order to find a best fit. It could be argued that some other combination of parameter values could as equally well fit the thermal strain measurements, since all parameters needed to be fit simultaneously to this one test. Moreover, no measurement data is presently available to attempt to calibrate the viscoplastic parameters for thermal loading rate-dependence. The plastic response is obviously expected to be sensitive to changes in the rate dependency parameters $\tau$ and n, where large n values approaches rate-independent plasticity. Future applications are planned to look at rate-effects (different thermo-mechanical loading rates), distinguishing between the effect of different choices of $\tau$ and n. Nevertheless, we emphasize that the calibrated parameters satisfy the purpose of suitably matching the measurements while demonstrating the effectiveness of the developed model and homogenization technique for modeling thermally induced plasticity (thermoplastic strain).

Lastly, with regard to the fitting of plastic parameters, we note that since average stress fields (and their tangent plane resolved tractions) drive the inter-granular plasticity in the model, the parameters are likely to be calibrated to effectually weaker interfaces than if the fluctuation of stress fields across SRG’s (i.e., the second moments of stress) were taken into account. Present work is aimed at including the second moments of stress in the model, where it is reasonable to expect that the highest fluctuations in stress will drive plasticity and calibrated interface yield parameters would thus be effectually stronger.

Conclusions

A novel thermo-elastoplastic self-consistent homogenization model for granular materials exhibiting inter-granular plasticity has been provided. The new model, TEPSCA, is made possible by the identification of a novel inter-granular plastic Eshelby-like tensor, not unsimilar to the elastic so-called modified Eshelby tensor of Qu [51] describing imperfect interface between grains (or equivalently an interphase layer responsible for the overall damaged state of the assembly as described in Bennettet al. [6]). A method for integrating the resolved irreversible inter-granular displacements to obtain plastic strains within the paradigm of self-consistent homogenization has been provided, making use of discretization, optimization, and averaging operations over the grain surface. Connection between the Eshelby-like plasticity tensor and Mohr–Coulomb type inter-granular contact yield and hardening laws are further provided through a novel second order viscoplastic interface constitutive tensor. The thermo-elastoplastic self-consistent homogenization equations inclusive of inter-granular plasticity are attained, where the contribution to the effective plastic strain from thermally induced inter-granular plasticity under globally stress-free conditions is identified, called here “thermoplastic strain” and “ratchet growth strain” interchangeably.

The micromechanical description of inter-granular plasticity as the deformation of interphase layers comprised of the local porosity surrounding each grain has been shown to capture well the phenomenon of thermal ratcheting, also called ratchet growth, of TATB pressed pellets by comparison of model simulations to the measurements of Woznick et al. [65]. Notably, TEPSCA is able to predict thermoplastic strains even in the absence of any global stress fields by modeling explicitly the locally heterogeneous and anisotropic stress/strain fields and resultant yielding between individual grains.

By explicitly modeling inter-granular slip (and dilation), the TEPSCA model accesses valuable local information about grain deformation fields (e.g., Figs. 5b, 7, 8). Present efforts to include second moments of the stress field, as mentioned in Sect. 6.1, may provide for interpreting more information from microscale measurements such as neutron diffraction (cf. [6]). Extension to coupled thermomechanical loading is also presently underway, where observed dependence of strength properties on temperature is expected to be able to be modeled with the provided framework with additional accounting for effects of the polymer binder.

Appendix A: Detailed derivation of modified thermoelastic and inter-granular plastic modified Eshelby tensors

The following simultaneous derivation of the thermoelastic modified Eshelby tensor, $\mathbb {S}^M$, and the inter-granular plastic modified Eshelby tensor, $\mathbb {T}$, has similarities to the derivation of the elastic modified Eshelby tensor provided in Qu and Cherkaoui [52], but with a significant difference in the treatment of the discontinuous surface displacements that allows for the identification of an inelastic part. The extension to thermo-elastoplasticity is obtained by introducing a more general perturbation in the displacement field, $\tilde{\varvec{u}}$. Although many of the details have been provided in Qu and Cherkaoui [52], they are repeated here to render a clear description of what is different when irreversible displacements are also included.

Consider a thermoelastic inclusion $\Omega$ as a subdomain within an infinite homogenous elastic matrix D. The equilibrium condition (satisfying balance of linear momentum) for a homogeneous body with an eigenstrain (cf. [44]) can be written in terms of the displacement field,

$$\begin{aligned} L_{ijkl} u_{k,lj} -L_{ijkl}e_{kl,j}=0, \end{aligned}$$

(67)

which we will call EC1.

Green’s function, $\phi _{ij}$, describes the displacement field in response to an arbitrary body force $F_j$ distribution, which can be calculated by the integral over the body $\Omega$,

$$\begin{aligned} u_i(\varvec{x})=\int _V \phi _{ij}(\varvec{x}-\varvec{\xi }) F_j(\varvec{\xi })\, dV \, ,\quad V\equiv \Omega . \end{aligned}$$

(68)

An additional equilibrium requirement can be derived from the integral of the traction over the body’s surface S, making use of the alternative definition for the stress field in terms of Green’s function, $\sigma _{kp}(\varvec{x}) = L_{kpim}\phi _{ij,m}(\varvec{x}-\varvec{\xi })F_j$, and the definition of the Dirac delta function, $\delta (\varvec{x}-\varvec{\xi })$. Equilibrium requires that if $\varvec{\xi } \in V$, then $F_j$ must be balanced by tractions over the volume’s surface S, i.e., since V and $F_j$ are arbitrary, this implies the strong form (cf. [44]),

$$\begin{aligned} L_{kpim}\phi _{ij,mp}(\varvec{x}-\varvec{\xi }) +\delta _{kj}\delta (\varvec{x}-\varvec{\xi }) =0, \end{aligned}$$

(69)

which we will call EC2.

We now introduce a perturbed field (fluctuation) $\tilde{\varvec{u}}\ne \varvec{u}$. More precisely, we will eventually define the fluctuation to contain both reversible and irreversible displacements, but will keep it more general for now. Multiplying Eq. (69) through by $\tilde{\varvec{u}}(\varvec{x})$, recalling that by the chain rule $\tilde{u}_{k,jl}\phi _{ij} =(\tilde{u}_{k,j} \phi _{ij})_{,l} - \tilde{u}_{k,j}\phi _{ji,l}$, and that the strong form of a field equation implies the weak (integral) form, we find

$$\begin{aligned}&\int _V \tilde{u}_{i}(\varvec{x}) L_{ijkl}\phi _{km,lj}(\varvec{x},\varvec{\xi })\, dV(\varvec{x}) + \int _V \tilde{u}_i (\varvec{x}) \delta _{im} \delta (\varvec{x}-\varvec{\xi })\, dV(\varvec{x})= 0\nonumber \\&\int _V L_{ijkl}\left[ \left( \tilde{u}_i(\varvec{x})\phi _{km,l}(\varvec{x}, \varvec{\xi })\right) _{,j} - \tilde{u}_{i,j}(\varvec{x}) \phi _{km,l}(\varvec{x}, \varvec{\xi }) \right] \, dV(\varvec{x}) \nonumber \\&\quad =-\underbrace{\int _V \tilde{u}_i (\varvec{x}) \delta _{im} \delta (\varvec{x}-\varvec{\xi })\, dV(\varvec{x})}_{=-\tilde{u}_m(\varvec{\xi }) \text { if } \varvec{\xi } \in V}\nonumber \\&\int _S L_{ijkl}\tilde{u}_i(\varvec{x})\phi _{km,l}(\varvec{x},\varvec{\xi })n_j(\varvec{x})\, dS(\varvec{x}) -\int _V L_{ijkl}\tilde{u}_{i,j}(\varvec{x}) \phi _{km,l}(\varvec{x},\varvec{\xi }) \, dV(\varvec{x}) \nonumber \\&\quad ={\left\{ \begin{array}{ll} -\tilde{u}_m(\varvec{\xi }), &{} \text {if } \varvec{\xi } \in V\\ 0, &{} \text {if } \varvec{\xi } \notin V. \end{array}\right. } \end{aligned}$$

(70)

Returning to EC1 (67), multiplying through by $\phi _{ij}$ provides

$$\begin{aligned}&L_{ijkl} u_{k,lj}(\varvec{x})\phi _{im}(\varvec{x},\varvec{\xi }) -L_{ijkl}e_{kl,j}(\varvec{x})\phi _{im}(\varvec{x},\varvec{\xi }) =0\nonumber \\&L_{ijkl} \left[ (u_{k,l}(\varvec{x})\phi _{im}(\varvec{x},\varvec{\xi }))_{,j} -u_{k,l}(\varvec{x})\phi _{im,j}(\varvec{x},\varvec{\xi }) \right] -L_{ijkl}e_{kl,j}\phi _{im}(\varvec{x},\varvec{\xi }) =0\nonumber \\&\int _V L_{ijkl} \left[ (u_{k,l}(\varvec{x})\phi _{im}(\varvec{x},\varvec{\xi }))_{,j} -u_{k,l}(\varvec{x})\phi _{im,j}(\varvec{x},\varvec{\xi }) \right] \, dV(\varvec{x})\nonumber \\&\quad -\int _V L_{ijkl}e_{kl,j}(\varvec{x})\phi _{im}(\varvec{x},\varvec{\xi })\, dV(\varvec{x}) =0\nonumber \\&\int _S L_{ijkl} u_{k,l}(\varvec{x})\phi _{im}(\varvec{x},\varvec{\xi })n_j(\varvec{x})\, dS(\varvec{x}) - \int _V L_{ijkl} u_{k,l}(\varvec{x})\phi _{im,j}(\varvec{x},\varvec{\xi }) \, dV(\varvec{x}) \nonumber \\&\quad -\int _VL_{ijkl}e_{kl,j}(\varvec{x})\phi _{im} (\varvec{x},\varvec{\xi })\, dV(\varvec{x}) =0 \end{aligned}$$

(71)

Subtracting the augmented forms of EC1 and EC2, that is the final forms of Eqs. (70) and (71), we obtain

$$\begin{aligned}&\int _S L_{ijkl}[u_{k,l}(\varvec{x})\phi _{im}(\varvec{x},\varvec{\xi }) n_j(\varvec{x}) - \tilde{u}_{i}(\varvec{x})\phi _{km,l}(\varvec{x},\varvec{\xi })n_j(\varvec{x})]\, dS(\varvec{x}) \nonumber \\&\quad - \int _V L_{ijkl}e_{kl,j}(\varvec{x})\phi _{im}(\varvec{x},\varvec{\xi })\, dV(\varvec{x}) = {\left\{ \begin{array}{ll} \tilde{u}_m(\varvec{\xi }), &{} \text {if } \varvec{\xi } \in V\\ 0, &{} \text {if } \varvec{\xi } \notin V, \end{array}\right. } \end{aligned}$$

(72)

and finally, by virtue of recognizing that $\varvec{x}$ in Eq. (72) is a bound (dummy) variable (and that Eq. (72) holds for all $\varvec{\xi } \in V$ and therefore for all $\varvec{x} \in V$), we equivalently obtain²:

$$\begin{aligned}&\int _S L_{ijkl}\left[ u_{k,l}(\varvec{\xi })\phi _{im}(\varvec{\xi },\varvec{x}) n_j(\varvec{\xi }) - \tilde{u}_{i}(\varvec{\xi })\frac{\partial \phi _{km}(\varvec{\xi },\varvec{x})}{\partial \xi _l}n_j(\varvec{\xi })\right] \, dS(\varvec{\xi }) \nonumber \\&\quad - \int _V L_{ijkl}e_{kl,j}(\varvec{\xi })\phi _{im}(\varvec{\xi },\varvec{x})\, dV(\varvec{\xi }) ={\left\{ \begin{array}{ll} \tilde{u}_m(\varvec{x}), &{} \text {if } \varvec{x} \in V\\ 0, &{} \text {if } \varvec{x} \notin V, \end{array}\right. } \end{aligned}$$

(73)

where to be clear, we have written out the long hand notation for the partial derivative of $\varvec{\phi }$ (because it depends on both $\varvec{x}$ and $\varvec{\xi })$, and note that comma subscript notation for functions of a single variable imply the partial derivative with respect to that variable, i.e., $a_{i,j}(\varvec{\xi }) =\partial _{\xi _j}a_i(\varvec{\xi })$. Note also that

$$\begin{aligned} \phi _{mk}(\varvec{x}-\varvec{\xi })=\phi _{km}(\varvec{\xi }-\varvec{x}), \end{aligned}$$

(74)

and

$$\begin{aligned} \phi _{mk,l}(\varvec{x}-\varvec{\xi })=\partial _{x_l}\phi _{mk}(\varvec{x}-\varvec{\xi }) =\partial _{\xi _l}\phi _{km}(\varvec{\xi }-\varvec{x}) =-\partial _{\xi _l}\phi _{mk}(\varvec{x}-\varvec{\xi }). \end{aligned}$$

(75)

Making use of the chain rule and divergence theorem, Eq. (73) becomes

$$\begin{aligned}&\int _S L_{ijkl}\big \{\phi _{im}(\varvec{\xi },\varvec{x})[u_{k,l}(\varvec{\xi })-e_{kl}(\varvec{\xi })] -\tilde{u}_{i}(\varvec{\xi })\frac{\partial \phi _{km}(\varvec{\xi },\varvec{x})}{\partial \xi _l} \big \} n_j(\varvec{\xi })\, dS(\varvec{\xi }) \nonumber \\&\quad + \int _V L_{ijkl}e_{kl}(\varvec{\xi })\phi _{im,j}(\varvec{\xi },\varvec{x})\, dV(\varvec{\xi }) = {\left\{ \begin{array}{ll} \tilde{u}_m(\varvec{x}), &{} \text {if } \varvec{x} \in V\\ 0, &{} \text {if } \varvec{x} \notin V. \end{array}\right. } \end{aligned}$$

(76)

The fields within and without the inclusion can be examined. The traction is assumed continuous across the interface, which in general may be imperfect in the sense of having non-vanishing thickness. The inclusion surface thus is distinguished as $S^-$ when approached from the interior and $S^+$ when approached from the exterior (cf. Fig. 1).

Firstly, consider $\varvec{x}\in \Omega$, and the volume explicitly taken as that of the inclusion ($V\equiv \Omega$). Eq. (76) then becomes

$$\begin{aligned} \tilde{u}_m(\varvec{x})|_{\varvec{x} \in \Omega }= & {} \int _{S^-} L_{ijkl}\big \{\phi _{im}(\varvec{\xi },\varvec{x})[u_{k,l}(\varvec{\xi })-e_{kl}(\varvec{\xi }] - \tilde{u}_{i}(\varvec{\xi })\frac{\partial \phi _{km}(\varvec{\xi },\varvec{x})}{\partial \xi _l} \big \} n_j(\varvec{\xi })\, dS(\varvec{\xi }) \nonumber \\&+ \int _V L_{ijkl}e_{kl}(\varvec{\xi })\phi _{im,j}(\varvec{\xi },\varvec{x})\, dV(\varvec{\xi }). \end{aligned}$$

(77)

For $\varvec{x} \in \Omega \supset S^-$ Hooke’s law given by

$$\begin{aligned} \sigma _{ij}(\varvec{x})= L_{ijkl}[u_{k,l}(\varvec{x})-e_{kl}(\varvec{x})], \end{aligned}$$

(78)

can be used to introduce the stress into the integrand,

$$\begin{aligned} \tilde{u}_m(\varvec{x})|_{\varvec{x} \in \Omega }= & {} \int _{S^-} \big \{\phi _{im}(\varvec{\xi },\varvec{x})\sigma _{ij}(\varvec{\xi }) -L_{ijkl}\tilde{u}_{i}(\varvec{\xi })\frac{\partial \phi _{km}(\varvec{\xi },\varvec{x})}{\partial \xi _l} \big \} n_j(\varvec{\xi })\, dS(\varvec{\xi }) \nonumber \\&+ \int _V L_{ijkl}e_{kl}(\varvec{\xi })\phi _{im,j}(\varvec{\xi },\varvec{x})\, dV(\varvec{\xi }). \end{aligned}$$

(79)

Next, consider the volume to to be the exterior of the inclusion, i.e., $V=D/ \Omega$. By definition, there is no eigenstrain outside of the inclusion, so

$$\begin{aligned} \int _{S^+}\big [\phi _{im}(\varvec{\xi },\varvec{x})\sigma _{ij}(\varvec{\xi }) - L_{ijkl}\tilde{u}_i(\varvec{\xi })\frac{\partial \phi _{km}(\varvec{\xi },\varvec{x})}{\partial \xi _l} \big ] n_j \, dS(\varvec{\xi }) = 0, \quad \forall \varvec{x} \, \text {if } V\nsupseteq \Omega , \end{aligned}$$

(80)

where the continuity of traction across the inclusion was used to keep the stress term (actually traction since it is contracted with the surface normal).

The exterior Eq. (80) can be subtracted from the interior (79) to yield that within the inclusion

$$\begin{aligned} \tilde{u}_m(\varvec{x})|_{\varvec{x} \in \Omega }= & {} \int _V \frac{\partial \phi _{im}}{\partial \xi _j} L_{ijkl} e_{kl}(\varvec{\xi }) \, dV(\varvec{\xi })\nonumber \\&- \int _S \llbracket \tilde{u}_i(\varvec{\xi }) \rrbracket L_{ijkl} \frac{\partial \phi _{km}}{\partial \xi _l} n_j \, dS(\varvec{\xi }), \end{aligned}$$

(81)

for $V=\Omega$.

Assuming the eigenstrain is constant within the inclusion, the non-symmetric total strain field within the inclusion is found

$$\begin{aligned} \tilde{u}_{m,n}(\varvec{x})|_{\varvec{x}\in \Omega }= & {} \int _V \frac{\partial ^2 \phi _{im}}{\partial \xi _j \partial x_n} L_{ijkl} e_{kl}(\varvec{\xi }) \, dV(\varvec{\xi })\nonumber \\&- \int _S \llbracket \tilde{u}_i(\varvec{\xi }) \rrbracket L_{ijkl} \frac{\partial ^2 \phi ^e_{km}}{\partial \xi _l \partial x_n} n_j \, dS(\varvec{\xi }), \end{aligned}$$

(82)

wherefore the tensor having major and minor symmetries is defined³

$$\begin{aligned} G_{ijnm}(\varvec{x},\varvec{\xi }) := \frac{1}{4} \left[ \frac{\partial ^2 \phi _{mi}(\varvec{x},\varvec{\xi })}{\partial x_j \partial \xi _n} +\frac{\partial ^2 \phi _{mj}(\varvec{x},\varvec{\xi })}{\partial x_i \partial \xi _n} + \frac{\partial ^2 \phi _{ni}(\varvec{x},\varvec{\xi })}{\partial x_j \partial \xi _m} +\frac{\partial ^2 \phi _{nj}(\varvec{x},\varvec{\xi })}{\partial x_i \partial \xi _m}\right] . \end{aligned}$$

(83)

We can then write the symmetric strain tensor within the inclusion as

$$\begin{aligned} \tilde{\varepsilon }_{ij}(\varvec{x})|_{\varvec{x}\in \Omega } =L_{mnkl}e_{kl}\int _\Omega G_{mnji}(\varvec{\xi },\varvec{x}) dV(\varvec{\xi }) -\int _S L_{mnkl} \llbracket \tilde{u}_m(\varvec{\xi }) \rrbracket G_{klji}(\varvec{\xi }, \varvec{x}) n_n\, dS(\varvec{\xi }), \end{aligned}$$

(84)

or equivalently making use of the major symmetries of $G_{ijkl}$ and $\mathbb {L}$,

$$\begin{aligned} \tilde{\varepsilon }_{ij}(\varvec{x})|_{\varvec{x}\in \Omega } =L_{klmn}e_{kl}\int _\Omega G_{ijmn}(\varvec{\xi },\varvec{x}) dV(\varvec{\xi }) -\int _S L_{klmn} \llbracket \tilde{u}_k(\varvec{\xi }) \rrbracket G_{ijmn}(\varvec{\xi }, \varvec{x}) n_l\, dS(\varvec{\xi }), \end{aligned}$$

(85)

which we see from the relation (cf. [44]) $S_{ijkl}=L_{mnkl}\int _\Omega G_{ijmn}dV$ is

$$\begin{aligned} \tilde{\varepsilon }_{ij}(\varvec{x})|_{\varvec{x}\in \Omega } = S_{ijkl}e_{kl} - \int _S L_{klmn} \llbracket \tilde{u}_k(\varvec{\xi }) \rrbracket G_{ijmn}(\varvec{\xi }, \varvec{x}) n_l\, dS(\varvec{\xi }). \end{aligned}$$

(86)

Since $S_{ijkl}$ and $G_{ijmn}$ do not change with time (vary with position but are constant $\forall \varvec{x}$), the rate form of Eq. (86) is simply (taking the time derivative),

$$\begin{aligned} \dot{\tilde{\varvec{\varepsilon }}}_{ij}(\varvec{x})|_{\varvec{x}\in \Omega } =S^e_{ijkl}\dot{e}_{kl} - \int _S L_{klmn} \llbracket \dot{\tilde{u}}_k(\varvec{\xi }) \rrbracket G^e_{ijmn} (\varvec{\xi }, \varvec{x}) n_l\, dS(\varvec{\xi }). \end{aligned}$$

(87)

Clearly, the surface integral in Eq. (86) involving the displacement jump quantifies the difference in the strain field due to imperfection of the interface. Making use of the assumption of (constant linear) proportionality between the displacement jump and the traction vector at the interface ($\llbracket \dot{\tilde{u}}_i \rrbracket = \eta _{ij}\dot{\sigma }_{jk}n_k + \bar{\zeta }_{ij}\sigma _{jk}n_k$), the integral can be written in terms of the traction,

$$\begin{aligned} \dot{\tilde{\varepsilon }}_{ij}(\varvec{x}) = S_{ijkl}\dot{e}_{kl} - \int _S L_{klmn}(\eta _{kp}\dot{\sigma }_{pq} +\bar{\zeta }_{kp}\sigma _{pq})G_{ijmn}(\varvec{\xi }, \varvec{x}) n_q n_l\, dS(\varvec{\xi }). \end{aligned}$$

(88)

Alternatively, again making use of Hooke’s law within the inclusion (with no super-imposed fields) $\varvec{\sigma } =\mathbb {L} ({{\,\mathrm{\varvec{\varepsilon }}\,}}-\varvec{e})$ (and $\dot{\varvec{\sigma }}=\mathbb {L}(\dot{{{\,\mathrm{\varvec{\varepsilon }}\,}}} -\dot{\varvec{e}})$), the traction term can be eliminated to obtain

$$\begin{aligned} \dot{\tilde{\varepsilon }}_{ij}(\varvec{x})= & {} S_{ijkl}\dot{e}_{kl} - \int _S L_{klmn}L_{pqst}(\eta _{kp}\dot{\varepsilon }_{st} + \bar{\zeta }_{kp}\varepsilon _{st})G_{ijmn}(\varvec{\xi }, \varvec{x}) n_q n_l\, dS(\varvec{\xi }) \nonumber \\&+ \int _S L_{klmn}L_{pqst}(\eta _{kp}\dot{e}_{st}+\bar{\zeta }_{kp}e_{st}) G_{ijmn}(\varvec{\xi }, \varvec{x}) n_q n_l\, dS(\varvec{\xi }). \end{aligned}$$

(89)

The asymptotic solution for small imperfection is considered, in which case $\dot{\tilde{\varvec{\varepsilon }}}^{(0)}=\mathbb {S}\dot{\varvec{e}}$, and consider a perturbation representing slight imperfection.

$$\begin{aligned} \dot{\tilde{\varepsilon }}_{ij}= & {} \dot{\tilde{\varepsilon }}_{ij}^{(0)} -\int _S L_{klmn}L_{pqst}(\eta _{kp}\dot{\tilde{\varepsilon }}_{st}^{(0)} + \bar{\zeta }_{kp}\tilde{\varepsilon }_{st}^{(0)})G_{ijmn}(\varvec{\xi }, \varvec{x}) n_q n_l\, dS(\varvec{\xi }) \nonumber \\&+ \int _S L_{klmn}L_{pqst}(\eta _{kp} \dot{\tilde{\varepsilon }}_{xz}^{(0)} + \bar{\zeta }_{kp}\tilde{\varepsilon }_{xz}^{(0)})S_{stxz}^{-1} G_{ijmn}(\varvec{\xi }, \varvec{x}) n_q n_l\, dS(\varvec{\xi })\nonumber \\= & {} \dot{\tilde{\varepsilon }}_{ij}^{(0)} - \int _S L_{klmn}L_{pqst}\eta _{kp} \left( \dot{\tilde{\varepsilon }}_{st}^{(0)} - \dot{\tilde{\varepsilon }}_{xz}^{(0)}S_{stxz}^{{-1}}\right) G_{ijmn}(\varvec{\xi }, \varvec{x}) n_q n_l\, dS(\varvec{\xi })\nonumber \\&- \int _S L_{klmn}L_{pqst}\bar{\zeta }_{kp} \left( \tilde{\varepsilon }_{st}^{(0)} - \tilde{\varepsilon }_{xz}^{(0)}S_{stxz}^{-1}\right) G_{ijmn}(\varvec{\xi }, \varvec{x}) n_q n_l\, dS(\varvec{\xi })\nonumber \\= & {} S_{ijkl}\dot{e}_{kl} - \int _S L_{klmn}L_{pqst}\eta _{kp} \left( S_{xzst}\dot{e}_{st}- \dot{e}_{st}\right) G_{ijmn}(\varvec{\xi }, \varvec{x}) n_q n_l\, dS(\varvec{\xi })\nonumber \\&- \int _S L_{klmn}L_{pqst}\bar{\zeta }_{kp} \left( S_{xzst}e_{st}- e_{st}\right) G_{ijmn}(\varvec{\xi }, \varvec{x}) n_q n_l\, dS(\varvec{\xi })\nonumber \\= & {} S_{ijkl}\dot{e}_{kl} +Z_{ijst}(\varvec{x})(I_{stkl}-S_{stkl})\dot{e}_{kl} + Y_{ijst}(\varvec{x})(I_{stkl}-S_{stkl})e_{kl} , \end{aligned}$$

(90)

where we have defined

$$\begin{aligned} Z_{ijst}:= L_{klmn}L_{stpq} \int _S \eta _{kp}G_{ijmn}(\varvec{\xi },\varvec{x})n_q n_l \, dS(\varvec{\xi }), \end{aligned}$$

(91)

and

$$\begin{aligned} Y_{ijst}:= L_{klmn}L_{stpq} \int _S \bar{\zeta }_{kp}G_{ijmn}(\varvec{\xi },\varvec{x})n_q n_l \, dS(\varvec{\xi }), \end{aligned}$$

(92)

The additive decomposition $\dot{\tilde{\varvec{\varepsilon }}} =\dot{\tilde{\varvec{\varepsilon }}}^{te} +\dot{\tilde{\varvec{\varepsilon }}}^{\delta p}$ suggests

$$\begin{aligned} \dot{\tilde{\varepsilon }}_{ij}^{te}(\varvec{x}) = S_{ijkl}\dot{e}_{kl} + Z_{ijst}(\varvec{x})(I_{stkl}-S_{stkl})\dot{e}_{kl} , \end{aligned}$$

(93)

and

$$\begin{aligned} \dot{\tilde{\varepsilon }}_{ij}^{\delta p}(\varvec{x}) =Y_{ijst}(\varvec{x})(I_{stkl}-S_{stkl})e_{kl}, \end{aligned}$$

(94)

such that the modified elastic Eshelby tensor is

$$\begin{aligned} \mathbb {S}^M(\varvec{x}):= \mathbb {S} +\mathbb {Z}(\varvec{x})\left( \mathbb {I} - \mathbb {S} \right) , \end{aligned}$$

(95)

and the inter-granular plastic Eshelby tensor is

$$\begin{aligned} \mathbb {T}(\varvec{x}) :=\mathbb {Y}(\varvec{x}) \left( \mathbb {I}-\mathbb {S}\right) . \end{aligned}$$

(96)

The elastic and plastic strain rates can then be written in terms of the eigenstrain,

(97)

or in total strain rate form,

$$\begin{aligned} \dot{\tilde{\varvec{\varepsilon }}}(\varvec{x}) = \mathbb {S}^M(\varvec{x}) \dot{\varvec{e}} + \mathbb {T}(\varvec{x})\varvec{e}. \end{aligned}$$

(98)