The frequency of cyclic loading has a strong effect on the fatigue lifetime of materials. High loading frequencies allow reaching very high cycle fatigue regimes. The effects of cyclic loading with different frequencies are studied using the finite element method within the framework of crystal plasticity. The numerical study represents copper polycrystal that undergoes tension -compression loading with frequencies 80 Hz and 20 kHz respectively. The spatial distributions o f stresses, strains, and plastic slip are used for estimation of volume fraction, morphology, and plasticity mechanism of persistent slip bands. The results show that at 80 Hz cyclic loading, plasticity in persistent slip band cells is dominated by the cre ation and annihilation of dislocations while at 20 kHz cyclic loading, plasticity is driven by the back -and-forth movement of dislocations inside the persistent slip band cells. The numerical study shows a correlation between amplitude and frequency of loa ding with spatial distribution and morphology of persistent slip bands that can be related to the fatigue lifetime. Keywords: Finite elements ; Crystal plasticity ; Copper ; VHCF ; Slip bands 1. Introduction The new applications require materials with enhanced fatigue lifetime which has to be extended to the range of 1010 - 1011 cycles. Such a large number of cycles enters the field of very high cycle fatigue (VHCF) which is usually combined with high -frequency loading. The combination of VHCF and high -frequency loading brings new phenomena as the material response is different at these conditions compared to low/middle -frequency loading conventionally used for low and high cycle fatigue regions. The extension of fatigue lifetime by several orders of mag nitude suggests a change in the damage mechanism [1][2][3]. Fatigue damage is initiated in the locations of localized plastic deformation [4]. In the case of FCC metals, plastic strain localization leads to the creation of persistent slip bands (PSBs) which have, in the case of copper, ladder -like dislocation struct ures with high dislocation density regions called walls and low dislocation density regions called channels [5][6][7]. The PSBs are the initiation sites for fatigue cracks. The most prominent feature of PSBs are slip marks (SMs) that are created at the free surface of PSBs [8]. Phung et al. [9] identified for copper three different types of SMs related to the different fatigue regions. Low cycle fatigue (below 106 cycles) was characterized by long and straight SMs occasionally occurring in families created by parallel markings (type I SMs). The Type I SMs can cross the grain boundaries. The high cycle fatigue region (106 - 108) was also characterized by long and straight SMs occasionally occurring in families, however, type II SMs do not cross the grain boundaries and are located mostly along the grain boundaries and even along the twin boundaries. The type III SMs preferentiall y occurring in the VHCF region (108 - 1011) had a totally different character. Type III SMs were created by a cluster of very fine markings formed on the grain boundary or in the grains. Phung et al. [10] performed numerical simulations using anisotropic elasticity and showed the stress concentrations at the grain boundaries of copper polycrystalline aggregate explaining the creation of SMs pr eferentially at the grain boundaries. Fintová et al. [11] performed a similar study on polycrystalline copper comparing the behavior of the material in the high -cycle and very high -cycle fatigue regions enhancing the knowledge about the description of the acting fatigue damage mechanism. The study confirmed the localization of cyclic plastic deformation and the creation of SMs on the sample surface, corresponding to the test frequency and applied stress amplitude. Three different types of SMs were observed to be created on samples during cycling, differing depending on the test conditions. Comparable results, however, in less extent can be found in the work of Weidne r et al. [12], showing the localization of the cyclic plastic deformation in the pure copper on the surface and below it, resulting in the formation of cavities an d primary fatigue cracks. The study was, however, limited only to the very high -frequency testing (19 kHz) and thus does not provide a comparison with the material damage under the low/middle -frequency loading. The analysis of PSBs and SMs evolution is cha llenging due to the complexity of their plastic processes. The first approach was done by the phenomenological models describing the evolution of plasticity and lattice defects inside PSBs using different assumptions [13]: i) surface stress -assisted models take into account stress created by dislocations dipoles in the longitudinal axis of PSBs [14][15]; ii) micromechanical models are based on the description of active slip systems under local loading conditions [16][17]; and iii) vacancy models describe the distribution of vacancies due to the creation of dislocation dipoles inside PSBs [18]. Polak’s [19] model belongs to the vacancy category and describes the transfer of mass by vacancy migration from PSB channels to PSB walls. Such transfer creates SM in the form of intrusions and extrusion at the surface where the fatigue crack initiates. The numerical models provide more detailed insight into the evolution of PSBs and SMs. The dislocation structures and their evolution are best described by discrete dislocation dynamics (DDD) models [20][21][22]. However, their computational costs l imit them to analysis of individual PSB under simplified loading conditions. The continuum models combined with the FEM method are able to describe the spatial evolution of stress and strains during cyclic loading and identify the stress/strain concentrati on locations in different materials such as copper [23][24], steel [25] or nickel superalloys [26][27]. The limitation of the continuum approach is in the inability of the inherent description of a ladder -like structure. Such a drawback can be overcome by a -priori defined PSB structures [28] or by a combination of continuum models with phenomenological ones. Sinha et al. [29] combined 3D FEM simulation with a ratcheting model to predict the fatigue lifetime of HSLA steel. Recently Fourel et al. [30] applied the Tanaka model [17] to predict the location and length of PSBs in steel polycrystals during contact fatigue. The presented s tudy uses copper as the model material, which is a typical representat ive of ductile, single phase material without extrinsic internal defects (type I material according to Mughrabi et al. [31]). Therefore, the development of plasticity is based on dislocations interactions. Such material can be very well described by the crystal plasticity model. The chosen approach in this study combines 3D FEM crystal plasticity simulations with a phenomenological model of PSB derived by Essmann and Mughrabi [32]. The crystal plasticity FE model is able to describe in much more details the stress distribution and activity of slip systems with respect to the microstructural features under cyclic loading with different amplitude and frequency compared to classical von Mises plasticity . The output of the FE model is consequently used in the PSB model to analyze the dislocation mechanisms and their relation to the experimentally observed SMs. Such a combination of models allows to analyze the evolution of plasticity and PSBs in polycrystalline aggregate with out excessive computational costs. The analysis reveal s the effects of loading amplitude and frequency on the evolution of plastic slip and consequently the relation between the loading conditions and the spatial distribution and morphology of PSBs in polycr ystal with a given texture. It can also help to identify the dislocation interaction mechanisms inside PSBs for given loading conditions. The simulated polycrystalline aggregate also allows the comparison with the experimental results. Such a comparison is done with the experimental work on copper polycrystals performed by Fintova et al. [11]. 2. Methods 2.1. Experiments The experiments directly relevant to the present ed numerical study include the description of material microstructure (grain size and crystallographic orient ation) and observation of slip marks at the fatigue test specimens ’ surface. These analyses were performed using the Tescan LYRA 3 XMU FEG/SEMxFIB scanning electron microscope (SEM) combined with AZtec software . The data of tensile properties at different strain rates were obtained using the universal Zwick Roell Z50 screw - driven testing machine (strain rates 4.5 ×10-4 - 8.9×10-1 s-1) and the Amsler RKP300 instrumented impact pendulum (strain rates 73.0 – 141.3 s-1). The complete and detailed description of all experimental methods can be found in [11]. 2.2. FE simulations The studied material is represented by the simulati ons as the cubic polycrystalline aggregate containing 50 grains. The FE simulations are performed with the finite element code Z -set developed by Mines Paris and ONERA. The aggregate was created by the V oronoi tessellation and it is shown in Figure 1a). Th e grain size and the grain orientation correspond to the experimental observations [11]. The representative cube with grains has a volume equal to one. The texture corresponds to the experimental measurements , therefore the simulated polycrystal contains 40 grains with (111) plane axis oriented perpendicular to the surface, 5 grains with (011) orientation and 5 grains oriented with (001) plane axis perpendicular to the surface . Five different representations of microstructure were created in the following way. The polycrystalline aggregate (e.g. grain shape) remains the same in all cases as well as the set of crystallographic orientations. However, the crystallographic orientation of each grain is different in each microstructure representation. A higher number of polycrystal representations provides better statistics for the results analysis. The polycrystalline aggregate is meshed by the 136137 tetrahedral quadratic elements (C3D10 – 3D element with 10 nodes ). Boundary conditions simulate the cubic volume element of material with one free surface side. The applied load is cyclic tension -compression with R = -1. The remaining 5 sides represent the interior of the material volume. Therefore, they are constricted to remain flat and parallelepiped, but they can move in the direction o f their normal according to the applied load (Ui = f(σ)). Such a condition simulates the effect of surrounding material on the analyzed volume. An additional layer of elements is added to the top and bottom surfaces in order to apply uniaxial stress loading at the outer surfaces of the extension layer. The grap hical representation of boundary and loading conditions is shown in Figure 1b). Two different loading levels are used in simulations: ±100 MPa and ±140 MPa respectively, which are taken according to the experimental measurements [11]. The 140 MPa is the maximal loading applied during the experiments. The ±100 MPa load is chosen, because the material has a similar fatigue lifetime for conditions ±100 MPa at ±80 Hz an d 140 MPa at 20 kHz. This load level is also very close to the fatigue limit for 20 kHz, which is 95 MPa [11]. At this loading (stress amplitude) only minor struct ural changes are expected to occur at the material at the very high frequency of 20 kHz, since the samples reached the VHCF region (above 107 cycles to failure). More pronounced structural changes are expected to occur under 80 Hz as the sample reached 106 cycles to the failure [11]. Each simulation is performed for 100 cycles . Material behavior is described by the constitutive model based on crystal plasticity for FCC (face -center cubic) lattice with active {111} slip systems family [33][34]. It combines cubic elasticity and plastic slip (γ) as described by the equation: 𝛾˙𝑠=⟨|𝜏𝑠−𝑥𝑠|−𝑟𝑠 𝐾⟩𝑛 𝑠𝑖𝑔𝑛(𝜏𝑠−𝑥𝑠) (1) where τs is the resolved shear stress in s-th slip system, xs is the shear stress resulting from kinematic hardening, rs is the critical resolved shear stress resulting from isotropic hardening and constants K and n are parameters of Norton’s viscoplastic law. The isotropic hardening defines the expansion of yield surface and it is described by the equation: 𝑟𝑠=𝑟0+𝑄∑ℎ𝑠𝑟(1−𝑒𝑥𝑝(−𝑏𝑣𝑟))𝑘 𝑟=1 , 𝑣˙𝑠=|𝛾˙𝑠| (2) where r0 is the initial value of critical resolved shear stress, b, Q and v are model parameters, k is the total number of slip systems and hsr is the interaction matrix between different slip systems. Kinematic hardening describes the shift of yield surface center during the cyclic loading and its equation is given as: 𝑥𝑠=𝐶𝛼𝑠, 𝛼˙𝑠=𝛾˙𝑠−𝐷|𝛾˙𝑠|𝛼𝑠 (3) with C, D being constant parameters. Two sets of parameters are used to describe the material’s behavior during loading with 80 Hz and 20 kHz, respectively. Parameters for kinematic hardening (C, D) and interaction matrix (hrs) are taken from literature where the applied model was fitted for cyclic loading of copper bi -crystal [35]. Parameters for isotropi c hardening (r0, Q, b) are obtained by fitting the experimental tensile curves for copper at different strain rates [11]. The tensile strain rates correspond to the rates which are expec ted in the material during cyclic loading with 80 Hz and 20 kHz, respectively. A comparison of experimental and fitted curves is shown in Figure 1c) . The most important region is the low plastic strain (< 2%) where the numerical curves correspond closely t o the experimental ones. The discrepancy between the curves at higher plastic strains is caused by the hardening law containing exponential function (see equation 2) while the experimental curves show nearly linear hardening. However, such a large plastic strain is not expected in simulations due to the applied load, therefore the differences between the fitted and experimental curves will not a have significant effect on simulations results. The parameters of the constitutive model together with copper ela stic constants are listed in Table 1. Table 1: Copper elastic constants and constitutive model parameters describing copper plastic behavior at different loading frequencies. Elasticity C1111 [GPa] C1122 [GPa] C1212 [GPa] 159.3 122.0 81.0 Plasticity Frequency [Hz] N K [MPa ∙s] r0 [MPa] Q [MPa] B C [MPa] D 80 10 5 19.35 150 30 4500 600 20×103 10 5 111 30 5 4500 600 Figure 1: a) Illustration of analyzed polycrystalline aggregate with FE mesh, b) Boundary conditions applied in FE simulations, c) Fitting of constitutive model using experimental tensile test curves for different strain rates. The stress/strain state in side the polycrystalline aggregate is complex and it is fully described by the tensorial variables, therefore it is useful to contract these tensors into scalar values. Such a contraction can be done by using the von Mises stress (σmises) and equivalent st rain (εeq) that combine both axial and shear components of stress/strain tensors. The effect of loading frequency is described by the equivalent strain rate ( dεeq) and the amount of plastic deformation is described by the total plastic slip (γ), which is t he sum of plastic slip on all active slip systems. These four variables are described by the following equations: 𝜎𝑚𝑖𝑠𝑒𝑠=√0.5[(𝜎11−𝜎22)2+(𝜎22−𝜎33)2+(𝜎33−𝜎11)2+6(𝜎122+𝜎232+𝜎312)] (4) 𝜀𝑒𝑞=√(23⁄)[(𝜀112+𝜀222+𝜀332)+2(𝜀122+𝜀232+𝜀312)] (5) 𝑑𝜀𝑒𝑞=(𝜀𝑒𝑞1−𝜀𝑒𝑞0)𝑡⁄ (6) 𝛾=∑𝛾𝑠12 𝑠=1 (7) where stress tensor components (σii) are taken at the point of maximal external load ( ±100 MPa, ±140 MPa), the equivalent strain rate is calculated as the difference of strains at maximal external load (εeq1) and zero external load (εeq0) divided by the time of ¼ of loading cycle ( t = 0.003125 s for 80 Hz and t = 1.25 ×10-5 s for 20 kHz). 2.3. Phenomenological model s The applied numerical model calculates plastic slip as the continuous variable. However, for the further analysis of persistent slip bands , it is useful to estimate the dislocation density in the polycrystal volume. Fan et al. [36] derived the relation between flow s tress, strain rate and dislocation density. This relation can be modified for polycrystalline materials into the form: 𝜎𝑚𝑖𝑠𝑒𝑠=𝑀2𝐵𝑑𝜀𝑒𝑞 𝑓𝑎𝜌𝑏2+𝑀𝛼𝐺𝑏√𝜌 (8) where M = 3.06 is the Taylor factor, B is the dislocation mobility factor ( B = 21×10-6 Pa∙s [37][38]), fa is the factor describing the ratio of mobile dislocation density, b is Burgers vector ( b = 2.556 ×10-10 m), dεeq is strain rate. Parameter α is the dislocations interactions factor and G is the shear modulus (G = 45 GPa). The equation contains two parameters fa and α that can be used for its fitting to the particular material. The equation (8) with fitted parameter is used for the estimation of local dislocation density by inserting the values of stress and strain rate. The presented s tudy of the effect of external loading stress and frequency is based on the comparison of the thickness of PSBs and their volumetric density inside the polycrystal. The stress/strain variables obtained from equations (4 -7) combined with dislocation densities obtained from equation (8) can be applied as input for the basic PSB model shown in Figure 2. The illustratio n in Figure 2a represents the ladder -like structure of PSB which is typical for copper. It consists of dislocations walls separated by cells with low dislocation density [41]. Figure 2b shows the detail of one cell bounded by two walls. The cell contains moving screw dislocations [32]. The dimensions of the individual PSB str ucture are walls height hw, cells height hc and the PSB thickness tPSB respectively . The geometrical assumptions are that cell height hc = 10hw and dislocation densities are ρw = 100ρc [42][43]. The analyzed parameter is the PSB thickness. Figure 2: a) Illustration of PSB model with ladder -like s tructure . b) Detail of one cell bounded by walls with t he screw dislocations in side the cell represented by red lines. Plastic slip inside PSB is realized by the motion of screw dislocations inside the cell and its magnitude can be estimated as: 𝑐𝑠𝛾=12⁄𝜌𝑐𝑏𝑡𝑃𝑆𝐵 (9) where cs = 1 is the ratio of screw dislocations. Therefore, taking into account average dislocation density ρ = 10ρc, the size of PSB can be calculated as: 𝑡𝑃𝑆𝐵=20𝑐𝑠𝛾𝑏𝜌⁄ (10) Dislocation density inside the PSB cell can be also described by the equation: 𝜌𝑐=𝑁𝐷𝑥𝑡𝑃𝑆𝐵⁄ (11) where ND is a number of dislocations and x is their mutual distance. Combining the equations ( 9) and (11) gives the relation between the amount of plastic slip and a number of dislocations: 𝛾=𝑁𝐷𝑏2𝑥⁄ (12) This equation sets the limit for a minimal amount of plastic slip in PSB for a given number of dislocations. Loading direction hw hc hw wall wall cell tPSB b x a) b) 3. Results All parameters are evalu ated after simulation of 100 cycles at loading stress amplitudes (±100 MPa, ±140 MPa) to analy ze the effect of cyclic loading. Values from all five representations of polycrystal are combined together into one set to get the statistical distributions of the analyzed variables. The cumulative distribution of the von Mises stress is shown in Figure 3a. The results show similar stress distribution for both loading frequencies . The steep slope of the curves is the sign , that stress values vary within the inte rval of about 200 MPa throughout the whole volume of the polycrystalline aggregate . An increase of applied stress amplitude causes the overall increase of stress level and also the values variation represented by the decrease of the curves ’ slopes. The ave rage values of von Mises stress for different loading cases are listed in Table 2. A similar picture shows the distribution of equivalent strain rate (see Figure 3b). The calculated values of equivalent strain rate correspond to the ones chosen for the constitutive model fitting (see Figure 1 c) as shown by the average values listed in Table 2. The distribution of plastic slip is shown in Figure 3c. There is three to four orders of magnitude difference between the volume that undergoes plastic deformation for different loading frequencies . The 80 Hz loading causes plasticity in whole volume while 20 kHz loading only in 0.01 % of volume for 100 MPa and 0.1 % of volume for 140 MPa respectively. Plastic slip values show also smaller dispersion for 80 Hz loading which decreases with increasing applied stress amplitude. The 20 kHz loading causes broader dispersion of plastic slip values from 10-5 to 10- 1. The average values are higher for 20 kHz compared to 80 Hz loading (see Table 2) . This means that plastic strain is localized in a small volume where it reaches higher values for 20 kHz loading. Table 2 : Average values of variables inside the polycrystalline aggregates Loading case 80 Hz ±100 MPa 80 Hz ±140 MPa 20 kHz ±100 MPa 20 kHz ±140 MPa σmises [MPa] 307.6 ± 73.46 419.1 ± 103.3 278.1 ± 65.80 385.6 ± 91.04 dεeq [s-1] 0.1753 ± 0.0325 0 0.2364 ± 0.0448 0 40.18 ± 7.426 55.71 ± 10.25 γ [×10-3] 4.811 ± 1.900 9.327 ± 3.100 6.580 ± 11.10 14.40 ± 30.70 ρ [×1015 m-2] 2.569 ± 1.341 5.554 ± 2.923 2.244 ± 1.041 4.093 ± 1.695 Figure 3: Cumulative distribution s of a) von Mises stress, b) equivalent strain rate, c) plastic slip in the total volume of five polycrystalline aggregate s for different loading conditions. The dislocation density is estima ted from equation (8), therefore the parameters fa and α need to be identified. The f itting of these parameters is based on the results of the uniaxial tensile test s of polycrystalline copper at different strain rates published in [11]. The initial dislocation density was estimated from the size of the dislocations’ cells observed by TEM. The diameter of the cells is dc = 0.8245 ± 0.1532 μm and its relation to the d islocation density is given by the relation [39][40]: 𝐾 𝛼=𝑑𝑐√𝜌, 𝛼=−0.034log(𝜌)+0.767 (13) where K = 7.5. Solving this relation gives the initial dislocation density ρ = 1.29×1015 m-2 which corresponds with the typical density of wrought materials. The best fitting of equation (8) is obtained when the data are divided in two intervals according to the strain rate, and the resulting equations can be written as: 𝑓𝑎=1.765×10−10(𝑑𝜀𝑒𝑞4.5×10−4⁄ )0.768 for dεeq < 75 s-1 (14a) 𝑓𝑎=1.091×10−6(𝑑𝜀𝑒𝑞4.5×10−4⁄ )0.040 for dεeq ≥ 75 s-1 (14b) The value 4.5 ×10-4 is the smallest experimental strain rate applied in [11] and here it is taken as the initial static case. The resulting cumulative distribution of dislocation density is shown in Figure 4a. The distributions have a similar shape for all loading cases. The increase of dislocation density is caused by the increase of loading stress. The average values of dislocation densities are listed in Table 2. Th e comparison of average values of dislocation density and von Mises stress with equation (8) for different equivalent strain rates is shown in Figure 4b. Figure 4: a) Cumulative distribution s of estimated dislocations density for different loading cases . b) The average values of von Mises stress and dislocation density for different loading cases compared with the equation (8) for different equivalent strain rates. The comparison shows that the significant effect of strain r ate on strength occurs for strain rates higher than 1 s-1 and polycrystals loaded at ±140 MPa are in the range of dislocation densities where the sensitivity to strain rate is diminished. 4. Discussion The fatigue damage mechanism is mostly based on the rearrangement of dislocations structures during cyclic loading. Such a rear rangement leads to the creation of PSB s which consequently creates SMs on the surface that consist of intrusions and extrusions which act as the initiation sites for fatigue cracks. Three types of SMs are usually observed on the surface of copper polycrystal under the cyclic loading [9][11][12]. The SMs of type I are usually long and perpendicular to the grain boundaries. Type II have similar shape but they are aligned along grain and twin boundaries. Finally , type III SMs are small, short and less pronounced. Type I and II are created at high er loading stresses and low loading frequencies, which are typical for high -cycle fatigue regime (HCF). Type III SM s are observed at low applied stresses and high loading frequencies and it can transform into type I or II with an increasing number of cycles. They are typical for VHCF. The different morphology of SM s types is explained by the differences in dislocations structures under the surface. The PSBs have a typical “ladder -like” structure consisting of dislocations walls separated by channels with low dislocation density as depicted in Figure 2a [41]. The PSB structures are created when sufficient driving force is applied on the dislocations and when the frequency of the loading is low enough that the dislocations ’ moment of inertia does not significant ly affect dislocations motion . There are two mechanisms of plasticity in the channel based on the dislocation density. The first one (mechanism MI) is base d on the back -and-forth dislocations motion which occurs when the dislocation density is lower than the saturation level and there is no multiplication and annihilation of dislocations. The second mechanism (mechanism MII) occurs when the dislocation densi ty reaches the saturation level and the mutual dislocation spacing decreases to the annihilation distance. Then the plasticity is driven by the dynamic equilibrium between the multiplication and annihilation of dislocations. In the following analysis the m echanism MI is represented by the limiting case of one dislocation inside the channel ( ND = 1) which moves back and forward. The typical number of dislocations inside the channel is between 1 and 10 [42][44]. Therefore, mechanism MII is described in the analysis as five dislocations ( ND = 5) inside the channel distanced by x = 50 nm which is the annihilation distance for screw dislocations in copper [45]. Inserting the respective numbers ND into equation (12) provides the minimal plastic slip values of γ = 0.002556 for ND = 1 (mechanism MI) and γ = 0.01278 for ND = 5 (mechanism MII) respectively. Therefore, the PSB is present in the given location of the polycrystal only when the local plasti c slip is higher than the corresponding limit. The thickness tPSB is then calculated by equation (10). Cumulative distributions of PSB thicknesses in polycrystal volume for different loading stresses and frequencies are shown in Figure 5a,b and correspondi ng histograms are illustrated in Figure 5c,d. Assuming the operation of mechanism MI (ND = 1) and loading frequency of 80 Hz, 88 % and 99 % of the polycrystal volume are suitable for PSB creation for loading ±100 MPa and ±140 MPa , respectively. The histograms show the same peak for PSB thickness (100 - 200 nm) and an increase in loading causes an increase of the volume fraction of thicker PSBs. The range of PSB thicknesses is large , from 50 to 1400 nm. Such large volume fraction of PSBs in the polycrystal which does not change with amplitude of loading stress , does not correspond to the experimental observations . Therefore, the mechanism MI is unlikely in the case of 80 Hz loading . Results for the mechanism MI in the case of 20 kHz give a different picture. The loading level influences the possible amount of PSB s which make 0.004 % and 0.1 % of polycrystal volume for ±100 MPa and ±140 MPa loading , respectively . The peak of PSB thickness is similarly within the range of 100 - 200 nm but the range of thicknesses is relatively narrow , from 50 to 400 nm. The small volume fraction , thickness range and strong loading effect indicate, that mechanism MI can probably operate during high -frequency loading. Assuming the mechanism MII (ND = 5) for 80 Hz loading cases , it provides different scenario . The volume of PSBs decreases to 0.12 % and 12.3 % for ±100 MPa and ±140 MPa , respectively. Most of the predicted PSBs have a thickness within 100 - 200 nm, but the overall range decreases to the interval 50 – 600 nm, which makes the mechanism MII more probable when compared with TEM experimental observations of PSBs in [11]. The results for high - frequency loading cases do not change significantly for the mechanism MII. The volume fraction of PSBs decreases to 0.0013 % (±100 MPa) and 0.049 % (±140 MPa) which is within the same order of magnitude as the case for mechanism MI. The peak and range of thicknesses also remain similar to mechanism MI which makes both mechanisms probable to be operating under 20 kHz loading . Fatigue lifetime is related to SMs at the surface which are the initiation points for fatigue cracks. There fore, the analysis of the PSBs at the surface is performed. Figure 6 shows the relative surface fraction of PSBs with respect to their average thickness for different loading cases and assumed dislocation motion mechanisms. The results show that the assump tion of mechanism MI (ND = 1) for 80 Hz loading would lead to the creation of PSBs at the majority of the surface without a significant difference for different loading stress es. The activation of multiplication and an nihilation of dislocations (mechanism MII, ND = 5) corresponds better with experimental observations . Figure 6 shows that for the case of mechanism MII, t he surface area of PSBs depends o n loading stress (1.8 % for ±100 MPa and 16 % for ±140 MPa) and the average thickness is about 160 nm with a small dispersion of values which predicts the similar morphology of PSBs. The high -frequency loading (20 kHz) shows a strong effect of loading stress on the surface area of PSBs for both dislocations ’ mechanisms (MI, MII) . The relative area increase s from about 0.02 - 0.05 % for ±100 MPa to 0.6 - 1 % for ±140 MPa. The average PSB thickness increases with increasing loading stress and it is higher for mechanism MII (ND = 5). Such larger thickness together with wider dispersion of values compare d to 80 Hz loading suggests that the PSB s have different morphology and are less localized when induced at 20 kHz . Such a morpholog y can be related to experimental ly observed type III SM s. The spatial distribution of von Mises stress, plastic slip and PSB thicknesses for mechanism MI (ND = 1) and mechanism MII (ND = 5) for all loading cases for one representation of polycrystal line aggregate free surface are shown in Figure 7. The stress and plastic strain maps show the concentrations at grain boundaries and junctions. There is a clear difference in plastic slip distribution for cases of 80 Hz and 20 kHz loading as the latter ones have the majority of the surface in an elastic state. The map of PSB thicknesses estimated according to the mechanism MI (ND = 1) shows that for 80 Hz loading the PSBs are located at the grain boundaries as well as in the grains ’ interior. They can be found also in the areas of low stress and plastic strain which doe s not correspond to the experimental observations. Distribution of mechanism MI ( ND = 1) PSBs for 20 kHz loading shows their location at the grain junction with the highest stress and plastic slip concentration which is a realistic prediction. The mechanis m II ( ND = 5) gives a more realistic map of PSBs for 80 Hz loading. The location of PSBs is at the grain boundaries and locations of stress and plastic slip concentrations and the overall amount of PSBs increases with loading stress. The map s of PSBs for b oth mechanism (MI and MII) are very similar for 20 kHz. The images of experimentally observed surface s for different loading conditions are illustrated in Figure 8. Comparison of the location of experimentally observed SMs with the predicted maps in Figure 7 can clarify the difference in the plasticity mechanism at different loading frequencies. The experimentally observed SMs are isolated features at the surface. Similar distribution can be seen in Figure 7 for mechani sm MII ( ND = 5) for 80 Hz loading and mechanisms MI and MII for 20 kHz loading. Therefore, t he low -frequency loading causes the creation of PSBs and SM s via the mechanism of dislocations multiplication and annihilation, while creation of PSBs and SMs at high- frequency loading can be possibly driven by both multiplication/annihilation and also by back and forth motion of dislocations. The estimated thicknesses of PSBs are within the range of 100 - 500 nm which is within the same order of magnitude as experim ental observation shown in Figure 8 and [11]. Despite the simplicity of the model, it provides a suitable prediction of PSBs locations and thicknesses. Furthermore, the model predicts similar PSBs distribution s for the loading cases 80 Hz ±100 MPa and 20 kHz ±140 MPa. These two cases show the same fat igue life [11]. Therefore, the model predictions can be also correlated with the experimentally measured fatigue life. Figure 5: Cumulative distributions of PSB thicknesses for different loading cases and dislocations movement mechanisms a) mechanism MI (ND = 1), b) mechanism MII (ND = 5). Histograms of PSB thicknesses for different loading cases and dislocations movement mechanisms c) mechanism MI (ND = 1), d) mechanism MII (ND = 5). Figure 6: Average PSB s thickness at the surface of polycrystals for different loading cases an d assumed dislocation movement mechanisms (MI and MII) . Figure 7: Maps of von Mises stress, plastic slip and PSB thicknesses for different dislocation movement mechanisms MI and MII for different loading cases . Comparison is made for one representation of polycrystalline aggregate. Figure 8: Experimentally observed slip marks on the surface of samples tested at analyzed loading conditions. 5. Conclusions Numerical simulations of copper polycrystalline aggregate subjected to cyclic loading with different frequencies and stress levels were performed using the finite element method combined with crystal plasticity. The analysis was focused on the development of persistent slip bands based on the localization of stress and plastic slip inside the polycrystal and its correlation with experimentally observed slip marks and fatigue lifetime. The findings can be summarized into the following points: • The areas of p lastic deformation develop preferentially at the grain boundaries near the free surface of the polycrystalline aggregate. • More pronounced deformation evolution and localization were observed for the lower te st frequency (80 Hz) when compared to high test f requency (20 kHz) showing only very limited regions undergoing cyclic plastic deformation. • The applied model of PSB s predicts that the most probable plastic ity mechanism in PSB for 80 Hz loading is based on the multiplication and annihilation of dislocatio ns during cyclic loading, while the plasticity in PSB at 20 kHz can be also driven by the back -and-forth motion of dislocations. • The applied model predicts thinner PSBs with more uniform thickness at the surface for 80 Hz loading which corresponds to the appearance of type I SMs, while the prediction of thicker and heterogen eous PSBs for 20 kHz loading corresponds with type III SMs. • Similar morphology and spatial distribution of PSBs at the surface predicted by the model for different loading condition c an be related to the similar fatigue lifetime that is experimentally measured . Therefore, the model predictions can be used for fatigue lifetime estimations. Acknowledgment: This publication was supported by the project “Mechanical Engineering of Biolog ical and Bio - inspired Systems”, funded as project No. CZ.02.01.01/00/22_008/0004634 by Programme Johannes Amos Commenius, call Excellent Research. Data availability: Data related to the paper are available in the Zenodo repository with DOI: 10.5281/zenodo.11183542. 6.