Crystal Plasticity Finite Element Method for Slip Systems Evolution Analysis of α / β Duplex Titanium Alloys during Quasi-Static Tensile Testing

: The crystal plasticity ﬁnite element method, modeled on a realistic microstructure image, was developed to investigate the evolution of slip systems in grains of α / β titanium alloys during quasi-static tensile testing. By analyzing the data of slip evolution of simulation during the overall plastic deformation process, it was found that the prismatic slip systems in the α phase and the {112} < 111 > slip systems in the β phase played a leading role. By calculating the Schmid factors, it was found that the values calculated from the local stress, which was represented by major principal stress, were larger than the values calculated from the primary uniaxial tensile direction, which was due to the deviation of the local stress direction from the primary uniaxial tensile direction. Furthermore, the deviation of local stress of α phase was di ﬀ erent from that of β phase, which was related to the deformation mechanism. During the deformation, the stress and strain were concentrated in the grains of the α phase, producing a driving e ﬀ ect on the neighboring grains of the β phase. Subsequently, the incompatible deformation produced the concentration of strain at the grain / interphase boundary, thus strengthening the grain interactions and leading to the deviation.


Introduction
It has been reported that the mechanical properties of titanium and titanium alloys are greatly influenced by their microstructure [1][2][3]. In particular, the effect of the phase ratio, size and morphology of the α and β phases makes it more difficult to accurately regulate their mechanical properties. By reviewing the literature, it was found that the plastic deformation of titanium alloys at room temperature included both dislocation slip and twinning [4][5][6][7][8][9]. However, the twinning mode of duplex titanium alloys is nearly completely suppressed because of the fine microstructure and high solute contents. Therefore, dislocation slip is the dominant factor that determines the plastic deformation of duplex titanium alloys at room temperature. Dislocation slip causes slip lines and it results in a height difference on the crystal surface, with a number of slip lines constituting a slip band. Extensive studies have been devoted to slip bands in crystalline materials. For example, some studies analyzed the activation of slip systems and the slip transfer of titanium alloys at different deformation amounts based on the Electron Backscatter Diffraction (EBSD) [10,11] method. Some works have been devoted to the understanding of the interaction between slip bands and grain boundaries in commercial pure titanium evolution, Zhang et al. [18,21] adopted a duplex titanium model; however, the authors analyzed the slip evolution of the α phase while the analysis of the β phase was lacking. Furthermore, for the relationship between slip system activation and the SFs, they reported a phenomenon where the basal slip system activation occurs preferentially with lower SFs [10]. However, their work lacked further discussion and quantitative analysis. In addition, for the grain interactions, Zhang et al. [21] found the initial inter-granular misorientations greater than about 5 • can influence the subsequent micromechanical grain behavior including slip, lattice rotation and GND density. However, in this study, the initial misorientation was measured by EBSD and the data during the deformation process was lacked. In addition, the author proposed a single-phase model, which resulted in the impossibility of analyzing the grain interactions of different phases.
In this study, which was based on in-situ tensile testing of duplex titanium (Ti-5Al-2.5Cr-0.5Fe -4.5Mo-1Sn-2Zr-3Zn), a CPFEM model, which was based on a realistic microstructure image, has been developed to investigate the activation of slip systems in the α and β phases, as well as the correlation between the grain interactions and slip system evolution during the plastic deformation at room temperature. In addition, we addressed the problem of slip system evolution in duplex titanium alloys. The influence of the local stress state of grains was explained by calculating the SFs. The related work was helpful to further understand the quasi-static process of duplex titanium alloys at room temperature.

Material and Characterization Techniques
The material used in this study with the nominal composition Ti5Al2.5Cr 0.5Fe4.5Mo1Sn2Zr3Zn was cut into dog-bone-shaped specimens, which were produced with a gauge length of 15 mm and a width of 4 mm (see Figure 1b). According to quantitative statistics of grain information, the ratio of volume fraction of the α and β phases in the region of interest (ROI) was~6:4. The quasi-static tensile experiment was performed on an INSTRON 5565 with a constant strain rate of~10 −3 /s. The tensile test was stopped at strains of 2.4%, 4%, 7.8%, and 15.5% in order to image with SEM before the final fracture.

Crystal Plasticity Models
The deformation of crystalline materials is a complicated process, which is not only related to macro factors, such as loading direction, rates, deformation amount, etc., but also to some micro The EBSD test had strict requirements on the surface of the sample and it needed to be polished before the test. In this study, the sample was polished to #1500 by SiC abrasive paper and etched in a solution of 2% HF, 10% HNO 3 , and 88% H 2 O. In addition, the sample was electropolished at a voltage of 25 V for 30 s at room temperature in an electrolyte of 6% HClO 4 , 34% CH 3 (CH 2 ) 3 OH, and 60% CH 3 OH to remove the residual stress and improve the surface quality.
The ROI was 90 µm × 90 µm (see Figure 1a), which was marked by micro indentation at the center of the gauge area for the convenience of subsequent observation. Detailed microstructures were analyzed by using a FEI Nanosem 430 field emission scanning electron microscope (SEM, HITACHI S-4800N, Tokyo, Japan) operating at 20 kV and the results were analyzed by using Oxford HKL Channel 5 software.

