Plastic Strain Estimation for Linear Elastic Simulations

####Matthew Priddy, ME graduate student @ GT

####Noah Paulson, ME graduate student @ GT

Introduction

In the high-cycle fatigue (HCF) regime, plastic deformation plays a small role in the overall response of the material. However, fatigue indicator parameters (FIPs) are estimated from plastic quantities such as plastic shear strain. This requires the use of a material model capable of capturing plastic deformation. In the past, our group has used CPFEM models to capture local deformation (both elastic and plastic), which are computationally inefficient and not easily translatable to the computational abilities of other research groups or industry partners. Therefore, we plan to explore the process of calculating FIPs soley from elastic deformation information.

For the purpose of this explanation, the HCF regime is defined as cyclic loading that does not exceed a maximum macroscopic strain of 75% of the strain to yield.

Plastic Strain Calculation in CPFEM

In crystal plasticity finite element method simulations (CPFEM), the FEM software passes information (typically the total deformation gradient, $ \mathbf{F} $) to the crystal plasticity model which determines the total stress and plastic strain, along with other model parameters. In this particular model, the total deformation gradient is written as a multiplicative elastic-plastic decomposition, $ \mathbf{F} = \mathbf{F}^\mathrm{e}\mathbf{F}^\mathrm{p} $.

The first step is to initialize the plastic (superscript p) deformation gradient $ \mathbf{F}^\mathrm{p} $. This is done by either assuming no plastic deformation has occurred (i.e. $ \mathbf{F}^\mathrm{p} = \mathbf{I} $) or using the plastic deformation gradient from the previous loading increment. Next, the elastic deformation gradient (superscript e) can be found through the multiplicative decomposition,

$\begin{equation} \mathbf{F}^\mathrm{e} = \mathbf{F} \cdot \left(\mathbf{F}^\mathrm{p}\right)^{-1} \;. \end{equation}$

From there, the elastic Green strain tensor can be determined,

$\begin{equation} \mathbf{E}^\mathrm{e} = \frac{1}{2} \left[ \left(\mathbf{F}^\mathrm{e}\right)^\text{T} \cdot \ \mathbf{F}^\mathrm{e} - \mathbf{I} \right] \;. \end{equation}$

As well as the second Piola-Kirchhoff stress tensor (symmetric)

$\begin{equation} \boldsymbol\sigma^\text{PK2} = \mathbb{C}_0 : \mathbf{E}^\mathrm{e} \;, \end{equation}$

where $\mathbb{C}_0$ is the fourth-order elastic stiffness tensor in the intermediate configuration. The Cauchy stress can be found by mapping the second Piola-Kirchhoff stress to the current configuration, as shown by

$\begin{equation} \boldsymbol\sigma = \frac{1}{\det\left(\mathbf{F}^\mathrm{e}\right)} \left[ \mathbf{F}^\mathrm{e} \cdot \boldsymbol\sigma^\text{PK2} \cdot \left(\mathbf{F}^\mathrm{e}\right)^{\text{T}} \right] \;. \end{equation}$

The resolved shear stress, $\boldsymbol\tau$, resides in the intermediate configuration and is therefore calculated from the second Piola-Kirchhoff. The resolved shear stress on each slip system is

$\begin{equation} \tau^{(\alpha)} = \boldsymbol\sigma^\text{PK2} : \left(\mathbf{s}_0^{(\alpha)} \otimes \mathbf{n}_0^{(\alpha)}\right) \;. \end{equation}$

With a first estimate of the resolved shear stress (RSS), the slip rate on each slip system can be evaluated with the flow rule. For this work, a power-law flow rule is employed in the following form,

$\begin{equation} \dot{\gamma}^\alpha = \dot{\gamma}_0 \left< \frac{\left| \tau^\alpha - \chi^\alpha \right| - \kappa^\alpha}{ D^\alpha} \right>^{M} sign\left( \tau^\alpha - \chi^\alpha \right) \;, \end{equation}$

where $\dot{\gamma}^\alpha$ is the shearing rate on the $ \alpha $ slip system, $ \dot{\gamma}_0 $ is the reference rate of shearing, $ \chi^\alpha $ is the back stress, $ \kappa^\alpha $ is the threshold stress, $ D^\alpha $ is the drag stress, and $M$ is the rate sensitivity parameter. With this formulation, a particular slip system is defined as active if the value inside the McCauley brackets, $\langle \rangle$, is greater than zero. As a side note, the rate independent limit occurs as $ M\rightarrow\infty $.

The shearing rate, $\dot{\gamma}^{\alpha}$, is related to the plastic velocity gradient through the relation

