A multiphase phase-field study of three-dimensional martensitic twinned microstructures at large strains

Abstract

The nanoscale multiphase phase-field model for stress and temperature-induced multivariant martensitic transformation under large strains developed by the authors in Basak and Levitas (J Mech Phys Solids 113:162–196, 2018) is revisited, the issues related to the gradient energy and coupled kinetic equations for the order parameters are resolved, and a thermodynamically consistent non-contradictory model for the same purpose is developed in this paper. The model considers $N+1$ order parameters to describe austenite and N martensitic variants. One of the order parameters describes austenite$\leftrightarrow$martensite transformations, and the remaining N order parameters, whose summation is constrained to the unity, describe the transformations between the variants. A non-contradictory gradient energy is used within the free energy of the system to account for the energies of the interfaces. In addition, a kinetic relationship for the rate of the order parameters versus thermodynamic driving forces is suggested, which leads to a system of consistent coupled Ginzburg–Landau equations for the order parameters. An approximate general crystallographic solution for twins within twins is presented, and the explicit solution for the cubic to tetragonal transformations is derived. A large strain-based finite element method is developed for solving the coupled Ginzburg–Landau and elasticity equations, and it is used to simulate a 3D complex twins within twins microstructure. A comparative study between the crystallographic solution and the simulation results is presented.

Appendices

Appendix

A Thermodynamic formalism and dissipation inequality

We use the general thermodynamically consistent nonlocal framework developed by the authors in [20, 73] for deriving the governing coupled mechanics phase-field equations. For the completeness of the presentation, we briefly show the derivation of the dissipation inequality, which is used for deriving the Ginzburg–Landau equations for the order parameters and the constitutive relations in Sect. 2. We begin with the following dissipation inequalities obtained using the first and second laws of thermodynamics in $\Omega _0$ [20, 73]:

$$\begin{aligned} \rho _0{{\mathcal {D}}}= & {} {\varvec{P}}:\dot{\varvec{F}}^T-\rho _0{\dot{\psi }}-\rho _0 s {\dot{\theta }}+\nabla _0\cdot (\varvec{Q}_0^\eta {\dot{\eta }}_0) +\sum _{i=1}^N\nabla _0\cdot (\varvec{Q}_i^\eta {\dot{\eta }}_i) \ge 0, \quad \text {and} -\frac{1}{\theta } \varvec{h}_0\cdot \nabla _0\theta \ge 0, \end{aligned}$$

(A.1)

where ${\mathcal {D}}$ is the power dissipation per unit mass. The symbol ${\varvec{P}}$ denotes the total first Piola–Kirchhoff stress tensor, $\psi$ is the specific (per unit mass) Helmholtz free energy, $\rho _0$ denotes the mass density of $V_0$, $\theta >0$ is the absolute temperature, s is the specific entropy, $\varvec{Q}_l^\eta$ ($l=1,\ldots ,N$) are the generalized force vectors in $\Omega _0$ introduced for the balance of some of the terms in the inequality (A.1) [20, 73], and $\varvec{h}_0$ is the heat flux vector. Based on inequality (A.1)$_2$, we write the Fourier’s law $\varvec{h}_0=-\varvec{K}_\theta \cdot \nabla _0\theta$, where $\varvec{K}_\theta$ is the heat conductivity tensor which is symmetric and positive semi-definite. If the body is at a uniform temperature, as assumed for the numerical simulation in Sect. 4.2, $\varvec{h}_0=\varvec{0}$.

The material time derivative of the system free energy given by in Eqs. (2.6) and (2.7) is obtained as