Crystal Plasticity Models
The deformation of crystalline materials is a complicated process, which is not only related to macro factors, such as loading direction, rates, deformation amount, etc., but also to some micro factors, such as crystal orientations, microstructure and deformation mechanism. The mechanical behaviors of materials are closely related to the loading direction due to the anisotropy of elastic tensor of the material. Besides, the microstructure of materials has crystallographic orientation, which results in the different properties of polycrystalline materials under different loading direction during the elastoplastic transformation process. In the numerical simulation research based on the phase structure of alloy materials, the mechanical constitutive phase of each constituent phase can be simplified into an isotropic constitutive model. In the CPFEM, the anisotropic constitutive parameters used for different crystal orientations and the factors such as the crystal structure type, slip deformation, hardening effect, and the change of crystal orientation are considered, which make the simulation calculations more accurate.
The model was developed from an EBSD image. The local area selected was 16.5 µm × 10.1 µm (see Figure 2a). The area consisted of 12 grains, including seven α-phase grains and five β-phase grains. The grains were numbered (i.e., 1 , 2 , 3 , and so on) and the excel file including crystal structure, centroid position, Euler angles (ϕ 1 , φ, ϕ 2 ) was output by channel5 software. Based on the RGB values of the image pixels, a 2D finite element model was developed. However, the 2D model used the plane strain or plane stress assumption, which may affect the accuracy of the simulation. Therefore, a unit thickness finite element model was established with 8-node 3D hexahedron solid elements, using the microstructure-based finite element model construction software developed by our group [27] (see Figure 2b), so it is called a 2.5D finite element model [19,28], which can reduce the calculation time and well reflect the actual situation. Each grain was considered as a single crystal, with a uniform crystallographic orientation, and its Euler angle was determined by centroid result. The model was meshed into 165 × 101 × 1 elements, with each element having dimensions of 0.1 µm × 0.1 µm × 0.1 µm. To simulate the quasi-static tensile experiment, boundary conditions were defined to obtain a uniaxial average stress state while upside and downside were a traction-free surface. Then, homogeneous normal displacement, i.e., in the loading direction X, was prescribed on left side in order to correspond to strains imposed during the experiments described in Figure 2d, which produced a displacement of 2.56 µm at 1.55 × 10 5 ms. Constraints on right side are defined in order to keep the face in stationary state. Finally, a total of 155 steps of loading were performed while the total deformation was 15.5%, and the commercial software Ansys/LS-Dyna was used to simulate the deformation process. Besides, we expanded the functions of Ansys/LS-Dyna [29] and the specific process is as follow: Appl. Sci. 2019, 9, x FOR PEER REVIEW 6 of 21

Constitutive Parameters of Single Phase
The constitutive parameters of the titanium in this study were determined by synchrotron radiation carried out at the Advanced Photon Source 11-ID-C station in the Argonne National Laboratory in the US [30]. The X-ray energy (E) and its wavelength (λ) function are derived from De Broglie relation and the formula for the calculation is as follows: The Fit2D software was used to analyze the data, mainly regarding the crystal plane spacing of each phase during the tensile process by Bragg's law (2d sinθ = Nλ). Subsequently, calculating the lattice strain ε: (6) and fitting to determine the constitutive parameters of each phase. Furthermore, dhkl is the crystal plane spacing of (h k l) at different deformations and d0,hkl is the initial crystal plane spacing before loading.
In order to determine the constitutive parameters of single phase, the first step was to obtain the elastic constants of the single crystal [30]. Comparing the elastic stiffness Ehkl of experiment and selfconsistent fitting by: the fitted single-crystal elastic constants Cij (GPa) of the α and β phases was shown in Table 1. (1) Read the parameter file containing the slip systems, elastic constants, critical resolved shear stress (CRSS), hardening coefficient, crystal orientation information, and stress-strain information in the sample coordinate system.
(2) Through the rotation matrix g, the mapping relationship between the sample coordinate system and the crystal coordinate system is established, so that the input constitutive parameters can be converted during the calculation.
cos ϕ 1 cos ϕ 2 − sin ϕ 1 sin ϕ 2 cos φ sin ϕ 1 cos ϕ 2 + cos ϕ 1 sin ϕ 2 cos φ sin ϕ 2 sin φ − cos ϕ 1 sin ϕ 2 − sin ϕ 1 cos ϕ 1 cos φ − sin ϕ 1 sin ϕ 2 + cos ϕ 1 cos ϕ 2 cos φ cos ϕ 2 sin φ sin (3) Start the simulation under the current loading step and the mechanical response behavior of the model is obtained. Based on the projection relationship of the strain increment of the grain, the change of the Euler angle is calculated, which is added to the original Euler angle to represent the rotation of the grain.
(4) In the crystal coordinate system, the total stress Σσ c,ij is obtained by accumulating the stress increments calculated from the previous steps. The τ RSS is defined as: Appl. Sci. 2020, 10, 7782 6 of 20 where ϕ 1 and ϕ 2 are respectively the angles between the normal direction and the tangential direction of the slip surface and the Σσ c,ij direction, c is the crystal coordinate system. The activation of slip systems is determined by comparing with the CRSS, τ 0 . If τ RSS > τ 0 , the slip system is activated and the grain enters the plastic deformation stage. With the activation of slip systems, the CRSS changed because of the strain hardening: where τ i is the increment of CRSS, γ j is the shear strain rate of the slip system j, Γ is the accumulating increment of shear strain of the slip system, h ij is the hardening coefficient. (5) In this simulation, it is considered that plastic deformation is caused by slip and the plastic strain ε p c is the cumulative value of resolved shear strain of all slip systems. According to the state of the grain, it is determined that the grain is in the elastic stage or plastic stage. After that, the instantaneous stiffness matrix is obtained, and the total stiffness coefficient matrix can be calculated: where L is the stiffness coefficient matrix, p and e represent respectively in the plastic stage and elastic stage. (6) Use the updated constitutive parameters to continue the next step of the finite element calculation until the last step.