$\begin{equation} \mathbf{L}^\mathrm{p} = \sum_{\alpha=1}^{\mathrm{N}}{\dot \gamma}^{(\alpha)}\left(\mathbf{s}_0^{(\alpha)} \otimes \mathbf{n}_0^{(\alpha)}\right) \;, \end{equation}$

where $ \hat{\mathbf{L}}^\mathrm{p} $ is the plastic velocity gradient in the intermediate configuration, $ \mathbf{s}_0^{(\alpha)} $ is the slip direction unit vector, and $ \mathbf{n}_0^{(\alpha)} $ is the slip plane normal unit vector. The plastic velocity gradient can also be defined as $ \mathbf{L}^\mathrm{p} = \mathbf{F}^\mathrm{p} \cdot \left( \dot{\mathbf{F}}^{p} \right)^{-1} $.

The rate of change of the plastic deformation gradient can be estimated from the numerical integration via

$\begin{equation} {\mathbf{F}_{t+\Delta t}^\mathrm{p} = \exp{\left( \mathbf{L}^\mathrm{p} \Delta t \right)} \cdot \mathbf{F}_{ t}^\mathrm{p}} \;. \end{equation}$

An updated estimate for the plastic velocity gradient has now been found and can be compared with the trial version. If the difference between the trial and updated is greater than some defined tolerance, further iterations are performed until said tolerance is achieved.

Estimation of Plastic Strain from Elastic Quantities

In crystal plasticity, the integration of the shearing rate, $\dot{\gamma}^{\alpha}$, is performed with an implicit or explicit integration scheme. Implicit integration schemes are typically used for quasi-static loading (i.e. small changes in deformation over a time increment) while explicit integration schemes are reserved for high-rate deformation. Further explanations can be found in (citations).

The purpose of this study is to ascertain the accuracy of the numerical integration of $\dot{\gamma}^{\alpha}$. As shown in a previous equation, it is a function of a variety of quantities. For now, lets assume $\dot{\gamma}^{\alpha} = f(\tau^{\alpha})$, meaning it is only a function of the resolved shear stress on a particular slip system. If that is the case, the numerical integration of $\dot{\gamma}^{\alpha}$ relies on providing accurate values of $\tau^{\alpha}$.

In the previous section, it is shown the resolved shear stress is a function of the second Piola-Kirchhoff stress and the slip direction and slip plane normal. In a linear-elastic simulation, the output information is the total (Cauchy) stress and total strain. However, if the deformation is minimal, the value for the Cauchy stress is approaches the value for the second Piola-Kirchhoff stress, so the resolved shear stress can be approximated as

$\begin{equation} \tau^{(\alpha)} = \boldsymbol\sigma : \left(\mathbf{s}_0^{(\alpha)} \otimes \mathbf{n}_0^{(\alpha)}\right) \;. \end{equation}$

Alternatively, the second Piola-Kirchhoff stress is related to the Cauchy stress through the elastic deformation gradient. If the plastic deformation is minimal, the elastic deformation gradient can be approximated as the total deformation gradient $(\mathbf{F}^\mathrm{e} \approx \mathbf{F})$. Also shown previously was the relation for the elastic Green strain, which is a function of the elastic deformation gradient. Unfortunately, the total deformation gradient has 9 independent components while the strain tensor only has 6 independent components, so the strain can be determined from the deformation gradient but the reverse is not possible. However, under small strain approximations, the certain properties hold true and allow for an approximation of the deformation gradient. For example, the Jacobian of the deformation gradient, $\mathrm{J}$, can be approximated as

$\begin{equation} \mathrm{J} = 1 + \epsilon_{11} + \epsilon_{22} + \epsilon_{33} \;. \end{equation}$

Additionally, under uniaxial deformation, the diagonal components of the deformation gradient can be approximated as

$\begin{equation} \mathrm{F}_{11} = 1 + \epsilon_{11}; \mathrm{F}_{22} = 1 + \epsilon_{22}; \mathrm{F}_{33} = 1 + \epsilon_{33} \;. \end{equation}$

There are still 6 unknown components of the deformation gradient, but there are also 6 independent components of the strain tensor. By writing the Green strain equation in component form, the relationship between the the deformation gradient and the strain total strain are