$$\begin{aligned} {\dot{\psi }}= & {} \frac{J_t}{\rho _0}\frac{\partial \psi ^e}{\partial \varvec{F}_e}\cdot \varvec{F}_t^{-T}:\dot{\varvec{F}}^T -\frac{J_t}{\rho _0}\varvec{F}_e^T\cdot \frac{\partial \psi ^e}{\partial \varvec{F}_e}\cdot \varvec{F}_t^{-T}:\dot{\varvec{F}}_t^T + \frac{J_t\psi ^e}{\rho _0}\varvec{F}_t^{-T}:\dot{\varvec{F}}_t^T +J(\breve{\psi }^{\theta }+\psi ^\nabla ){\varvec{F}}^{-T}:\dot{\varvec{F}}^T \nonumber \\{} & {} +\sum _{i=0}^N\left( J\frac{\partial (\breve{\psi }^{\theta }+\psi ^\nabla )}{\partial \eta _i}+\frac{\partial ({\tilde{\psi }}^\theta +\psi ^p)}{\partial \eta _i}+\frac{J_t}{\rho _0}\left. \frac{\partial \psi ^e}{\partial \eta _i}\right| _{\varvec{F}_e}\right) {\dot{\eta }}_i +J\sum _{i=0}^N\frac{\partial \psi ^\nabla }{\partial \nabla \eta _i}\cdot \dot{\overline{\nabla \eta _i}} + \frac{\partial \psi }{\partial \theta }{\dot{\theta }}, \nonumber \\ \end{aligned}$$

(A.2)

where we have used $\dot{\overline{det\,\varvec{T}}}=({det\,\varvec{T}}) \varvec{T}^{-1}:\dot{\varvec{T}}$, and $\dot{\overline{\varvec{T}^{-1}}}=-\varvec{T}^{-1}\cdot \dot{\varvec{T}}\cdot \varvec{T}^{-1}$ for an arbitrary invertible second-order tensor $\varvec{T}(t)$. Using the relation

$$\begin{aligned} \sum _{l=0}^N\frac{\partial \psi ^\nabla }{\partial \nabla \eta _l}\cdot \dot{\overline{\nabla \eta _l}} = \sum _{l=0}^N\frac{\partial \psi ^\nabla }{\partial \nabla \eta _l}\cdot \dot{\overline{(\varvec{F}^{-T}\cdot {\nabla _0\eta _l})}} =\sum _{l=0}^N\left( \varvec{F}^{-1}\cdot \frac{\partial \psi ^\nabla }{\partial \nabla \eta _l}\cdot \dot{\overline{ \nabla _0\eta _l}} -\nabla \eta _l\otimes \varvec{F}^{-1}\cdot \frac{\partial \psi ^\nabla }{\partial \nabla \eta _l}:\dot{\varvec{F}}^T \right) , \nonumber \\ \end{aligned}$$

(A.3)

where we have used $\nabla \eta _l=\varvec{F}^{-T}\cdot \nabla _0\eta _l$, and noticing that $\varvec{F}_t$ is a function of all the order parameters [see Eq. (2.5)], we rewrite Eq. (A.2) as

$$\begin{aligned} {\dot{\psi }}= & {} \left( \frac{J_t}{\rho _0}\frac{\partial \psi ^e}{\partial \varvec{F}_e}\cdot \varvec{F}_t^{-T} +J(\breve{\psi }^{\theta }+\psi ^\nabla ){\varvec{F}}^{-T} -J \sum _{l=0}^N\nabla \eta _l \otimes \varvec{F}^{-1}\cdot \frac{\partial \psi ^\nabla }{\partial \nabla \eta _l}\right) : \dot{\varvec{F}}^T \nonumber \\{} & {} \quad +\sum _{l=0}^N J\varvec{F}^{-1}\cdot \frac{\partial \psi ^\nabla }{\partial \nabla \eta _l} \cdot \dot{\overline{ \nabla _0\eta }}+\frac{\partial \psi }{\partial \theta }{\dot{\theta }} + \nonumber \\ {}{} & {} \sum _{l=0}^N\left( J\frac{\partial (\breve{\psi }^{\theta }+\psi ^\nabla )}{\partial \eta _l}+ \frac{\partial ({{\tilde{\psi }}}^\theta +\psi ^p)}{\partial \eta _l} +\frac{J_t}{\rho _0}\left. \frac{\partial \psi ^e}{\partial \eta _l}\right| _{\varvec{F}_e}\right. \nonumber \\{} & {} \left. + \left( \frac{J_t\psi ^e}{\rho _0}\varvec{F}_t^{-T} - \frac{J_t}{\rho _0}\varvec{F}_e^T\cdot \frac{\partial \psi ^e}{\partial \varvec{F}_e}\cdot \varvec{F}_t^{-T} \right) :\frac{\partial {\varvec{F}}_t^T}{\partial \eta _l}\right) {\dot{\eta }}_i. \end{aligned}$$