Constitutive Parameters of Single Phase
The constitutive parameters of the titanium in this study were determined by synchrotron radiation carried out at the Advanced Photon Source 11-ID-C station in the Argonne National Laboratory in the US [30]. The X-ray energy (E) and its wavelength (λ) function are derived from De Broglie relation and the formula for the calculation is as follows: The Fit2D software was used to analyze the data, mainly regarding the crystal plane spacing of each phase during the tensile process by Bragg's law (2d sinθ = Nλ). Subsequently, calculating the lattice strain ε: and fitting to determine the constitutive parameters of each phase. Furthermore, d hkl is the crystal plane spacing of (h k l) at different deformations and d 0,hkl is the initial crystal plane spacing before loading. In order to determine the constitutive parameters of single phase, the first step was to obtain the elastic constants of the single crystal [30]. Comparing the elastic stiffness E hkl of experiment and self-consistent fitting by: the fitted single-crystal elastic constants C ij (GPa) of the α and β phases was shown in Table 1. Table 1. Single-crystal elastic constants C ij (GPa). Furthermore, the modified Vocé equation [30] is used to calculate the strain hardening effect of the slip system, i, in a single crystal: where Γ is the accumulating increment of shear strain of all slip systems during deformation, τ is the transient flow stress. The deformation at early stage is characterized by CRSS (τ 0 ) and hardening coefficient θ 0 . The deformation at later stage is characterized by (τ 0 + τ 1 ) and θ 1 . The values of CRSS and hardening coefficient used in this study were shown in Table 2. Table 2. The critical resolved shear stress (CRSS) and hardening coefficient of slip systems [30].

Misorientation Angle
The misorientation angle between two objects is defined as the smallest rotation angles among equivalent rotations relating two given orientations of the objects, which is generally used to characterize the rotation of the grains. For duplex titanium alloys, because the α phase has a hexagonal lattice structure, it is necessary to convert the three-axis coordinate system (hkil)[uvtw] to the four-axis coordinate system (HKL) [UVW]: crystal f ace index : (HKL) → (hkil), and The misorientation matrix M is defined as: where g 1 is the orientation matrix before deformation and g 2 is the orientation matrix after deformation. Furthermore, the real misorientation matrix M' is calculated by M and the symmetry operator T: Appl. Sci. 2020, 10, 7782 8 of 20 There are 12 symmetry operator matrices for the α phase and 24 symmetry operator matrices for the β phase. Therefore, the M' has a number of equivalent solutions and the smallest one was chosen in this study. Finally, the misorientation angle θ was calculated: For the misorientation angles between the α and β phases, the (0001) plane of hexagonal crystal structure is parallel to the (011) plane of the cubic crystal structure and the 121 0] crystal direction of the hexagonal crystal structure is parallel to the 111 crystal direction of the cubic crystal structure. As a result, the (0001) plane of the α phase and the (011) plane of the β phase were determined in this study to calculate the misorientation angles.

Resolved Shear Strain
Based on the crystal plasticity theory [31,32] and the Schmid's law [33], it is considered that when the shear stress of the slip system reaches the critical resolved shear stress (CRSS), the slip system becomes active: where m is the SF, τ s rss and γ s rss are the shear stress and the shear strain, respectively, σ and ε are the axial stress and the axial strain in the crystal coordinate system, respectively. In this simulation, a total of 155 steps of loading were performed and each step provided strain increment, . ε (the rotations were disregarded): where . γ s rss is the resolved shear strain increment. Therefore, the shear strain of slip systems at step i, γ s i , can be calculated by accumulating . γ s i (see Figure 3): Appl. Sci. 2019, 9, x FOR PEER REVIEW 9 of 21

