Effect of 3D Representative Volume Element (RVE) Thickness on Stress and Strain Partitioning in Crystal Plasticity Simulations of Multi-Phase Materials

: Crystal plasticity simulations help to understand the local deformation behavior of multi-phase materials based on the microstructural attributes. The results of such simulations are mainly dependent on the Representative Volume Element (RVE) size and composition. The effect of RVE thickness on the changing global and local stress and strain is analyzed in this work for a test case of dual-phase steels in order to identify the minimal RVE thickness for obtaining consistent results. 100 × 100 × 100 voxel representative volume elements are constructed by varying grain size and random orientation distribution in DREAM-3D. The constructed RVEs are sliced in depth up to 1, 5, 10, 15, 20, 25, 30, 40, and 50 layers to construct different geometries with increasing thickness. Crystal plasticity model parameters for ferrite and martensite are taken from already published data and assigned to respective phases. Although the global stress/strain behavior of different RVEs is similar (<5% divergence), the local stress/strain partitioning in RVEs with varying thickness and grain size shows a considerable variation when statistically compared. It is concluded that two-dimensional (2D) RVEs can be used for crystal plasticity simulations when global deformation behavior is of interest. Whereas, it is necessary to consider three-dimensional (3D) RVEs, which have a specific thickness and number of grains for determining stabilized and more accurate local deformation behavior. This estimation will help researchers in optimizing the computation time for accurate mesoscale simulations.


Introduction
The micro-structure of a material plays an important role in defining the mechanical properties [1] and service life of a component [2,3]. Numerical models can help to understand and improve the component's life by providing detailed insight into the local [4,5] and global [6][7][8] deformation behaviors. A concept of Representative Volume Element (RVE) is used in order to numerically simulate continuous yet locally heterogeneous materials [9,10]. A lot of work in the recent past has been carried out to estimate the local deformation behavior of single and multi-phase materials based on two-dimensional (2D) and three-dimensional (3D) RVEs [11][12][13][14]. 2D and 3D RVEs can be constructed using single [15] or multilayer [16] Voronoi tessellation using measured or virtual local grain size, phase, without compromising the accuracy of results is still missing. Different researchers in the past have approached this problem differently. The method that is presented in the current work is unique, fast, and effective in establishing the desired output.
In this work, the effect of RVE thickness-with the intention of identifying the amount of finite thickness required for consistent results-on the changing global and local stress and strain is analyzed for a test case of dual-phase steels. RVEs were constructed by varying mean feature ESD (D f ) between 3 to 18 µm using DREAM-3D with grains having randomly assigned orientation distributions. The constructed RVEs were sliced in up to 1,5,10,15,20,25,30,40, and 50 layers to produce different geometries comprising the same microstructure-with increasing thickness. Crystal plasticity model parameters for ferrite and martensite are taken from already published data and assigned to respective phases. Probability Distribution Functions (PDF) and Cumulative Distribution Functions (CDF) of all simulation results are compared in order to estimate the solution convergence with changing grain size statistically. In the end, a simple function is proposed for calculating the sufficient RVE thickness that is necessary for obtaining a converged solution.
Section 2 provides the details of material data, RVE construction, simulation scheme, and CP material model parameters. Section 3 presents the results that were obtained in this study. In Section 4, the results are discussed in comparison with state of the art, and insight into the outlook is provided. Eventually, the study is concluded in Section 5.

Numerical Simulation Model Development
For the numerical simulation modeling, RVEs were virtually constructed while using open source tool DREAM-3D [29]. The case of dual-phase steel is considered as a case study in the current work. Micro-structural attributes were adopted from the literature [45]. 3D RVEs were constructed-with varying grain size-using open source tool DREAM-3D [29]. These constructed RVEs were sliced into varying thickness ranges from one layer to 50 layers and used as input geometries for numerical simulations. Crystal plasticity based open source tool DAMASK was used to model and solve numerical simulations by adopting already published methodology [22].