(A.4)

We apply the Coleman–Noll procedure [77] (see [78] for its generalization to a nonlocal theory) to obtain the constitute relations. Using Eq. (A.4) in the inequality (A.1)$_1$ and assuming that the dissipation rate is independent of ${\dot{\theta }}$, and $\dot{\overline{ \nabla _0\eta _l}}$ (for all $l=0,1,\ldots ,N$), and neglecting the viscous dissipation within the solid, we get the constitutive relations for entropy $s=-\displaystyle \frac{\partial \psi }{\partial \theta }$, generalized force vector $\varvec{Q}_{l}^\eta = \rho _0J\varvec{F}^{-1}\cdot \displaystyle \frac{\partial \psi ^\nabla }{\partial \nabla \eta _l}$ for all $l=0,1,\ldots ,N$ [77, 78], and the total first Piola–Kirchhoff stress tensor ${\varvec{P}}$ is given by Eq. (2.10)$_1$. Considering these relations, the dissipation inequality (A.1)$_1$ finally simplifies to (2.17) given in Sect. 2.4.2.

B Finite element procedure and numerical implementation

In this appendix, we briefly describe our FE formulation and derive the system of algebraic equations which are solved to compute the nodal displacements and order parameters in our 3D domain. We discretize the volume in the reference body into $n_{el}$ finite elements, i.e., $V_0\approx \cup _{el=1}^{n_{el}}V_0^{el}$. The total number of grid points is $n_\textrm{grid}$. The boundary $S_0$ consists of the edges and areas of some elements. We have developed the FE procedure for solving the couple mechanics, and phase-field equations using a non-monolithic scheme, where the order parameters are assumed to remain fixed while solving the mechanics equations and the displacements are assumed to remain fixed while solving the phase-field equations (see [82] for details).

FE procedure for mechanics problem. Writing the weak form of the equilibrium equation Eq. (2.9)$_1$, we linearize it and discretize the linearized equation using the standard procedure [91] to obtain the following system of algebraic equations for computing the increment in the nodal displacements (see [82] for derivation)

$$\begin{aligned} \varvec{K}\cdot \Delta \varvec{u}^p = -{\varvec{r}}_u, \quad \text {where} \end{aligned}$$

(B.1)

$\varvec{K}$ is the $3\,n_\textrm{grid}\times 3\,n_\textrm{grid}$ global tangent stiffness matrix which is symmetric, $\varvec{r}_u$ is the $3\,n_\textrm{grid}\times 1$ global residual matrix, and $\Delta \varvec{u}^p$ is the $3\,n_\textrm{grid}\times 1$ incremental displacement matrix at the p-th iteration. The global matrix $\varvec{K}$ is obtained by using the standard assembly operation on the elemental tangent matrix, which is obtained using the following $3\times 3$ matrix related to the nodes $\iota$ and $\kappa$ of each element (index as el)

$$\begin{aligned} {\varvec{ K}}^{el}_{\iota \kappa }= \int _{V_0^{el}}\left( {{\varvec{B}}_{\iota }^{el}}^T\cdot {{{\textbf {{\textsf { C}}}}}}^{el}\cdot {\varvec{B}}_{\kappa }^{el}+G_{\iota \kappa }{\varvec{I}}\right) \, dV_0^{el}. \end{aligned}$$