Determination of Slip Trace
The slip trace was the intersecting line between the slip plane and the grain surface (the X-Y plane in this study), as shown in Figure 4. Fitting the slip systems and the rotation matrix g can transfer the slip plane from the crystal coordinate system to the sample coordinate system. A slip

Determination of Slip Trace
The slip trace was the intersecting line between the slip plane and the grain surface (the X-Y plane in this study), as shown in Figure 4. Fitting the slip systems and the rotation matrix g can transfer the slip plane from the crystal coordinate system to the sample coordinate system. A slip trace can be determined when the slip plane intersects with the grain surface and the angle between the simulation slip trace and X-axis is calculated, which should be compared with the experiment result.

Determination of Slip Trace
The slip trace was the intersecting line between the slip plane and the grain surface (the X-Y plane in this study), as shown in Figure 4. Fitting the slip systems and the rotation matrix g can transfer the slip plane from the crystal coordinate system to the sample coordinate system. A slip trace can be determined when the slip plane intersects with the grain surface and the angle between the simulation slip trace and X-axis is calculated, which should be compared with the experiment result.

Reliability Verification
To confirm the reliability of the quasi-static tensile simulation of duplex titanium alloys, the true stress-strain curves obtained by simulation and experiment are shown in Figure 5. The values of relative error are calculated by (Sim. − Exp.)/Exp. × 100%. When the macro strain was 2.4%, 4%, 7.8%, and 15.5%, the values of relative error are −2.63%, 0.25%, 1.39%, and 4.92%, respectively. It was found that the results of the CPFEM were in good agreement with the experimental results.

Reliability Verification
To confirm the reliability of the quasi-static tensile simulation of duplex titanium alloys, the true stress-strain curves obtained by simulation and experiment are shown in Figure 5. The values of relative error are calculated by (Sim. − Exp.)/Exp. × 100%. When the macro strain was 2.4%, 4%, 7.8%, and 15.5%, the values of relative error are −2.63%, 0.25%, 1.39%, and 4.92%, respectively. It was found that the results of the CPFEM were in good agreement with the experimental results. The grains of No.2, No.5, No.6, and No.9 (see Figure 2b) with obvious slip traces were selected and the activation of slip systems was analyzed statistically at the deformation of 15.5%. The gray figures were SEM images, and the color figures are simulation results. The information provided on sample surface by slip traces may be misleading if the "visibility coefficient", that is the inclination of the Burgers vectors on that surface, is not carefully examined. Therefore, errors may exist in the Figure 6. In spite of that, we believe that the results are still acceptable, because the slip traces of experiment and simulation agree very well with each other. Subsequently, the angles of slip traces surface by slip traces may be misleading if the "visibility coefficient", that is the inclination of the Burgers vectors on that surface, is not carefully examined. Therefore, errors may exist in the Figure 6. In spite of that, we believe that the results are still acceptable, because the slip traces of experiment and simulation agree very well with each other. Subsequently, the angles of slip traces were calculated by the method described in Section 2.2.5, which were compared with experiment (see Figure 6). It was found that the maximum error was 9.825 • , as determined by (0001)[1210] of grain No. 6, while its simulation angle was 30.175 • and the experimental value was 40 • . In addition, the minimum error was 0.446 • , which was determined by (1010)[1210] of grain No. 9, while the simulation angle was 109.654 • and the experiment was 110 • . Furthermore, by calculating statistically, the average error was only 4.667 • . In conclusion, the simulation results were in good agreement with the experimental results.

Evolution of Stress and Strain
To further analyze the evolution of stress and strain, the Von Mises effective stress and effective strain field at macro strains of 0.1%, 0.7%, 1.2%, 4%, 7.8%, and 15.5% are shown in Figure 7. Considering the needs of comparison and analysis, the figures were scaled to the same size. It can be found that the tendencies of Von Mises effective stress and effective strain are somewhat different. At the beginning, the initial stress and strain were relatively small. As the strain increases, the stress also increases. When entering the plastic stage, the stress value increases to a certain value and remains unchanged, while the accumulated strain continues to increase.
As can be seen from Figure 7a-d, the inhomogeneous distribution of stress became more evident with the development of deformation. It was found that the maximum stress was 333 MPa and the minimum was 162 MPa at a macro strain of 0.7%. However, at the macro strain of 4%, there were some grains whose internal stress was higher than 1000 MPa and some grains whose internal stress was lower than 300 MPa. Furthermore, a high stress concentration was observed at the grain/interphase boundary, which may be related to the incompatible deformation of grains. As can be seen from Figure 7e,f, the grains at low stress maintained the same status.