RVE Construction and Geometry Files Production
Dual-phase steel is chosen as a case study in the current work, because, on the one hand, it is a simple material that consists of ferrite matrix and embedded martensite islands (consisting of laths with different crystallographic orientations) in it, while, on the other hand, it is an excellent example of multi-phase materials with two phases comprising of drastically different mechanical properties. The microstructural attributes i.e., grain size and martensite percentage, were adopted from the previous work of Jiang et al. [46]. It has been reported that, in DP steels, the grain size of ferrite ranges from 5 to 25 µm, whereas martensite grain size ranges from 3 to 18 µm. Considering this upper and lower grain size limit, DREAM-3D was used for the construction of virtual RVEs with varying grain sizes, while keeping the ferrite and martensite grain size ratio to be 1:1 intentionally to keep the problem at hand simple. Table 1 provides the total number of grain in each 100 3 RVE for a corresponding ESD given as a dimensionless entity. A pipeline was built with the initialization of virtual data generation while using a stats generator filter. The ratio of 90 vol. % primary phase(ferrite), 10 vol. % secondary phase(martensite) was chosen for all cases to simplify the simulation model results and generalize the comparison of statistical data. As the crystal structure of both Ferrite and Martensite is body-centered cubic. Therefore, they were generated as "cubic equiaxed" phases. Different RVEs were generated by varying the ESD statistical distribution data. The grain size distribution of ferrite and martensite for each RVE taken is given in Table 1.
The synthetic volume size of all RVEs was kept 100 × 100 × 100 voxels (where the size of each voxel is 1 µm 3 ) and specified number of particles could pack in the defined volume, which means that more grains packed in the defined volume when D f was small and lesser when D f was large as given in last column of Table 1. For both phases, grain shape was stated to be "ellipsoid". The maximum number of iterations (swaps) allowed was 100,000. Using already published methodology [30], the generated RVEs were recorded as .xdmf file, which is readable by Paraview for visualization, and .geom files were saved, which are readable by DAMASK. Figure 1 presents the general flow of information in the pipeline with associated inputs and outputs.   Table 1 provides details of statistical input data for A, B, C, D, and E. Figure 2 shows a schematic diagram representing the generated 3D RVEs and their slicing sequence for the construction of simulation geometries. RVE-A and E are shown for representation, whereas B, C, and D were modeled in a similar fashion. Ferrite (F) and Martensite (M) are shown separately in the sliced images for better visualization of phase distributions. This research's primary objective is to identify the optimal RVE thickness with respect to the material grain size for obtaining converged local results. For this reason, each constructed 3D RVE with varying grain size was sliced into geometries with increasing thickness from one to 50 layers. Slicing sequence is also shown in Figure 2 for each RVE with yellow dotted lines, whereas green dotted lines in the schematic diagram represent more sliced geometries that are not shown in the current figure, but were constructed for simulations. There are 100,000 gaussian mesh elements in one layer, which means that each geometry contains 100,000 gaussian mesh elements multiplied by the number of layers considered in each case. A complete grain comes into account with a lesser number of slices in the case of small D f , whereas more slices are needed in order to represent full-grain with bigger D f .

Material Properties and Loading Conditions
Although there is a large difference in the mechanical properties of ferrite and martensite, in current research, both of the phases were assigned elastic-viscoplastic properties by adopting already developed phenomenological power law available in DAMASK [47]. Details of the hardening model with equations is provided in Appendix A for the readers who are not familiar with the constitutive equations used in DAMASK. The elastic coefficients, variables defining plastic flow behavior and fitting parameters for both phases were adopted from the already published literature [45]. Table 2 presents the adopted parameters. In the spectral method each grid point is considered as a computation point which is assigned a phase and initial crystallographic orientation before the simulations start. In the current work, mixed boundary conditions were applied uni-axially along x-direction as follows, while keeping the out-of-plain surfaces in all geometries stress-free: Periodicity of the solution is inherent to the spectral method due to the FOURIER approximation of the deformation gradient field [47]. Yet, it is known that the immediate neighborhood mainly influences the strain heterogeneity, and the influence of artificial periodicity introduced by the boundary description is confined to a narrow zone [22]. Therefore, in the current study, this effect is ignored. For strongly different phases i.e., ferrite and martensite, the reference stiffness has a strong influence on stability and convergence rate, as was shown earlier by Michel et al. [48]. The readers are encouraged to refer to the work by Diehl et al. [19] for the explanation of the assumptions in the FFT regarding the reference stiffness that are used within the framework of the crystal plasticity provided by DAMASK. The simulation results were post-processed by using already available subroutines in the DAMASK installation module, and data were further statistically analyzed using the Seaborn library in Python. The local stress-strain distributions were visualized using open source tool Paraview [49].

