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

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 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., 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 = b·c is the scalar product of vectors b and c and a = 1:A is the trace of second order tensor A with 1 the second order identity.
Symmetric dyadic products are specified by the s ⊗ symbol, e.g., A = a s ⊗ b implies that a and b are first order tensors (vectors) and A is a symmetric second order tensor. Fourth order tensors are indicated with double-bold font, e.g., I is the fourth order symmetric identity, such that I:A = A.
Average quantities are denoted with over-bars. For example, the ensemble average of A over volume V is A := 1/V V A dV ≈ A := α=N α=1 φ α A α , where the discretized form is exact if A is uniform over each of the α th (for α = 1 to N ) φ α 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 := σ , i.e., the ensemble average of the mean fields over each fraction. Similarly, the global strain is defined E := ε . 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., = L e :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.,Ė is the global strain rate, and fluctuations are denoted by tildes, e.g.,ε is the local fluctuation in the strain field andε is the local fluctuation in the average strain rate of a grain. 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.
Composite grain in SCH 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. Adapted from Bennett et al. [6] The local and global stress constitutive laws are described in terms of generally anisotropic linear elastic stiffness tensors in rate form bẏ respectively. The plastic strain rate, alternatively, is related to the stress bẏ where we use the general notation M • = L • −1 . For example, The local and global plastic compliances, M p and 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, for x on S (see Fig. 2). The surface traction rateṫ is assumed continuous across the interface, and resolved in the standard way such thatṫ =ṫ n +ṫ r , whereṫ n is the normal component andṫ r is the tangential (see Fig. 3). The velocity jump vector is similarly resolved to its a b 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  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ϕ (dotted lines); b surface patch of ellipsoid with associated tangent plane and surface traction vector, t, with normal and tangent parts also shown normal u n and tangential u r components. Moreover, the velocity jump is additively decomposed into elastic, u e , and inelastic (plastic), u p , parts, i.e., 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 = t (force per unit area) andδ = u (speed), where the various components are distinguished with subscripts, e.g., f = f 2 n + f 2 r andδ p = (δ p n ) 2 + (δ p r ) 2 , etc. 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 whereσ·n =ṫ is the traction increment and the non-evolving in time second order tensor η represents the elastic interface compliance [50,51] such that η → 0 recovers an elastically perfect interface, and η → ∞ represents complete debonding. The tractionseparation 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 where strictly non-negative α and β represent tangential and normal compliance, respectively. Bennett et al. [6] described the interphase effective shear modulus G int and bulk modulus K int conjugate to the interphase thickness, h, such that where G int and K 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 where c r is the cohesive frictional strength between particles and ξ (q) is a hardening function of the stress-like internal state variable (ISV) q. A hyperbolic hardening law is posed, where φ 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 ξ asymptotically approaches that of the apparent inter-particle friction angle (i.e., tan[φ]). 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 ≡ z 0 + δ p r , 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, u p =γ ∂F ∂t r + λγ n, (12) which in implicit integrated (incremental) form is, where γ is the strictly positive incremental plastic multiplier over an increment of time, t = t (n+1) − t (n) , and λ 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., ∂ t r F = 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, where, explicitly noting the dependence on (surface) position and dropping the n + 1 subscript for brevity, for x ∈ S. The point on the surface where inter-particle slip is predicted to initiate, say x * , is the tangent plane inclined at the angle ψ * where the so-called obliquity is greatest, i.e., at ψ * = max[tan −1 [f r /f n ]]. For example, the approximation is often made in two dimensions that ψ * ≈ 45 • + φ/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 x * . In our experience, we have found that 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, 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 int , with v m the mobilized volume. The mobilized volume is obtained by identifying the mobilized area, a m , of yielding such that v m = a m h, where h is the interphase thickness. A characteristic particle interaction 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 n * is described relative to the local grain axes aligned with crystals' basal plane normal by angle θ angle ϕ describes a m as a function of its mean normal n = n(x * ), i.e., the normal to the tangent plane at x * (also called the standard normal), depicted in Fig. 4b. The interaction angle ϕ is the angle between the mean normal n and normals along the perimeter of a m (see Fig. 3). The mobilized area is related to the discretized tangent plane area δ A by where a m 0 describes phenomenologically excess area greater than the discretized tangent plane being mobilized. 1 For a sphere, δ A is constant, corresponding to half the angle uniformly discretizing the surface (see Fig. 3). However, for ellipsoidal shaped particles, δ A varies with x * with constant ϕ. The relation between the mobilized area on the surface of the ellipsoid δ A and the characteristic angle ϕ is provided by the Gauss map between differential areas at n and that corresponding to a uniformly discretized sphere in the standard way (cf. [49]), written in terms of Gaussian curvature, K , where A σ 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 1 For example, a m 0 may describe multi-grain interactions or other more-complex micro-structural mechanisms that effectually mean that a greater area becomes mobilized than that associated with a certain interface discretization. An alternative approach would be to let a m = δ A ; however, this would impose more of a restriction on the mobilized area, since the concept of a surface patch would become nonsensical at large a m .
where m := η + ζ is the (nonlinear) elastoplastic compliance, and which implies The evolution of the hardening ISV is given bẏ where

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, is obtained for each particle through their respective trial concentration tensors, {B (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 laẇ whereγ 0 is a normalization factor with dimensions of 1 over time, τ is the so-called fluidity with dimensions of stress, and n is the rate-sensitivity exponent. The (implicit) integrated form of Eq. (24) is This closes the microscale constitutive and evolution equations. The connection to local stress/strain fields is made possible by recognizing from Eq. (14) that m of Eq. (19) can be written where the coefficients of m are given bỹ 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], where it is clear from Eq. (14) that R(m) = R e (η) + R p (ζ). A constitutive update equation for the traction can then be derived from Eq. (14), 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], The plastic part is given bẏ

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], where S M := S + S:R e :L e :(I − 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,ε δ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, where 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 (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), u p = ζ·t = ζ·σ·n, allows the integral to be written in terms of traction, i.e.,ε The where The asymptotic solution for small imperfection of Eq. (35) similar to that provided by Qu [50,51] where the definition is made, again making explicit note of the dependence on position, Equation (36) gives rise to the notion of an Eshelby-like tensor, defining Integration provides the relation of the average fields related by T, i.e., where T := S:R p :L e : (I − S) .
Taking the time derivative of the (average) plastic Eshelby tensor provides, whereṘ p can be evaluated in (implicit) integrated form, i.e., by solving R p (n + 1) and storing R 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.

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 S M and Eshelby-like tensor 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, where the M • are the interaction tensors (also called the overall constraint tensors by Hill [27]),ε The te superscript is used to denote thermoelastic, e.g.,ε te =ε e +ε th . The fluctuation fields are defined aṡ The interaction tensors are solved making use of the Eshelby tensor [21] relations between fluctuation in strain and eigenstrain, e, e.g., byε te = S:e. For example, the thermoelastic interaction tensor for perfect inter-particle interface (cf. [27]) is given by and the more general thermoelastic interaction tensor allowing also for imperfect interface (cf. [6]) is given by where 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ε te = S M :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, 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: which is rearranged to providė defining, such thaṫ 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, 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 B and b for the case of imperfect interfaces: It is found thatB ≡ B, butb = b (cf. [69]). Rather, For the self-consistent solution allowing for different grain shapes and orientations, we obtain where and

Plastic part
Multiplying the stress concentration Eq. (56) by R p , we obtain the expression for the global plastic strain rate, where the ensemble average was used, and the global plastic compliance is identified 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,Ė 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 elastoviscoplastic 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),Ė te\ te , and the contribution from the interphase compliance,Ė te , i.e., respectively. Obviously,Ė te\ te =Ė e\ e +Ė th\ th andĖ te =Ė e +Ė th .
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 R e , R e :σ = R e :B:˙ + R e :ḃ The thermoelastic strain rate exclusive of the contribution from the interphase compliance is given bẏ We note that the expressions for M e\ e andĖ th\ th have essentially the same form as the ones proposed by Zecevic and Lebensohn [69] in the context of elasto-viscoplastic selfconsistent 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, where from Eqs. (62) and (63) The thermal strain rate is identified from Eqs. (62) and (63), completely specifying the self consistent equations. Table 1 summarizes the SCH equations. Significantly, the thermoplastic strain rate,Ė p 0 , (also called the "ratchet growth rate") has been identified, which, as can be seen by its dependence on R p and 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 R p → 0, thenĖ p → 0 and M p → 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 R e → 0 (and R p → 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. 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 and temperature θ are given. At t = t (n) , the initial global state θ (n) , E p (n) , E e (n) , E th (n) , (n) , and the set of local states {σ,ε p ,ε th , z} (n) (for all grains) are known. The global thermo-elastoplastic constitutive response E e (n+1) , E th (n+1) , E p (n+1) are required. In solving the constitutive response, we also determine for numerical analysis the global tangents L e and L p and the updated local states {σ,ε p ,ε th , z} (n+1) . On first call, a purely thermoelastic increment is assumed with { γ = 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 Algorithm 1 TEPSCA stress integration algorithm.  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. Algorithm 2 Iterative self-consistent solution of global tangents and local concentration tensors. The (n + 1) subscript to denote current quantities is dropped here for brevity, with the previous (initial) values denoted by the (n) subscript.
1 Given a timestep t with an initial guess for global tangents L e , L p , known/computed local tangents {L e , R e , R p ( γ )}, initial global (n) , E th (n) , and local {z, σ, ε th , T} (n) states, and the shape, volume fraction and orientations of all grains,  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  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 • C/min of an unconfined and stress-free specimen from room temperature (23 • C) up to a maximum temperature of 153 • 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 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 (L e and ε th (θ )) are direct inputs into the model without any need for calibration. Elastic interface-damage parameters α and β 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. 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 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 • C at a rate of 1 • 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 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.  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 θ = 0) Table 3 Polynomial fitting coefficients of the form c 0 + c 1 T + c 2 T 2 , where T is temperature, for single crystal TATB lattice parameter variation fit to the measurements of [31]

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ḣ = 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 int and K 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 τ 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 τ 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 intergranular 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. which we will call EC1.
Green's function, φ 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 , 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, σ kp (x) = L kpim φ ij,m (x − ξ)F j , and the definition of the Dirac delta function, δ(x − ξ). Equilibrium requires that if ξ ∈ 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]), which we will call EC2. We now introduce a perturbed field (fluctuation)ũ = 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ũ(x), recalling that by the chain ruleũ k,jl φ ij = (ũ k,j φ ij ) ,l −ũ k,j φ ji,l , and that the strong form of a field equation implies the weak (integral) form, we find Returning to EC1 (67), multiplying through by φ ij provides Subtracting the augmented forms of EC1 and EC2, that is the final forms of Eqs. (70) and (71), we obtain and finally, by virtue of recognizing that x in Eq. (72) is a bound (dummy) variable (and that Eq. (72) holds for all ξ ∈ V and therefore for all x ∈ V ), we equivalently obtain 2 : where to be clear, we have written out the long hand notation for the partial derivative of φ (because it depends on both x and ξ), 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 (ξ) = ∂ ξ j a i (ξ). Note also that and Making use of the chain rule and divergence theorem, Eq. (73) becomes 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 x ∈ , and the volume explicitly taken as that of the inclusion (V ≡ ). Eq. (76) then becomes ∂φ km (ξ, x) ∂ξ l n j (ξ) dS(ξ) For x ∈ ⊃ S − Hooke's law given by can be used to introduce the stress into the integrand, Next, consider the volume to to be the exterior of the inclusion, i.e., V = D/ . By definition, there is no eigenstrain outside of the inclusion, so 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 inclusioñ for V = .
Assuming the eigenstrain is constant within the inclusion, the non-symmetric total strain field within the inclusion is found wherefore the tensor having major and minor symmetries is defined 3 G ijnm (x, ξ) := 1 4 We can then write the symmetric strain tensor within the inclusion as 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 ( u i = η ijσjk n k +ζ ij σ jk n k ), the integral can be written in terms of the traction, ε ij (x) = S ijklėkl − S L klmn (η kpσpq +ζ kp σ pq )G ijmn (ξ, x)n q n l dS(ξ).
The asymptotic solution for small imperfection is considered, in which caseε (0) = Sė, and consider a perturbation representing slight imperfection.