(B.2)

The size of the elemental matrix $\varvec{K}^{el}$ is $3\,n_g\times 3\,n_g$, where $n_g$ is the number of nodes in an element. The global $\varvec{r}_u$ matrix is obtained by assembling the elemental residual matrix obtained using the following $3\times 1$ matrix related to node $\iota$ of each element

$$\begin{aligned} ({\varvec{r}_\iota ^{el}})_u = \int _{V_0^{el}}J^{el}{\varvec{B}_{\iota }^{el}}^T\cdot \varvec{\sigma }^{el}\, dV_0^{el}. \end{aligned}$$

(B.3)

The size of the elemental ${\varvec{r}}^{el}$ matrix is $3\,n_g\times 1$. Here, we enlist the other symbols in Eq. (B.2). ${{{\textbf {{\textsf {C}}}}}}^{el}= {{{\textbf {{\textsf {C}}}}}}^{el}_e+{{{\textbf {{\textsf {C}}}}}}^{el}_{st}$ is the $6\times 6$ elemental stiffness matrix and ${{\textbf {{\textsf { C}}}}}^{el}_e$ and ${{{\textbf {{\textsf { C}}}}}}^{el}_{st}$ are the $6\times 6$ elemental elastic and structural stiffness matrices, respectively. Their components are obtained from the following two fourth-order tensors in $\Omega$ (see [82] for derivation):

$$\begin{aligned} \mathsf{{C}}_{(e)abcd}= & {} J_e^{-1} F_{(e)ap}F_{(e)bq}F_{(e)cr}F_{(e)ds}\hat{{\mathcal {C}}}_{(e)pqrs}, \quad \text {and} \nonumber \\ \mathsf{{C}}_{(st)abcd}= & {} \sigma _{(st)ab}\delta _{cd}+\delta _{ab}\sigma _{(st)cd}-(\sigma _{(st)ac}\delta _{bd}+\sigma _{(st)ad} \delta _{bc}+\delta _{ac}\sigma _{(st)bd}+\delta _{ad}\sigma _{(st)bc})\nonumber \\{} & {} +\rho _0J(\breve{\psi }^\theta +\psi ^\nabla )(\delta _{ac}\delta _{bd}+\delta _{ad}\delta _{bc}-\delta _{ab}\delta _{cd}), \end{aligned}$$

(B.4)

respectively, and using the Voigt representation. In Eq. (B.4)$_1$, $\displaystyle \hat{\varvec{{\mathcal {C}}}}_{(e)}:=\frac{\partial ^2\psi ^e}{\partial \varvec{E}_e\partial \varvec{E}_e}$ is the fourth-order elasticity tensor defined in $\Omega _t$ [82]. We denote the shape functions for the $\iota$th grid point in element $V_0^{el}$ as $N_\iota$ for $\iota =1,\ldots , n_g$. The standard $\nabla N_\iota$ and $\varvec{B}_\iota ^{el}$ matrices of sizes $3\times 1$ and $6\times 3$, respectively, are given by [91]

$$\begin{aligned} \nabla N_\iota= & {} \begin{bmatrix} N_{\iota , 1} \\ N_{\iota , 2} \\ N_{\iota , 3} \\ \end{bmatrix}; \quad \text {and} \quad \varvec{B}_{\iota }^{el} = \begin{bmatrix} N_{\iota , 1}&{} 0&{} 0 \\ 0 &{} N_{\iota , 2} &{} 0 \\ 0 &{} 0 &{} N_{\iota , 3} \\ N_{\iota , 2} &{} N_{\iota , 1} &{} 0 \\ 0 &{} N_{\iota , 3} &{} N_{\iota , 2} \\ N_{\iota , 3} &{} 0 &{} N_{\iota , 1} \\ \end{bmatrix}, \end{aligned}$$