$\begin{equation} \begin{split} \epsilon_{11} &= \frac{1}{2} \left( \mathrm{F}^2_{11} + \mathrm{F}^2_{21} + \mathrm{F}^2_{31} - 1.0 \right) \\ \epsilon_{22} &= \frac{1}{2} \left( \mathrm{F}^2_{12} + \mathrm{F}^2_{22} + \mathrm{F}^2_{32} - 1.0 \right) \\ \epsilon_{33} &= \frac{1}{2} \left( \mathrm{F}^2_{13} + \mathrm{F}^2_{23} + \mathrm{F}^2_{33} - 1.0 \right) \\ \epsilon_{12} &= \frac{1}{2} \left( \mathrm{F}_{11} \cdot \mathrm{F}_{12} + \mathrm{F}_{21} \cdot \mathrm{F}_{22} + \mathrm{F}_{31} \cdot \mathrm{F}_{32} \right) \\ \epsilon_{13} &= \frac{1}{2} \left( \mathrm{F}_{11} \cdot \mathrm{F}_{13} + \mathrm{F}_{21} \cdot \mathrm{F}_{23} + \mathrm{F}_{31} \cdot \mathrm{F}_{33} \right) \\ \epsilon_{23} &= \frac{1}{2} \left( \mathrm{F}_{12} \cdot \mathrm{F}_{13} + \mathrm{F}_{22} \cdot \mathrm{F}_{23} + \mathrm{F}_{32} \cdot \mathrm{F}_{33} \right) \;. \end{split} \end{equation}$

Then, the plastic shear strain rate ($ \dot{\gamma}^\alpha $) is calculated from the equation below. Finally, the plastic strain can be determined from

$\begin{equation} {\boldsymbol\epsilon^\mathrm{p} = \sum_{\alpha=1}^{\mathrm{N}}{\gamma}^{(\alpha)}\left(\mathbf{s}_0^{(\alpha)} \otimes \mathbf{n}_0^{(\alpha)}\right)_{sym}} \;, \end{equation}$

where $ \gamma^{(\alpha)} = \dot \gamma^{(\alpha)} \Delta t $.

From linear-elastic FEM simulations, the total stress and total strain fields for each element (voxel) are known quantities for the entirety of the deformation. Therefore, the only unknown in this set of equations is $ \dot \gamma^{(\alpha)} $, which can be estimated through a numerical integration scheme.

The numerical integration scheme selected for this particular problem is the Simpson’s $ \frac{1}{3} $ rule,

$\begin{equation} \int_{a}^{b} f(x)dx \approx \frac{h}{6} \left( f(x_0) + 4f(x_1) + f(x_2) \right) - \frac{1}{90}h^5 f^{(4)} (\xi) \;, \end{equation}$

where $ f(x) = \dot \gamma^{(\alpha)}(\tau^{(\alpha)}) $, $ h=\Delta t $, and $ \left[a,b\right] = \left[0,t\right] $. In this form, $ \gamma^{(\alpha)} $ would be estimated over $ \Delta t $ with only three data points. Therefore, the integration interval can be divided into subincrements to reduce the error of the integration scheme. The Simpson’s $ \frac{1}{3} $ rule was selected over the Simpson’s $ \frac{3}{8} $ rule because they both attain third-order accuracy but the Simpson’s $ \frac{1}{3} $ uses one less point.

In order of completeness, the back stress and threshold stress also evolve with changes in $ \dot \gamma^{(\alpha)} $. Therefore, their evolution also needs to be captured in this numerical integration scheme.

Preliminary Results

Comparison of Linear Elastic (LE) and CPFEM Stress and Strain Fields for 0.5% Strain Amplitude

The first assumption to check is the comparison of the total stress from linear elastic simulations and CPFEM simulations, assuming the microstructure is identical between both. These simulations were performed for $\alpha$-Ti with single element grains with random orientations. Fully-reversed (R=-1) Cyclic loading was performed for 3 cycles with a strain amplitude of 0.5%. A comparison of the von-Mises total stress is shown below with very good agreement between the simulations.

Additionally, the von Mises total strain is also nearly identical between the two simulations, as shown below.

It is also important to note that the spatial variation of the stress (and strain) values also seem to be nearly identical between the simulations.

Comparison of CPFEM and LE Strain Fields for Varied Strain Amplitudes and Material Phases

The assumption of identical strains between LE and CPFEM simulations was further examined through the variation of the applied strain amplitude and the comparison of the total strain fields. Simulations were run with strain amplitudes varying between 0.25% and 2.00% with increments of 0.25% strain. Two sets of CPFEM simulations were performed with single element grains in random orientations. The first set was comprised of only $\alpha$-Ti grains while the second was completely comprised of colony grains, which mix the mechanical properties of both $\alpha$ and $\beta$ Titanium phases. While the elastic properties of each set were identical, the plastic parameters were different. All other parameters were the same as in the previous study.

The plot below compares the maximum von - Mises strain error between the CPFEM and LE results for each strain amplitude and material phase.

From this plot it is clear that the assumption of identical strains falls apart somewhere between 0.50% and 0.75% strain amplitude. The following von - Mises strain fields at 1.00% strain amplitude for colony grains show a significant difference between the CPFEM and LE simulations.

