Numerically efficient microstructure-based calculation of internal stresses in superalloys

Ni-base and Co-base superalloys are known for their outstanding creep resistance at high temperatures. The creep behavior and corresponding mechanisms related to the dislocation motion in the microstructure of these alloys have been well investigated experimentally by different research groups, such as the general introduction of creep and dislocation structures in some typical superalloys [1–3], the formation of dislocation networks on γ/γ' interfaces at high temperatures [4], the orientation dependence of creep under different loading conditions [5, 6], as well as the shearing of $\gamma ^{\prime}$ precipitates with stacking faults [7]. In addition, through comprehensive dislocation dynamics modeling, the important creep mechanisms were clarified and validated [8]. Since the creep behavior is attributed to the particular microstructure where softer face centered cubic γ matrix is strengthened by stronger cuboidal or spherical L1₂-ordered $\gamma ^{\prime}$ precipitates [9], the crystal plasticity modeling of creep deformation should be also based on this microstructure. Due to the γ/γ' lattice mismatch and the inhomogeneous deformation in γ matrix and $\gamma ^{\prime}$ precipitate during creep, internal stresses are generated, which influences further deformation and rafting of precipitates [10–14]. Although several dislocation-based constitutive models have taken the importance of internal stresses into account [15–17], their calculations mainly depend on the common microstructure with narrow γ channels and cubic $\gamma ^{\prime}$ precipitate. To overcome this restriction, a new method is developed here to calculate internal stresses in γ/γ' microstructure of superalloys numerically, by generalization of Eshelby inclusion problem [18, 19] and identification of regions with uniform deformation behavior. This method can be used not only for cubic $\gamma ^{\prime}$ precipitates, but also for different morphologies, such as spherical $\gamma ^{\prime}$ precipitates. It is worthy to note that the calculated internal stresses in the present work result from the eigenstrains in different regions which are sum of the accumulated plastic strains and misfit strains at a certain time, where the microstructure and the deformation state are known. During the deformation, the internal stresses can be updated by using this method in every time step. This method can be easily implemented in crystal plasticity finite element (FE) models to consider the evolution of internal stresses while simulating creep or plastic deformation of larger structures or components. In our work we only solve the linear elastic problem relating the eigenstrains to internal stresses, but we note for the sake of completeness that recent developments on spectral solvers also allow them to be used for finite deformations and nonlinear materials [20–22].

Under external load the lattice mismatch in superalloys leads to the inhomogeneous deformation in different channels [13]. The internal stresses caused by the misfit strain and the deformation heterogeneity play an important role on the deformation. If the elastic and plastic strains are known, the total stress ${{\boldsymbol{\sigma }}}_{\mathrm{tot}}$ including the external stress ${{\boldsymbol{\sigma }}}_{\mathrm{ext}}$ and the internal stress ${{\boldsymbol{\sigma }}}_{\mathrm{int}}$ are given as

$$\begin{eqnarray}&&{{\boldsymbol{\sigma }}}_{\mathrm{tot}}={{\boldsymbol{\sigma }}}_{\mathrm{ext}}+{{\boldsymbol{\sigma }}}_{\mathrm{int}}={\mathbb{C}}:{{\boldsymbol{\varepsilon }}}_{\mathrm{el}}+{{\mathbb{C}}}^{{\prime\prime} }:({{\boldsymbol{\varepsilon }}}_{\mathrm{pl}}+{{\boldsymbol{\varepsilon }}}_{\mathrm{mis}}),\end{eqnarray} \tag{ 1 }$$

where ${\mathbb{C}}$ is the elastic stiffness matrix, ${{\boldsymbol{\varepsilon }}}_{\mathrm{el}}$ is the elastic strain, ${{\boldsymbol{\varepsilon }}}_{\mathrm{pl}}$ is the plastic strain, ${{\boldsymbol{\varepsilon }}}_{\mathrm{mis}}$ is the misfit strain, and ${{\mathbb{C}}}^{{\prime\prime} }$ is an assumed effective stiffness matrix connecting the internal stress with plastic strain and misfit strain. In this paper, we aim to determine ${{\mathbb{C}}}^{{\prime\prime} }$, so as to calculate the internal stress through it.

In view of the characteristic γ/γ' microstructure of superalloy single crystals, a cubic representative volume element (RVE), similar to that used in [13, 23, 24], containing one cubic $\gamma ^{\prime}$ precipitate surrounded by narrow γ matrix channels, is utilized to determine internal stresses. With known eigenstrains in the RVE with periodic boundary conditions, the internal stresses can be calculated by the fast Fourier transformation (FFT) method [13, 25–27]. The RVE is treated as a material with a reference stiffness $\bar{{\mathbb{C}}}$, which is a volume average of local elastic stiffnesses. Analogous to the solution of the Eshelby inclusion problem [18, 19], by assigning the eigenstrain $\bar{{\boldsymbol{\varepsilon }}}$ to the local strain, the internal stress on each FFT grid can be obtained by

$$\begin{eqnarray}&&{{\boldsymbol{\sigma }}}_{\mathrm{int}}={\bar{{\boldsymbol{\sigma }}}}_{\mathrm{int}}+{\tilde{{\boldsymbol{\sigma }}}}_{\mathrm{int}},\end{eqnarray} \tag{ 2 }$$

