Temperature-Dependent Misfit Stress in Gamma Double Prime Strengthened Ni-Base Superalloys

Abstract

Misfit stresses in the γ″-strengthened Ni-base superalloy Inconel 718 were calculated from the measured constrained misfit strain using Eshelby’s inclusion method. The constrained misfit strains of the γ″ precipitates were measured using neutron diffraction at various temperatures with the aid of the stress-induced variant selection method. Eshelby tensor was calculated using the expressions for the case of anisotropic matrix given by Mura. Results show the presence of significant compressive misfit stresses in the γ″ precipitates with an anisotropic distribution, namely 3.0 GPa along the habit plane and 1.7 GPa along the plane normal direction at room temperature, and 2.0 and 1.2 GPa at 664 °C. The decrease in misfit stresses was due to the decrease in stiffness and the different coefficients of thermal expansion of the γ and γ″ phases. The average internal stresses in the γ matrix due to lattice misfit were found to be ~ 329 MPa at room temperature and ~ 186 MPa at 664 °C in tension. The possibility of relieving such high levels of misfit stresses in precipitates by loss of coherency during continued growth of precipitates is also discussed.

1 Introduction

Ni-base superalloys derive strengthening from a coherent precipitate phase embedded in a matrix phase.[1] The precipitate phase would normally possess lattice parameter(s) different from that of the matrix phase which leads to a lattice misfit between the two phases.[2] The lattice misfit generates misfit stresses in both the precipitates and the matrix. Misfit stresses will influence microstructure evolution including precipitate shape required to accommodate the misfit stresses, precipitate coarsening kinetic and rafting phenomenon at elevated temperatures during processing or servicing of the alloys.[2] Such misfit stresses will also interact with dislocations during plastic deformation and influence the strength of the alloys.[3] Misfit stresses are determined by the misfit strains and elastic moduli which are both temperature dependent. Therefore, misfit stresses are temperature dependent and their influence on the alloys will vary with temperature. Quantitative analysis of misfit stresses in the system would thus help deepen the understanding of the materials in different aspects.

Inconel 718 (IN718) is a widely used Ni-base superalloys in turbine engines. IN718 derives its strengthening mainly from the γ″ precipitates which is body-centered tetragonal.[3] The γ″ precipitates are coherently embedded in the γ matrix with different misfit strains along a and c-axes of the γ″ precipitates. Due to the anisotropic misfit strains that lead to anisotropic misfit stresses, the γ″ precipitates grow into oblate ellipsoids.[4] Lattice misfit between the two phases can be up to 3 pct along the c-axis,[5] which is much larger than the misfits in γ′-strengthened Ni-base superalloys where the misfits between γ and γ′ phases are normally less than 1 pct.[6] The large lattice misfits lead to large misfit stresses. Phase field simulations have shown the misfit stresses in the γ″ precipitates at the level of 1 GPa in compression at 790 °C.[7] Misfit stresses can vary significantly at different temperatures due to the different thermal expansion of the strengthening and matrix phases and the temperature-dependent elastic moduli.

The misfit stresses arising from the lattice misfit between the strengthening phase and matrix phase can be calculated using the Eshelby inclusion method.[8] The dilute dispersion of ellipsoidal γ″ precipitates in this study finely meet the case of Eshelby’s thought experiment described in his paper, which provides analytical solutions to the elastic strain and stress field in an ellipsoidal inclusion.[8] Eshelby’s method has substantial applications in research of solid mechanics and details of the method can be found in literature.[8,9,10,11,12,13]

2 The Eshelby Method

2.1 The Transformed Inclusion

Starting from a situation where inclusion and matrix are elastically homogeneous, in Eshelby’s method, the problem is solved by a series of cutting and welding exercises as shown in Figure 1[8]: (i) a region (the inclusion) is cut from a material (the matrix), (ii) the region undergoes a shape transformation strain (ε^T) to reach its unconstrained state, (iii) apply an appropriate surface traction on the inclusion to restore it to the original shape, (iv) then put it back to the hole where it was cut off from the matrix and weld the interface between the inclusion and the matrix, (v) remove the surface traction and the equilibrium between the inclusion and matrix is reached at a constrained strain (ε^C) in the inclusion relative to its original shape in stage (i).

Fig. 1

Schematic illustration of Eshelby’s cutting and welding processes

Now the elastic strain within the inclusion at stage (v) is $\varepsilon^{\text{e}} = \varepsilon^{\text{C}} - \varepsilon^{\text{T}}$, and the stress in the inclusion σ_I can be expressed by Hooke’s law:

$$ \sigma_{\text{I}} = C_{\text{I}} \left( {\varepsilon^{\text{C}} - \varepsilon^{\text{T}} } \right) $$

(1)

where C_I is the elastic stiffness of the inclusion.

In precipitate-strengthening Ni-base superalloys, the precipitate phase (γ′ and γ″) and the γ matrix phase are elastically inhomogeneous. The inhomogeneity has a stress-free transformation strain (SFTS) $\varepsilon^{{{\text{T}}*}}$ which differs from $\varepsilon^{\text{T}}$ of a homogeneous inclusion. Such inhomogeneous problems can be solved using the equivalent inclusion method[8, 9] that replaces the inhomogeneity with elastic stiffness C_I by a homogeneous inclusion with elastic stiffness C_M, maintaining identical stresses in the inclusion:

$$ \sigma_{\text{I}} = C_{\text{M}} \left( {\varepsilon^{\text{C}} - \varepsilon^{\text{T}} } \right) = C_{\text{I}} \left( {\varepsilon^{\text{C}} - \varepsilon^{{{\text{T}}*}} } \right) $$

(2)

Eshelby found that the stress in an ellipsoidal inclusion is uniform and the constrained strain ε^C can be related to the stress-free transformation strain (SFTS) ε^T by the so-called Eshelby tensor S[8]:

$$ \varepsilon^{\text{C}} = S\varepsilon^{\text{T}} $$

(3)

The Eshelby tensor depends solely upon the inclusion shape and the elastic properties of the matrix, such as the Poisson’s ratio and Young’s modulus for an isotropic matrix or the elastic constants for an anisotropic matrix.[8, 10] The calculation of the Eshelby tensor is lengthy and here we adopt the expressions given by Mura for the case of anisotropic matrix[10] which are shown in Appendix.

2.2 Determination of constrained misfit strain

