Integrating phase field and crystal plasticity finite element models for simulations of titanium alloy Ti-5553

The ability to simultaneously predict the microstructure and bulk material properties of 3D printed (additively manufactured or AM) metals is critical to the development of process intelligence that can be used by a digital-twin for forecasting and optimising alloy composition and fabrication parameters. This study proposes a simulation framework for predicting the microstructure and corresponding meso- and macro-scale properties of AM materials. This is achieved by integrating phase-field and crystal plasticity modelling techniques, whereby the phase field model predicts the microstructure and the crystal plasticity constitutive model computes the stress–strain evolution using the microstructure as the input. The simulation of multiple microstructures demonstrates that this integrated approach can be used to test the influence of different microstructures on the mechanical properties of titanium alloy Ti-5553. This includes the influence of grain size and grain orientation on both the meso- and macro-scale behaviour.


Introduction
A digital twin (DT) of additive manufacturing (AM) can be a powerful contribution to the development of AM. A DT has the potential to allow for improved understanding of the influence of processing parameters on process-structure-property relationships via bi-directional data connections between the physical and virtual AM processes. One such contribution to the virtual replica of the AM process is validated computational models to calculate the structure and properties of a material from its thermal history. This process-structure-property information is a description of the material, such as grain morphology or dendrite arm spacing given some process or thermal history, and the relationship between microstructure and properties. In order to optimise the properties of a build, integrated computational models can be used to determine how that microstructure affects the elasto-plastic behaviour of the material under various loading conditions. A multi-scale, 'multi-physics' approach is required to accurately compute the microstructure and the resulting mechanical properties.
There is a growing body of work in the development of integrated computational tools to enhance the design process of materials. Collins et al [1], developed a predominately continuum-scale approach to assess the influence of processing parameters on the properties and performance of the material. Keller et al [2] performed simulation analysis of a Ni-based superalloy using a combination of finite element (FE) and phase-field (PF) modelling. In that paper, the thermal history was estimated by FE and in situ thermographic measurements, extracting geometric and thermal features of the meltpool, which were then used as input to PF simulations. A multi-stage integration tool was developed by Yan et al [3] to simulate different stages in the manufacture of Ti-6Al-4V by electron beam melting (EBM). Their integrated process used thermal-fluid flow models, a cellular automaton, and self-consistent clustering analysis to predict mechanical properties and fatigue-life, with results comparable to the reported experimental findings. Liu et al [4] also investigated a modelling framework for Ti-6Al-4V. Their tool used FE to predict thermal fields based on the scanning strategy and temperature, which was then passed to PF simulations to predict different microstructures. These microstructures were then analysed using fast Fourier transform-based elasto-viscoplastic (EVP-FFT) models to predict the mechanical behaviour. Simulating AM of Ti-6Al-4V was also investigated by Pinomaa et al [5] who developed an approach to predict the stress-strain evolution and fatigue-life of the material based on a PF simulated microstructure, using information from thermal history simulations. This approach was then utilised to predict the link between processing parameters and mechanical proprieties of AM 316 L [6].
Many of the contributions to integrated computational models of AM processes has been primarily focused on the predicted macro-scale properties of the material. However, these computational tools offer a significant advantage in the development of understanding how features at the microstructural level influence local deformation. This information is vital in assessing meso-scale failures arising from cracks or voids in the material. To investigate how integrated computational models can be used to gain deeper insight into not only macro-scale behaviours but also meso-scale deformation, an integrated PF and CP model is used to simulate monotonic tensile deformation of titanium alloy Ti-5553 manufactured using the selective laser melting (SLM) technique. Due to difficulties in modelling the process-structure relationship under AM conditions [7] this study focuses on the structure-property aspect. The developed PF model is used to generate trial data. Five microstructures with different grain environments are used in CP simulations to investigate the influence of these environments on macro-and meso-scale deformation. Phase-field simulations are conducted within the open-source phase-field package OpenPhase [8], while CP simulations are conducted using the fast Fourier transform (FFT) based solver AMITEX [9]. Ti-5553 was chosen because it is a β alloy and the composition can be approximated as a binary Ti-Mo alloy through the equivalents heuristic [10], both of which simplify modelling in a commercially available alloy. The results highlight the advantages of integrated computational models to provide greater insight into the influence of microstructural features on the meso-and macro-scale behaviour of AM materials.