where the mean internal stress ${\bar{{\boldsymbol{\sigma }}}}_{\mathrm{int}}$ is the volume average of ${{\boldsymbol{\sigma }}}_{\mathrm{int}}$. The local internal stress fluctuation ${\tilde{{\boldsymbol{\sigma }}}}_{\mathrm{int}}$ is identical to the constraint stress caused by the corresponding strain fluctuation $\tilde{{\boldsymbol{\varepsilon }}}$:

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{\sigma }}}}_{\mathrm{int}}=\bar{{\mathbb{C}}}({\bf{E}}-{\bf{I}})\tilde{{\boldsymbol{\varepsilon }}},\end{eqnarray} \tag{ 3 }$$

where ${\bf{E}}$ is the Eshelby tensor and ${\bf{I}}$ is the unit tensor. We extend equation (3) to consider the interaction among FFT grids

$$\begin{eqnarray}&&{\tilde{\sigma }}_{\mathrm{int}{ij},j}={[\bar{{\mathbb{C}}}({\bf{E}}-{\bf{I}})]}_{{ijkl}}{\tilde{\varepsilon }}_{{kl},j}={[\bar{{\mathbb{C}}}({\bf{E}}-{\bf{I}})]}_{{ijkl}}{\tilde{u}}_{k,{lj}},\end{eqnarray} \tag{ 4 }$$

where ${\tilde{u}}_{k,{lj}}$ is the displacement to produce ${\tilde{\varepsilon }}_{{kl},j}$. We use the same Fourier transformation method shown in [13] to solve equation (4), and finally to gain the internal stress.

Based on this FFT method and assuming the deformation of the RVE follows a simple constitutive equation as

$$\begin{eqnarray}&&\dot{\varepsilon }=A{\sigma }_{\mathrm{tot}}^{n}=A{({\sigma }_{\mathrm{ext}}+{\sigma }_{\mathrm{int}})}^{n},\end{eqnarray} \tag{ 5 }$$

where $\dot{\varepsilon }$ is the strain rate, A is the reference parameter and n is the stress exponent, the full fields of plastic strain and internal stress in the RVE can be calculated. Before the deformation starts, the misfit strain determines the internal stress. Here, in $\gamma ^{\prime}$ precipitate ${{\boldsymbol{\varepsilon }}}_{\mathrm{mis}}=\delta {\bf{I}}$ and in γ channels ${{\boldsymbol{\varepsilon }}}_{\mathrm{mis}}=0$, where δ is the lattice mismatch ratio. As the plastic strain is generated, we combine the misfit strain and the plastic strain as the general eigenstrain to determine the internal stress in each time step. The new internal stress influences the plastic strain in the next time step. From an example for a RVE with approx. 66% volume fraction of a cube $\gamma ^{\prime}$ precipitate in figure 1, we can see that the plastic strains and the internal stresses in some regions are similar under a constant uniaxial tensile load of 300 MPa along [001] direction. Meanwhile, the internal stress varies with the plastic strain gradually during the deformation. Furthermore, the regions with similar plastic strains are also found in figure 2 for other two loading directions. These results depend on the model parameters in table 1 and the elastic constants at room temperature of a Ni-base single crystal superalloy [28]. Since the $\gamma ^{\prime}$ phase is stronger than γ phase, we use a smaller A for $\gamma ^{\prime}$ phase.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Model parameters.

n	3
A in γ (MPa−n)	$1.0\times {10}^{-12}$
A in $\gamma ^{\prime}$ (MPa−n)	$1.0\times {10}^{-14}$
δ	−0.0015

If this RVE is only considered as an integration point in a FE method, the detailed plastic strains and internal stresses on each FFT node are not required, and it is sufficient to know the solution in the regions of homogeneous deformation behavior. According to the observation and the analysis of plastic strains in the RVE, as well as the approach proposed by Fedelich [16], the plastic deformation can be taken approximately uniform in each channel and in the $\gamma ^{\prime}$ precipitate. Therefore, the RVE generally consists of four main regions which are one $\gamma ^{\prime}$ phase and three γ channels in three given orthogonal directions of cartesian coordinate system, respectively. Based on equation (8) in the literature [16], the uniform plastic strain and the misfit strain in these four regions can be expressed as

$$\begin{eqnarray}\bar{{\boldsymbol{\varepsilon }}}\left\{\begin{array}{ll}{\bar{{\boldsymbol{\varepsilon }}}}^{{x}}\,: & {\rm{plastic}}\,{\rm{strain}}\,{\rm{in}}\,{\rm{channels}}\,{\rm{normal}}\,{\rm{to}}\,{x} \mbox{-} {\rm{d}}{\rm{i}}{\rm{r}}{\rm{e}}{\rm{c}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}};\\ {\bar{{\boldsymbol{\varepsilon }}}}^{{y}}\,: & {\rm{plastic}}\,{\rm{strain}}\,{\rm{in}}\,{\rm{channels}}\,{\rm{normal}}\,{\rm{to}}\,{y} \mbox{-} {\rm{d}}{\rm{i}}{\rm{r}}{\rm{e}}{\rm{c}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}};\\ {\bar{{\boldsymbol{\varepsilon }}}}^{{z}}\,: & {\rm{plastic}}\,{\rm{strain}}\,{\rm{in}}\,{\rm{channels}}\,{\rm{normal}}\,{\rm{to}}\,{z} \mbox{-} {\rm{d}}{\rm{i}}{\rm{r}}{\rm{e}}{\rm{c}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}};\\ {\bar{{\boldsymbol{\varepsilon }}}}^{{p}}\,: & {\rm{plastic}}\,{\rm{strain}}\,{\rm{and}}\,{\rm{misfit}}\,{\rm{strain}}\,{\rm{in}}\,{\rm{the}}\,{\rm{precipitate}},\end{array}\right.\end{eqnarray} \tag{ 6 }$$

