Dissolution of the Primary γ′ Precipitates and Grain Growth during Solution Treatment of Three Nickel Base Superalloys

Second phase particles (SPP) play an essential role in controlling grain size and properties of polycrystalline nickel base superalloys. The understanding of the behavior of these precipitates is of prime importance in predicting microstructure evolutions. The dissolution kinetics of the primary γ′ precipitates during subsolvus solution treatments were investigated for three nickel base superalloys (René 65, AD730 and N19). A temperature-time codependency equation was established to describe the evolution of primary γ′ precipitates of each material using experimental data, the Thermo-Calc software and the Johnson–Mehl–Avrami–Kolmogorov (JMAK) model. The dissolution kinetics of precipitates was also simulated using the level-set (LS) method and the former phenomenological model. The precipitates are represented using an additional LS function and a numerical treatment around grain boundaries in the vicinity of the precipitates is applied to reproduce their pinning pressure correctly. Thus, considering the actual precipitate dissolution, these simulations aim to predict grain size evolution in the transient and stable states. Furthermore, it is illustrated how a population of Prior Particle Boundaries (PPB) particles can be considered in the numerical framework in order to reproduce the grain size evolution in the powder metallurgy N19 superalloy. The proposed full-field strategy is validated and the obtained results are in good agreement with experimental data regarding the precipitates and grain size.


Introduction
In polycrystalline nickel base superalloys, the γ phase is found in several types of precipitates that can basically be distinguished based on their size and specific role in the microstructure. The so-called primary precipitates are the largest ones, typically a few microns (µm). Secondary and tertiary ones are much finer and essentially control the mechanical behavior [1][2][3]. After forging operations, polycrystalline γ-γ superalloys are submitted to partial solution treatment, thus below the γ solvus temperature (T solvus ), where primary γ precipitates are kept to avoid excessive grain growth (GG), while the fine hardening precipitates dissolved to reprecipitate again with an optimized size distribution during the subsequent controlled cooling.
Primary γ precipitates control the grain size by hindering grain boundary (GB) motion through the Smith-Zener pinning (SZP) mechanism [4][5][6][7]. In order to properly control the grain size during a specific subsolvus solution treatment, the evolution of these primary precipitates as a function of time and temperature must be understood.
Such simulations may help clarify, predict, and control microstructure evolution when precipitate dissolution occurs during heat treatments (HT), at subsolvus and near-solvus temperatures for each considered material. In the present work, only the effect of primary γ precipitates will be considered, as secondary and tertiary ones are dissolved during subsolvus solution treatments.
All these full-field approaches can integrate SPP into their formulation. A spherical shape to describe them is commonly adopted, although the influence of more complex shapes has also been investigated [28,31,34]. In some cases, real precipitate morphology, obtained from Electron BackScatter Diffraction (EBSD) maps, secondary electron (SE) images [31,32] or backscattered electron (BSE) images [31], has been used. However, all these works consider the precipitates as static objects. This reasoning is valid at low temperature (relatively to T solvus ) where precipitates do not evolve (or very slowly) or at high temperature where precipitates dissolve very fast, in which case SZP does not need to be modeled at all. However, at a temperature just below T solvus , precipitate dissolution and grain growth can be concomitant. Moreover, the Ostwald ripening mechanism can also take place at long annealing times.
For high temperatures close to T solvus , the precipitate volume fraction and size distribution can evolve significantly at the same timescale of recrystallization and GG phenomena, which motivates this work focused on GG, where the SZP plays a role of prime importance.
The SZP-GG coupling exhibits complex behaviors [9,31,35] that include precipitate evolution. To our knowledge, few numerical works consider the precipitate evolution, some of them using the PF method [30,36] and other using the LS approach [13]. This work adopts the LS approach presented in [13], where a special treatment at the precipitate/GB interface is made to model incoherent precipitates (i.e., assuming the precipitate/GB contact angle to be 90 • ).
The studied materials are presented in Section 2 as well as the experimental methods used to get data regarding the evolution of precipitates and matrix grains during annealing at different temperatures and times. A phenomenological model is then established in Section 3 for describing the evolution of precipitate area fraction and size as a function of time and temperature by combining Thermo-Calc simulations and a JMAK model [11,12]. Section 3 also presents how the parameters which will be used for describing GB motion in the subsequent full-field simulations have been obtained from the experimental grain growth kinetics. Both the phenomenological model for precipitate evolution and GB motion parameters are then used in full-field simulations conducted within an LS framework. Section 4 introduces the LS method and the required numerical treatment for considering precipitate evolution as a function of the dissolution velocity. The results of the LS fullfield simulations are presented in Section 5 and compared with the experimental data for validation. The last section summarizes the conclusions and opens perspectives for further works.

Materials and Methods
Three nickel base superalloys, René 65, AD730 and N19, were submitted to heat treatments to study the dissolution kinetics of primary γ precipitates. Table 1 shows their chemical compositions. These are materials of high technological interest as they are candidates for high temperature applications in aeronautical components, where strict grain size criteria must be fulfilled.

