*2.1. Model Setup*

Figure 2 illustrates the geometries of the simulation model. To investigate the characteristics of the interaction between PSBs and surface coatings, only the tensile loading portion of one fatigue cycle is applied. The intended boundary conditions (BCs) on the lateral coating surfaces are free—i.e., the surface tractions are zero. A rate-controlled uniaxial tension was applied by enforcing only the zz component of the deformation gradient rate tensor (*F*˙ zz), while ensuring that σxx = 0 and σyy = 0. Since this work employed the spectral solver implemented in PETSc, which imposed a full periodic BC on the computational cells [30], the presence of free surfaces was mimicked by adding two soft ~1-μm thick buffers layers on both sides of the sample (see Figure 2a), which is similar to the approach by [4,33,34]. A strong contrast in elastic moduli and strengths existed between the buffer layers and the samples (the elastic constants were at least one order of magnitude lower and the strengths were at least three orders of magnitudes lower for the buffer layer). The effective BCs are therefore free, periodic, and periodic in the x-, y-, and z-directions, respectively.

**Figure 2.** Schematic representation of model geometries: (**a**) a three-dimensional view of the overall geometry of the models omitting the details of the sample, and (**b**) the front view of the simulation cells highlighting the interior geometry and dimensions of the sample.

To conserve computational resources, the y-dimension of the simulation cells was minimized and kept constant at ~0.2 μm. Thus, the models were thin slab-shaped. The width of the substrate along the x-direction was also kept constant at ~34 μm. The PSBs were modeled to be 45◦ off the loading axis (z) and to have a constant thickness of tPSB = 1 μm in accordance with direct experimental observations in the open literature. Indeed, for metallic materials such as Cu, Ni, and SS 316, the thickness of the PSBs was ~1 μm [12,35]. Constrained by the periodic BC, the height of the models along the z-direction was dictated by the thickness (tPSB) of and spacing (dPSB) between the PSBs—i.e., h <sup>=</sup> <sup>√</sup>2 (tPSB <sup>+</sup> dPSB). The spacing dPSB varied between 1 and 8 μm, which corresponds to a PSB volume fraction of 50%~11%, and an overall shear plastic strain range of approximately Δγpl = 0.01~0.002 assuming a 0.01 plastic strain amplitude in the PSBs. Three different coating thicknesses (tcoat), namely 0.5, 1.0, and 2.0 μm, were considered here. These parameters of the models are listed in Table 1.

**Table 1.** Design of simulations performed in the current study. The meaning of tPSB and dPSB are shown in Figure 2b.


Both isotropic and anisotropic plastic flow rules have been utilized in the models. The buffer layers and the coatings were treated by isotropic plasticity. This assumption is sound since the buffer layers only have marginal resistance to deformation and the nanocrystalline Cr coatings exhibit isotropic mechanical behavior. Correspondingly, the plastic flow rule is written as [30]

$$
\dot{\gamma}\_p = \dot{\gamma}\_0 \left(\frac{\sqrt{3f\_2}}{M\xi}\right)^n,\tag{1}
$$

where . <sup>γ</sup>*<sup>p</sup>* is the plastic shear strain rate, . γ<sup>0</sup> is a reference strain rate, *n* is the stress exponent, *J*<sup>2</sup> is the second invariant of the deviatoric stress tensor, and *M* is the Taylor factor. The ξ term in the denominator is the resistance to plastic flow. The rate of ξ is given as

$$\dot{\xi} = \dot{\nu}\_p h\_0 \left| 1 - \frac{\dot{\xi}}{\dot{\xi}\_{\infty}} \right|^a \text{sgn} \left( 1 - \frac{\dot{\xi}}{\dot{\xi}\_{\infty}} \right) \tag{2}$$

where . γ*<sup>p</sup>* is the plastic shear strain rate, *h*<sup>0</sup> is the strain hardening coefficient, ξ<sup>∞</sup> is the saturation resistance to plastic flow, and *a* is a material-dependent exponent.

Anisotropic plasticity was used for the substrate (including both the matrix and the PSBs) and a phenomenological hardening law was used. The flow rule is written as

$$\dot{\gamma}\_p^a = \dot{\gamma}\_0^a \Big| \frac{\tau^a}{\xi^a} \Big|^n \text{sgn}(\tau^a),\tag{3}$$

where <sup>τ</sup><sup>α</sup> is the resolved shear stress on the slip system <sup>α</sup>, ξα is the slip resistance on the slip system, . γα 0 is the reference strain rate, and *n* is the stress exponent. Since only the loading portion of a cyclic loading period was modeled, the back-stress term—which is necessary to capture the kinematic hardening effect in cyclic loading—is not included in the present study. The rate form of the resistance ξα is given as