Results
If global/averaged stress-strain behavior of all the simulations is plotted and compared for multiple layers of a specific RVE, as shown in Figure 3, similar results with slight variation of slope are observed for varying thickness of geometry (Figure 3a,b) or variation of grain size (Figure 3c). It can be clearly observed in these figures that the trend of stress-strain curve is the same for all grain sizes with slight variation (higher stress response to same strain with increasing grain size). It is observed that there is not much difference in the results and almost any geometry yields very similar results. Some numerical simulation models crashed due to the non convergence of global stress/strain equilibrium at the spectral code level, and therefore are presented up to that certain level of strain.
Although stress-strain curves are used by engineers and scientists to understand the overall material behavior, they can be very misleading in the case of such full phase simulations where the local results may vary drastically. Still, the averages' response remains the same. To make this point, during the post-processing of the data, the RVEs were sectioned as schematically represented in Figure 4. This scheme was adopted to expose the top surface and middle section of the RVE-E as an example case. Local stress and strain values are represented with the same scale to show how they change with varying RVE thickness at 25% of true strain. Figure 5 shows the local von Mises true stress distribution in all geometries of RVE-E at 25% of true strain. Only the ferrite phase is shown here for better visualization by filtering out the martensite phase, which, due to very high stresses (≈2.5 GPa), distorts the scale. It is observed that there is a high contrast of stress distribution in the 01-layer RVE with some areas of very high stresses and others with very low stresses. As the thickness of the RVE is increased from 01 layers to 50 layers, the stresses on the surface diminish and they are relatively more homogeneously distributed within the matrix, and the high contrast for stresses diminishes. There is very less difference in the local stress distribution of 40-layer and 50-layer simulations. In these 3D simulations, it is observed that, although the phase interface is more prone to higher stresses, it is not always the case. In the middle section, it is observed that the local stress distribution becomes consistent with similar areas of high and low stresses. This similarity in obtained results-with increased RVE thickness-represents the convergence of the point to point local solutions.   Figure 6 shows local von Mises true strain distribution in all geometries of RVE-E at 25% of global true strain. In this figure, both-ferrite and martensite-phases are shown. Embedded martensite grains undergo negligible plastic strain during overall deformation and, hence, exhibit almost zero strain in Figure 6 (pointed out by green arrows). In 01-layer simulation results, it is observed that the local strain contrast is quite large with sharp strain channels around martensite grains oriented 45 • to the applied load direction. As the RVE thickness increases from 01-layer to 50-layers, it is observed that this sharp strain contrast on the top surface diminishes due to strain distribution in the third dimension. In the middle section of simulated geometry, it is observed that the local strain distribution converges with similar solution output in case of 40-layer and 50-layer geometries, respectively (pointed out by red arrows). The magnitude and position of local strain distribution in these cases are identical. The visual comparison of local stress and strain distribution in multiple varying geometries to observe the convergence of results-as shown in Figures 5 and 6-is a very challenging task. Visual inspections are primarily dependent on subjective choices; therefore, statistical data analysis tools are adopted in the current research in order to work out the convergence of the observed results. For statistical analysis, PDFs and CDFs of true local stress and strain distributions in each phase are constructed for each simulated geometry at the maximum global strain. Local stress and strain distributions in each phase of each geometry are compared. This detailed comparison is shown in Figure 7 by intelligently grouping data in different subplots.
It is observed that, with small D f , i.e., in the case of RVE-A, as shown in Figure 7a, the local stress, and strain distribution in both phases is quite different for 01-layer geometry as compared with thicker geometries. It is observed that, with increasing geometry thickness, the PDFs and CDFs become similar after more than 10-layers.
When D f increases, i.e., in case of RVE-C as shown in Figure 7c, similar trend of convergence with increasing geometry thickness is observed. The distribution varies up to 20-layer geometry for the current case and it does not change with a further increase in geometry thickness. It is observed that, with increasing D f in all RVEs i.e., in Figure 7a-e, the PDFs and CDFs converge with increasing geometry thicknesses, but more geometry thickness is needed when D f is large. It is an expected response because for RVEs with large D f more geometry thickness is needed to define a grain completely and hence the flow of stresses and strains around it becomes possible.
The stress and strain distribution behavior of martensite in thin geometries (one-layer to 15-layers) is different from thick geometries (20-layers to 50-layers), as observed in Figure 7(iii,iv). One-layer simulations in Figure 7(iii) represent a large and packed strain distribution profile compared to higher thickness results. Distribution is relatively more dispersed over a broad strain range, independent of D f . The stress distribution profile of one-layer simulations in Figure 7(iv) for all RVEs shows two peaks and a wider dispersion, whereas the distribution in close to bell shape when the RVE thickness is increased.
In Figure 7a-e, it is observed that at 50-layers geometry-due to very less change in the local stress and strain distribution-a converged solution for all RVEs is obtained. The PDF and CDF plots for both phases and all geometries with varying D f are compared in Figure 8. It is observed that the curves accurately match with a slight difference in the peak values. When considering the same composition and material properties in all cases, this comparison confirms the convergence of the obtained results. From these data, one can interpret that with increasing layers in the RVE, the material volume and number of mesh points increase, or more specifically, a total number of randomly oriented grains increase. Therefore, the statistical behavior converges towards an average and, hence, produces a false notion of local convergence. This interpretation is not correct, as it can be verified by comparing local plots in Figures 5 and 6 that the local point-to-point convergence of the results happens, which is captured by the statistical comparison presented in Figure 7.  Although with increasing geometry thickness, the stress and strain distributions for both phases in Figure 7 are observed to move in multiple dimensions with varying shape. The maxima were noted and normalized against 50-layer geometries to simplify the convergence criteria. Peak values of PDF normalized against the 50-layer thickness simulation for ferrite strain are shown in Figure 9a and for ferrite stress in Figure 9b. Here, it is important to mention that the convergence of a result was analyzed against 50-layer thickness simulations, assuming them as perfectly converged, which might not be the case.