As-Received States
The initial ("as-forged") state of the samples studied in this work come either from a conventional cast-and-wrought process (for AD730 and René 65) or from a powder metallurgy route and isothermal forging (for the N19). Figure 1 illustrates the initial microstructure of each alloy. Images on the left with low magnification show the distribution of the primary γ precipitates in the microstructure. The images on the right present the morphology of primary and secondary γ populations. The initial microstructures show noticeable differences, even between AD730 and René 65 alloys, which have quite similar chemical compositions and are both cast-and-wrought alloys. This highlights the importance of the thermomechanical background on the fraction and distribution of γ precipitates in the microstructure. The presence of secondary and tertiary precipitates in the microstructure is notably determined by the last cooling undergone by the material, but they will be dissolved again at the targeted temperature in this work.

Thermo-Calc Simulations and Experimental Data
Thermo-Calc software (TCNI9 database, Thermo-Calc, Stockholm, Sweden) [40] was used to estimate the γ phase volume fraction of each superalloy as a function of temperature by using the TCNI9 database and the compositions given in Table 1. Figure 2 presents the volume fraction of γ phase as a function of temperature for the three alloys. The γ volume fraction estimated at room temperature was calculated as 37% for René 65 and AD730 and 43% for N19, which are in agreement with the values presented in the literature [38,41,42]. According to the curves in Figure 2, the full γ phase dissolution temperatures (T solvus ) estimated by Thermo-Calc correspond to 1093 • C, 1108 • C and 1138 • C for René 65, AD730 and N19, respectively.