(B.5)

where the comma followed by the number designates the derivative with respect to $r_i$ (for $i=1,2,3$) in $\Omega$. In Eq. (B.3), $\varvec{\sigma }^{el}=\{\sigma _{11},\,\sigma _{22},\,\sigma _{33},\,\sigma _{12},\,\sigma _{23},\,\sigma _{13} \}^T$ is the $6\times 1$ matrix whose elements are the components of the total stress $\varvec{\sigma }$ given by Eq. (2.13)$_{1}$. In Eq. (B.2)$_1$, $G_{\iota \kappa }{\varvec{I}} = J^{el}\,(\nabla N_\iota ^T\cdot \varvec{\sigma }\cdot \nabla N_\kappa )\,{\varvec{I}}$ is the geometric part of the tangent matrix $\varvec{K}$, and here note that $\varvec{\sigma }$ is a $3\times 3$ matrix (see Chapter 3 of [91]). Finally, solving Eq. (B.1) iteratively using Newton’s method, the $3n_{grid}\times 1$ nodal displacement matrix is updated after the $p^{th}$ iteration using

$$\begin{aligned} \varvec{u}^p = \varvec{u}^{p-1}+\Delta \varvec{u}^p. \end{aligned}$$

(B.6)

FE procedure for Ginzburg–Landau equations. We discretize the time rate of the order parameters as

$$\begin{aligned} {\dot{\eta }}_l = (c_1\eta _l^n+c_2\eta _l^{n-1}+c_3\eta _l^{n-2}) /\Delta t^n \quad \text { for all }l=0,1,2,\ldots ,N, \end{aligned}$$

(B.7)

where the constants $c_1$, $c_2$, and $c_3$ take the values 1, $-1$, and 0, respectively, for BDF scheme of order one, and their values are 1.5, $-2$, and 0.5, respectively, for BDF scheme of order two. The time step size $\Delta t^n$ at the nth iteration is chosen which yields a converged solution for the N order parameters. The superscript n designates the time index, and the time instance after $(n-1)$th iteration is given by $t^n=t^{n-1}+\Delta t^n$. Using Eq. (B.7) in the Ginzburg–Landau equations (2.23) and (2.29), we write down their weak forms, linearize the weak forms and discretize them and finally, do the assembly operation using the standard procedure [91] to obtain the system of algebraic equations for the N independent order parameter $\eta _l$ (say $\eta _0,\eta _1,\ldots ,\eta _{N-1})$:

$$\begin{aligned} {\varvec{T}_l} \cdot \Delta \varvec{\eta }_l^{n,q}=-\varvec{r}_l\quad \text {for }l=0,1,\ldots , N-1. \end{aligned}$$

(B.8)

In Eq. (B.8), $\varvec{T}_l$ is $n_\textrm{grid}\times n_{grid}$ symmetric global matrix related to $\eta _l$ and it is given by

$$\begin{aligned} \varvec{T}_l(\eta ^{n,q-1}_k) = c_1\varvec{M}_l +\Delta t^n\varvec{H}_l+\Delta t^n\varvec{G}_l. \end{aligned}$$

(B.9)

The global mass matrix $\varvec{M}_l$, global Laplace matrix $\varvec{H}_l$, and global nonlinear matrix $\varvec{G}_l$ are obtained by applying the standard assembly operation of the corresponding $n_g\times n_g$ elemental mass, Laplace, and nonlinear matrices, respectively. The $\iota \kappa$ components of the elemental mass, Laplace, and nonlinear matrices are given by