Evolution of Stress and Strain
To further analyze the evolution of stress and strain, the Von Mises effective stress and effective strain field at macro strains of 0.1%, 0.7%, 1.2%, 4%, 7.8%, and 15.5% are shown in Figure 7. Considering the needs of comparison and analysis, the figures were scaled to the same size. It can be found that the tendencies of Von Mises effective stress and effective strain are somewhat different. At the beginning, the initial stress and strain were relatively small. As the strain increases, the stress also increases. When entering the plastic stage, the stress value increases to a certain value and remains unchanged, while the accumulated strain continues to increase. According to the statistics, when the macro strain was less than 0.7%, the average stress of the α phase increased by ~120 MPa for every 0.1% macro-strain increment, while the average stress of the β phase increased by ~104 MPa. When the macro strain was 0.7%, the average stress of the α phase was 720 MPa and the average stress of the β phase was 624 MPa. When the macro strain was 1.2%, the average stress of the α phase was 1035 MPa and the average stress of the β phase was 1055 MPa, which meant the β phase became the main stress carrier. This was mainly due to the fact that when the macro strain was over 0.7%, the slip systems of the α phase were activated and the simulation entered the stage of plastic deformation. As a result, the ability of the α phase to bear stress decreased and the stress transferred from the α phase into the β phase gradually.
After reaching the macro strain of 0.7%, as shown in Figure 7g,h, the simulation was in the stage of elastic deformation, where the strain within the grains is a constant of 0.7%. After reaching the macro strain of 1.2%, as shown in Figure 7i, the material began to deform by slip and the strain is As can be seen from Figure 7a-d, the inhomogeneous distribution of stress became more evident with the development of deformation. It was found that the maximum stress was 333 MPa and the minimum was 162 MPa at a macro strain of 0.7%. However, at the macro strain of 4%, there were some grains whose internal stress was higher than 1000 MPa and some grains whose internal stress was lower than 300 MPa. Furthermore, a high stress concentration was observed at the grain/interphase boundary, which may be related to the incompatible deformation of grains. As can be seen from Figure 7e,f, the grains at low stress maintained the same status.
According to the statistics, when the macro strain was less than 0.7%, the average stress of the α phase increased by~120 MPa for every 0.1% macro-strain increment, while the average stress of the β phase increased by~104 MPa. When the macro strain was 0.7%, the average stress of the α phase was 720 MPa and the average stress of the β phase was 624 MPa. When the macro strain was 1.2%, the average stress of the α phase was 1035 MPa and the average stress of the β phase was 1055 MPa, which meant the β phase became the main stress carrier. This was mainly due to the fact that when the macro strain was over 0.7%, the slip systems of the α phase were activated and the simulation entered the stage of plastic deformation. As a result, the ability of the α phase to bear stress decreased and the stress transferred from the α phase into the β phase gradually.
After reaching the macro strain of 0.7%, as shown in Figure 7g,h, the simulation was in the stage of elastic deformation, where the strain within the grains is a constant of 0.7%. After reaching the macro strain of 1.2%, as shown in Figure 7i, the material began to deform by slip and the strain is locally concentrated. In addition, as shown in Figure 7j-l, it can be seen that the localized strain concentration developed with increasing macro strain. Firstly, after reaching the macro strain of 4%, the maximum strain is 11% and the minimum is 2%. At the macro strain of 7.8%, there are regions where the strain is less than 4% and regions where the strain is more than 16%. Finally, when the macro strain was 15.5%, there was no region where the strain was less than 4% and the region where the strain was greater than 16% increases further. The localized concentration mainly related to the grain orientation changes and the grain interactions, which will be specifically analyzed in Sections 3.3.2 and 3.3.3.