Discussion
Full phase simulation models are extensively used to study microstructural attributes' effects on local and global material deformation behaviors. The validation of such simulation models simultaneously on a global and local scale is a challenge that has not been addressed in the existing literature [32,41,42], but with its limitations. Crystal plasticity-based full-phase simulation models rely upon many constitutive and fitting parameters identified by comparing the experimental evolution of deformation mechanisms with the simulation results. These averaged stress-strain curves are compared to represent the accuracy of simulation results [17,22]. Based on such-global results based calibrated-models, the local material behavior is analyzed and studied. This methodology is useful for the inexpensive employment of such models. From the results of the current work, it is now clear that such methodology can be misleading as a wide variety of microstructures, which might result in very different local results, can yield the same global results.
It is shown in Figures 5 and 6 that the problems arise when the local deformation behaviors are compared. Owing to time and capital-intensive critical task, not many researchers in the past have carried out such analysis and observed that results only match qualitatively [19,31]. It has been reported by earlier researchers [19] that, in 2D DAMASK simulations, there is a high local stress and strain contrast due to the non-availability of the third direction. Experimental analysis of stresses and strains in 3D is almost impossible due to which the 3D simulation results cannot be validated. Additionally, paradoxically, if 3D EBSD data are measured, then there is no sample left to perform experimentation. On the other hand, if in-situ tests are executed on a sample to analyze deformation behavior, 3D EBSD data cannot be collected.
In the past, it has been elaborated that the 3D results are different from the 2D results, and realistic 3D geometries are better than the generalized 3D geometries [19]. How much of the 3D dimension is required for a converged solution was reported by Diehl et al. [32] by running simulations with fixed grain size RVEs. In this study, the effect of varying grain size of ferrite and martensite in DP steel is analyzed in order to study its effect. RVEs with varying mean feature size were synthetically constructed for dual-phase steel case in the current work. The RVEs were sliced to construct geometries with varying thicknesses from one to 50 layers. In 2D simulations with actual EBSD data or synthetically constructed RVEs, non-realistic high contrast stress and strain distribution are observed. It is evident that there is a drastic difference in the local results in 2D geometry when compared with 3D geometry.
The current work has not analyzed the effect of free surface on the convergence of the results, which can be important when comparing the converged local solutions with in-situ test observations. This should be kept in consideration while adopting the current methodology for such analysis. The readers are encouraged to refer to [32,50] for such modelling methodology.
In current work, statistical probability and cumulative distribution curves were compared for each phase in order to comprehensively analyze the local stress and strain distributions. The following relationship can be drawn for the obtained results, which is in accordance with the previous publications [32,42]: An empirical function can be drawn and it generalizes the convergence results trend by normalizing the peak of stress and strain PDFs in the ferrite phase. Figure 10 shows the constructed empirical convergence criteria from the results of RVE-A, B, and C converged simulations. Only the first three points were considered in the development of the proposed empirical model. This is because the convergence of results here is used as a relative term against 50-layer thickness simulations, which might be misleading in the case of larger grain size RVEs. More simulations are needed with bigger RVEs in the future to be sure about their convergence. It is understood from previous work [18] that a large number of grains is vital for a solution convergence and, therefore, it is not good practice to reduce the total number of grains below 500 in an RVE. This criterion helps in setting the upper bound limit of D f to 15. The number of layers required for a specific D f can be calculated using the following derived empirical function: this equation is very specific and it is only valid for the similar grain sizes of two phases having 100 3 µm 3 RVE size. To generalize the conclusions for a more general use, the equation can be modified as: This given criterion is the new finding of the current research. It should be adopted and fulfilled while carrying out full phase simulations and discussing the local stress and strain evolution in a given microstructure.