$$\begin{aligned} (\varvec{M}_l^{el})_{\iota \kappa }= & {} \int _{V_0^{el}} N_\iota N_\kappa \,dV_0^{el}, \nonumber \\ (\varvec{H}_l^{el})_{\iota \kappa }= & {} \int _{V_0^{el}}h_l^{n,q-1} {L_l^\beta }J^n(\nabla N_\iota \cdot \nabla N_\kappa ) dV_0^{el},\nonumber \\ (\varvec{G}_l^{el})_{\iota \kappa }= & {} \int _{V_0^{el}}\left. \frac{\partial f_l^n(\eta _l^{n,q-1})}{\partial \eta _l^n}\right| _{\varvec{F}} N_\iota N_\kappa dV_0^{el}, \end{aligned}$$

(B.10)

respectively, where $h_0^{n,q-1}=1$, $L^\beta _0=L_{0M}\beta _{0M}$, $h_l^{n,q-1}={\tilde{\varphi }}(a_\beta ,a_c,\eta _0^{n,q-1})$, and $L^\beta _l=\sum _{k=1,\ne l}^NL_{lk}\beta _{lk}$ for $l=1,\ldots ,N$. The expressions for $f_l^n$ for $l=0$ and $l=1,\ldots ,N$ are given by Eqs. (B.13) and (B.14), respectively, and $n_\textrm{grid}$ is the total number of degrees of freedom for $\eta _l$. In Eq. (B.8), $\Delta \varvec{\eta }_l^{n,q}$ is $n_\textrm{grid}\times 1$ matrix for the increment of $\eta _l$ at $q^{th}$ Newton iteration in $n^{th}$ time step, and $\varvec{r}_l$ is $n_\textrm{grid}\times 1$ residual matrix given by

$$\begin{aligned} \varvec{r}_l = (c_1\varvec{M}_l+\Delta t^n\varvec{H}_l)\cdot \varvec{\eta }_l^{n,q-1} + c_2\varvec{M}_l\cdot \varvec{\eta }_l^{n-1} + c_3\varvec{M}_l\cdot \varvec{\eta }_l^{ n-2}+\Delta t^n\varvec{f}_l, \end{aligned}$$

(B.11)

where $n_\textrm{grid}\times 1$ global matrix $\varvec{f}_l$ is obtained by assembling the $n_g\times 1$ elemental matrices whose $\iota$ component is

$$\begin{aligned} (\varvec{f}_l^{el})_\iota = \int _{V_0^{el}} f_l^n N_\iota dV_0^{el} \end{aligned}$$

(B.12)

for all $l=0,1,\ldots ,N-1$. The order parameters after each iteration are updated using

$$\begin{aligned} {\varvec{\eta }_l^{n,q}}={\varvec{\eta }_l^{n,{q-1}}}+\Delta {\varvec{\eta }_l^{n,q}} \quad \text {for }l=0,1,\ldots , N-1. \end{aligned}$$

The nodal matrix $\varvec{\eta }_N^{n,q}$ is then obtained using the constraint Eq. (2.1).

Following the procedure of [82], the expressions for $f_0$ and $f_i$ for $i=1,2,\ldots ,N$ appearing in Eqs. (B.10)$_3$ and (B.12) are obtained using Eqs. (2.54), (2.55), (2.56), and (2.57) as