Activation of Slip Systems
In order to further study the grain interactions of different phases during the plastic deformation, it was necessary to analyze the evolution of slip systems. In this simulation, a total of 155 load steps (macro strain = 15.5%) were carried out to analyze the data of the slip systems, as shown in Figure 8, where the relative frequency meant the ratio of the number of the activated slip systems to the total number activated. In this figure, the relative frequencies showed a tendency of saturation after the macro strain of 2.4%, but their absolute quantity increased, which was similar to the experimental results reported by Liu et al. [34]. On the basis of the data of the α phase, it is found that the deformation was mainly determined by basal and prismatic slip systems, which was because of the low critical resolved shear stress. With the development of deformation, the pyramidal slip systems activated and the number of pyramidal <c + a> slip systems significantly increased. According to the statistics, after reaching the macro strains of 2.4%, 4%, 7.8%, and 15.5%, the relative frequencies of pyramidal <c + a> slip systems were 31%, 32%, 32%, and 32%, respectively, which was the highest proportion of the corresponding deformation of the α phase.
By analyzing the slip systems of the β phase, it was found that the slip systems became activated when the step was nine (macro strain = 0.9%), where the relative frequency of {110} <111> was 17% and that of {110} <111> was 8%. With the development of deformation, the relative frequency of {112} <111> gradually increased. When the macro strain was 15.5%, the relative frequency of {110} <111> was 18% and that of {112} <111> was 15%. Therefore, the {110} <111> was the highest proportion during plastic deformation of the β phase. In addition, the proportion of the relative frequency of the β phase increased with increasing macro strain, and it was concluded that the grains of the β phase distorted themselves to adapt to the overall plastic deformation.
In order to determine the slip systems that played a leading role during the plastic deformation of this simulation, the contribution of each slip system to the plastic deformation was further analyzed. When the macro strain was 15.5%, the γ 155 s of the α phase is shown in the bar chart of Figure 9 In order to determine the slip systems that played a leading role during the plastic deformation of this simulation, the contribution of each slip system to the plastic deformation was further analyzed. When the macro strain was 15.5%, the  According to the data in Figure 9, when the macro strain was 15.5%, the prismatic and basal slip systems accounted for 86% and 14%, respectively, of the slip systems with the highest while the relative frequency in Figure 8 was 32%, which was more than others. In order to determine  In order to determine the slip systems that played a leading role during the plastic deformation of this simulation, the contribution of each slip system to the plastic deformation was further analyzed. When the macro strain was 15.5%, the  According to the data in Figure 9, when the macro strain was 15.5%, the prismatic and basal slip systems accounted for 86% and 14%, respectively, of the slip systems with the highest  while the relative frequency in Figure 8 was 32%, which was more than others. In order to determine According to the data in Figure 9, when the macro strain was 15.5%, the prismatic and basal slip systems accounted for 86% and 14%, respectively, of the slip systems with the highest γ 155 s of each grain in the α phase. In addition, the γ 155 s of the 011 0 2110 slip system was 38.97 × 10 −3 , which was much higher than others in grain No. 1. This suggested that the 0110 2110 slip system was suitable for slip for grain No. 1. For pyramidal<c + a> slip systems, they had lower γ 155 s , while the relative frequency in Figure 8 was 32%, which was more than others. In order to determine the leading role of slip systems during the plastic deformation, the relative frequency and γ 155 s were combined to analyze the contribution of slip systems: where γ cum is the sum of the γ 155 s of the slip system, s. According to the line chart of Figure 9, the γ cum values of the prismatic slip systems were higher than others, especially the γ cum that was 40.15 × 10 −3 in grain No. 1, which means the prismatic slip systems play a leading role during the plastic deformation. Furthermore, the highest γ cum of the pyramidal <a> slip system was 1.466 × 10 −3 in grain No. 9, which was significantly lower than other slip systems. Meanwhile, it was found that the number of activated slip systems was also less than others. Therefore, the contribution of pyramidal<a> to the plastic deformation was the least.
When the macro strain was 15.5%, the γ 155 s and γ cum of the β phase are shown in Figure 10. As we can see, some individual{112} <111> slip systems, such as (211) [111] and 121 111 in grain No. 2, the γ 155 s values were 12.523 × 10 −3 and 13.271 × 10 −3 , respectively, which were significantly higher than others in the same grain. According to the diagram, it was found that all γ cum values of the {112} <111> slip systems were higher than the γ cum of {110} <111>, which meant the contribution of {112} <111> to the plastic deformation was much larger than {110} <111>. Therefore, we can draw the conclusion that the {110} <111> slip systems were easily activated because of the lower critical resolved shear stress, while the contribution of them to the plastic deformation is lower than {112} <111>. On the basis that some γ 155 s of the {112} <111> slip systems are much higher than γ 155 s of {110} <111> slip systems, it was considered that the orientation of some of {112} <111> slip systems was conducive to the plastic deformation, as also reported in the literature [14].
higher than others, especially the cum γ that was 40.15 × 10 −3 in grain No. 1, which means the prismatic slip systems play a leading role during the plastic deformation. Furthermore, the highest cum γ of the pyramidal <a> slip system was 1.466 × 10 −3 in grain No. 9, which was significantly lower than other slip systems. Meanwhile, it was found that the number of activated slip systems was also less than others. Therefore, the contribution of pyramidal<a> to the plastic deformation was the least.
When the macro strain was 15.5%, the 155 s γ and cum γ of the β phase are shown in Figure 10.
As we can see, some individual{112} <111> slip systems, such as (21 1) γ values were 12.523 × 10 −3 and 13.271 × 10 −3 , respectively, which were significantly higher than others in the same grain. According to the diagram, it was found that all cum γ values of the {112} <111> slip systems were higher than the cum γ of {110} <111>, which meant the contribution of {112} <111> to the plastic deformation was much larger than {110} <111>. Therefore, we can draw the conclusion that the {110} <111> slip systems were easily activated because of the lower critical resolved shear stress, while the contribution of them to the plastic deformation is lower than {112} <111>. On the basis that some

Analysis of Schmid Factors
According to the Schmid's law [35,36], it was believed that when the CRSS was reached, some specific crystal planes would slip along a specific direction, that is, when the shear stress on the slip system reached a critical value, the plastic deformation began. In this study, we adopted the uniaxial tensile simulation, and the SFs were analyzed. For grains with different orientations, the SFs were calculated in the crystal coordinate system and the primary uniaxial tensile direction should be transferred from sample coordinate system to the crystal coordinate system by rotation matrix g (see