The γ phase has a face-centered cubic crystallographic structure, and the γ″ phase possesses an ordered D0₂₂ body-centered tetragonal crystallographic structure. The γ″ phase maintains a γ″//{100}γ and γ″//<100>γ crystallographic relationship with the γ matrix, and three possible variants of γ″ would exist in approximate equal quantities.[7] The normal of γ″ precipitate habit plane is parallel to the c-axis and therefore, the constrained strain $\varepsilon_{{ {\text{a}}}}^{\text{C}}$ along a-axis and $\varepsilon_{ c}^{\text{C}}$ along c-axis for γ″ can be determined by

$$ \varepsilon_{\gamma^{\prime\prime}a}^{\text{C}} = \frac{{d_{\gamma^{\prime\prime}200}^{\text{C}} - d_{\gamma 200}^{\text{C}} }}{{d_{\gamma 200}^{\text{C}} }} $$

(4a)

$$ \varepsilon_{\gamma^{\prime\prime}c}^{\text{C}} = \frac{{d_{\gamma^{\prime\prime}004}^{\text{C}} - d_{\gamma 200}^{\text{C}} }}{{d_{\gamma 200}^{\text{C}} }} $$

(4b)

where $d_{200}^{\text{C}}$, $d_{004}^{\text{C}}$ and $d_{200}^{\text{C}}$ are spacings of hkl lattice planes of the γ″ and γ phases measured in bulk IN718, respectively. The superscript capital C denotes the constrained state. Neutron or X-ray diffraction techniques are desirable in measuring lattice parameters of constituent phases in multi-phase alloys. Lattice parameters of γ′ and γ phases in γ′-strengthened Ni-base superalloys such as Li720, CMSX-4 are very similar (within 1 pct) and lead to overlapping diffraction peaks. Such overlapping problems are commonly solved using the deconvolution method developed by Stone,[14] and thereby lattice spacings of both the γ′ and γ phases are obtained. In the case of IN718 where γ″ and γ phases exist, the deconvolution method to separate all three peaks becomes unsuitable. Mukherji et al. showed that a direct separation of the overlapping peak was feasible if the volume fraction of each phase was known.[15] In our previous study, we have developed a method using stress-induced variant selection (SIVS) effect to facilitate the deconvolution of overlapping peaks from IN718.[16] Using such a method, lattice spacings of the γ and γ″ phases can be determined from diffraction patterns.

2.3 Elastic Constants

In the calculation of the Eshelby tensor, a knowledge of elastic constants of the γ matrix phase is required. In order to calculate the SFTS $\varepsilon^{{{\text{T}}*}}$ of the inhomogeneity, the elastic constants of the γ″ phase also need to be known. In the current study, tensile tests on IN718 were performed at room temperature (RT) and at elevated temperature (664 °C). Young’s modulus E and DECs of the γ phase at these two temperatures were obtained and used to derive elastic constants at temperatures in between. Here we assume the Young’s modulus E follows a linear relationship with temperature in the current temperature range according to the measurements in literature,[17] and the elastic constants of the γ phase at various temperatures follow the same linear relationship with temperature. For the γ″ phase, elastic constants are adopted from first principle calculations by Lin et al.[18] for various temperatures.

3 Experimental Procedures

3.1 Sample Preparation

Polycrystalline IN718 with nominal chemical composition listed in Table I was used in the current study. Specimens were extracted from forged IN718, which was solution heat treated at 1000 °C for 1 hour before machining into tensile bar shape with gauge diameter of 8 mm and length of 42 mm. In order to make the deconvolution feasible for the subsequent neutron diffraction experiment, the SIVS method was employed, which applied 300 MPa tensile stress to the specimens during aging heat treatment at 790 °C for 5 hours. Depending on the grain orientation, the SIVS led to a reduction in the number of γ″ variant in grains. For example, in {100} oriented grains, only one γ″ variant with its habit plane normal parallel to the loading axis exist. A detail description of the SIVS method can be found in our previous study.[16] After aging heat treatment, the sample were air cooled to RT to minimize the generation of residual stresses.

Table I Chemical Composition (Wt Pct)

3.2 Microstructural Characterization

As aforementioned, the misfit strains were measured along the principal axes of the crystals. The microstructural observation by scanning electron microscopy (SEM) was performed on a [100] oriented grain with respect to the axial direction of the specimen. In order to measure the aspect ratio of the γ″ precipitates, thin foils were prepared by focus-ion beam (FIB) and subjected to transmission electron microscopy (TEM) characterization.

3.3 Misfit Strain Measurement Using Neutron Diffraction

Diffraction measurements were performed in the ENGIN-X neutron diffractometer at ISIS neutron source in the UK. The specimens were mounted horizontally on a stress rig (INSTRON, 100 KN) with the axial direction 45 deg to the incident beam. Bank 1 and bank 2 detectors were fixed at 2θ = ±90 deg to the incident beam, so that the d-spacing parallel and perpendicular to the axial direction were measured by the bank 1 and bank 2 detectors, respectively (Figure 2). One specimen was subjected to tensile loading test at RT to obtain the Young’s modulus E and DECs. The other specimen was subjected to heating from RT to 600 °C in a stepwise manner under a negligible load (5 MPa). During heating, the step interval was set to 100 °C with a heating time of 5 minutes, and each temperature was held for 20 minutes for neutron scanning. After this, a tensile loading test was performed on the sample at the 600 °C. Three K-type thermocouples were connected to the surface of the rod specimens by twisting the thermocouple wires to control and monitor the temperature. The heating was carried out by a radiative furnace open to air. Neutron diffraction acquisition time was 20 minutes.

Fig. 2

Schematic illustration of neutron experiment setup at ENGIN-X, ISIS, UK

3.4 Temperature Calibration

After the experiments, a tiny gap (~ < 0.5 mm) between the thermocouple contact pin and the sample surface was noticed. In order to rectify the measured temperatures, a post-test temperature calibration was performed. The calibration was performed using two thermocouples, with one directly in contact with the sample surface while the other was twisted around the sample. The sample was heated by the radiative furnace with temperature increased stepwise. At each step, the temperature was held for about 5 minutes for temperature homogenization. The plots of measured temperatures by the two thermocouples are plotted in Figure 3. A polynomial relationship (Eq. [5]) was obtained between the two sets of temperature. The derived relationship was then used to calibrate the temperature data measured during the previous heating experiment. The set temperature of 600 °C was calibrated to be 664 °C.

Fig. 3

Temperature calibration plot