Conclusions
In this research, the effect of RVE thickness on the changing global and local results is analyzed for an exemplary case of dual-phase steels. This research's primary objective is to identify the optimal RVE thickness with respect to the defined feature size for obtaining converged local results. RVEs were constructed by varying the mean feature sizes between 3 to 18 µm using DREAM-3D with grains having randomly assigned orientation distributions. The constructed RVEs were sliced into 1,5,10,15,20,25,30,40, and 50 layers to construct different geometries with increasing thickness. Crystal plasticity model parameters for ferrite and martensite phases were adopted from already published work. The global and local simulation results were directly and statistically compared in order to draw the following conclusions from the study: 1. As long as the orientation distribution and composition of the material is the same-even with changing mean feature size, the total number of grains and geometry thickness-the global deformation and stress-strain behavior of the material does not change drastically (<5% divergence) 2. Although the global stress-strain behavior of different RVEs is similar, a large variation in the local stress-strain of RVEs with varying thickness and the feature size is observed during the visual and statistical comparison. 3. Stress and strain probability distribution in 2D RVEs is drastically different from 3D RVEs. In 2D RVEs, the strain prediction in the soft phase is higher, in the hard phase is lower, whereas the stress prediction in the soft phase is lower and the hard phase is higher. 4. It is necessary to consider 3D RVEs, which are at least five times larger than the average grain size for stabilized and more accurate local deformation behavior determination. This estimation helps in optimizing the computation time for accurate mesoscale simulations. 5. Figure 10 and Equations (4) & (5) provide a criterion for choosing the mean feature ESD and effective RVE thickness. This criterion should be fulfilled for full phase simulations carried out using DAMASK in the future for obtaining converged simulation results for multi-phase materials i.e., DP steel.

Appendix A. Phenomenological Crystal Plasticity Model
The model is adopted for the body-centered cubic (bcc) crystals of the phenomenological crystal plasticity description by Peirce et al. [51]. The microstrcuture is parametrized in terms of a slip resistance S α {011} on each of the 12 {011} 111 slip systems, and S α {211} on each of the 12 {211} 111 slip systems which are indexed by α = 1, . . . , 24. These resistances increase asymptotically towards S α ∞ with shear γ according to the relationshipṠ with interaction (h αβ ) and fitting (w, h 0 ) parameters. Given a set of current slip resistances, shear on each system evolves at a rate ofγ with τ α = S.(b α ⊗ n α ), a reference shear rateγ 0 and a stress exponent n. The superposition of shear on all slip systems in turn determines the plastic velocity gradient: where b α and n α are unit vectors along the slip direction and slip plane normal respectively. Using the quasi-static strain rate of 1 × 10 −4 all RVEs were loaded in tension for 2500 s (load direction with reference to the top surface is shown in Figure 2 with red arrows) to reach a total strain of average 25%.