where ${\bar{{\boldsymbol{\varepsilon }}}}^{{x}}$, ${\bar{{\boldsymbol{\varepsilon }}}}^{{y}}$, ${\bar{{\boldsymbol{\varepsilon }}}}^{{z}}$, and ${\bar{{\boldsymbol{\varepsilon }}}}^{{p}}$ are deemed as eigenstrains for the calculation of internal stresses.

As shown in figure 3, when ${\bar{{\boldsymbol{\varepsilon }}}}^{{x}}$, ${\bar{{\boldsymbol{\varepsilon }}}}^{{y}}$, ${\bar{{\boldsymbol{\varepsilon }}}}^{{z}}$, and ${\bar{{\boldsymbol{\varepsilon }}}}^{{p}}$ are assigned to the corresponding FFT grids in x-, y-, z-channels and precipitate, respectively, and the FFT grids in the edge and corner areas take the average values of eigenstrains in adjacent channels, the full field of internal stresses in the RVE can be obtained by the introduced FFT method as well. By averaging the values in the respective regions, the ${\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{x}}$, ${\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{y}}$, ${\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{z}}$, and ${\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{p}}$ are adopted to represent the internal stresses there. For a large precipitate and narrow channels, the areas in the edge and corner can be neglected. The general subdivision of the RVE is displayed in figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Even though the FFT method is an efficient way to calculate internal stresses, it is still time-consuming for long-time deformation of a large component with lots of material points. We introduce a mathematical linear relation between internal stresses and eigenstrains with an effective stiffness matrix as

$$\begin{eqnarray}&&{\bar{{\bf{S}}}}^{\mathrm{INT}}={{\mathbb{C}}}^{{\prime\prime} }:\bar{{\boldsymbol{\varepsilon }}}.\end{eqnarray} \tag{ 7 }$$

Extending ${\bar{{\bf{S}}}}^{\mathrm{INT}}$ and $\bar{{\boldsymbol{\varepsilon }}}$ to four sections respectively, equation (7) becomes

$$\begin{eqnarray}&&\left(\begin{array}{c}{\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{x}}\\ {\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{y}}\\ {\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{z}}\\ {\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{p}}\end{array}\right)=\left(\begin{array}{cccc}{{\mathbb{C}}}_{{xx}}^{{\prime\prime} } & {{\mathbb{C}}}_{{xy}}^{{\prime\prime} } & {{\mathbb{C}}}_{{xz}}^{{\prime\prime} } & {{\mathbb{C}}}_{{xp}}^{{\prime\prime} }\\ {{\mathbb{C}}}_{{yx}}^{{\prime\prime} } & {{\mathbb{C}}}_{{yy}}^{{\prime\prime} } & {{\mathbb{C}}}_{{yz}}^{{\prime\prime} } & {{\mathbb{C}}}_{{yp}}^{{\prime\prime} }\\ {{\mathbb{C}}}_{{zx}}^{{\prime\prime} } & {{\mathbb{C}}}_{{zy}}^{{\prime\prime} } & {{\mathbb{C}}}_{{zz}}^{{\prime\prime} } & {{\mathbb{C}}}_{{zp}}^{{\prime\prime} }\\ {{\mathbb{C}}}_{{px}}^{{\prime\prime} } & {{\mathbb{C}}}_{{py}}^{{\prime\prime} } & {{\mathbb{C}}}_{{pz}}^{{\prime\prime} } & {{\mathbb{C}}}_{{pp}}^{{\prime\prime} }\end{array}\right)\cdot \left(\begin{array}{c}{\bar{{\boldsymbol{\varepsilon }}}}^{{x}}\\ {\bar{{\boldsymbol{\varepsilon }}}}^{{y}}\\ {\bar{{\boldsymbol{\varepsilon }}}}^{{z}}\\ {\bar{{\boldsymbol{\varepsilon }}}}^{{p}}\end{array}\right).\end{eqnarray} \tag{ 8 }$$

Taking six components of each stress (strain) tensor into account, there are 24 independent stress (strain) components totally. Hence, the entire effective stiffness ${{\mathbb{C}}}^{{\prime\prime} }$ should be a 24 × 24 matrix. In order to determine the effective stiffness matrix, we assign the value of '1' to each of strain component (e.g., ${\bar{\varepsilon }}_{11}^{{x}}$, ${\bar{\varepsilon }}_{22}^{{x}}$, ${\bar{\varepsilon }}_{33}^{{x}}$, ...), keeping other strain components to '0'. By means of the method mentioned above, the obtained stress components could be one column values in the effective stiffness matrix. For instance, if the first strain component is '1' and the rest of them are '0', the calculated stress components are the first column values of matrix. By repeating this calculation for 24 times, the completed effective stiffness matrix can be obtained. As long as the geometry of RVE and elastic constants are unaltered, this effective stiffness matrix is fixed and can be simply applied to determine the internal stresses.

Our new method is derived from the RVE with a cubic $\gamma ^{\prime}$ precipitate, and it can extend to the case for a spherical $\gamma ^{\prime}$ precipitate. However, the effective stiffness matrix has to be modified, because the different morphology of precipitate leads to a different stress state in the RVE. For the approx. 40% volume fraction of precipitate and the lattice mismatch ratio of −0.0015, under uniaxial tensile stress of 300 MPa along x-direction, the stress is more concentrated in small areas of the γ matrix for the spherical $\gamma ^{\prime}$ precipitate than the cubic one. Furthermore, the values are even higher, as displayed in figure 5. This will cause more inhomogeneous deformation in the microstructure. In addition, the stress field in large corner areas around the spherical $\gamma ^{\prime}$ precipitate is considerable. Therefore, the stress state in the whole RVE with a spherical precipitate can not be represented by four regions as simple as the cubic case.