$$ T_{\text{Cali}} \; = \; - 0.229\; \times \;10^{ - 3} \; \times \;T_{\text{twist}}^{2} \; + \;1.25\; \times \;T_{\text{twist}} \; - \;5.25 $$

(5)

4 Experimental Results

4.1 Microstructural Characterization

Figure 4(b) shows the micrograph of the cross-section of a grain with [100] orientation with respect to the sample axial direction, selected from an electron backscattering diffraction (EBSD) map as shown in Figure 4(a). The γ matrix was eliminated by chemical etching with a solution (5 g CuCl2 + 100 ml hydrochloric acid + 100 ml ethanol) in order to reveal the γ′ and γ″ precipitates. As shown in Figure 4(b), the smaller round particles are the γ′ precipitates, whilst the larger particles are the γ″ precipitates. As mentioned above, SIVS allowed only the γ″ variant with its habit plane normal to the stress axis (along the viewing direction in this micrograph) to exist (Figure 4(b)). Although the volume fraction of both precipitates in the bulk sample could not be estimated from the micrograph, the γ′-to-γ″ volume fraction ratio was estimated by counting the number of particles and their average particle size. In the current study, the volume ratio was estimated to be 0.36.

Fig. 4

(a) EBSD map on cross-section normal to the axial axis of the specimen, (b) SEM micrograph on a (100) oriented grain selected from the EBSD map as shown

4.2 Aspect Ratio

Figure 5(a) shows the [100]-oriented bright field TEM image of a FIB foil from the aged IN718. Only one variant of the γ″ precipitates displaying in ellipsoidal shape was observed in the TEM image. The inset selected area electron diffraction (SAED) shows only one set of the superlattice diffraction pattern confirming the existence of only one γ″ variant in the micrograph. The major and minor axes of the ellipse were counted and the aspect ratio (minor-to-major axis) was obtained for each particle in the TEM image and plotted against the major axis in Figure 5(b). The aspect ratio, α, decreases with increasing precipitate size. This relationship has been reported previously in literature.[19,20,21] The volume-averaged major and minor axes were determined to be 55.8 nm and 12.1 nm, respectively, and hence the overall volume-averaged aspect ratio was determined to be 0.22, which is consistent with the reported values.[20, 21]

Fig. 5

(a) [100]-oriented bright field TEM image of γ″ precipitates, inset is the SAED pattern, (b) aspect ratio vs major axis

4.3 Peak Deconvolution

Since the three major phases (γ′, γ″ and γ) have similar lattice parameters, the diffraction peaks from the three phases overlap. In order to obtain isolated d-spacing of the γ and γ″ phases, deconvolution of the overlapping peaks is required. Mukherji et al. showed the feasibility to deconvolute the overlapping peaks from three constituent phases if the volume fraction of each phase is known a priori.[22] However, the volume fraction of each phase in IN718 is difficult to determine experimentally and the reported volume fractions vary in literature. Thermodynamic simulations provide the equilibrium volume fractions based on a pseudo-ternary database, but the actual sample would have not reached the equilibrium state in the absence of sufficient aging time.

In the current study, the samples were tensile-stress aging heat treated before the neutron diffraction experiment. The γ″ variants were selected by the applied stress[16] and for the {100} oriented grains, only the γ″ variant with disc-plane normal to the axial direction is present. Such a selected variant distribution will facilitate the peak deconvolution which is described below. Figure 6(a) shows diffraction peaks for {200} grains obtained from bank 1 detector at RT. One major peak (left) and one minor peak (right) were identified. The minor peak is the {004} γ″ peak, since only one γ″ variant exist, all the volume of γ″ phase in these oriented grains contributes to the {004} γ″ peak, and thereby, the major peak consists of only the γ′ and γ phase. Since the volume fraction of the γ′ phase is very low and has a lattice parameter very close to that of the γ phase,[5] the major peak can be treated as a single γ peak. As shown in Figure 6(a), the major and minor peaks can be fitted by two Pseudo-Voigt functions. Using this assumption, the γ-to-γ″ peak intensity ratio was obtained to be around 8. This intensity ratio was used in the subsequent deconvolution of the diffraction peak obtained from bank 2 detector.

Fig. 6

(a) Fitting of peaks obtained from bank 1 detector and (b) deconvolution of overlapping peaks obtained from bank 2 detector

Figure 6(b) shows the deconvolution of the overlapping peak obtained from bank 2 detector which consists of {200} peaks from the γ″ and γ phases. The deconvolution was performed by fitting the peak with two Pseudo-Voigt functions representing the two phases. The intensity of the γ peak was set to be 8 times of that of the γ″ peak. Such a set of constraints were used in a trial fitting to each overlapping peak obtained at different temperatures.

The d-spacing of {200} γ and {004} γ″ peaks were obtained by fitting the peaks from the bank 1 detector, and d-spacing of {200} γ and {200} γ″ peaks from the bank 2 detector, as shown in Figure 6. Such a strategy was employed to fitting peaks obtained at different temperatures, and lattice strain of both the γ and γ″ phases $\varepsilon_{hkl}^{i,t}$ due to heating was calculated by

$$ \varepsilon_{hkl}^{i,T} = \frac{{d_{hkl}^{i,T} - d_{hkl}^{i,0} }}{{d_{hkl}^{i,0} }} $$

(6)

where $d_{hkl}^{i, T}$ and $d_{hkl}^{i, 0}$ are the d-spacing of hkl planes for a phase measured at a temperature T and at RT, respectively. The lattice-strain plots against temperature are shown in Figure 7. All the plots show a polynomial relationship with temperature which agrees with the general polynomial thermal expansion of the material. The lattice strain along the c-axis of γ″ is apparently smaller than that of the γ phase, while lattice strain along the a-axis is very similar to that of γ phase.

Fig. 7

Lattice-strain evolution with temperature for both phases obtained from (a) bank 1 detector, (b) bank 2 detector. The lattice strain along the c-axis of γ″ is apparently smaller than that of the γ phase, while lattice strain along the a-axis is very similar to that of γ phase. Errors of the lattice strains associated with the peak fitting are about 20 microstrains for the γ phase and 150 microstrains for the γ″ phase. For the reason of clarity, error bars are not shown

4.4 Constrained Misfit Strain