$$\begin{aligned} f_0^n= & {} L_{0M}\left[ -\left( J^n {\varvec{F}^n}^{-1}\cdot \varvec{\sigma }_e^n\cdot \varvec{F}^n-J_t^n{\psi ^e}^n{\varvec{I}}\right) :\varvec{Y}_0+ {J}_t^n{\left. \frac{\partial {\psi ^e}^n}{\partial \eta _0}\right| _{\varvec{F}_e}}+ \rho _0\frac{\partial \varphi (a_\theta ,\eta _0^n)}{\partial \eta _0}\Delta \psi ^\theta + \right. J^n\rho _0{\bar{A}} \nonumber \\{} & {} \times \sum _{i=1}^{N-1}\sum _{j=i+1}^{N}{\eta _i^n}^2 {\eta _j^n}^2 \frac{\partial \varphi (a_b,\eta _0^n)}{\partial \eta _0} +J^n\rho _0 A_{0M}(\theta ) (2\eta _0^n- 6{\eta _0^n}^2+4{\eta _0^n}^3) + \frac{J^n}{2}\frac{\partial {\tilde{\varphi }}(a_\beta ,a_c,\eta _0^n)}{\partial \eta _0} \nonumber \\{} & {} \times \sum _{i=1}^{N-1}\sum _{j=i+1}^{N} \beta _{ij}\nabla \eta _i^n\cdot \nabla \eta _j^n +\rho _0\left( \sum _{i=1}^{N-1}\sum _{j=i+1}^N K_{0ij}{\eta _i^n}^2{\eta _j^n}^2+ \sum _{i=1}^{N-2}\sum _{j=i+1}^{N-1}\sum _{k=j+1}^N K_{0ijk}{\eta _i^n}^2{\eta _j^n}^2{\eta _k^n}^2\right) \nonumber \\{} & {} \times \left( 2(1-\varphi (a_K,\eta _0^n))\eta _0^n-\frac{\partial \varphi (a_K,\eta _0^n)}{\partial \eta _0}{\eta _0^n}^2\right) + \rho _0\frac{\partial \varphi (a_K, \eta _0^n)}{\partial \eta _0}\left( \sum _{i=1}^{N-1}\sum _{j=i+1}^N K_{ij}( \eta _i^n+\eta _j^n-1)^2 {\eta _i^n}^2{\eta _j^n}^2 \right. \nonumber \\{} & {} \left. \left. + \sum _{i=1}^{N-2}\sum _{j=i+1}^{N-1}\sum _{k=j+1}^{N}K_{ijk} {\eta _i^n}^2{\eta _j^n}^2{\eta _k^n}^2 + \sum _{i=1}^{N-3}\sum _{j=i+1}^{N-2}\sum _{k=j+1}^{N-1}\sum _{l=k+1}^N K_{ijkl} {\eta _i^n}^2{\eta _j^n}^2{\eta _k^n}^2{\eta _l^n}^2\right) \right] , \quad \text {and} \end{aligned}$$

(B.13)

$$\begin{aligned}{} & {} f_i^n({{\tilde{\eta }}}^n,({{\tilde{\eta }}}^\nabla )^n) = -\sum _{m=1,\ne i}^{N} L_{im} \left( {X}_{im}^{loc}+ {X}_{im}^{\nabla (1)}\right) \quad \text {for all }i=1,\ldots ,N, \quad \text {where} \end{aligned}$$

(B.14)