Analysis of Schmid Factors
According to the Schmid's law [35,36], it was believed that when the CRSS was reached, some specific crystal planes would slip along a specific direction, that is, when the shear stress on the slip system reached a critical value, the plastic deformation began. In this study, we adopted the uniaxial tensile simulation, and the SFs were analyzed. For grains with different orientations, the SFs were calculated in the crystal coordinate system and the primary uniaxial tensile direction should be transferred from sample coordinate system to the crystal coordinate system by rotation matrix g (see Equation (1)). After that, the SFs were calculated by the slip system and the transferred loading direction. Therefore, the activation of slip systems in different grains was determined by the SFs.
The distribution of SFs, corresponding to the macro strains of 0.7%, 0.9%, 2.4%, 4%, 7.8%, and 15.5%, is shown in Figure 11. According to the data of α phase, it was shown that the prismatic slip systems with the SFs in a range of 0.1-0.2 activated when the macro strain = 0.7% and 0.9%. When the macro strain reached 2.4%, the SFs lower than 0.1 of pyramidal<a> slip systems activated. In comparison, the SFs of β phase were much larger, especially at the macro strain = 0.9% and 2. systems with the SFs in a range of 0.1-0.2 activated when the macro strain = 0.7% and 0.9%. When the macro strain reached 2.4%, the SFs lower than 0.1 of pyramidal<a> slip systems activated. In comparison, the SFs of β phase were much larger, especially at the macro strain = 0.9% and 2.4%, it was shown that all the SFs were in the range of 0.4-0.5 when the slip systems activated. When the macro strain reached 4%, there existed SFs lower than 0.1 of {110} <111> and SFs in the range of 0.1-0.2 of {112} <111>. By comparing the difference of the SFs of the α and β phases, it was found that the SFs of the α phase were lower than the SFs of the β phase on the whole, and there was a situation that the slip systems, with lower SFs in the range of 0.1-0.2, became activated preferentially, which was also described in the literature [10]. However, according to the Schmid's law, the slip system with the higher SF can be activated easily, contrary to the result above. On this basis, an assumption was proposed: the grains of α phase were deviated from uniaxial stress due to grain interactions, which meant the local realistic loading direction of grains of α phase, deviated from the predominant loading direction.
Generally, a resolved shear stress is determined from the full Schmid tensors unless the stress state is uniaxial, what is unlikely for a local stress state. In this study, we analyze the major principal stress to verify the assumption and prove the influence of grain interactions on the SFs quantitatively. For example, when the macro strain reached 0.7%, the prismatic slip systems of grain No. 3 and 11 By comparing the difference of the SFs of the α and β phases, it was found that the SFs of the α phase were lower than the SFs of the β phase on the whole, and there was a situation that the slip systems, with lower SFs in the range of 0.1-0.2, became activated preferentially, which was also described in the literature [10]. However, according to the Schmid's law, the slip system with the higher SF can be activated easily, contrary to the result above. On this basis, an assumption was proposed: the grains of α phase were deviated from uniaxial stress due to grain interactions, which meant the local realistic loading direction of grains of α phase, deviated from the predominant loading direction.
Generally, a resolved shear stress is determined from the full Schmid tensors unless the stress state is uniaxial, what is unlikely for a local stress state. In this study, we analyze the major principal stress to verify the assumption and prove the influence of grain interactions on the SFs quantitatively. For example, when the macro strain reached 0.7%, the prismatic slip systems of grain No. 3 and 11 were activated, while the SF of 0110 2110 in grain No. 3 was 0.136 and the SF of 0110 2110 in grain No. 11 was 0.321. It was found that the 0110 2110 in grain No. 3 should not be the primary slip system. Therefore, as shown in Figure 12a, the macro deformation direction was parallel to X axis, and the direction of major principal stress δ 1 deviated from the X-axis and the deviation angle was θ (ignoring the deviation angle between the direction of δ 1 and Z-axis), which was used to fit with the Euler angle and recalculate the SFs. As shown in Figure 12b, the θ of grain No. 3 was 30 • and the θ of grain No. 11 was close to 0 • . After recalculation, the SF of 0110 2110 in grain No. 3 was 0.333, which meant that the 0110 2110 slip system was easily activated. direction was parallel to X axis, and the direction of major principal stress δ1 deviated from the X-axis and the deviation angle was θ (ignoring the deviation angle between the direction of δ1 and Z-axis), which was used to fit with the Euler angle and recalculate the SFs. As shown in Figure 12b, the θ of grain No. 3 was 30° and the θ of grain No. 11 was close to 0°.  Therefore, the recalculated SFs were shown in Figure 13. By comparing Figure 11 and Figure 13, it was found that SFs overall were higher in Figure 13 than that in Figure 11, especially for the values of prismatic slip systems and pyramidal<a> slip systems. Therefore, the recalculated SFs were shown in Figure 13. By comparing Figures 11 and 13, it was found that SFs overall were higher in Figure 13 than that in Figure 11, especially for the values of prismatic slip systems and pyramidal<a> slip systems. The results proved the assumption that some grains of the α phase were deviated from uniaxial stress at the early stage of deformation, which resulted in the deviation of the local stress from the primary uniaxial tensile direction. Figure 12b-d shows the direction of major principal stress and the effective strain increment The results proved the assumption that some grains of the α phase were deviated from uniaxial stress at the early stage of deformation, which resulted in the deviation of the local stress from the primary uniaxial tensile direction. Figure 12b-d shows the direction of major principal stress and the effective strain increment fields at macro strains of 0.7%, 2.4%, and 7.8% (the effective strain increment was defined as the increment of effective strain of elements for every increase of 0.1% of macro strain and to show the effect of current step). As shown in Figure 12b, when the macro strain reached 0.7%, the grains of the α phase deviated from uniaxial stress, while the grains of the β phase kept the uniaxial stress. This was mainly due to the fact that the effective strain increment was concentrated on the grains of the α phase, resulting in the inconsistent deformation and the grain interactions, which forced the grains to deviate from the uniaxial stress. Then, as shown in Figure 12c, the effective strain increment of the β phase increased and the direction of major principal stress slightly changed, while it was found that the SFs of the activated slip systems in the β phase were in the range of 0.4-0.5 (see Figure 11). As shown in Figure 12d, the deviation of the β phase further increased, while SFs smaller than 0.1 appeared (see Figure 11), which meant the deviation influenced the activation of slip systems in the β phase.