Microstructure model
This section describes the details of the PF approach, the relevant material details, and other information relevant to the generation of trial microstructures to be integrated into the CP model. Most approaches to PF modelling of the additive manufacturing process focus on the formation of the fine cellular-dendritic solidification microstructure [2,[11][12][13]. These approaches address fine scale features of the solidification microstructure, such as the cell spacing, dendrite arm scale, and solute partition, but neglect details at the scale of multiple grains. The present approach is to use a PF model at the scale of the melt track, using an approach to capture coarse grained grain-scale information that can be passed to the CP model described later. Here, 'coarse-grained' implies that smaller scale solidification details are neglected, such as cell formation, and the focus is on the formation of grains. Neglecting the solidification structure allows the sampling of the evolution of the grain in the absence of strong interface complexity. However, microsegregation cannot be modelled, as the larger curvature of the solidification microstructure cannot be resolved. Despite neglecting this complexity, broad details of the grains can be resolved by connection to the approximate thermodynamics. This is achieved through the use of a linearised phase diagram for the free-energy functional and control of the initial distribution of grains. The microstructure is modelled using a multi-grain PF model [14] with the equilibrium partition model implemented in OpenPhase [8], an open-source phase field solver. While it is possible to develop quantitative phase field simulations for up to four diffusing species, curating thermo-kinetic data of a quinternary system is difficult, and no public domain database currently exists for Ti-5553. Therefore, a conventional binary approximation is applied in order to use a single diffusing species.
Ti-5Al-5Mo-5V-3Cr, or Ti-5553, is a metastable β-titanium alloy. When processed by SLM using parameters optimised for density, the grain structure of Ti-5553 has distinct regions of columnar and equiaxed grains [15]. Smaller equiaxed grains form at melt pool edges, where the initial undercooling in the melt is sufficient for nucleation, T < T Nu ; and large columnar grains at larger undercooling that can persist through several build layers, where the dominant mode of solidification is epitaxial. At a smaller scale, internal to the grain, a cellular solidification structure exists with cell diameter on the order of d cell ≲ 10 −6 m.

Binary analogue
Solute atoms are represented by a single solute species, molybdenum (Mo), via the molybdenum equivalents heuristic [10], denoted Mo Eq . Using the nominal composition of Ti-5553 the system composition is computed to be 9.25 wt% (≈ 4.84 at%) Mo Eq . At the liquidus, T L , a liquid with composition 4.84 atom % Mo Eq (T L ≈ 1977.48K) coexists with a solid with composition 6.54 atom % Mo Eq . A binary phase diagram, computed using ThermoCalc, is presented in figure 1. At this concentration, the slopes of the liquidus and solidus are well approximated as linear. In the case of a lens shaped solid-liquid region, the solid is solute rich with respect to the liquid, and thus at the expected high rate of solidification, the depletion of the solute, the analogous phenomena to solute rejection for solute lean solids, is expected to be small. However, the concentration difference between the solid and the liquid is small, and as a result the barrier to formation of the solid phase with the concentration of the liquid requires only a small undercooling. Therefore, at the large thermal undercooling expected, significant solvent trapping will lead to essentially flat concentration profiles across the interface [15,16], although the presence of a cellular solidification pattern observed in experiments may imply a small degree of micro-segregation it could also be related to the production of latent heat [15].

Driving force and free-energy related quantities
The PF evolution equations require an estimate of the driving force to model the motion of the interface and quantities controlling the diffusive behaviour of any solutes. Free-energy related quantities are derived from a linear fit to the solidus and liquidus on the phase-diagram, and the driving force, ∆g SL , is computed by [14].
where T ′ is the intersection of the liquidus and solidus with the c = 0 composition, i.e. the melting temperature of pure Ti (1941 K), m L and m S are the slope of the liquidus and solidus respectively, c i is the current composition, c L and c S are the composition of the corresponding liquid and solid phases, T i is the current temperature, and ∆S SL is the entropy of fusion. Their corresponding values are presented in table 1.

Thermal field
There is no general agreement about the evolution of the thermal profile in the melt [11,[17][18][19][20][21]. Convention is to use FE simulations of the laser-metal interaction and extract linear approximations of the thermal history, and apply these quantities to the temperature history via the frozen temperature approximation.
Since there is no data available specific to Ti-5553 four assumptions were made about the initial conditions.  [24] (a) At the interface between melt pool and substrate, or the previous build layer, the temperature, T, is equal to or less than the solidus, T S in the liquid directly ahead of the interface, T ⩽ T S . Therefore the initial temperature is set at the starting point of a cell as T S , the corresponding value of the solidus directly below the liquidus. (b) For a new layer, at small undercooling before epitaxial growth is established, the thermal undercooling may decrease below the critical nucleation temperature, and a group of new grains may form [15]. Therefore, at the beginning of a layer's growth a new a set of grain nuclei is generated in the liquid at a small distance from the growing solid. (c) It is expected that the release of latent heat may alter the local temperature of the liquid at the solidification front, but it was assumed that the heat extraction rate is sufficiently large that the latent heat produced will be dominated by heat conduction through the substrate. Therefore the latent heat generation is neglected, greatly reducing the computational cost of our simulations. (d) It is assumed that varying the laser power and scanning speed will affect the overall temperature at the solidification front such that an increase in laser scanning speed may not necessarily lead to an increase in solidification rate. Therefore a set of combinations of cooling rate and temperature gradient conditions is simulated.
The following presents the details of the thermal history as they are applied to the simulation. Firstly, the cooling rate, C, temperature gradient, G, and solidification rate, R, at the steady state condition are related by where G is the temperature gradient in units of K m −1 , C is the cooling rate in units of K s −1 , and R is the solidification rate in units of m s −1 . In practice it is difficult to impose a value for R since it is a complicated function of the thermo-kinetics, interface mobility and stiffness, and thermal history. Studies that apply a pulling velocity tend to expect the velocity to approach a constant value at steady state [2,12]. Values of G in the range 10 3 ⩽ G ⩽ 10 6 and C in the range 10 2 ⩽ C ⩽ 3 × 10 6 were used. It is noted that with the phase-diagram model and fast cooling rates, the initial transient modifies the solidification rate term such that aR = C/G, where a > 1. The evolution of the thermal field is modelled by cooling, C, which is a linear function of time, as where T 0 is the initial system temperature, set to the value of the solidus, z 0 is the z coordinate of the temperature gradient zero point, where T(z 0 ) = T 0 (K) at the beginning of a simulation, t is the current time in s. Table 1 tabulates the parameters for the PF. In the absence of data for the solidification velocity under AM conditions, the measurement of a related system [21] is used.

Crystal plasticity model
In the following sections, the governing classical crystal plasticity kinematics are defined followed by the constitutive relations which control the evolution of slip.

Kinematics
The kinematics of finite deformation describe the evolution of deformation which occurs from a reference configuration to a current state based on the deformation gradient [25]. However, the plastic and elastic deformations can be separated into components by assuming the deformation gradient obeys a multiplicative decomposition, where F p is the plastic part, and F e the elastic part of the deformation gradient. The total velocity gradient tensor (L) has the following relationship, where L e and L p are the elastic and plastic parts of the velocity gradient, which themselves can be divided in symmetric parts (deformation rate tensors (D e , D p )) and asymmetric parts (spin tensors (Ω e ,Ω p )). L p can be calculated from the summation of crystallographic slip overall slip systems, where N s is the number of slip systems,γ α is the slip rate on slip systems α, s * α and m * α represent the vectors along the slip direction and normal to the slip plane of system α in the deformed configuration respectively. These vectors are also related to the reference configuration based on the following relationships, where s α and m α are the unit vectors in the slip direction and normal to the slip plane in the reference configuration respectively. An increment of shear stress can then be calculated using the following relationship, where C is the elastic stiffness tensor and σ is the Cauchy stress tensor.

Constitutive models
In this study, slip was modelled using a visco-plastic power law, whereγ 0 is a reference rate, τ α is the resolved shear stress on the slip system α, g α is the material's resistance to slip on the slip system α and m is a strain rate sensitivity parameter. Slip system hardening evolves according to the following relationship, where H αβ are the slip hardening moduli, where for latent hardening, where q is a constant which allows for the inclusion of interactions between slip systems and H αα are the self hardening moduli representing the hardening evolution on a single slip system. The initial hardness of each slip system (g α (0)) is assumed to be the value of the initial critical resolved shear stress for each slip system. The evolution of the slip hardening moduli adopted was that proposed by Peirce et al [26] and Asaro [27], where h 0 is the initial hardening modulus, τ 0 is the critical resolved shear stress and τ s represents the saturation value of flow stress. γ is the Taylor cumulative shear strain on all slip systems which evolves according to equation (13), Since this initial study applies classical crystal plasticity constitutive models, the length scale is not inherently considered. Consequently, length scale dependent deformation is not recognised, resulting in different sized grains of the same grain environment and orientation experiencing the same level of deformation. This has led to the development of other crystal plasticity theories, most notably strain gradient crystal plasticity, which has undergone significant development as reviewed by Voyiadjis and Song [28]. A length scale dependence on deformation is introduced through the association with strain gradients, which is based on the influence of strain gradients on mechanical behaviour proposed by Aifantis [29]. However, since classical crystal plasticity is being applied in this study, the modification to introduce a length scale is based on a similar approach by Weng [30] and subsequently implemented by a number of researchers [31][32][33][34][35]. In this approach, the influence of grain size on yield stress was included using the Hall-Petch relation [36,37]. This relation was included in the definition of the initial value of the critical resolved shear stress τ 0 , since this parameter has a significant influence on the predicted macro-scale yield stress. The formulation of the modification of the value of τ 0 is provided in the following equation, where τ 0d represents the critical resolved shear stress modified according to the diameter of each grain, K y is the Hall-Petch slope, and d is the equivalent spherical/circular diameter for each grain.

Simulation framework
The constitutive relations outlined in section 2.2.2 were implemented using a specifically modified version of the user-defined material subroutine (UMAT) developed by Huang [38]. The nonlinear equations were solved using the Newton-Raphson iterative method. All simulations were conducted using the fast Fourier transform (FFT) based solver AMITEX [9]. The metastable β-titanium alloy Ti-5553 has a body centred-cubic (BCC) crystal structure, which for these simulations the most active slip systems were assumed to be the slip families {110} < 111 > and {211} < 111 > with both containing twelve slip systems each. Strain controlled tensile simulations were conducted using a strain rate of 3.33 × 10 −4 s −1 . A maximum load of 10% strain was applied. For FFT simulations, periodic boundary conditions are applied. Therefore, since the predicted structure of the PF microstructure is a 2-dimensional periodic unit cell (of one voxel thickness in the y-direction), the simulations are equivalent to a two-dimensional model under a generalised plane strain hypothesis (an infinitely thick two-dimensional section).
Since the proposed approach used to simulate an AM material requires a smooth transition from PF to CP simulations, previously calibrated parameters for Ti-5553 from the literature [39] were used for the flow rule and hardening evolution. For the critical resolved shear stress, the Hall-Petch modification was adopted to introduce a length scale dependence. This introduced the Hall-Petch constant from [40]. The parameters are listed in table 2.

Integrating phase field and crystal plasticity models
The proposed approach to integrating PF and CP-FFT models makes use of the voxel geometry created by the PF model simulations and needed for FFT based solvers, which is often given in Visualization Toolkit (VTK) format.
Firstly, the resolution (the volume of voxels per grain) of the PF predicted grains was altered slightly to improve the computational efficiency of the CP simulations. This was done by extracting all voxel data from the VTK file outputted by the PF model. This data formed a matrix (A = (a ij ) ∈ R m×n ) which was then vectorised, vec(A) = [a 1,1 , ..., a m,1 , a 1,2 , ..., a m,2 , ..., a 1,n , ..., a m,n ] T (15) where T is the transpose. The vector was then reshaped to form a square matrix, where k = √ mn, and A ′ = (a ′ ij ) ∈ R k×k . The resolution of the matrix was then altered by taking every u th element (where u | k) of the row and column of the square matrix (A ′ ), where s = k/u and A r is the matrix formed from the change in resolution.
The equivalent circular diameter used in the constitutive models was calculated for each grain using the the total voxels per grain. This was achieved using the following equation, where Area was calculated from the voxels which define a grain. This information was then combined with the orientation data for each grain predicted by the PF model and used to develop an input file for CP-FFT simulations.
The procedure to alter the resolution and calculate the equivalent circular diameter of the microstructure was developed in Python [41] using the libraries NumPy [42] and Pandas [43]. The advantage of the proposed integration approach is that it allows for a seamless transition from PF to CP simulations since the process utilises the VTK format of the outputted PF files, with all internal manipulation of the data to generate input files for the CP simulations conducted through Python. This ensures a closed loop tool.
A schematic of the process used to take the PF microstructure and develop an input for CP simulations is provided in figure 2.

Generating microstructures using phase field
Using trial and error reductions of the time step, ∆t, it was found that ∆t < 10 −6 was sufficient for the given discretisation; a value of ∆t = 10 −8 was used in this study. Two dimensional grain microstructures were generated for a set of different thermal history conditions. The simulations were organised on a grid of G z vs. C to compute the solidification rate. The process was refined to address single and multilayer structures to emulate the presence of successive build layers. Combinations of G z and C were selected from best estimates of the solidification conditions. Large 2D simulations were then performed using the parameters contained in table 1. In general, all structures obtained had large columnar regions, with small regions of approximately equiaxed grains between layers, or in layers that did not have the preferred growth orientation and thus did not grow quickly enough. In the following sections, single, double, and triple layer simulations are presented. All extracted figures of the solidification do not distinguish between liquid and solid, instead show the interface, given that the liquid is at the top of the simulation boxes where the temperature is comparatively high, and the solid is growing into it from below.
All nucleation steps are restricted to the time where the solidus isotherm is moved back to the 1/2 box lengths 1/3 and 2/3 location, in that sense they are only allowed to take place at the beginning of a simulation, or if a successive layer is deposited. The creation of new grains is artificially controlled. The new grains are generated in a band in z, around the critical undercooling (approximately 2 K). However, because the grid size is larger than the critical nucleus size, the growth of new 'nuclei' depends only on their relative thermal and curvature undercooling. New grains are assigned a random orientation. This angle determines the misorientation with the thermal gradient. More details about the specific phase-field implementation can be found in [44].

Single layer solidification
Single layer simulation begins by generating a series of nucleates placed at random at the bottom of the cell (see the bottom of the left most configuration in figure 3). Each nucleus is assigned a random orientation defined by the angle of favourable growth orientation, as defined by the interface energy anisotropy, with the vertical coordinate (z). The thermal gradient is placed such that T 0 lies at the z offset, allowing a portion of the box to be below the T S , thus allowing for nucleation to occur. The system is allowed to evolve until it is fully solid. Figure 3 shows the evolution of a single layer of solid undergoing columnar growth in the direction of the thermal gradient. Large anisotropic grains dominate, and are selected by their preferred growth orientation relative to G.

Double layer solidification
Double layer solidification begins in the same way as single layer. When the total solid fraction reaches 0.5, the liquidus isotherm is reset to the simulation domain midpoint. The layer size is relatively small, approximately half the actual melt pool depth of ≈100 µm obtained in [15]. Figure 4 shows snapshots of the simulation.

Large scale triple layer microstructures
Using the data collected from the previous simulations, a set of simulations were performed on domains of size 300 × 300 µm 2 to produce layers of height ≈ 100 µm. The final simulations are collected in figure 5. Additional layers are added in the manner described above, and the system evolved until solidification was complete. While all grains obtained are generally columnar, decorated with smaller grains at layer interfaces, grains are generally more angular under high solidification velocities.
From figure 5, differences in morphology of the elongated and equiaxed grains between PF microstructures are apparent. Additionally, there are different grain environments, particularly evident from the changes in equiaxed grain clustering between layers. The morphology of the grains is dictated primarily by their relative interaction in addition to their orientation. More favourably orientated grains will grow faster than those which are not, resulting in a variety of competing interactions which contributes to the observed differences in morphology. This contributes to the spread of grain sizes which can be observed from  the grain diameter histograms in figure 6. Since each nucleus is assigned a random orientation at the beginning of layer deposition, this leads to differences in the combination of grain morphologies between the five microstructures. Additionally, this also results in slightly different textures between the simulated microstructures, which is subtly evident from the pole figures in figure 7.

Microstructure-property simulations
Using the microstructures generated by PF in figure 5, the integration approach in section 3 was applied to create the input file for CP simulations. The microstructure was discretised into 250 000 voxels. The final simulations used the parameters in table 2. The result for each microstructure is compared to the experimental monotonic tensile data from [15] in figure 8.
No noticeable difference in predicted macro-scale stress-strain evolution is apparent for each of the PF microstructures. This suggests that although the predicted microstructures contain slightly different grain environments, the average deformation is very similar. The Cauchy stress (σ 33 ) contour plot for the different PF microstructures are provided in figure 9. Regions of high stress are noticeably different between the simulated microstructures. This is the consequence of the different grain environments which exist between the predicted microstructures.
The noticeable differences in the meso-scale results can be explained by recognising the influence of differences in grain orientations between microstructures. Although the overall spread of grain sizes is very similar between the PF microstructures, the simulated orientations for each of these grains are different. This leads to slight changes in the influence of local deformation.
Deeper insight into regions of localised damage can be developed based the accumulated effective slip, which can be used to understand potential crack initiation. This is because effective plastic strain has been  well reported as a precursor for crack initiation [45,46], an approach utilised in [47,48]. The effective plastic slip rate (ṗ) can be calculated from the plastic velocity gradient (L p ), The effective plastic slip can then be calculated using the following, Application of this approach provides a means of calculating the effective plastic slip across the grains within the RVE which results in the deformation contour plots in figure 10.    We note that the magnitude of p is dependent on the parameters being applied as well as the spatial discretisation. A sensitivity analysis was conducted to understand the influence of the rate sensitivity parameter (m) and the spatial discretisation (detailed in appendix), which showed that the magnitude of p within each grain changed with the size of m. However, the locations of the maximum accumulation remained consistent. It is therefore important that such parameters which do not have a significant influence on the homogenised results are still considered for their influence on meso-scale deformation. The results also showed that increasing the spatial discretisation from that employed in this study (250 000 voxels) did not have any discernible impact on the evolution of p.
From figure 10, localised regions of p accumulate as slip bands. Additionally, the highest magnitude of p occurs at the grain boundaries of the small equiaxed grains. Comparing the five simulations, it is also evident how areas of localisation change. Although the PF simulated microstructures are almost identical in the  distribution of grain size and overall texture, what these deformation plots highlight is how slight changes in the microstructure through the combination of morphology and grain environments can produce potentially more favourable combinations to prevent crack initiation.
The statistics of the average p for each grain extracted from each microstructure is provided as histograms in figure 11. These histograms provide some further insight on the influence of the differences in grain environments within each microstructure. What is immediately evident is that there is a lognormal distribution of p in microstructure (d), while p in the other microstructures is normally distributed. The lognormal distribution of p presents a more favourable distribution across the microstructure. This is due to the peak of the distribution being located at a smaller value of p, which suggests that there is a lower magnitude of effective plastic slip. This combination of grain environments has the potential to be less susceptible to crack initiation.
The histogram for microstructure (e) is also noticeably different to the other distributions, with it being narrower with a higher peak and shorter tails. These extracted statistics suggest that p in this microstructure is significantly more localised to a select few grains. Consequently, p is not as diffuse across the microstructure.
Further validation of the findings can be acquired by considering additional statistics. The histograms of the standard deviation of p across each individual grain is provided in figure 12. Unlike figure 11 which considered the average p of each grain, these statistics consider the range of p across the grain.
Microstructures (a), (d), and (e) have similar shaped histograms. All three peaks of the histogram are at lower standard deviations of p than microstructures (b) and (c). This suggests that in these microstructures, the distribution of p is more diffuse across the grains. For microstructures (b) and (c) more localisation of p is expected with the statistics suggesting that there are more grains containing higher values in standard deviation of p. These results further support microstructure (d) as the most favourable microstructure against crack initiation.

Conclusion
This work has successfully demonstrated a means of linking PF and CP-FFT modelling to predict the microstructure and tensile behaviour of an SLM titanium alloy. This computational approach will support the development of a DT, which can be used for material design. A method for creating a microstructure indicative of that seen in SLM material has been discussed in detail. Additionally, an approach to simulate the mechanical response of a PF generated microstructure using CP-FFT was presented. From this work, the following conclusions and observations can be drawn: • The developed PF model can accurately simulate a layered microstructure representative of what has been observed experimentally in SLM Ti-5553 [15]. This included equiaxed grains between layers of large columnar grains. • An approach to linking PF and CP-FFT predictions has been developed. This method allows for the flexibility of altering the number of cells per grain (resolution) to improve computational efficiency. Additionally, this approach can be automated allowing for a resolution change in addition to CP input file creation to automatically occur after the completion of PF simulations. This ensures a seamless transition from PF to CP simulations and enables the potential of future close-loop simulations including a feedback loop were information gained from CP simulations can be fed back into PF simulations. • The simulation of multiple microstructures demonstrates the possibility of using an integrated approach to test the influence of different microstructures on the mechanical properties of the material. This includes the influence of grain sizes in addition to the combination of grain size and orientation on both the mesoand macro-scale behaviour. • By extracting and analysing meso-scale data, an approach to investigating the influence of different features of the microstructure was presented. This was done by considering how the different grain environments can influence the accumulation of effective plastic strain.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).  Additionally, there is the potential for the level of spatial discretisation to influence the evolution of p. This was also investigated by increasing the extent of spatial discretisation from 250 000 voxels to 390 625 and 693 889 voxels. A comparison of the deformation contour plots of these three different spatial discretisation levels is provided in figure A3.
From the deformation contours there seems to be no discernible difference between them. The cdf plots were also made of the average p for each grain and the standard deviation of p across each grain within each of the tested levels of spatial discretisation. The results of which are provided in figure A4.
The results of the cdf plots also confirm that increasing the spatial discretisation from that used in the analysis would not influence the distribution of p across the microstructure. Therefore, these results support the use of 250 000 voxels.