Zoom In Zoom Out Reset image size

Making use of the same FFT method and equation (5) with the same model parameters, the full fields of plastic strains in the RVE with approx. 40% volume fraction of a spherical $\gamma ^{\prime}$ precipitate for different loading conditions in a certain time step are obtained, which are shown in figure 6. By analyzing the geometrical symmetry and the similarity of deformation behavior, as well as considering the efficiency and accuracy of the calculation of internal stresses, the RVE with a spherical $\gamma ^{\prime}$ precipitate is divided into 10 regions as shown in figure 7. Thus, internal stresses can be calculated by equation (9) as

$$\begin{eqnarray}&&\left(\begin{array}{c}{\bar{{\bf{S}}}}^{\mathrm{INT}1}\\ {\bar{{\bf{S}}}}^{\mathrm{INT}2}\\ \vdots \\ {\bar{{\bf{S}}}}^{\mathrm{INT}10}\end{array}\right)=\left(\begin{array}{cccc}{{\mathbb{C}}}_{11}^{{\prime\prime} } & {{\mathbb{C}}}_{12}^{{\prime\prime} } & \cdots & {{\mathbb{C}}}_{110}^{{\prime\prime} }\\ {{\mathbb{C}}}_{21}^{{\prime\prime} } & {{\mathbb{C}}}_{22}^{{\prime\prime} } & \cdots & {{\mathbb{C}}}_{210}^{{\prime\prime} }\\ \vdots & \vdots & \ddots & \vdots \\ {{\mathbb{C}}}_{101}^{{\prime\prime} } & {{\mathbb{C}}}_{102}^{{\prime\prime} } & \cdots & {{\mathbb{C}}}_{1010}^{{\prime\prime} }\end{array}\right)\cdot \left(\begin{array}{c}{\bar{{\boldsymbol{\varepsilon }}}}^{1}\\ {\bar{{\boldsymbol{\varepsilon }}}}^{2}\\ \vdots \\ {\bar{{\boldsymbol{\varepsilon }}}}^{10}\end{array}\right).\end{eqnarray} \tag{ 9 }$$

This is a simplified equation similar to equation (7). In fact, the whole effective stiffness matrix contains 60 × 60 components.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As the expense of a certain accuracy, the subdivision method overcomes the limit of the computer memory for a long-time full field calculation and increases the efficiency. The subdivisions for two shapes of $\gamma ^{\prime}$ phase proposed in our paper are only two of feasible ways. If the more accurate solution is required, a more complex subdivision can be made by the analysis of specific local stresses and deformations.

3.1. Comparison of internal stresses

We assume that the volume fraction of $\gamma ^{\prime}$ precipitate is approx. 66% and the total number of FFT grids in the RVE with a cubic $\gamma ^{\prime}$ precipitate is 32 × 32 × 32. Thus, the precipitate occupies 28 × 28 × 28 grids and each channel occupies 28 × 28 × 4 grids. Based on this microstructure and elastic constants [28], the effective stiffness matrix with 24 × 24 components can be determined. Since various loading conditions give rise to distinct eigenstrains in the RVE, we assigned many different values to eigenstrains to validate the reliability of equation (7) with the effective stiffness matrix. Here, we purposefully take arbitrary large eigenstrains shown in table 2 as an example to demonstrate the generality of our method. The calculated internal stress fields by FFT method are shown in figure 8. The corresponding effective values in the four sections of the RVE determined by averaging the FFT results within these regions are listed in table 3. The calculated internal stresses by effective stiffness method are listed in table 4. These internal stresses calculated by two methods are almost equivalent. Table 5 shows that the deviation of each stress component in the RVE between two methods, which is the absolute percentage of difference of individual stress component between two methods accounting for the corresponding stress component determined by FFT method. It can be seen that the maximum and average deviations are approx. 5.9% and 0.7%, respectively. Although the maximum difference of individual stress component could be in the order of 10³ MPa for these given large eigenstrains, the deviation is small. If the magnitude of eigenstrains decreases, the value of internal stresses will decrease and the maximum difference of internal stresses would be smaller than 10 MPa which is negligible.

Zoom In Zoom Out Reset image size

Table 2.  Assigned eigenstrains in the four sections of RVE with a cubic $\gamma ^{\prime}$ precipitate.

Strain component	11	22	33	12	13	23
${\bar{{\boldsymbol{\varepsilon }}}}^{{x}}$	0.7	0.8	0.9	0.8	0.7	0.6
${\bar{{\boldsymbol{\varepsilon }}}}^{{y}}$	0.5	0.4	0.3	0.2	0.1	0.0
${\bar{{\boldsymbol{\varepsilon }}}}^{{z}}$	1.0	0.1	0.2	0.3	0.4	0.5
${\bar{{\boldsymbol{\varepsilon }}}}^{{p}}$	0.1	0.2	0.3	0.4	0.5	0.6

Table 3.  Effective internal stresses (MPa) in the four sections of RVE with a cubic $\gamma ^{\prime}$ precipitate obtained by averaging FFT results within these regions.