Analysis of Grain Interactions
In order to study the grain interactions during deformation, the average grain misorientation angles were statistically analyzed (the calculation method is described in Section 2.2.3. As shown in Figure 14a,b, the curves were divided into three stages: macro strain less than 2.4% is stage I; macro strain between 2.4% and 4% is stage II; macro strain more than 4% is stage III. In stage I, the average misorientation angles of the α/α phase decreased, while those of the β/β and α/β phases remain invariant. In stage II, the average misorientation angles of the α/α and α/β phases decreased, while those of the β/β phase remain invariant. In stage III, the average misorientation angles of the α/α and α/β phases decreased, while those of the β/β phase began to increase. Therefore, we can draw the conclusion that the α phase produced plastic deformation first in stage I. When the deformation reached a certain degree (entering stage II), the grains of the α phase forced the neighboring grains of the β phase to deform and rotate. After that, when the deformation caused by "the driving effect" reached a certain degree (entering stage III), the average misorientation angles of the β/β phase began to increase to adapt to the overall deformation.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 19 of 21 invariant. In stage II, the average misorientation angles of the α/α and α/β phases decreased, while those of the β/β phase remain invariant. In stage III, the average misorientation angles of the α/α and α/β phases decreased, while those of the β/β phase began to increase. Therefore, we can draw the conclusion that the α phase produced plastic deformation first in stage I. When the deformation reached a certain degree (entering stage II), the grains of the α phase forced the neighboring grains of the β phase to deform and rotate. After that, when the deformation caused by "the driving effect" reached a certain degree (entering stage III), the average misorientation angles of the β/β phase began to increase to adapt to the overall deformation.

Conclusions
In this study, the crystal plasticity finite element method (CPFEM), modeled on a α/β realistic microstructure image, was proposed to simulate the quasi-static tensile process to investigate the evolution of slip systems. The following conclusions can be drawn: (1) A method of slip trace determination based on the CPFEM was proposed. By comparing the angles between the slip trace and the X-axis of the simulation and experiment, the average error was only 4.667°, indicating good reliability of the method.
(2) The mechanism of the slip deformation of an α/β duplex titanium alloy was revealed by analyzing the evolution of slip systems. In the α phase, the prismatic slip systems were activated preferentially with the highest contribution to the plastic deformation. In the β phase, by comparing

Conclusions
In this study, the crystal plasticity finite element method (CPFEM), modeled on a α/β realistic microstructure image, was proposed to simulate the quasi-static tensile process to investigate the evolution of slip systems. The following conclusions can be drawn: (1) A method of slip trace determination based on the CPFEM was proposed. By comparing the angles between the slip trace and the X-axis of the simulation and experiment, the average error was only 4.667 • , indicating good reliability of the method.
(2) The mechanism of the slip deformation of an α/β duplex titanium alloy was revealed by analyzing the evolution of slip systems. In the α phase, the prismatic slip systems were activated preferentially with the highest contribution to the plastic deformation. In the β phase, by comparing the sum of contribution to the deformation, it was found that the values of {112} <111> slip systems were much larger than the values of {110} <111> slip systems.
(3) By calculating the Schmid factors, it was found that the values calculated from local stress were larger than the values calculated from primary uniaxial tensile direction, which was due to the deviation of the local stress from the primary uniaxial tensile direction, while the local stress of the α phase more easily deviated than that of the β phase.