Initial Attempt at the Estimation of Plastic Strains from LE Simulations

Given the results of the previous studies the components of the total stress values were extracted at each peak and valley of the cyclic loading curve and those values, along with CPFEM model parameters, were used to estimate gamma with the numerical integration scheme outlined above. A subincrement value of 10 was used between each peak and valley in order to minimize error. Subincrement values of larger than 10 were also investigated with very minimal changes in the results. The results below indicate the variation in the spatially resolved effective plastic strain values follow the same trend for specific voxels.

However, there is a difference in the CPFEM effective plastic strain and the estimated values from numerical integration of the LE results.

Comparison of FIP Distributions for ‘Realistic’ Microstructures

The previous results have focused on the comparison of the stress, total strain, or plastic strain values between the CPFEM and LE+PY simulation results. However, the purpose for performing this type of analysis is to extend the plastic strain calculations to determine FIP values for multiple microstructure instantiations. That is exactly what we have done in this particular portion of the study.

In this section, we used Dream3D and fictitious microstructure statistics to generate 100 digital microstructures. The microstructure generation was automated through Dream3D’s new command line executable (Pipeline Runner) on MAC OSX. One example of the digital microstructures is shown below, along with some statistics (pole figure, grain size distribution, and misorientation distribution) associated with this particular digital microstrucutre.

Fully-reversed cyclic loading with a strain-amplitude of 0.5% was performed for 3 cycles on the 100 microstructure instantiations. These simulations were performed with CPFEM and LE material models. The LE simulation results were extended to determine plastic strains via the numerical integration scheme previously discussed. A comparison of the von-Mises stress and strain for the microstructure shown above indicates the CPFEM and LE material responses are identical. A histogram of the von-Mises strain values is shown below.

However, the plastic strain distribution shows a difference between the CPFEM and the numerical integration scheme.

The extreme value Fatemi-Socie FIPs for each simulation were combined to generate a Gumbel distribution of the data. From the semi-log plot, the CPFEM and LE plue numerical integration scheme results follow the same trend, but the results are quite distinct. For the purposes of this study, the data has been normalized. The actual magnitude of the values isn’t important; we are focusing on the relationship between the 2 datasets.

Discussion

From these results, it is important to note that the plastic strain values are multiple orders of magnitude smaller than the total strain, so it is believable that they would have little impact on the local total stress (and strain). However, their utility lies in the calculation of FIP values.

Additionally, we hypothesize that the difference in plastic strain values between CPFEM and the numerical integration is mainly due to slight variations in the RSS for each slip system. We are currently investigating this matter and hope to resolve it soon (no pun intended).

The next planned step is to perform the FIP analysis with Ti-64 digital microstructures, generated with Dream3D and EBSD scan data collected by Jordan. We are currently working on inputting the EBSD scan data into Dream3D for microstructure generation. Once that step is complete, the generation, preparation, simulation, and post-processing of the simulation data will be mostly automated. The automation of this process will allow us to analyze more sets of data over the coming weeks.

Appendix: General Overview of fatigue indicator parameter (FIP) calculation

The order of operations for calculating FIPs is as follows:

Calculate the difference in the plastic strain (2nd rank tensor) over a single cycle, $\Delta \epsilon^{p} = \epsilon^{p}_{max} - \epsilon^{p}_{min}$ .
Determine the eignevalues and eigenvectors for $\Delta \epsilon^{p}$ , $\Delta \epsilon^{p} \nu^{\epsilon^{p}} = \lambda^{\epsilon^{p}} \nu^{\epsilon^{p}}$ .
Calculate the maximum plastic shear strain, $\Delta \gamma^{p}_{max} = \frac{1}{2} \left| \lambda^{\epsilon^{p}}_{1} - \lambda^{\epsilon^{p}}_{3} \right|$ .
Find the critical plane for the maximum plastic shear strain, $\nu^{\epsilon^{p}}_{max} = \frac{\nu^{\epsilon^{p}}_{1} - \nu^{\epsilon^{p}}_{3}}{\| \nu^{\epsilon^{p}}_{1} - \nu^{\epsilon^{p}}_{3} \|}$ .
Calculate the maximum normal stress on the critical plane, $\sigma^{n}_{max} = ( \nu^{\epsilon^{p}}_{max} )^{T} \sigma_{max} ( \nu^{\epsilon^{p}}_{max} )$ .
Fatemi-Socie FIP, $FIPS_{FS} = \Delta \gamma^{p}_{max} ( 1 + K \frac{\sigma^{n}_{max}}{\sigma_{y}} )$ .

References

Bob Mcginty’s PhD thesis, 2001. Georgia Tech
Craig Przybyla’s PhD thesis, 2010. Georgia Tech
Numerical Methods for Engineers