$$\dot{\boldsymbol{\xi}}^{\alpha} = h\_0 \prescript{\boldsymbol{s} - \boldsymbol{\kappa}}{\boldsymbol{\kappa}} \sum\_{\alpha'=1}^{N\_s} \left| \dot{\boldsymbol{\gamma}}^{\alpha'} \right| \left| 1 - \frac{\boldsymbol{\xi}^{\alpha'}}{\boldsymbol{\xi}\_{\alpha \alpha}^{\alpha'}} \right|^a \text{sgn} \left( 1 - \frac{\boldsymbol{\xi}^{\alpha'}}{\boldsymbol{\xi}\_{\alpha \alpha}^{\alpha'}} \right) \mathbf{h}^{\alpha \alpha'}, \tag{4}$$

where . γα is shear strain rate on the slip system α , <sup>ξ</sup><sup>∞</sup> is the resistance saturation value, *<sup>h</sup>*αα' is the slip hardening matrix (including both self- and latent hardening), and *a* is a material-dependent exponent.

The elastic and plastic flow constants used for the SS 316, PSB, coatings, and buffer layers are summarized in Tables 2 and 3. The stress–strain behaviors produced by the elastic and plastic constants are shown in Figure 3. The elastic constants of both the Cr coating and the SS 316 substrate (including both the PSBs and the matrix) were obtained from the open literature [36–38]. To obtain the plastic flow constants of SS 316, the stress–strain response of a "virtual single crystal"—which was an average of many tensile tests (>50) on single crystals of randomized orientations under the isostrain assumption—was fitted to an experimental curve [39–41]. This technique, also referred to as the "material point" simulation, is a standardized practice to establish flow constants for crystal plasticity simulations [30,42,43].

**Table 2.** Anisotropic elastic and plastic material constants used for SS316 substrate (including PSB).


<sup>1</sup> These values are for the primary slip system in the PSB. The slip activities in the secondary slip systems were suppressed by using much higher slip resistances (at least 1000 times higher).


**Table 3.** Isotropic elastic and plastic material constants used for the buffer layers and coatings.

<sup>1</sup> Although the <sup>ξ</sup>**<sup>0</sup>** and <sup>ξ</sup><sup>∞</sup> appear to be quite low, when combined with the *M* and *n*, they produced a true yield strength of 750 MPa and a strength of 2.5 GPa at 0.5 strain (see Figure 3).

**Figure 3.** The stress–strain responses of the four material types considered in this study.

PSBs within the substrate were modeled as different materials with identical elastic constants and crystallographic orientations but lower shear resistances (~100 MPa in the PSBs compared to ~300 MPa in the matrix) on the primary slip system. As discussed in the introduction, due to the wider dislocation channels within the PSB compared to the matrix, the PSBs have substantially lower critical resolved shear stress (CRSS). The applied overall strains are therefore localized within the PSBs. For pure Cu, Ni, and Ag, the respective CRSSs are ~30, 50, and 20 MPa [44]. As for SS316, the CRSS of PSBs is not known to the authors' best knowledge and must be assumed. Considering the solid solution strengthening effect in SS316, the highest known CRSSs among the three aforementioned elemental metals, i.e., 50 MPa of Ni, was used. Assuming a Schmid factor of 0.5 (which applies for the current model geometries), the corresponding tensile yield strength is 100 MPa. As will be shown in Section 3, this choice of the CRSS appears to be sufficient to capture the strain localization within the PSBs. The slip activity on the secondary systems was completely suppressed by applying much higher critical resolved shear stresses (~40 GPa). The plastic flow constants of the coating and the buffer layer were calibrated so that they reproduce yield strengths of 750 and ~0 MPa, respectively.

Note that the crystallographic orientation of the substrate must be defined carefully, so that the primary slip system in the PSB experiences maximum shear stress under the tensile loading applied (see Figure 2b). In other words, the primary slip direction and slip plane should both be 45◦ off the loading axis. A cubic grid with the characteristic size ~0.09 × ~0.09 × ~0.09 μm has been chosen for all the models, leading to 2 Fourier points (FPs) along the y-direction, 400 to 440 (FP) along the x-direction, and 30 to 135 (FP) along the z-direction. The one-FP per voxel configuration is comparable to the one-integration point, C3D8R finite element type (according to ABAQUS), which was emulated in the model setup by DAMASK [45].