$$\begin{aligned} {X}_{im}^{loc}= & {} (J^n {\varvec{F}^n}^{-1}\cdot \varvec{\sigma }_e^n\cdot \varvec{F}^n-J_t^n{\psi ^e}^n{\varvec{I}}):(\varvec{Y}_i^n-\varvec{Y}_m^n)- {J}_t^n\left( {\left. \frac{\partial {\psi ^e}^n}{\partial \eta _i}\right| _{{\varvec{F}}_e}} - {\left. \frac{\partial {\psi ^e}^n}{\partial \eta _m}\right| _{{\varvec{F}}_e}}\right) -2J^n\rho _0{\bar{A}}\left( \sum _{j=1,\ne i}^{N}{\eta _j^n}^2\eta _i^n \right. \nonumber \\{} & {} \left. -\sum _{j=1,\ne m}^{N}{\eta _j^n}^2\eta _m^n\right) \varphi (a_b,\eta _0^n)- 2\rho _0\varphi (a_K,\eta _0^n) \left[ \sum _{j=1,\ne i}^{N} K_{ij}(\eta _i^n+\eta _j^n-1) (2\eta _i^n+\eta _j^n-1)\right. \eta _i^n{\eta _j^n}^2 \nonumber \\{} & {} \left. -\sum _{j=1,\ne m}^{N} K_{mj} (\eta _m^n+\eta _j^n-1)(2\eta _m^n+\eta _j^n-1)\eta _m^n{\eta _j^n}^2 \right] -2\rho _0\left[ \sum _{j=1,\ne i}^{N}K_{0ij}\eta _i^n{\eta _j^n}^2-\sum _{j=1,\ne m}^{N}K_{0mj}\eta _m^n {\eta _j^n}^2 \right. \nonumber \\{} & {} \left. + \sum _{j=1, \ne i}^{N-1}\sum _{k=j+1, \ne i}^{N} K_{0ijk}\eta _i^n {\eta _j^n}^2{\eta _k^n}^2- \sum _{j=1, \ne m}^{N-1}\sum _{k=j+1, \ne m}^{N} K_{0mjk}\eta _m^n{\eta _j^n}^2 {\eta _k^n}^2 \right] {\eta _0^n}^2(1-\varphi (a_K,\eta _0^n)) \nonumber \\{} & {} -2\rho _0\varphi (a_K,\eta _0^n) \left[ \sum _{j=1,\ne i}^{N-1}\sum _{k=j+1, \ne i}^{N}K_{ijk}\eta _i^n{\eta _j^n}^2{\eta _k^n}^2 + \sum _{j=1, \ne i}^{N-2}\sum _{k=j+1, \ne i}^{N-1}\sum _{l=k+1, \ne i}^{N} K_{ijkl}\eta _i^n{\eta _j^n}^2{\eta _k^n}^2{\eta _l^n}^2 \right. \nonumber \\{} & {} -\left. \sum _{j=1,\ne m}^{N-1}\sum _{k=j+1, \ne m}^{N} K_{mjk}\eta _m^n{\eta _j^n}^2{\eta _k^n}^2-\sum _{j=1, \ne m}^{N-2}\sum _{k=j+1, \ne m}^{N-1}\sum _{l=k+1, \ne m}^{N} K_{mjkl}\eta _m^n{\eta _j^n}^2{\eta _k^n}^2{\eta _l^n}^2 \right] \nonumber \\{} & {} \quad \text {for all }i,m=1,2,\ldots ,N, \end{aligned}$$

(B.15)

$$\begin{aligned} {X}_{im}^\nabla= & {} \nabla \cdot \left[ \frac{{\tilde{\varphi }}}{2} \left( \sum _{j=1,\ne i}^{N}\beta _{ij}\nabla \eta _j-\sum _{j=1,\ne m}^{N}\beta _{mj}\nabla \eta _j\right) \right] ={X}_{im}^{\nabla (1)}+{X}_{im}^{\nabla (2)}, \end{aligned}$$

(B.16)

$$\begin{aligned} {X}_{im}^{\nabla (2)}= & {} \nabla \cdot \left[ \frac{{\tilde{\varphi }}}{2} \left( \sum _{j=1,\ne i}^{N}\beta _{ij}\nabla \eta _j-\sum _{j=1,\ne i,m}^{N}\beta _{mj}\nabla \eta _j\right) \right] , \quad \text {and} \nonumber \\ {X}_{im}^{\nabla (1)}= & {} - \nabla \cdot \left[ \frac{{\tilde{\varphi }}}{2} \beta _{im} \nabla \eta _i \right] . \end{aligned}$$

(B.17)

In Eqs. (B.13) and (B.15), $\varvec{Y}_l^n=\displaystyle {{{\varvec{F}_t^n}}^{-1}}\cdot {\frac{\partial \varvec{F}_t^n}{\partial \eta _l^n}}$ for all $l=0,1,\ldots ,N$. Following the procedure of [91], the expressions for $\partial f_0^n/\partial \eta _0^n$ and $\partial f_i^n/\partial \eta _i^n$ appearing in Eq. (B.10)$_3$ are derived using Eqs. (B.13) and (B.14).