The constrained misfit strain for the γ″ phase along a- and c-axes was calculated using Eq. [4], and is plotted against temperature in Figure 8. The calculated constrained misfit strains for the γ″ phase were positive along both a- and c-axes. The magnitude of the constrained misfit strain along the a-axis was much smaller compared to c-axis, which agrees with the general knowledge of lattice parameter of the γ″ phase.[5] The constrained misfit strain along a-axis remained relatively stable (0.10 to 0.14 pct) with temperature, while the constrained misfit strain along c-axis showed an apparent decline from 3.30 to 2.95 pct. Such a reduction in constrained misfit strain indicates an apparent drop in misfit stresses with increasing temperature.

Fig. 8

Misfit strain evolution with temperature for the γ″ phase along a- and c-axes. The constrained misfit strain along a-axis remained relatively stable (0.10 to 0.14 pct) with temperature, while the constrained misfit strain along c-axis showed an apparent decline from 3.30 to 2.95 pct. Error bars are shown for the uncertainties associated with fitting

4.5 Estimation of γ Phase Single Crystal Elastic Constants

Since the FCC γ matrix phase has anisotropic elastic properties, the single crystal elastic constant (SCEC) of the γ phase is indispensable for the calculation of the Eshelby tensor S. The SCEC of the γ phase has only been reported at RT, but not at elevated temperatures. Therefore, in order to calculate the Eshelby tensor S at different temperatures, the variation in SCEC with temperature needs to be estimated.

Young’s moduli at RT and 664 °C were obtained from the tensile experiments. Lattice parameter at each stress step was obtained by whole pattern fitting using the Pawley refinement.[23] Young’s moduli were derived from the slopes of the stress–lattice-strain plots as shown in Figure 9. Similarly, the hkl-specific DEC of the γ phase at RT and 664 °C were measured. hkl-specific lattice strain was obtained by single peak fitting and the DECs were derived from the slopes of the stress–lattice-strain plots in Figure 10.

Fig. 9

Young’s modulus measured by neutron diffraction at RT and 664 °C. Error bars associated with fitting uncertainties are smaller than the symbol

Fig. 10

Stress–lattice strain plots of the γ phase for aged IN718 at (a) RT and (b) 664 °C. hkl-specific DEC of the γ phase was obtained from the slope of each linear fit. Error bars associated with fitting uncertainties are smaller than the symbol

According to the elastic properties study of IN718 by Aba-Perea et al.,[17] it is reasonable to assume that the Young’s modulus decreases linearly with increasing temperature. As shown in Figure 9, the Young’s modulus dropped by 22.3 pct from RT (25 °C) to 664 °C; thereby, it is assumed that the SCEC dropped by the same ratio from RT to 664 °C.

The SCEC of the γ phase at RT reported by Dye et al.[24] was adopted in the current study (Table II). On the basis of the adopted RT SCEC and the linear relationship with temperature, the SCEC at 664 °C was calculated and listed in Table II.

Table II SCECs of the γ Phase in Aged IN718 at RT and 664 °C

The SCECs were subsequently justified by comparing the measured DECs to the calculated DECs from the SCECs using a program ‘DECcalc’.[25] The measured and calculated DECs are listed in Table III, showing relatively good agreements in DECs between the measurements and calculations. This proves that the estimated SCECs are appropriate. Therefore, the SCECs at various temperatures were estimated and plotted in Figure 11.

Table III Measured and Estimated DECs of the γ Phase at RT and 664 °C

Fig. 11

Estimated temperature-dependent single crystal elastic constants of the γ phase

4.6 Calculation of Eshelby tensor and misfit stress

According to Eshelby’s method, the Eshelby tensor in Eq. [3] can be obtained by solving a surface integral, although the calculation is lengthy and challenging if the matrix is not isotropic.[8] For the case of an ellipsoidal inclusion embedded in an anisotropic matrix, Kinoshita and Mura[26] and Lin and Mura[27] have derived expressions for the Eshelby tensor S by using the theory of residues to reduce the surface integral to a line integral.[10, 27] With the measured aspect ratio and estimated SCECs of the γ matrix, the Eshelby tensor S was calculated. The detailed calculation of the Eshelby tensor S is described in Appendix.

$$ S_{ijmn} = \left[ {\begin{array}{*{20}c} {0.1599} & {0.0201} & { - 0.0194} \\ {0.0201} & {0.1599} & { - 0.0194} \\ {0.3278} & {0.3278} & {0.8531} \\ \end{array} } \right] $$

(7)

The Eq. [3] can be rewritten as