Stress component	11	22	33	12	13	23
${\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{x}}$	0.351e + 04	−0.108e + 06	−0.111e + 06	−0.657e + 04	−0.114e + 05	−0.508e + 04
${\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{y}}$	−0.372e + 05	0.780e + 04	−0.714e + 04	−0.162e + 04	0.876e + 05	0.827e + 04
${\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{z}}$	−0.103e + 06	−0.172e + 05	0.745e + 04	0.232e + 05	−0.893e + 04	0.415e + 04
${\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{p}}$	0.212e + 05	0.191e + 05	0.178e + 05	−0.318e + 04	−0.124e + 05	0.282e + 03

Table 4.  Calculated internal stresses (MPa) in the four sections of RVE with a cubic $\gamma ^{\prime}$ precipitate by effective stiffness method.

Stress component	11	22	33	12	13	23
${\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{x}}$	0.359e + 04	−0.108e + 06	−0.110e + 06	−0.657e + 04	−0.114e + 05	−0.478e + 04
${\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{y}}$	−0.371e + 05	0.786e + 04	−0.697e + 04	−0.162e + 04	0.878e + 05	0.827e + 04
${\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{z}}$	−0.103e + 06	−0.171e + 05	0.750e + 04	0.234e + 05	−0.893e + 04	0.414e + 04
${\bar{{\bf{S}}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{p}}$	0.213e + 05	0.191e + 05	0.178e + 05	−0.318e + 04	−0.124e + 05	0.279e + 03

Table 5.  The deviation (%) of the internal stresses in the four sections of RVE with a cubic $\gamma ^{\prime}$ precipitate between the results obtained by FFT method and effective stiffness method.

Stress component	11	22	33	12	13	23
${{\bf{D}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{x}}$	2.3	0.0	0.9	0.0	0.0	5.9
${{\bf{D}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{y}}$	0.3	0.8	2.4	0.0	0.2	0.0
${{\bf{D}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{z}}$	0.0	0.6	0.7	0.9	0.0	0.2
${{\bf{D}}}^{{\rm{I}}{\rm{N}}{\rm{T}}{p}}$	0.5	0.0	0.0	0.0	0.0	1.0

When the volume fraction of $\gamma ^{\prime}$ precipitate changes, the average deviation varies a little, but the fluctuation of maximum deviation is larger, as shown in figure 9. There is no monotonic regulation between deviations and volume fractions. Nevertheless, these deviations are generally acceptable. Furthermore, if the precipitate is a cuboid, which occupies 28 × 30 × 26 grids, the maximum and average deviations are approx. 6.8% and 0.6%, respectively. Therefore, the simple effective stiffness method is valid for calculation of internal stresses in both cubic and cuboid cases with different volume fractions of $\gamma ^{\prime}$ precipitates, and the accuracy increases with decreasing magnitude of eigenstrains.

Zoom In Zoom Out Reset image size

When we use different elastic constants of γ and $\gamma ^{\prime}$ phases [13] to redo the calculation for cubic precipitate with 66% volume fraction, the maximum and average deviations are approx. 2.9% and 0.4%, respectively. This demonstrates that the unequal elastic constants for two phases do not worsen the accuracy of this new method, and the same indication is also found for other volume fractions.

For given arbitrary values of eigenstrains in the 10 sections of the RVE with a spherical $\gamma ^{\prime}$ precipitate as shown in table 6, the maximum and average deviation between FFT method and effective stiffness method are approx. 4.7% and 0.4%, respectively, as displayed in table 7. The influence of volume fractions is shown in figure 10. It is found that for volume fraction of 20%, one stress component has a large deviation, however, its absolute difference is only 25 MPa. Furthermore, the assigned eigenstrains in this paper are complex, which may increase the deviation. The calculated internal stresses would be more accurate in the simple deformation case. Hence, the effective stiffness method is also available to obtain the reasonable results for spherical $\gamma ^{\prime}$ precipitates.

Zoom In Zoom Out Reset image size

Table 6.  Assigned eigenstrains in the 10 sections of RVE with a spherical $\gamma ^{\prime}$ precipitate.

Strain component	11	22	33	12	13	23
${\bar{{\boldsymbol{\varepsilon }}}}^{1}$	0.5	0.1	0.5	0.7	0.04	0.05
${\bar{{\boldsymbol{\varepsilon }}}}^{2}$	−0.1	0.0	0.3	0.8	0.8	0.0
${\bar{{\boldsymbol{\varepsilon }}}}^{3}$	0.3	0.4	−0.4	0.1	0.2	0.07
${\bar{{\boldsymbol{\varepsilon }}}}^{4}$	0.2	−0.2	0.0	0.0	0.1	0.01
${\bar{{\boldsymbol{\varepsilon }}}}^{5}$	0.5	−0.2	0.0	0.1	0.0	0.08
${\bar{{\boldsymbol{\varepsilon }}}}^{6}$	0.2	0.0	0.0	−0.2	0.7	0.0
${\bar{{\boldsymbol{\varepsilon }}}}^{7}$	0.3	0.0	−0.2	0.4	0.0	0.0
${\bar{{\boldsymbol{\varepsilon }}}}^{8}$	0.0	0.2	0.5	0.05	−0.6	0.0
${\bar{{\boldsymbol{\varepsilon }}}}^{9}$	0.4	0.3	0.4	0.7	0.0	0.0
${\bar{{\boldsymbol{\varepsilon }}}}^{10}$	0.4	0.4	0.0	0.007	0.3	1.0

Table 7.  The deviation (%) of the internal stresses in the 10 sections of RVE with a spherical $\gamma ^{\prime}$ precipitate between the results obtained by FFT method and effective stiffness method.

Stress component	11	22	33	12	13	23
${{\bf{D}}}^{\mathrm{INT}1}$	0.0	0.4	0.0	0.5	0.2	0.2
${{\bf{D}}}^{\mathrm{INT}2}$	0.6	0.2	0.8	0.2	0.2	0.3
${{\bf{D}}}^{\mathrm{INT}3}$	4.7	0.5	0.8	0.0	0.0	0.0
${{\bf{D}}}^{\mathrm{INT}4}$	0.2	0.0	0.0	0.5	0.3	1.8
${{\bf{D}}}^{\mathrm{INT}5}$	0.0	0.0	0.0	0.5	0.0	0.0
${{\bf{D}}}^{\mathrm{INT}6}$	0.0	0.0	0.0	0.4	0.1	0.0
${{\bf{D}}}^{\mathrm{INT}7}$	0.4	0.7	0.2	0.0	0.0	0.3
${{\bf{D}}}^{\mathrm{INT}8}$	0.4	0.0	1.5	0.9	0.0	0.2
${{\bf{D}}}^{\mathrm{INT}9}$	0.0	0.2	0.0	0.8	0.0	0.2
${{\bf{D}}}^{\mathrm{INT}10}$	2.3	0.0	0.7	0.2	0.0	0.1

3.2. Comparison of plastic strains

Following the present subdivision and assigning the internal stresses in different regions obtained by the FFT method or the effective stiffness method to equation (5), the plastic strain in each region in every time step can be calculated. Thus, the general equivalent plastic strain of RVE is given by

$$\begin{eqnarray}&&{\varepsilon }_{{\rm{g}}}=\displaystyle \sum _{i=1}^{m}{f}_{i}{\varepsilon }_{i},\end{eqnarray} \tag{ 10 }$$

where f_i and ${\varepsilon }_{i}$ are the volume fraction and the equivalent plastic strain of each region i, respectively, and m is the total number of regions. The used model parameters are the same as that in table 1. Figure 11 shows that the comparisons of the equivalent plastic strains for different methods in some typical regions of the RVEs with a cubic $\gamma ^{\prime}$ precipitate and a spherical $\gamma ^{\prime}$ precipitate, respectively, which demonstrates that the FFT method and the effective stiffness method almost lead to the same results, where small deviations occur at rather large strains. For [110] and [111] loading conditions, the general equivalent plastic strains obtained by two methods are nearly identical as shown in figure 12. These results denote that, in general, the deviation increases with increasing strain. This is consistent with the previous indication that the accuracy of internal stresses calculated by the effective stiffness method decreases with increasing magnitude of eigenstrains. Moreover, the deviation in the case of a cubic precipitate is larger than that in the case of a spherical precipitate for the same strain. Due to the more subsections in the RVE with a spherical precipitate, the description of the stress states is improved, so that the precision is enhanced.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1. Crystal plasticity constitutive equations

Based on the dislocation slip controlled phenomenological crystal plasticity model [29], a new flow rule for $\langle 110\rangle \{111\}$ slip in each γ matrix channel is

$$\begin{eqnarray}&&{\dot{\gamma }}_{\alpha }={\dot{\gamma }}_{0}{\left|\displaystyle \frac{{\tau }_{\alpha }+{\tau }_{\alpha }^{\mathrm{INT}}}{{\hat{\tau }}_{\alpha }^{\mathrm{slip}}+{\tau }^{\mathrm{oro}}}\right|}^{{p}_{1}}\mathrm{sign}({\tau }_{\alpha }+{\tau }_{\alpha }^{\mathrm{INT}}),\end{eqnarray} \tag{ 11 }$$

where α represents the slip system, ${\dot{\gamma }}_{\alpha }$ is the shear rate in a certain slip system, ${\dot{\gamma }}_{0}$ is the reference shear rate, ${\tau }^{\mathrm{oro}}$ is the Orowan stress, and p₁ is the inverse value of the strain rate sensitivity. The resolved shear stress ${\tau }_{\alpha }$ is given by

$$\begin{eqnarray}&&{\tau }_{\alpha }=\displaystyle \frac{{\mathbb{C}}}{2}({{\bf{F}}}_{{\rm{e}}}^{{\rm{T}}}{{\bf{F}}}_{{\rm{e}}}-{\bf{I}})\cdot {{\bf{M}}}_{\alpha },\end{eqnarray} \tag{ 12 }$$

where ${{\bf{F}}}_{{\rm{e}}}$ is the elastic deformation gradient, and ${{\bf{M}}}_{\alpha }$ is the Schmid matrix. The internal stress gives rise to an additional resolved shear stress as

$$\begin{eqnarray}&&{\tau }_{\alpha }^{\mathrm{INT}}={\bar{{\bf{S}}}}^{\mathrm{INT}}\cdot {{\bf{M}}}_{\alpha }.\end{eqnarray} \tag{ 13 }$$

${\hat{\tau }}_{\alpha }^{\mathrm{slip}}$ is the slip resistance determined by

$$\begin{eqnarray}&&{\dot{\hat{\tau }}}_{\alpha }^{\mathrm{slip}}=\displaystyle \sum _{\beta =1}^{12}{h}_{0}{\chi }_{\alpha \beta }{\left(1-\displaystyle \frac{{\hat{\tau }}_{\beta }^{\mathrm{slip}}}{{\hat{\tau }}^{\mathrm{sat}}}\right)}^{{p}_{2}}| {\dot{\gamma }}_{\beta }| ,\end{eqnarray} \tag{ 14 }$$

where β is the slip system, h₀ is the reference hardening parameter, ${\chi }_{\alpha \beta }$ is the cross hardening matrix, ${\hat{\tau }}^{\mathrm{sat}}$ is the saturated slip resistance due to dislocation density accumulation, and p₂ is a fitting parameter. Due to the shearing of $\gamma ^{\prime}$ precipitate by stacking faults [7, 16], the flow rule for $\langle 112\rangle \{111\}$ slip in $\gamma ^{\prime}$ precipitate is

$$\begin{eqnarray}&&{\dot{\gamma }}_{\alpha }={\dot{\gamma }}_{0}{\left|\displaystyle \frac{{\tau }_{\alpha }+{\tau }_{\alpha }^{\mathrm{INT}}}{{\hat{\tau }}_{\alpha }^{\mathrm{slip}}}\right|}^{{p}_{1}}\mathrm{sign}({\tau }_{\alpha }+{\tau }_{\alpha }^{\mathrm{INT}}).\end{eqnarray} \tag{ 15 }$$

Taking the RVE with a cubic $\gamma ^{\prime}$ precipitate as an example, the plastic velocity gradients in the four subsections as discussed above are given by

$$\begin{eqnarray}&&{{\bf{L}}}_{{\rm{p}}}^{{x}}=\displaystyle \sum _{\alpha =1}^{12}{\dot{\gamma }}_{\alpha }^{{x}}{{\bf{M}}}_{\alpha },\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{{\bf{L}}}_{{\rm{p}}}^{{y}}=\displaystyle \sum _{\alpha =1}^{12}{\dot{\gamma }}_{\alpha }^{{y}}{{\bf{M}}}_{\alpha },\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{{\bf{L}}}_{{\rm{p}}}^{{z}}=\displaystyle \sum _{\alpha =1}^{12}{\dot{\gamma }}_{\alpha }^{{z}}{{\bf{M}}}_{\alpha },\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{{\bf{L}}}_{{\rm{p}}}^{{p}}=\displaystyle \sum _{\alpha =1}^{12}{\dot{\gamma }}_{\alpha }^{{p}}{{\bf{M}}}_{\alpha }.\end{eqnarray} \tag{ 19 }$$

With the known plastic deformation gradient in the previous time step, the current plastic deformation gradients for different regions are

$$\begin{eqnarray}&&{{\bf{F}}}_{{\rm{p}}}^{{x}(t)}=({\bf{I}}+{{\bf{L}}}_{{\rm{p}}}^{{x}}\cdot {\rm{\Delta }}t){{\bf{F}}}_{{\rm{p}}}^{{x}(t-{\rm{\Delta }}t)},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{{\bf{F}}}_{{\rm{p}}}^{{y}(t)}=({\bf{I}}+{{\bf{L}}}_{{\rm{p}}}^{{y}}\cdot {\rm{\Delta }}t){{\bf{F}}}_{{\rm{p}}}^{{y}(t-{\rm{\Delta }}t)},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{{\bf{F}}}_{{\rm{p}}}^{{z}(t)}=({\bf{I}}+{{\bf{L}}}_{{\rm{p}}}^{{z}}\cdot {\rm{\Delta }}t){{\bf{F}}}_{{\rm{p}}}^{{z}(t-{\rm{\Delta }}t)},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{{\bf{F}}}_{{\rm{p}}}^{{p}(t)}=({\bf{I}}+{{\bf{L}}}_{{\rm{p}}}^{{p}}\cdot {\rm{\Delta }}t){{\bf{F}}}_{{\rm{p}}}^{{p}(t-{\rm{\Delta }}t)},\end{eqnarray} \tag{ 23 }$$

and the total plastic deformation gradient ${{\bf{F}}}_{{\rm{p}}}$ is calculated by

$$\begin{eqnarray}&&{{\bf{F}}}_{{\rm{p}}}^{(t)}=({\bf{I}}+({f}^{{x}}\cdot {{\bf{L}}}_{{\rm{p}}}^{{x}}+{f}^{{y}}\cdot {{\bf{L}}}_{{\rm{p}}}^{{y}}+{f}^{{z}}\cdot {{\bf{L}}}_{{\rm{p}}}^{{z}}+{f}^{{p}}\cdot {{\bf{L}}}_{{\rm{p}}}^{{p}})\cdot {\rm{\Delta }}t){{\bf{F}}}_{{\rm{p}}}^{(t-{\rm{\Delta }}t)},\end{eqnarray} \tag{ 24 }$$

where ${f}^{{x}}$, ${f}^{{y}}$, ${f}^{{z}}$, and ${f}^{{p}}$ are volume fractions of each γ channel and precipitate, respectively. ${\rm{\Delta }}t$ is the time step. The eigenstrains in each subsection originated from the plastic deformation gradients are

$$\begin{eqnarray}&&{\bar{{\boldsymbol{\varepsilon }}}}^{{x}}=\displaystyle \frac{1}{2}({{\bf{F}}}_{{\rm{p}}}^{{x}{\rm{T}}}{{\bf{F}}}_{{\rm{p}}}^{{x}}-{\bf{I}}),\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{\bar{{\boldsymbol{\varepsilon }}}}^{{y}}=\displaystyle \frac{1}{2}({{\bf{F}}}_{{\rm{p}}}^{{y}{\rm{T}}}{{\bf{F}}}_{{\rm{p}}}^{{y}}-{\bf{I}}),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{\bar{{\boldsymbol{\varepsilon }}}}^{{z}}=\displaystyle \frac{1}{2}({{\bf{F}}}_{{\rm{p}}}^{{z}{\rm{T}}}{{\bf{F}}}_{{\rm{p}}}^{{z}}-{\bf{I}}),\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{\bar{{\boldsymbol{\varepsilon }}}}^{{p}}=\displaystyle \frac{1}{2}({{\bf{F}}}_{{\rm{p}}}^{{p}{\rm{T}}}{{\bf{F}}}_{{\rm{p}}}^{{p}}-{\bf{I}})+\delta {\bf{I}}.\end{eqnarray} \tag{ 28 }$$

The internal stresses can be determined from these eigenstrains at each time step by equation (7) efficiently. The same method is able to be used for the RVE with a spherical $\gamma ^{\prime}$ precipitate.

4.2. FE model

The crystal plasticity constitutive equations are able to be implemented in the ABAQUS [30] platform by using a user-defined subroutine UMAT. A simple 3D FE model with dimensions of 80 mm × 40 mm × 40 mm is set up to represent a sample of single crystal superalloys. As shown in figure 13, the 1/8 complete sample is discretized with 2000 regular C3D8 elements, in which each integration point is represented by the previous RVE. The symmetric boundary conditions are applied in x-, y-, z-direction, respectively. The common model parameters are shown in table 8. According to different shapes and volume fractions of $\gamma ^{\prime}$ precipitate, as well as different γ channel widths, ${\tau }^{\mathrm{oro}}$ can be determined by the Orowan threshold equations proposed in [13, 31]. Assuming the volume fraction of $\gamma ^{\prime}$ precipitate is 40% and the length of the RVE is 500 nm, the Orowan stresses are 71 MPa and 22 MPa for cubic $\gamma ^{\prime}$ and spherical $\gamma ^{\prime}$ precipitates, respectively.

Zoom In Zoom Out Reset image size

Table 8.  Crystal plasticity model parameters. ${\hat{\tau }}_{0}^{\mathrm{slip}}$ is the initial slip resistance.

${\dot{\gamma }}_{0}$ (s−1)	$4.0\times {10}^{-6}$
${\hat{\tau }}_{0}^{\mathrm{slip}}\langle 110\rangle \{111\}$ (MPa)	60
${\hat{\tau }}_{0}^{\mathrm{slip}}\langle 112\rangle \{111\}$ (MPa)	500
p1	10
p2	0.05
h0 (MPa)	60
${\hat{\tau }}^{\mathrm{sat}}$ (MPa)	800
${\chi }_{\alpha \beta }$ coplanar{111}	1.0
${\chi }_{\alpha \beta }$ non-coplanar{111}	1.4
δ	−0.0015

Loading in z-direction (i.e., [001] crystallographic orientation) with a constant stress of 300 MPa for $1.0\times {10}^{5}$ s, the sample deforms almost homogeneously and the plastic strain evolving with time is displayed in figure 14. Since the internal stresses in the whole sample are nearly same, we select one arbitrary element to show the evolution of internal stresses during the deformation. For the cubic $\gamma ^{\prime}$ precipitate, from figure 15(a), it can be seen that the internal stresses in γ channels decrease in the beginning and increase gradually afterwards, but the internal stress in the precipitate increases continuously. Initially, the internal stresses are determined by the lattice mismatch. The following deformation accommodates the misfit strain first, so that the high internal stresses in channels decrease. This is also observed in the full field calculation shown in figure 1(b). However, when they reach certain values, they will increase again, because the further accumulation of deformation in the microstructure will cause the increase of internal stresses. In addition, due to the high strength of $\gamma ^{\prime}$ precipitate, the inhomogeneous deformation in the channels and the precipitate leads to increasing internal stresses. The previous deformation state determines the current internal stresses, the current internal stresses influence the further deformation as well. The similar phenomenon of internal stresses during the creep in the superalloy has been discussed in [32]. For the spherical $\gamma ^{\prime}$ precipitate, the evolution of internal stresses is shown in figure 15(b). It is found that the different shapes of $\gamma ^{\prime}$ precipitates result in different internal stresses which give rise to different plastic strains. Our method can discriminate the influence of the various shapes of $\gamma ^{\prime}$ precipitate on deformation behavior for the same volume fraction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The present FE model is one simple application of our effective stiffness method. The detailed analysis of deformation in this model is not the focus of this paper. In the future, the effective stiffness method can be used in more complex and more realistic models, in which the simulation results are able to rationalize the experimental observations.

On the basis of typical γ/γ' microstructures of superalloys, the internal stresses due to eigenstrains resulting from lattice mismatch and deformation heterogeneity can be calculated by an efficient linear relation with an effective stiffness matrix. The effective stiffness method can be considered as a generalization of the Eshelby inclusion problem for more complex microstructures. It relies on a subdivision of the microstructure into regions with uniform deformation behavior. The feasibility of such a subdivision has been demonstrated for γ/γ' microstructures with cubic and spherical precipitates. Compared with direct simulations with the FFT method, which determines the internal stress fields with the full spatial resolution, this simplified method is numerically much more efficient and thus well suitable for applications in macroscopic simulations of deformation and creep with FE method or other mechanics solvers. The effective stiffness matrix is calculated by the FFT method and its dimension varies with different morphologies of $\gamma ^{\prime}$ precipitates. The accuracy of this method increases with decreasing magnitude of eigenstrains, but it is not significantly influenced by the volume fraction of $\gamma ^{\prime}$ precipitate and elastic constants. Furthermore, if the volume fraction and shape of $\gamma ^{\prime}$ phase in superalloys remain unchanged during the deformation, the effective stiffness matrix will be constant throughout the simulation.

The authors are grateful for funding by the Deutsche Forschungsgemeinschaft (DFG) through Project C4 of the collaborative research center SFB/Transregio 103 superalloy single crystals under grant number INST 213/747-2.