Methods
Heat treatments were performed at different times and temperatures for all three materials (summarized in Table 2), followed by cooling at 100 • C/min or water quenching. The evolution of primary γ precipitates was considered to be negligible during cooling at 100 • C/min, as confirmed by comparison with the primary γ evolution obtained after water quenching.
Sample preparation for microstructure analysis was carried out in two different ways to obtain high-contrast images between the γ matrix and the γ precipitates. The first preparation protocol was a mechanical polishing down to 1 µm grid and then electrolytic polishing with a solution of 10% perchloric acid and 90% methanol at 10 • C, under 35 V for 5-8 s. The second was a mechanical polishing down to 0.25 µm grid and etching (with citric acid or HCl/CH 3 COOH).   Microstructural analysis was performed with two microscopes. The samples that received the first protocol were examined using the secondary electron detector (SE) of a Zeiss Supra40 Scanning Electron Microscope (SEM) (Zeiss, Oberkochen, Germany). SE images were acquired with an acceleration voltage of 3 kV, at 3-5 mm working distance and with 60 µm diaphragm size. The samples prepared with the second protocol were examined by optical microscopy. Ten images of 100 µm × 75 µm (for AD730 and N19) and 227 µm × 170 µm (for René 65) were taken for each sample for the sake of statistical representativity. In addition, 350-400 precipitates/image were found in René 65 microstructure after annealing at temperatures below 1050 • C, i.e., 60 • C below T solvus . The higher the annealing temperature, the lower the number of precipitates found, reaching 140-230 precipitates/image at 1090 • C. The same behavior was observed in AD730 and N19 alloys images, where the number of precipitates decreases as the temperature increases. ImageJ software (1.52a, LOCI, University of Wisconsin, Madison, WI, USA) [43] was used for image binarization and quantitative analysis of the γ phase area fraction and average precipitate size.
EBSD maps were acquired with an acceleration voltage of 20 kV, at 15-16 mm working distance, and 120 µm diaphragm size. They were analyzed with MTEX [44], an open MATLAB toolbox (R2017a, MathWorks, Natick, MA, USA). A GB was defined when the disorientation between two neighboring points exceeded 10 • . Regarding the γ phase identification, it cannot be done with EBSD data because this phase is indistinguishable from the γ matrix [45] because of their crystallographic structure being too close to each other. Therefore, in the case of AD730, a size criterion has been applied to separate the particles from the grains. The particles with 0.3 µm ≤ D spp ≤ 4 µm were considered as γ primary precipitates. The step size of the maps was adapted; as the annealing temperature increased, the grains became larger and larger. For instance, for the 4 h HT at 1070 • C, an EBSD map of size 1130 µm × 840 µm with a step size of 1.13 µm was analyzed, giving a total number of 5120 grains. For the 4 h HT at 1100 • C, two EBSD maps were taken to reach a good statistical representativity of the sample; the size of the maps was 2180 µm × 1700 µm with a step size of 4 µm, and they included a total number of 2272 grains.
For the N19 material, it was not possible to distinguish the primary γ precipitates from the grains (γ phase) with the same strategy because the primary γ precipitates have a size similar to that of the grains. Therefore, EBSD analysis was coupled with energy dispersive X-ray (EDS) analysis [46] to distinguish the two phases based on their chemical compositions. However, one drawback of this technique is that the spatial resolution of the EDS is limited to about 0.5 µm under the current acquisition setup. For these EBSD analyses, two maps of size 120 µm × 90 µm with a step = 0.08 µm were analyzed, giving a total number of around 1600 grains for the samples annealed at 1110 • C and 1200 grains at 1130 • C.  Figure 3 shows micrographs after the four-hour isothermal HT carried out at different temperatures for each alloy. Note that the γ precipitates of each alloy have different morphologies, being directly related to the previous thermomechanical treatments. However, a common point is the dissolution process itself, which is the reduction of the volume fraction of primary γ precipitates with increasing temperature. In addition, as is visible from the crystallographic contrast of some images, the primary γ dissolution leads to an increase in the grain size (Figure 3b,c, for instance).  The smaller the primary γ fraction in the microstructure, the larger the average grain size. This result has already been proven numerous times in different alloys [2,6,41,[47][48][49] and can be explained by the SZP phenomenon. Figure 4 shows the evolution of γ phase fraction as a function of T. The primary γ fraction significantly evolves (dissolution mechanism) when the temperature is close to the solvus temperature (50 • C −T solvus ≤ T ≤ T solvus ), and remains constant when the temperature is far from the solvus point, as illustrated in Figure 4a,c. The presence of secondary and tertiary γ precipitates was not considered in the measurement of the γ phase fraction. However, they are given reminders on these graphs with labels γ I I and γ I I I , respectively. An experimental or effective T solvus could be determined as the temperature at which the area fraction of primary γ is equal to zero after 4 h isothermal holding. Figure 4 also exhibits a comparison between the Thermo-Calc curves and the experimental data.

Area Fraction of Primary γ Precipitates after Annealing for 4 h
Note that thermodynamic simulation values do not match exactly with the experimental T solvus obtained for the René 65 and N19 alloys. For temperatures close to T solvus , the Thermo-Calc curve and the experimental points seem to follow the same evolution and shape. For the AD730 alloy (Figure 4b), there is a much better fit between experimental points and the Thermo-Calc results; however, for the other two alloys (Figure 4a,c), there is a slight offset. A shift of the Thermo-Calc curve was done to find a better fit with the experimental points, an offset of 17 • C and 12 • C were applied to the René 65 and N19 alloys, respectively. From these results, it was defined that the effective T solvus values for 4 h annealing corresponds to 1110 • C, 1108 • C and 1150 • C for the René 65, AD730 and N19 superalloys, respectively. The shifted curves can be described as a function of T − T solvus difference, as presented in Figure 5. The same curve can then be used to reproduce the experimental behavior of the René 65 and AD730, for temperatures in the 50 • C range below the T solvus . This shows that the overall area fraction evolution is similar for both conventional forged materials at the equilibrium state. Figure 5. Area fraction of primary γ precipitates after a four-hour isothermal treatment and the respective curves obtained through phenomenological fits established on the shape of the Thermo-Calc curves (described in Equations (1) and (2)).
The dissolution kinetics, i.e., the slope of the curve, for the N19 is steeper compared to the René 65 and AD730 alloys. This result is physically relevant as the annealing temperatures and T solvus are higher for this alloy. The equations, i.e., phenomenological models that describe the behavior presented in Figure 5 are Equation (1) for René 65 and AD730 alloys and Equation (2) for N19; they were obtained as second order polynomial approximations of the shifted curves: where f spp eq is the equilibrium volume fraction considered to be reached after 4 h and actually corresponds to the measured primary γ precipitate fraction as secondary and tertiary precipitates are already dissolved at about 50 • C below the solvus point. The shape of the shifted curves, as well as the values of the obtained second order coefficients, illustrates that linear approximations, although less precise, would be acceptable.

Area Fraction of Primary γ Precipitates as a Function of Time and Temperature
The main purpose of this part is to establish a model that can predict the area fraction of the primary γ phase as a function of time and temperature within the 50 • C window below T solvus . Through experimental data and Thermo-Calc, two equations that reproduce the fraction of primary γ as a function of temperature were obtained in Section 3.1 (Equations (1) and (2)). A JMAK model was then used to describe the diffusion-controlled dissolution of spherical precipitates. This relation has already been used in several works to represent the dissolution of γ phase for polycrystalline nickel base superalloys [11,12,50], as follows: The area fraction of precipitates at thermodynamic equilibrium is represented by f eq , f 0 corresponds to the primary γ area fraction at t = 0 s (as-received material) and t 1 is a time-dependent constant.
Assuming that Equation (3) is valid to describe the γ phase dissolution of the three alloys, and that Equations (1) and (2) are valid over the range 50 − T solvus ≤ T ≤ T solvus to provide the equilibrium area fraction of primary γ precipitates ( f spp eq alloy = f eq ), the following expression is obtained: This time dependent model was validated through several heat treatments ranging from 10 min to 4 h at different temperatures for each of the three studied alloys. The initial primary precipitates area fraction f 0 and the identified value of the t 1 constant were respectively 12.5% and 1300 s for René 65 and AD730 and 18.5% and 250 s for N19. It is worth mentioning that the constant t 1 , once identified for the best fit of one curve (one temperature), is in fact quite precise for all temperatures in the analyzed range, contrary to the work presented in [11], where the parameter t 1 needs to be identified for each temperature. Figure 6 shows a good agreement between the experimental results and Equation (4) model predictions. The evolution of the γ phase area fraction could then be described in the considered range of time and temperature. A few data available for the René 65 alloy have also been added on the plot of Figure 6a.  Figure 7 illustrates the evolution of the mean primary precipitate size (D spp ) for the three alloys. During short annealing times, the average precipitate size decreases, as could be expected from the dissolution process. However, after 1 h of heat treatment, a slight increase in the mean size of the precipitate can be seen. This behavior suggests the presence of the Ostwald ripening phenomenon, for which large precipitates grow and the small ones shrink and disappear, while keeping the particle fraction constant. The driving force of this mechanism is the minimization of interfacial energy (Gibbs Thomson effect). A composition gradient in the matrix between precipitates of different curvature radii induces a diffusion flux from the small precipitates towards the large ones leading to the shrinkage and disappearance of small precipitates and eventually to a decrease of the free energy associated with the overall matrix-precipitate interface area. .

Precipitate Dissolution Rate Derived from the JMAK Model
The precipitate area fraction as a function of time and temperature has already been established using a JMAK model given by Equation (3) where the precipitate area fraction is described as a function of the initial and the equilibrium fractions.
The precipitate dissolution rate v spp (defined here as the rate at which the precipitate interface moves) can then be described using the following relation: with S the total considered area, r spp (t, T) the arithmetic mean precipitate radius and N spp (t, T) the remaining number of precipitates (at each step of the simulations). This velocity expression will be introduced later in the full-field simulations to reproduce the precipitate interface velocity observed experimentally. Figure 8a illustrates the evolution of the experimental arithmetic mean grain size for the AD730 and René 65 alloys. The grain size does not change significantly when the annealing is performed at a temperature far from T solvus , as the grains are blocked by the rather high γ phase fraction present in the initial material (12.5%). However, for higher temperatures (>1080 • C), where the dissolution of the precipitates is more significant, the grains may grow more and more, especially when the microstructure is free from precipitates (at T solvus = 1110 • C).

Grain Size Evolution
For the N19 alloy (Figure 8b), a similar behavior is observed, with only slight grain growth at temperatures far below T solvus . However, for temperatures close to T solvus , a very significant evolution is noted in the first minutes, but after 1 h heat-treatment, the average grain size tends to keep the same value. The N19 superalloy behavior is thus different from that of the AD730 for temperatures close to the solvus domain. It will be necessary to understand why the mean grain size of this PM material behaves differently and then considers this effect in the simulations to compare the experimental data with the numerical results.

Experimental Data Used for the Identification of GB Migration Parameters
The procedure used to obtain the values of the GB migration parameters consistent with the experimental data for each material is presented below. This identification is mandatory for the forthcoming comparison between the experimental data and the numerical results.
The main parameter controlling the migration of a GB submitted to a driving force is its mobility, which can be approximated through an Arrhenius law as follows: with M 0 a pre-exponential factor considered to be constant, Q m the activation energy for GB migration, T the absolute temperature, and R the gas constant. Here, the mobility is considered to be the same for all GB, independently of crystallographic properties. The values of M, M 0 and Q m can be obtained from the experimental GG data when there are no precipitates present in the material to hinder GB motion, i.e., at the supersolvus domain.
The procedure used to calculate M 0 and Q m (Equation (6)) can be summarized by the following steps [51]:

1.
Supersolvus heat treatments were performed at different temperatures for several holding times for the AD730, up to 50 • C above the T solvus and for the N19 material, up to 30 • C above T solvus as presented in Figure 9. It can be noticed again here in the supersolvus domain that the grain size of the N19 alloy stagnates at about 20 µm, whereas the grain size of the two other alloys keeps increasing with time and temperature.

2.
A first approximation of the reduced mobility (Mγ BT ) is obtained for each temperature by the best fit of the experimental data with a Burke & Turnbull law (B&T) [52]: with γ the grain boundary energy set here and in the following to 0.6 J/m 2 , which corresponds to GB energy in pure nickel at 1000 • C [53].

3.
As the reduced mobility (Mγ) is actually model-dependent [51], these values will have to be further refined based on full-field simulations as described in Sections 5.1.1 and 5.2.1.

Full-Field Model Description
In the following, the LS method is used to perform full-field GG simulations under the influence of the primary γ precipitates submitted to dissolution as described with the formerly established phenomenological model. First, the global numerical framework is introduced below.

Level-Set Approach to Simulate GG
In this approach, each grain G (sub-domain) of the calculation domain Ω (polycrystal) can be defined by computing a signed Euclidean distance function (ϕ) to its interface. This function is considered to be positive inside the grain and negative outside, and the interface (i.e., grain boundary) is then described by the zero-isovalue as follows: where d(x, Γ(t)) is the distance between a point x and the boundary Γ(t) of the considered grain G. The interface migration can be deduced by the resolution of a set of convective equations, with N G the number of grains: When GG is studied, it is classically accepted that, at the mesoscopic scale, the behavior of ϕ is only driven by the curvature flow [54]. Thus, the grain boundary velocity can be approximated by: where M i is the local grain boundary mobility, γ i the local grain boundary energy, κ i = ∇ · n i the local mean curvature (i.e., the curvature in 2D and the sum of the main curvatures in 3D) and n i = − ∇ϕ i ∇ϕ i the outward unitary normal vector to Γ i = ∂G i . Classically, by using Equation (10) and if the LS functions remain distance functions at each time step ( ∇ϕ i = 1), Equation (9) can be rewritten as a set of diffusive equations [55] avoiding the exact calculation of the GB curvature: For this work, the material is considered as isotropic, i.e., M i = M(T) is assumed to be only temperature-dependent and γ i = γ is considered as constant [56]. However, LS can also be applied to more realistic microstructures considering the heterogeneity of Mγ, taking into consideration the misorientation angle [57][58][59][60][61], the inclination of the grain boundaries [60,62,63] and the influence of solute segregation [64].
In the considered isotropic assumption context and at a high temperature, the mobility can be approximated considering the Arrhenius law presented in Equation (6).
To summarize, the GG simulation under the effect of capillary pressure in the proposed finite element framework (FE) coupled to the LS method involves the following steps at each time increment: • Grain boundary migration is calculated by solving Equation (11) for each active LS function. In order to limit the computational cost, a Graph coloration/recoloration technique is used [65] allowing for reducing the number of LS functions by grouping non-neighboring grains in common LS functions, called global Level-Set (GLS) functions. Each GLS contains numerous grains that are separated by a minimum distance, hence N GLS N G [65]. • After the resolution of Equation (11), voids or overlaps can appeared at the GB and their junctions. In order to treat this issue and remove such non-physical regions, the methodology proposed in [57] and described by Equation (12) is used for each GLS function:φ whereφ i is used as the new corrected GLS function. • All active GLS functions are reinitialized ( ∇ϕ = 1) with the direct reinitialization method proposed in [66]. Indeed, one of the weaknesses of the LS formulation is to not naturally conserve the metric property of an LS function initially defined as a distance function. To ensure the validity of the diffusive form of Equation (11) and to describe curvature flow migration, distance functions are required at least near the interfaces. Several approaches have been developed to treat the reinitialization procedure in regular grids, or unstructured FE meshes [67][68][69]. Here, the parallel and direct reinitialization algorithm proposed by Shakoor et al. [66] was adopted as a fast and accurate method. Discussions concerning the residual error of this approach are detailed in [70]. • The negative GLS functions (grains which have disappeared) are removed from the system of equations. • To avoid numerical coalescence, the re-coloring technique presented in [32,65] is applied. • Interface remeshing operations following the methodologies proposed in [31,71] are performed if required.

Description of SPP
The multiple junctions treatment based on [57] has been applied both for 2D and 3D simulations of GG in single phase materials or with static SPP using the LS approach [32,55,56,72]. A Modified Multiple Junctions Treatment (MMJT) has recently been proposed in [13] to consider the precipitates as an LS function in FE-LS GG modeling. This new description of the particles opens the possibility of making SPP evolve. The MMJT can be summarized as follows: Without SPP, the MMJT is equivalent to the classical numerical treatment (Equation (12)). When SPP are present, it enables by successive iterations to impose the Young-Herring equilibrium for incoherent SPP (see [13] for further details).
At the beginning of the simulation, an initial mesh is generated, followed by the polycrystal creation. The latter must be statistically representative of the initial state of the studied material. Then, the grain distance functions ϕ i are adjusted to introduce the new particle distance function ϕ spp .

SPP Interface Velocity
A convective equation is considered for ϕ spp to make the precipitate interfaces evolve. A SPP interface velocity, v spp , is introduced at the SPP interfaces in order to reproduce the precipitate dissolution given by the JMAK model previously developed. More precisely, a continuous velocity field, v, is applied through the Laplacian equation Equation (14) considering Dirichlet boundary conditions [73] at the precipitate interfaces. In the end, this smoothed velocity field v is implemented to calculate the velocity v oriented towards the center of each precipitate (see Figure 10). and with n the unitary inside normal vector to the SPP, v spp the velocity imposed to the SPP interfaces (Equation (5)) and v the resulting velocity field imposed through a convection equation to ϕ spp (see [13]).

Full-Field Simulation for the Identification of GB Parameters and GG
Before developing the GG simulations, the parameters describing the GB migration of the considered material must be identified for each alloy to provide the best agreement between the numerical results and the experimental data.

Parameter Identification and Validation
Section 3.6 describes the procedure used to calculate M 0 and Q m , the steps are more detailed below.
Supersolvus heat treatments were performed between 1120 • C and 1160 • C for several holding times (see Figure 9a). A first approximation of Mγ BT was obtained through Equation (7) and then used to run the full-field simulations, where a large number of grains needs to be considered in order to be as representative as possible. Consequently, a highly performant Lagrangian model called ToRealMotion (TRM) [74] was preferred to the LS methodology in this identification part to improve the "Number o f grains/calculation time" ratio. In the case of the AD730 alloy, 180,000 initial grains were simulated, with an arithmetic mean diameter D = 9.2 µm and a domain size of 4.3 mm × 4.3 mm. Then, a L 2 error minimization was made between the numerical results and experimental solvus data to obtain the fitted Mγ value. Finally, the Q m and M 0 parameters were refined using the new Mγ values and the Arrhenius law, leading to M 0 = 2.9 · 10 25 m 4 /J· s, Q m = 9.8 · 10 5 J/mol. Supersolvus GG simulations were again performed with the identified GB mobility parameters, then compared with the experimental data of the AD730 material. Figure 11 presents the results of supersolvus GG simulations at supersolvus temperatures for the AD730 material. The results are in good agreement with the experimental points. We can conclude, therefore, that the approach used to identify the properties gave a good approximation of M 0 and Q m values at the supersolvus domain, and they will be used in the subsolvus simulations.
The René 65 alloy has similar behavior compared to the AD730 alloy concerning the γ area fraction evolution, due to the similarity of their chemical compositions and their respective T solvus . It will also be considered that their properties, such as the GB energy (γ), the pre-exponential factor (M 0 ) and the activation energy (Q m ), are similar for both materials. Given such similarity, the full-field simulations were performed only for the AD730 material, assuming that the simulation results will be very close for the René 65 superalloy (as confirmed in Figure 11). Once the material parameters have been validated, we can proceed with the full-field GG simulations taking into account the identified parameters in this section and the γ phase dissolution described in Section 3.2.

GG with Primary γ Precipitate Dissolution
A pre-simulation was made in order to generate an initial numerical microstructure similar to the one observed experimentally (were the SPP are mainly located at the GB); thus, a short heat treatment of 400 s ≈ 6 min at 1070 • C was simulated starting from an initial smaller mean grain size than the one measured in order to get the boundaries stopping at the SPP.
The resulting initial state (after 400 s) consists of about 13,000 grains generated with an LS Laguerre-Voronoi tessellation [75]. The arithmetic mean grain radius is equal to D = 9.4 µm, and an initial polydisperse and spherical precipitate population was generated by using the experimental γ precipitate size distribution with an initial area fraction f spp = 12.5% and a mean size D spp = 2.1 µm (about 27,000 SPP were generated); all SPP are considered to be incoherent.
The size of the domain is 1 mm × 1 mm, and a time step of 5 s was considered to describe the precipitate/GB interaction correctly. The simulation parameters are M 0 = 2.9 · 10 25 m 4 /J · s, Q m = 9.8 · 10 5 J/mol as calculated in Section 5.1.1, γ = 0.6 J/m 2 [53], and T solvus = 1108 • C, as determined in Section 3.1.
Different isothermal heat treatments of 4 h were considered from T = 1060 • C (i.e., T solvus − 50 • C) up to the solvus temperature. The precipitate interface velocity v spp is calculated using Equation (5).
The evolution of the primary γ precipitate area fraction for the AD730 alloy is given by Equations (1) and (4). Figure 12 presents the results obtained from the LS simulations. Regarding the precipitate area fraction evolution (Figure 12a), the full-field simulations reproduce, as expected, the phenomenological model. Most importantly, these results agree with the experimental grain size evolution, even if the mobility was extrapolated from the supersolvus domain to the subsolvus domain. Figure 12b illustrates the mean precipitate size evolution, for the numerical treatment at 1060 • C, where the primary γ phase fraction does not change, neither does the size of the precipitates; for temperatures between 1070 • C and 1090 • C, the size of the precipitate decreases during the first minutes due to the dissolution mechanism, but then the average size increases as the smallest precipitates disappear and eventually reaches a stable value. Finally, for the temperatures 1100 • C and 1110 • C, the dissolution is more significant and the mean size of the precipitates decreases until the γ phase fraction stops evolving, but for T solvus = 1108 • C, the precipitates were mostly dissolved after 1 h of annealing. As discussed in the experimental Section 3.3, the Ostwald ripening mechanism is likely to occur after 1 h when the γ phase fraction remains constant. However, this phenomenon is not taken into account in the mean-field model used to describe the evolution of the precipitates and thus also not in the full-field simulations. The evolution of the numbers of precipitates and grains along the simulation is shown in Figure 12c,d, respectively.
Finally, Figure 12e illustrates the arithmetic mean grain size evolution for all tests. When the heat treatment is far from the T solvus (i.e., 1060 • C, 1070 • C), the SPP completely blocks the grains boundaries. However, for temperatures between 1080-1100 • C where the precipitates dissolve, the grains can grow until they find their new equilibrium state, which is reached when precipitate dissolution ends (if the Ostwald ripening mechanism was taken into account, the grains would be expected to keep evolving slowly, as the increase of mean precipitate radius at constant fraction reduces the pinning pressure and releases the GB). For the HT at the solvus temperature (1110 • C), the precipitates dissolve rapidly and the grains can grow freely when no precipitate remains. These tests show that the used numerical framework can consider the evolution of SPP and predict the grain size not only at long annealing times when the SPP have reached their equilibrium fraction but also in the transient regime. These results validate the possibility of reproducing the precipitate dissolution with a phenomenological model used in LS full-field simulations to predict the material behavior with a good compromise between accuracy and simplicity.

Parameter Identification and Validation
Since the N19 is a PM superalloy, the identification procedure was slightly different. This type of alloy has indeed a specific GG evolution. At supersolvus temperatures, the mean grain size evolves up to a maximum size (about 20 µm) and then remains constant. Thus, the identification process was only carried out during the first 20 min of heat treatments, where GG was observed to occur. The supersolvus treatments were made at 1160 • C and 1180 • C, once the first approximation of Mγ BT was obtained through Equation (7). The first sets of TRM simulations were run with the precalculated parameters, and 40,000 initial grains were considered, with arithmetic mean diameter D = 3 µm and a domain size of 0.5 mm × 0.5 mm. Then, a L 2 error minimization was made using the numerical results and the experimental supersolvus data to obtain a fitted Mγ value. Finally, the Arrhenius law was used to calculate the parameter values of the N19 material, which are M 0 = 0.51 m 4 /J · s and Q m = 3.3 · 10 5 J/mol. Figure 13 presents the experimental data and the supersolvus simulations (dashed lines) obtained with the identified properties. These simulations only match the experimental points at the HT onset. These results show a lack in the material description that we need to consider to reproduce and predict the real microstructure evolution. The reason for the grain size stagnation is likely to be that Prior Particle Boundary (PPB) particles block the grains boundaries as often observed in PM materials. Then, to correctly predict GG, it is necessary to consider this additional type of particle population in the simulations. •

Simulation of Prior Particles Boundary (PPB) particles
Such kind of additional precipitate population is common in PM materials, and they have been largely studied experimentally [76][77][78]. The PPB form at the surface of the initial powder grains, and they usually consist of oxides, carbides, or oxy-carbides. It is commonly accepted that PPB particles result from segregations or surface contamination. A slight dissolution [79] and coarsening [80,81] of these precipitates have been observed in materials submitted to high temperatures. Some numerical works have been developed to account for their presence, for instance, using a particle packing in [82] and an FE framework in [83].
Here, the PPB particles will be treated as a population of oxides that will not evolve; thus, they will be considered as static particles located in the interfaces of the initial powder grains. Therefore, to generate the PPB network, a Laguerre-Voronoi tessellation was used considering the experimental size distribution of the initial powders (mean D powder = 25 µm), followed by the generation of the PPB particles at the boundaries of the tessellation (see Figure 14). The PPB particles are considered as a new distance function in our LS approach. Once the PPB particles were generated, the actual initial polycrystal was statistically generated (several grains/crystals per powder grain), and, finally, the GG simulation was performed. The validation of the identified properties for the N19 material was made once again with the new considered topology with f PPB = 2% and a homogeneous precipitate size D PPB = 1 µm. Although these values are large concerning the actual characteristics of the PPB [79,80], they enable, in Smith-Zener meaning (i.e., by respecting the ratio f PPB /D PPB while increasing by around a factor 10 the respective values of D PPB and f PPB ), to respect a representative pinning pressure resulting from this PPB particles population. Indeed, in terms of numerical cost of the meshing adaptation and subsequent simulations, it will not be realistic to consider a large number of nanoparticles.
The results are presented in Figure 13 (solid lines). The mean grain size prediction is in very good agreement with the experimental data, where the PPB particles block the grain boundaries. The predicted final mean grain size (about 20 µm) corresponds approximately to the mean size of the initial powders as a result of the strong pinning pressure exerted by the small PPB particles.
Once the material parameter values have been established, the full-field GG simulations can be performed, taking into account the identified GB migration parameters, the γ precipitate dissolution described in Section 3.2, and by considering or not PPB particles.

GG with Primary γ Precipitate Dissolution
A pre-simulation was also made as for the AD730, but for a HT of 100 s at 1110 • C. The initial state (after 100 s) consists of around 33,000 grains generated with an LS Laguerre-Voronoi tessellation [75]. The arithmetic mean grain diameter D = 3 µm and an initial spherical primary γ precipitate population with a area fraction f spp = 18.5% and a mean size D spp = 1.3 µm (about 32,000 incoherent SPP) are considered.
The size of the domain is 0.5 mm × 0.5 mm (see Figure 15 but without the red particles), and the considered time step was 5 s. The simulation parameters were M 0 = 0.51 m 4 /J · s, Q m = 3.3 · 10 5 J/mol as calculated in Section 5.2.1, γ = 0.6 J/m 2 [53] and T solvus = 1150 • C as determined in Section 2.1.1.
Simulations considering or not the PPB particles, in addition to the evolving γ precipitates, were carried out to understand how they affect the GG behavior.

•
Simulations without considering the PPB particles Different isothermal heat treatments (4 h) around T solvus were simulated. The primary γ precipitate interface velocity was calculated according to Equation (5).
The primary γ phase area fraction evolution for the N19 was obtained by Equations (2) and (4). Figure 16a shows that the γ phase area fraction evolution matches the phenomenological model and the experimental data. The mean primary γ precipitate size evolution is presented in Figure 16b where there is a slight decrease at the beginning of the simulation; then, as the small precipitates continue to dissolve, the mean grain size tends to increase. Figure 16c,d illustrate the number evolution of the γ precipitates and grains, respectively. The GB kinetics is so fast that it makes the smallest grains disappear during the first minutes of the simulations.
Finally, Figure 16e shows the arithmetic mean grain size evolution for all the simulations. As discussed before, the mean grain size evolves during the γ precipitate dissolution. However, the grain size reaches a stable value once the dissolution stops. As expected for this case, a grain size larger than the one observed experimentally is reached. Indeed, for the HT at 1150 • C, the predicted grain size continues to increase contrary to the experimental points that stagnate at about 20 µm.

•
Simulations considering the PPB particles The PPB particles were generated at the Laguerre-Voronoi tessellation boundaries based on the experimental granulometry of the initial powder of this material (see Section 3.2). The initial state is the same used previously for the N19 (see Figure 15). However, for this case, an additional static particle population was introduced ( f PPB = 2% of size D PPB = 1 µm, red particles). Figure 17 illustrates the microstructure evolution for some HT at different times and temperatures along the simulation and the results are summarized in Figure 18.  . .
. : : : : : : .          The γ area fraction evolution (Figure 18a) is similar to one of ( Figure 16a) because the primary γ interface velocity is the same. Considering the static PPB particles, the predicted arithmetic mean grain size is in better agreement with the experimental data as illustrated in Figure 18c. The LS simulations describe very well the limiting grain size independently of the temperature. Figure 19 illustrates the microstructural evolution for both the considered cases (with or without PPB particles) at 1150 • C, and allows us to visualize how the grains continue to grow when no γ precipitate remains and PPB particles are not taken into account. On the contrary, when static PPB particles are simulated, the grains' boundaries are blocked by their network, limiting the mean grain size slightly below the initial powder size of the material. Figure 19. Comparison of the simulated microstructure evolution of the N19 superalloy considering or not the PPB particles (red) at 1150 • C for a zoomed area of the 0.5 mm × 0.5 mm simulation domain.

Conclusions
A new approach using a mean-field model to describe the dissolution of second phase particles was introduced in a finite element framework using the Level-Set model to simulate grain growth for different γ-γ nickel base superalloys as primary γ precipitates dissolve.
The dissolution kinetics of primary γ precipitates was studied in the AD730, René 65 and N19 alloys. Thermo-Calc curves and a JMAK approximation were used to establish an equation for describing the primary γ phase fraction evolution as a function of temperature and time, which is valid from T solvus to 50 • C below for each alloy (T solvus − 50 ≤ T ≤ T solvus ). It was shown that the particle evolution for the two conventional cast and wrought alloys (AD730, René 65) was very similar. Therefore, the same equation was used to reproduce the primary γ dissolution in both alloys. For the powder metallurgy alloy (N19), the dissolution kinetics was faster.
The LS simulation results respect the imposed phenomenological equation and so there is a good agreement with the experimental data concerning the second phase precipitate evolution. These results confirm the possibility of considering a simple phenomenological model to describe the precipitate evolution to get good accuracy in terms of the predicted precipitate and grain sizes, including both the transient regime at short times and the steady-state at longer annealing times.
Moreover, the possibility of considering the Prior Particles Boundaries (PPB) particles has been introduced in the simulations of the powder metallurgy N19 alloy. The numerical results are consistent with the observations, where these PPB particles limit the mean grain size to about the mean size of the initial alloy powder.
The current perspectives of this work concern the use, in 2D and 3D, of the proposed methodology to discuss other mechanisms, such as the Ostwald ripening and precipitation.

Data Availability Statement:
The experimental data necessary to reproduce these findings are available from the corresponding author on request.