$$ \left[ {\begin{array}{*{20}c} {\varepsilon_{11}^{\text{C}} } \\ {\varepsilon_{22}^{\text{C}} } \\ {\varepsilon_{33}^{\text{C}} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {S_{1111} } & {S_{1122} } & {S_{1133} } \\ {S_{2211} } & {S_{2222} } & {S_{2233} } \\ {S_{3311} } & {S_{3322} } & {S_{3333} } \\ \end{array} } \right] \times \left[ {\begin{array}{*{20}c} {\varepsilon_{11}^{T} } \\ {\varepsilon_{22}^{T} } \\ {\varepsilon_{33}^{T} } \\ \end{array} } \right] $$

(8)

where $\varepsilon_{11}^{\text{C}}$, $\varepsilon_{22}^{\text{C}}$, and $\varepsilon_{33}^{\text{C}}$ are constrained misfit strains measured at a specific temperature, subscripts 11 and 22 denote direction along a-axis, 33 denotes direction along c-axis. At RT, they were measured to be $\varepsilon_{11}^{\text{C}} = \varepsilon_{ 22}^{\text{C}} = 0.0014$, $\varepsilon_{33}^{\text{C}} = 0.033$. The SFTS $\varepsilon_{11}^{T}$, $\varepsilon_{22}^{T}$ and $\varepsilon_{33}^{T}$ at RT were derived by solving Eq. [8]:

$$ \left[ {\begin{array}{*{20}c} {\varepsilon_{11}^{T} } \\ {\varepsilon_{22}^{T} } \\ {\varepsilon_{33}^{T} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {0.0110} \\ {0.0110} \\ {0.0302} \\ \end{array} } \right] $$

(9)

The misfit stresses at RT were then calculated by Eq. [2]

$$ \left[ {\begin{array}{*{20}c} {\sigma_{11} } \\ {\sigma_{22} } \\ {\sigma_{33} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {C_{11} } & {C_{12} } & {C_{12} } \\ {C_{12} } & {C_{11} } & {C_{12} } \\ {C_{12} } & {C_{12} } & {C_{11} } \\ \end{array} } \right] \times \left[ {\begin{array}{*{20}c} {\varepsilon_{11}^{\text{C}} - \varepsilon_{11}^{T} } \\ {\varepsilon_{22}^{\text{C}} - \varepsilon_{22}^{T} } \\ {\varepsilon_{33}^{\text{C}} - \varepsilon_{33}^{T} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} { - 3.0\;{\text{GPa}}} \\ { - 3.0\;{\text{GPa}}} \\ { - 1.7\;{\text{GPa}}} \\ \end{array} } \right]. $$

(10)

Misfit stresses at different temperatures were calculated in the same manner with the temperature-dependent constrained misfit strains (Figure 8) and elastic constants (Figure 11). The results are plotted in Figure 12.

Fig. 12

Calculated misfit stresses along a- and c-axes of the γ″ phase at various temperatures. It is shown a larger compressive stress existed along a-axis compared to c-axis, both stresses declined with increasing temperature. Uncertainties were estimated from the uncertainties of measured constrained misfit strain in Figure 8

The results show that the γ″ precipitates were in a high level of compressive stress state: 3.0 GPa along a-axis and 1.7 GPa along c-axis at RT, which dropped with increasing temperature to 2.0 GPa and 1.2 GPa at 664 °C along a- and c-axes, respectively. Results also show that the magnitude of misfit stress was much higher along the a-axis than along the c-axis, or equivalently, along the in-plane direction and the normal-to-plane direction, respectively. Phase field simulations by Zhou et al.[7] have shown that compressive misfit stresses of 1 GPa level exist in γ″ phase in IN718 at 790 °C. The discrepancy between results of the current work and their phase field simulations could be due to the use of a smaller SFTS, and the assumption of identical elastic constants of the γ and γ″ phases in the simulations. However, the simulation by Lin et al.[18] found that the BCT-γ″ phase is stiffer than the FCC-γ phase, for example, the elastic modulus along the c-axis of the γ″ phase is 30 pct higher than the γ phase.

4.7 Misfit Stress in Matrix

Since macro residual stress in the specimen is assumed to be negligible, the equilibrium condition of average stresses in the two phases within the sampling volume requires[28]

$$ \sigma ^{{\gamma ^{\prime\prime}}} f^{{\gamma ^{\prime\prime}}} \; + \;\sigma ^{\gamma } f^{\gamma } \; = \;0 $$

(11)

where f^γ″ and f^γ are the volume fraction of the γ″ and γ phases, respectively. The volume fraction can be estimated from the {004} γ″/{200} γ peak intensity ratio $\frac{{I_{{004}}^{{\gamma ^{{\prime \prime }} }} }}{{I_{{200}}^{\gamma } }}$ in Figure 6(a) via the equation[16]:

$$ \frac{{I_{{004}}^{{\gamma ^{\prime\prime}}} }}{{I_{{200}}^{\gamma } }} = \frac{{\left( {\lambda _{{004}} } \right)^{4} \left( {\eta _{{\gamma ^{\prime\prime}}} } \right)^{2} f_{{\gamma ^{\prime\prime}}} }}{{\left( {\lambda _{{200}} } \right)^{4} \left( {\eta _{\gamma } } \right)^{2} f_{\gamma } }} $$

(12)

where $\lambda$ is the diffracted neutron beam wave length, $\eta$ the scattering length density. The intensity ratio was measured to be 1/8, the scattering length densities of the two phases, $\eta_\gamma \; = \;7.206\; \times \;10^{10} \;{\text{cm}}^{ - 2}$ and ${\eta_{\gamma ^{\prime\prime}}}\; = \;7.328\; \times \;10^{10} \;{\text{cm}}^{ - 2}$ were adopted from Reference 29. Therefore, the volume fraction of the γ″ phase was estimated to be 10 pct, with an uncertainty of ±1 pct due to the uncertainty of the measured intensity ratio. The average stresses in the γ matrix at various temperatures were calculated from Eq. [11] and the results are plotted in Figure 13. Misfit stresses of 329 MPa and 186 MPa in tension were found in the matrix along the a-axis and c-axis directions, respectively, at RT. As expected, the tensile misfit stresses decreased with increasing temperature, which were 224 MPa and 132 MPa along a-axis and c-axis at 664 °C, respectively.

Fig. 13

Average internal stresses in the γ matrix. It is shown a larger tensile stress existed along a-axis compared to c-axis, both stresses declined with increasing temperature. Uncertainties were estimated using Eq. [11] based on the uncertainties of misfit stresses in the γ″ phase, the uncertainty of estimated volume fraction would give a systematic discrepancy of 10 pct in stress level

4.8 Stress-Free Transformation Strain of the γ″ Phase

As mentioned in Section II–A, the SFTS $\varepsilon^{T*}$ for an inhomogeneity is different from the SFTS ε^T for an equivalent inclusion. The SFTS ε^T in Eq. [9] was calculated for the equivalent inclusion, and the SFTS $\varepsilon^{T *}$ of the γ″ phase can be calculated using Eq. [2]. Since no experimentally determined values can be found from literature, the SCECs of the γ″ phase is adopted from the literature which were derived by first principle calculations[18] and listed in Table IV. Similar to the justification of the SCECs of the γ phase in Section IV–E, the SCECs of the γ″ can be justified by comparing the calculated DECs to the measured DECs of the γ″. DECs of the γ″ have been measured using neutron diffraction in a previous study[30] as listed in Table V, along with the calculated DECs from the adopted SCECs using the program ‘DECcalc’.[25] The discrepancies between the calculated and measured DECs are generally acceptable and thereby, the adopted SCECs of the γ″ seem appropriate.

Table IV SCECs of the γ″ Phase for RT Adopted from Ref. [18]

Table V DECs of the γ″ Phase at RT Obtained by Experimental Measurements[30] and Calculations from the Adopted SCECs[18] Using the Matlab Program ‘DECcalc’ Applying the Kröner Model[25]

Using the elastic constants for the γ and γ″ phases in Tables II and V, and the equivalent SFTS ε^T in Eq. [9], the SFTS $\varepsilon^{T *}$ at RT was calculated using Eq. [2] to be

$$ \left[ {\begin{array}{*{20}c} {\varepsilon_{11}^{T*} } \\ {\varepsilon_{22}^{T*} } \\ {\varepsilon_{33}^{T*} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {0.0086} \\ {0.0086} \\ {0.0313} \\ \end{array} } \right] $$

(13)

These SFTS $\varepsilon^{T *}$ agree well with those values reported in literature for IN718 aged at 750 °C for 4 hours,[5] wherein the γ″ precipitates were believed to maintain high coherency with the matrix.

4.9 Thermal Expansion

The linear coefficient of thermal expansion (CTE) of the γ and γ″ phases can be obtained by the deviation of the polynomial plots in Figure 7. The variation of CTEs with temperature is shown in Figure 14, CTE of the γ″ phase along a-axis is similar to that of the γ phase, while CTE along c-axis of the γ″ phase is apparently smaller. This is consistent with Figure 7 which shows that the lattice strain along c-axis is less than that along a-axis and the γ phase. However, one should note that these CTEs were obtained in the constrained state, and thereby could be different from those measured in stress-free state when one phase is isolated from the other.

Fig. 14

Coefficients of thermal expansion of the γ and γ″ phases obtained from the deviation of the polynomial relationships in Figure 7. The CTE of the γ phase was the average value obtained from bank 1 and bank 2

The misfit stresses decreased by 30 pct from RT to 664 °C, meanwhile, the stiffness decreased by 22 pct, which indicates that about 70 pct of the decrease in misfit stresses is due to the decrease in stiffness and about 30 pct of the decrease is due to the different CTEs of the two phases.

5 Discussion

5.1 Coherent γ″ Precipitate

The use of Eshelby’s method requires the γ″ precipitate is fully coherent with the γ matrix. It is generally believed that the γ″ is highly coherent with the γ matrix when the precipitate size is relatively small. This can be revealed by the high level of coherency strain field around the precipitates that have been shown in TEM images in the literature,[5, 19, 31] on contrast, no misfit/interfacial dislocation around the precipitates have been reported to the authors’ best knowledge. Two situations may reduce the coherency: (i) the misfit is relatively large. For plate-shaped precipitate, when the misfit $\varepsilon_{33}^{T} > 0.1$, suggested by Weatherly,[32] a dislocation loop may form at the interface. In the current material, the misfit $\varepsilon_{33}^{T} \sim 0.03$, which is not large enough for the formation of dislocation loop; (ii) the precipitate grows beyond a critical size. It is logical to assume that the gamma double prime precipitates form initially fully coherent with the gamma matrix. When the γ″ precipitate grows beyond a critical size, formation of stacking faults within the precipitates may occur and trigger the loss of coherency.[33, 34] The critical diameter has been suggested to be > 90 nm by a few studies.[5, 34] The diameter of γ″ in the studied material is around 55 nm, which is much less than the critical diameter for loss of coherency. Therefore, it is reasonable to consider the γ″ precipitates maintaining a fully coherent bonding with the γ matrix, or perfectly glued with the matrix in the context of Eshelby’s theory.

5.2 High Level of Misfit Stress

The results show that the misfit stresses in the γ″ phase are very high, since they are hydrostatic in compression. Such high level of stresses would not lead to failure of the precipitates. The misfit stresses in the γ″ precipitate are highly anisotropic. The magnitude of misfit stress along the a-axis is much larger than that along the c-axis due to the anisotropic misfit strain. Eshelby’s theory[35] predicts that the misfit stresses in a plate-like precipitate decrease faster along the normal direction than that along the parallel direction of the habit plane when the precipitate grows.

Many studies have shown that large misfit stresses (or called coherency stresses or internal stresses) can exist in multi-phase materials when lattice misfit occurs due to phase transformation or thermal contraction. Thermal misfit stress with a magnitude of 800 MPa in compression was reported by Harjo et al. in the ferrite phase of a Fe-Cr-Ni alloy which consisted of 32 pct volume fraction of the ferrite phase and 68 pct of the austenite phase.[36] In their study, the diffraction measured stresses agreed fairly well with the stresses predicted using Eshelby’s theory. In the cementite with a volume fraction of 11 pct in the as-drawn steel wire, residual tensile stresses as high as 2000 MPa has been reported.[37] Withers et al.[9] also pointed out that the misfit stresses in the fiber for an aligned Al/SiC composite are reversely proportional to the volume fraction of the fiber. The internal stress in the fiber can be up to 1400 MPa compressive along the longitudinal direction when the volume fraction is 10 pct, while the compressive stress was only 600 MPa when the volume fraction was 40 pct.[9] Meanwhile, the average tensile stresses in the matrix decrease with decreased fiber volume fraction.[9] Recently, phase field simulations predicted that the misfit stresses in γ″ precipitates with a volume fraction of 16 pct in IN718 can reach 1 GP in compression at 790 °C.[7] All these studies showed that high level of misfit stresses can exist in the second phase with an apparently lower volume fraction compared to the matrix phase.

In IN718, although the misfit stresses in the γ″ precipitates are at a high level, the low volume fraction of the γ″ precipitates allows the average stress in the γ matrix to be kept at a relatively low level. The misfit strains were calculated from the lattice spacings which sampled from a relatively large sample volume. The γ lattice spacing in the vicinity of the γ″ precipitate can be much different from the measured average lattice spacing. The magnitude of local misfit stress is therefore much higher than that at a position in the matrix far from the precipitate. Lattice strains at different distances from the surface of a coherent precipitate measured by convergent beam electron diffraction showed that most of the strain is relieved in the matrix at a distance comparable to the size of the precipitate.[38] The presence of a high level of coherency strain around the γ″ precipitates is often observed in TEM bright field images[31, 39, 40] indicating large misfit stress in the vicinity of the γ/γ″ interface. Withers et al.[9] proposed that rather than the local fluctuations in stress, it is the mean stresses in the two phases that determines much of the mechanical behavior. The local stress will not promote the motion of dislocations, but can encourage the occurrence of microplasticity which is often responsible for the onset of plasticity.[9]

5.3 Loss of Coherency

It has been seen that a large magnitude of internal stress is established in the γ″ phase. Such high stresses can be partly relieved by relaxation mechanisms such as cracking, loss of coherency, interface gliding, local dislocation glide and diffusion.[9] The relaxation of internal stresses may have effects on the thermal stability of IN718, which is an important concern for the applications of IN718 at elevated temperatures. Thermal stability of IN718 largely depends on the resistance of phase transformation from the meta-stable γ″ phase to the equilibrium δ phase. Sundararaman et al.[33] suggested that such a phase transformation could be trigger by the formation of stacking faults in the γ″ precipitates, which is correlated to the loss of coherency.

Loss of coherency has been observed in the γ″ precipitates when they grow beyond a critical size. The critical size to coherency loss was proposed to be about 90-100 nm.[5, 34] The precipitates in this study with an average diameter about 56 nm, were believed to maintain full coherency with the matrix. The large misfit stress is expected to be relieved by the loss of coherency through different mechanisms correlated to the dimensions of the precipitates and the amount of dislocations at the interface.[32, 41] Sundararaman et al.[33] observed that it is the generation of stacking faults in the γ″ that leads to coherency loss. The experiment by Slama [5] showed that after aging for 50 hours, the $\varepsilon_{11}^{T}$ reduced significantly compared to that only aged for 4 hours, while the $\varepsilon_{33}^{T}$ showed little change, indicating that coherency loss happens along the habit plane. This would be expected since the magnitude of misfit stress $\sigma_{11}$ is much larger than $\sigma_{33}$ when the γ″ precipitates grow to a few tens of nanometers in diameter, the misfit stress $\sigma_{11}$ is more favorable to be relaxed. On the other hand, for a given particle size and aspect ratio, such findings indicate that a smaller level of internal stresses will be less favored for the loss of coherency, and thereby achieving a better resistance to phase transformation and thermal stability. According to Withers et al.,[9] a smaller level of internal stresses may be achieved by having a larger fraction of the γ″ phase; however, this is limited by the solubility of Nb in IN718. Another way to adjust the internal stresses is by changing the chemical composition, particularly the content of Nb, of which a higher content leads to a larger lattice misfit and a lower content to a less misfit along c-axis.[42] A profound thermal stability of IN718 can be achieved by forming the compact morphology which consist of a cube-shaped γ′ core coated by γ″ shell,[43] or by increasing the volume fraction of the γ′ phase comparable to that of the γ″ phase.[44]

6 Summary

In this study, the misfit stresses in IN718 as a function of temperature have been characterized using Eshelby’s inclusion theory. The key findings are:

• 1. The constrained misfit strains in a stress-induced variant selection heat treated IN718 are characterized using neutron diffraction at various temperatures, showing a very small positive constrained misfit strain along the a-axis (~0.1 pct) and a much larger constrained misfit strain along the c-axis (2.95 to 3.30 pct).

• 2. The coefficient of thermal expansion of the γ phase is found to be very similar to that of the γ″ phase along a-axis, while the coefficient of thermal expansion of the γ″ phase along c-axis is smaller compared to that along a-axis.

• 3. The calculated misfit stresses showed large compressive and anisotropic misfit stresses in the γ″ precipitates: 3.0 GPa along the a-axis and 1.7 GPa along the c-axis at RT, which reduced to 2.0 GPa and 1.2 GPa at 664 °C, respectively.

• 4. Tensile stresses calculated in the matrix showed a magnitude of ~ 329 MPa at RT and ~ 186 MPa at 664 °C.

This study reports a simple and reasonably accurate method to calculate misfit stresses from the constrained misfit strain obtained by in-situ diffraction experiment showing a large level of compressive misfit stresses in the highly coherent precipitates. This work is expected to deepen the understanding of misfit stress which plays an important role in the precipitation strengthening alloys. The proposed method can be extended to in-situ heat treatment or creep studies.

Appendix

For the case of inclusions embedded in an anisotropic matrix, explicit expressions for anisotropic Green’s functions are not available, the technique used in calculations of Eshelby tensor $S$ in isotropic matrix is not applicable. Mura[10,27,28] gave explicit expressions for the cases of spheroid inclusion in cubic crystals. When the inclusion has an ellipsoidal shape, the stress field inside the inclusion is constant and can be expressed in the form of a surface integral on a unit sphere. The surface integral can be reduced to a line integral from the surface integral and is ready for numerical calculations.

Eshelby tensor S can be expressed

$$ S_{ijmn} = \left( {\frac{1}{8\pi }} \right)C_{pqmn} \left( {\overline{G}_{ipjq} + \overline{G}_{jpiq} } \right) $$

(A1)

For the spheroid inclusion which principal directions are coincident with the matrix crystalline directions, the nonzero components of $\overline{G}_{ijkl}$ are:

$$ \overline{G}_{1111} \; = \;\overline{G}_{2222} \; = \;\frac{2\pi }{a}\mathop \smallint \limits_{0}^{1} \frac{{1 - x^{2} }}{pq}\left( {1 - x^{2} + \rho^{2} x^{2} } \right)\left[ {\mu^{2} (1 - x^{2} + \rho^{2} x^{2} ) + \beta \rho^{2} x^{2} } \right]{\text{d}}x + \frac{\pi }{a}\mathop \smallint \limits_{0}^{1} \frac{{\left( {1 - x^{2} } \right)^{2} }}{{p\left( {p + q} \right)}}\left[ {\beta \left( {1 - x^{2} + \rho^{2} x^{2} } \right) + \gamma \rho^{2} x^{2} } \right]{\text{d}}x $$

(A2a)

$$ \bar{G}_{3333} = \frac{4\pi }{a}\mathop \smallint \limits_{0}^{1} \frac{{\rho^{2} x^{2} }}{pq}\left( {1 - x^{2} + \rho^{2} x^{2} } \right)\left[ {\mu^{2} (1 - x^{2} + \rho^{2} x^{2} ) + \beta (1 - x^{2} )} \right]{\text{d}}x + \frac{\pi \gamma }{a}\mathop \smallint \limits_{0}^{1} \frac{{\rho^{2} x^{2} \left( {1 - x^{2} } \right)}}{{p\left( {p + q} \right)}}{\text{d}}x $$

(A2b)

$$ \bar{G}_{1122} \; = \;\bar{G}_{2211} \; = \;\frac{2\pi }{a}\mathop \smallint \limits_{0}^{1} \frac{{1 - x^{2} }}{pq}\left\{ {\left( {1 - x^{2} + \rho^{2} x^{2} } \right)\left[ {\mu^{2} (1 - x^{2} + \rho^{2} x^{2} ) + \beta \rho^{2} x^{2} } \right] + \left( {1 - x^{2} } \right)\left[ {\beta \left( {1 - x^{2} + \rho^{2} x^{2} } \right) + \gamma \rho^{2} x^{2} } \right]} \right\}{\text{d}}x - \frac{\pi }{a}\mathop \smallint \limits_{0}^{1} \frac{{\left( {1 - x^{2} } \right)^{2} }}{{p\left( {p + q} \right)}}\left[ {\beta \left( {1 - x^{2} + \rho^{2} x^{2} } \right) + \gamma \rho^{2} x^{2} } \right]{\text{d}}x $$

(A2c)

$$ \overline{G}_{1133} = \overline{G}_{2233} = \frac{2\pi }{a}\mathop \smallint \limits_{0}^{1} \frac{{\rho^{2} x^{2} }}{pq}\left\{ {2\left( {1 - x^{2} + \rho^{2} x^{2} } \right)\left[ {\mu^{2} (1 - x^{2} + \rho^{2} x^{2} ) + \beta \rho^{2} x^{2} } \right] + \left( {1 - x^{2} } \right)\left[ {\beta \left( {1 - x^{2} + \rho^{2} x^{2} } \right) + \gamma \rho^{2} x^{2} } \right]} \right\}{\text{d}}x $$

(A2d)

$$ \overline{G}_{1212} = - \frac{{\pi \left( {\lambda + \mu } \right)}}{a}\mathop \smallint \limits_{0}^{1} \frac{{\left( {1 - x^{2} } \right)^{2} }}{{p\left( {p + q} \right)}}\left\{ {\left[ {\mu (1 - x^{2} + \rho^{2} x^{2} ) + \beta \rho^{2} x^{2} } \right] + \mu^{\prime}\rho^{2} x^{2} } \right\}{\text{d}}x $$

(A2e)

$$ \bar{G}_{1313} = \bar{G}_{2323} = - \frac{{2\pi \mu \left( {\lambda + \mu } \right)}}{a}\mathop \smallint \limits_{0}^{1} \frac{{\rho^{2} x^{2} \left( {1 - x^{2} } \right)\left( {1 - x^{2} + \rho^{2} x^{2} } \right)}}{pq}{\text{d}}x - \frac{{\pi \mu^{\prime}\left( {\lambda + \mu } \right)}}{a}\mathop \smallint \limits_{0}^{1} \frac{{\rho^{2} x^{2} \left( {1 - x^{2} } \right)^{2} }}{{p\left( {p + q} \right)}}{\text{d}}x $$

(A2f)

$$ \bar{G}_{3311} \; = \;\bar{G}_{3322} \; = \;\frac{2\pi }{a}\mathop \smallint \limits_{0}^{1} \frac{{1 - x^{2} }}{pq}\left( {1 - x^{2} + \rho^{2} x^{2} } \right)\left[ {\mu^{2} (1 - x^{2} + \rho^{2} x^{2} ) + \beta \left( {1 - x^{2} } \right)} \right]{\text{d}}x\; + \;\frac{\pi \gamma }{2a}\mathop \smallint \limits_{0}^{1} \frac{{\left( {1 - x^{2} } \right)^{3} }}{{p\left( {p + q} \right)}}{\text{d}}x $$

(A2g)

where:

$$ \lambda = C_{12} $$

(A3a)

$$ \mu = C_{44} $$

(A3b)

$$ \mu^{\prime} = C_{11} - C_{12} - 2C_{44} $$

(A3c)

$$ \rho = a_{1} /a_{3} $$

(A3d)

$$ a = \mu^{2} \left( {\lambda + 2\mu + \mu^{\prime}} \right) $$

(A3e)

$$ b = a^{ - 1} \mu \mu^{\prime}\left( {2\lambda + 2\mu + \mu^{\prime}} \right) $$

(A3f)

$$ c = a^{ - 1} \mu^{\prime2} \left( {3\lambda + 3\mu + \mu^{\prime}} \right) $$

(A3g)

$$ \beta = \mu \left( {\lambda + \mu + \mu^{\prime}} \right) $$

(A3h)

$$ \gamma = \mu^{\prime}\left( {2\lambda + 2\mu + \mu^{\prime}} \right) $$

(A3i)

$$ p = \left\{ {\left( {1 - x^{2} + \rho^{2} x^{2} } \right)^{3} + b\rho^{2} x^{2} \left( {1 - x^{2} } \right)\left( {1 - x^{2} + \rho^{2} x^{2} } \right) + \frac{1}{4}\left( {1 - x^{2} } \right)^{2} \left[ {b\left( {1 - x^{2} + \rho^{2} x^{2} } \right) + c\rho^{2} x^{2} } \right]} \right\}^{1/2} $$

(A3j)

$$ q = \left\{ {\left( {1 - x^{2} + \rho^{2} x^{2} } \right)^{3} + b\rho^{2} x^{2} \left( {1 - x^{2} } \right)\left( {1 - x^{2} + \rho^{2} x^{2} } \right)} \right\}^{1/2} $$

(A3k)

here $C_{11}$, $C_{12}$ and $C_{44}$ are the Voigt elastic constants of the γ matrix, which are equivalent to $C_{1111}$, $C_{1122}$ and $C_{2323}$, respectively. $a_{1}$ is the average major axis and $a_{3}$ the average minor axis of the γ″ precipitates. The elastic constants are estimated in Section IV–E, the average dimension is measured in Section IV–A. Therefore, the nonzero components of $\overline{G}_{ijkl}$ are calculated as:

$$ \overline{G}_{1111} = \overline{G}_{2222} = 3.1784 $$

(A4a)

$$ \overline{G}_{3333} = 12.4511 $$

(A4b)

$$ \overline{G}_{1122} \; = \;\overline{G}_{2211} \; = \;3.0154 $$

(A4c)

$$ \overline{G}_{1133} \; = \;\overline{G}_{2233} \; = \;0.9799 $$

(A4d)

$$ \overline{G}_{1212} \; = \; - \;0.5709 $$

(A4e)

$$ \overline{G}_{1313} \; = \;\overline{G}_{2323} \; = \; - 1.6127 $$

(A4f)

$$ \overline{G}_{3311} \; = \;\overline{G}_{3322} \; = \;3.0804. $$

(A4g)

From Eq. [A.1], with the estimated elastic constants $C_{pqmn}$ and calculated $\overline{G}_{ijkl},$ the Eshelby tensor for the case of anisotropic matrix at RT is obtained:

$$ S_{ijmn} = \left[ {\begin{array}{*{20}c} {0.1599} & {0.0201} & { - 0.0194} \\ {0.0201} & {0.1599} & { - 0.0194} \\ {0.3278} & {0.3278} & {0.8531} \\ \end{array} } \right]. $$

(A5)