Simulation of Diffusion Bonding of Different Heat Resistant Nickel-Base Alloys

: Currently, an important fundamental problem of practical importance is the production of high-quality solid-phase compounds of various metals. This paper presents a theoretical model that allows one to study the diffusion process in nickel-base refractory alloys. As an example, a two-dimensional model of ternary alloy is considered to model diffusion bonding of the alloys with different compositions. The main idea is to divide the alloy components into three groups: (i) the base element Ni, (ii) the intermetallic forming elements Al and Ti and (iii) the alloying elements. This approach allows one to consider multi-component alloys as ternary alloys, which greatly simpliﬁes the analysis. The calculations are carried out within the framework of the hard sphere model when describing interatomic interactions by pair potentials. The energy of any conﬁguration of a given system is written in terms of order parameters and ordering energies. A vacancy diffusion model is described, which takes into account the gain/loss of potential energy due to a vacancy jump and temperature. Diffusion bonding of two dissimilar refractory alloys is modeled. The concentration proﬁles of the components and order parameters are analyzed at different times. The results obtained indicate that the ternary alloy model is efﬁcient in modeling the diffusion bonding of dissimilar Ni-base refractory alloys.


Introduction
The theory and practice of diffusion welding of dissimilar metals has been actively developed over the past decades due to the undoubted importance of this process for various technologies [1,2]. Some alloys can be successfully bonded without the use of an interlayer [3][4][5]. However, in most cases, to avoid the formation of brittle intermetallic inclusions, an interlayer metal is used. Examples are bonding of titanium to steel [6][7][8][9], nickel [10], copper [11] and the use of copper interlayer for bonding stainless steel to Zircaloy 4 [12]. Heat-resistant intermetallic alloys can have unique physical and mechanical properties (heat resistance, resistance to oxidation, corrosion, creep performance, etc.) that makes them ideal for a number of applications in aerospace industry [13][14][15][16][17]. Structure and properties of Ni-base alloys can be improved by laser treatment [18,19] and thermomechanical treatment [20,21]. Electro-discharge machining process of Inconel-718 superalloy has been optimized in [22].
Interest in these alloys is due to the possibility of achieving higher alloy characteristics. Intermetallic alloys based on Ni 3 Al are used in parts of combustion chambers of aircraft engines with operating temperatures up to 1300 • C. Methods for joining metals and superalloys include diffusion welding [15,23], fusion welding, friction welding [24][25][26][27][28], diffusion welding in a liquid phase [29], friction stir welding [28], etc. When dissimilar metals are joined together in solid phases, undesirable brittle intermetallic phases are often formed, which reduce the mechanical properties of parts [3,6,7,12,30].
One of the ways to predict the properties of intermetallic phases is computer modeling of the structure of alloys, which allows to reveal undesirable structural states for a given alloy prior to experiment. Various computational methods can be used, for example, molecular dynamics (MD) [31], phase-field [32], finite element [33], Monte Carlo [34][35][36][37] and others methods. Diffusion processes have also been modeled in frame of two-dimensional continuum model [38]. It is difficult to apply MD to simulate the diffusion process, since the time scale achieved by this method is on the order of 1 ns, and diffusion bonding takes about 15 min. The time difference is huge, twelve orders of magnitude. This is why the development of the Monte Carlo method is very important. The use of continuum models for diffusion in multi-component media is not straightforward, since it requires knowledge of the matrix of diffusion coefficients [39]. Various simulation methods work at different scales and they make it possible to establish ways to achieve the required alloy structure, to trace the formation of intermetallic inclusions and their distribution over the volume and to study the effect of alloying and temperature changes. Modeling the ordering process of alloys makes it possible to establish an equilibrium structure and calculate its parameters, such as energy, degree of short-range and long-range order, depending on the concentration of components.
The main difficulty in atomistic simulation of heat-resistant Ni-base alloys is that they include many alloying elements. This prevents the use of MD modeling due to the complexity of constructing interatomic potentials. It is necessary to adopt additional simplifying assumptions to create models describing diffusion processes.
In this work, within the framework of a two-dimensional model, the structure of heat-resistant VKNA-25 and EP975 alloys is investigated. These alloys are based on the Ni 3 Al intermetallic compound containing alloying elements with different concentrations. Alloy components are divided into three groups depending on the role they play in the alloy. The Monte Carlo algorithm is used to simulate the diffusion process by the vacancy mechanism in a three-component metallic systems.
In Section 2, we describe the experimental results on diffusion bonding of VKNA-25 and EP975 alloys, present the derivation of order parameters for ternary alloys following [37], describe the three-component alloy model and the vacancy diffusion model. Then, in Section 3, the numerical results and their discussion are presented. Section 4 concludes this work.

Experimental Results
The chemical composition of VKNA-25 and EP975 alloys is shown in Table 1 in weight percent [40]. Both alloys are based on nickel, while aluminum and titanium are elements that form the intermetallic phase Ni 3 Al, where aluminum atoms can be replaced by titanium atoms. Many alloying elements are also present. Figure 1 shows the microstructures of (a) VKNA-25 and (b) EP975 alloys. In (c), the microstructure of the solid phase joint of the two alloys is presented. In the images, the γ' phase, which is the Ni 3 Al intermetallic compound, has a dark color, and the solid solution of alloying elements in Ni 3 Al is the γ phase, which has a light color. The microstructures of the alloys and the joint zone were studied by scanning electron microscope Mira 3LMH (TESCAN, Brno , Czech Republic). EP975 alloy is a wrought polycrystalline alloy containing 55% of the hardening γ' phase, and VKNA-25 alloy is a cast intermetallic alloy with a single-crystal structure based on Ni 3 Al, in which the content of the intermetallic γ' phase is 85-90% [15,16]. Cylindrical samples of VKNA-25 and EP975 alloys were 20 mm in diameter and 30 mm in height. Pressure welding of VKNA-25 and EP975 alloys was performed in vacuum and welding temperature was 1125 • C. The two cylindrical samples were put in contact and compressed with the strain rate of˙ = 10 −4 s −1 . Note that VKNA-25 is a hard-to-deform single-crystal alloy, while EP975 is polycrystalline alloy that flows superplastically at the selected temperature and strain rate. The superplasticity of EP975 alloy is important because it helps to achieve a better physical contact during pressure welding and to minimize the formation of voids. Strain control loading is performed so that pressure is not a controlled parameter. A compressive strain is applied normal to the interface to achieve better physical contact between the parts to be welded.
The microstructure of the joint is presented in Figure 1c, where VKNA-25 and EP975 alloys are above and below the contact plane, respectively. A high-quality joint is obtained with no voids or other imperfections. Close to the interface, the γ phase in VKNA-25 alloy has lighter color than in the bulk. This indicates a change in the composition of the alloys in the vicinity of the contact plane. This is confirmed by the components concentration profiles in the vicinity of the contact plane presented in Figure 2a. The distribution of elements in the solid-phase joint zone was determined by the X-ray energy-dispersive spectroscopy (EDS) on a Tescan VEGA 3SBH scanning electron microscope with an EDX attachment. The measurements were carried out only from the γ-phase, moving from the EP975 alloy to the VKNA-25 alloy with a step of 10 µm. The diffusion zone was found to be about 35 µm wide. To plot the distribution of alloying elements, the EDS data were averaged over at least three measurements per point. In Figure 2b, the concentrations of three groups of alloy components are shown: the base element Ni, the intermetallic-forming Al and Ti are combined and named Al and all other elements are combined and named Cr. The merit of combining alloy elements into these three groups is explained in Section 2.4.

Order Parameters for Ternary Alloy
Multi-component alloys can contain domains of ordered and disordered phases, as well as segregation of pure components. It is important to quantify the degree of ordering and the content of segregations. In this section, the order parameters for ternary alloys derived in [37] are reproduced for the sake of the readers. In addition, a new result is presented-the derivation of the parameters that determine the proximity of the alloy to decomposition into pure components.
All the expressions presented below are valid for any crystal lattice of any dimension. The information about the type and dimensionality of the lattice is introduced through the coordination numbers N i , i.e., the numbers of atoms on the i-th coordination sphere.
Let us consider a three-component alloy of composition A m B n C k , in which atoms of types A, B and C are located at lattice sites. The lattice has N i atoms on the i-th coordination sphere. The concentrations of atoms of types A, B and C are equal to, respectively, We denote by p KL the probability that the i-th coordination sphere of an atom of the type K contains an atom of the type L, where K, L = {A, B, C}. In a ternary alloy, the following relationships exist between the nine probabilities p Let us denote by ϕ KL the binding energy of a pair of atoms of types K and L located at a distance equal to the radius of the i-th coordination sphere.
The potential energy of the structure per atom, taking into account the interaction of atoms in the first I coordination spheres, can be written in the form The energy of the completely disordered state of the structure is determined by Equation (6) for The energy of decomposition into pure components is determined by Equation (6) for p Let us choose the E disord energy as a reference point, and characterize the energy of any structure by the difference where we have introduced the order parameters and ordering energies When deriving Equation (9), Equations (2)-(5) were used.
On the other hand, if the E decomp energy is used as a reference point, then the energy difference is where a set of parameters defining the proximity of the alloy to decomposition into pure components is introduced When deriving Equation (12), we used Equations (2)-(5).
As we can see from Equation (9), the deviation of the energy of any structure from the energy of a disordered alloy is uniquely determined by the coordination numbers N i , ordering energies ω For a two-component alloy (c C = 0), the above equations are greatly simplified. For example, Equations (2)-(5) reduce to and Equation (9) assumes the form On the other hand, the complexity of the relationships between the probabilities p (i) KL increases very rapidly with an increase in the number of alloy components. For this reason, in the following section we propose a ternary alloy model that can be used to describe diffusion processes in multi-component refractory nickel-base alloys.

Vacancy Diffusion Model
Let us describe the mathematical model of the diffusion process in the alloy according to the vacancy mechanism in the model of hard spheres, which can be applied to a multi-component alloy of the composition, specified on a lattice of any type and any dimension. In this model, it is assumed that an elementary act of diffusion is the transition of one of the atoms surrounding a vacancy to its place. It is assumed that any atom from the first K coordination spheres can jump into the vacant lattice site. The number of such atoms is M = ∑ K i=1 N i , where N i is the coordination number for the i-th coordination sphere. Each of the M atoms is assigned the probability p m , m = 1, ..., M, to occupy the place of the vacancy in an elementary act of diffusion, so that ∑ M m=1 p m = 1. For this purpose, the change in the energy of the alloy, ∆E m , associated with the transition of the m-th atom to a vacant place for a given temperature of the alloy T is calculated. The required probabilities are determined as follows where and k = 8.61733 × 10 −5 eV/K is the Boltzmann constant. Note that for increasing temperature the probabilities p m will tend to become equal, regardless the values of the energies ∆E m , and this will lead to the disordered state. On the other hand, for relatively small temperatures, the probabilities p m will depend on ∆E m , leading to the ordering process. Pressure is not taken into account in this model, because we use the rigid lattice assumption without taking into account the effect of atomic relaxation.

Ternary Alloy Model and Simulation Setup
For modeling multi-component refractory alloys, we divide all elements into three groups. The first group includes a single element, nickel, which is the base element. The second group combines the intermetallic-forming elements aluminum and titanium. Finally, the third group includes all alloying elements. We will call these groups according to their main elements, that is, the first group is A = Ni, the second is B = Al and the third one is C = Cr. It should be noted that the ternary alloy model can describe the considered multi-component superalloys only at a qualitative level.
In Table 2, we present the content of the three components in EP975 and VKNA-25 alloys in weight percent and atomic percent. In our experiments and in our model, we do not have solidification because joining takes place in the solid state without melting. The high temperature activates the diffusion and superplastic flow of the EP975 alloy, and these two factors play a positive role in achieving a high-quality joint.
The main assumptions made in the model are the approximation of pair interatomic interactions and the approximation of a rigid crystal lattice, i.e., the assumption that the atoms occupy the sites of undeformed crystal lattice, or, in other words, the effects of atomic relaxation are not taken into account.
Our simulation will be carried out with the use of the toy model. We consider two-dimensional square lattice with the interatomic distance equal to unity, which can be always achieved by a proper choice of the unit of length. The atoms in the computational cell are numbered by the indices 1 ≤ n x ≤ N x and 1 ≤ n y ≤ N y in the horizontal and vertical directions, respectively. The size of the computational cell is N x × N y = 200 × 400 and periodic boundary conditions are applied. The typical interatomic distance in metals is about 2.5 Å, then the size of the computational cell is about 0.1 µm. This is two orders of magnitude smaller than the experimentally observed diffusion zone of 35 µm, see Figure 2. An increase in the size of the computational cell requires an increase in the simulation time.
To describe interatomic interactions, for simplicity, the Morse pair potential is chosen where D KL , θ KL and R KL are the parameters of the potential for the atoms of sorts K and L, and r (i) is the radius of the i-th coordination sphere. The following values of the potential parameters are chosen: As it can be seen from Equation (19), for simplicity, we take equal values of the equilibrium distance R KL , which reflects the fact that we do not take into account atomic relaxation effects. The equilibrium distance is slightly (5%) greater than the interatomic distance, which is always so when long-range interactions are considered. In our simulations, three coordination shells are taken into account. We also take equal values for the parameters θ KL , that define the rigidity of bonds, and the reason is the same, the lattice is assumed to be non-deformable. The depth of the potentials (binding energies) are chosen in a way to obtain the desired structure of the alloys. In particular, we take D NiAl larger than D NiNi and D AlAl to get negative ordering energy ω NiAl = ϕ NiNi + ϕ AlAl − 2ϕ AB , so that the formation of the intermetallic compound Ni 3 Al is energetically favorable. The Cr-Cr potential is the shallowest one in order to prevent segregation of the alloying elements because we want them to be dissolved in the intermetallic matrix.

Results and Discussion
Firstly, we create in the computational cell the completely disordered state by assigning the sorts of the atoms Ni, Al or Cr randomly with the probability equal to their concentrations. We do this for EP975 and VKNA-25 alloys taking the concentrations of the elements from Table 2 (in at.%). Then one vacancy is created in the computational cell and the vacancy diffusion algorithm described in Section 2.3 is used to model the equilibration of the alloy structure. The simulation is conducted at the temperature T = 1200 • C until the total energy of the alloy reaches the minimal value. The minimal energy value was achieved in 10 10 jumps of the vacancy. This way, the equilibrium structures of EP975 and VKNA-25 alloys were obtained.
Secondly, in the half of the computational cell, 1 ≤ n y ≤ N y /2 equilibrium EP975 alloy is placed, while in the second half, (N y /2) + 1 ≤ n y ≤ N y , we place the equilibrium VKNA-25 alloy. Again, single vacancy is placed in the computational cell and the interdiffusion of the components between the two alloys is simulated at the temperature T = 1200 • C. The duration of this simulation is 3 × 10 11 vacancy jumps.
In Figure 3, we present the evolution of the structure of the alloys during diffusion bonding. In (a), the initial structure at t = 0 is shown. In the lower (upper) half of the computational cell, EP975 (VKNA-25) alloy is found. The inset shows the atomic resolution of the structure. Dark blue and light blue colors are used for Ni and Al atoms, respectively. They form ordered structure with Ni 3 Al stoichiometry. Alloying elements (Cr) appear in yellow color. They are dissolved in the intermetallic matrix forming small clusters. In Figure 3b,c structure is shown at t = 4.5 × 10 10 and 3.0 × 10 11 (note that the number of vacancy jumps is the measure of time in our simulations). Gradual redistribution of the alloy elements can be observed.
Time evolution of the concentration profiles is presented in Figure 4a-c for Ni, Al and Cr, respectively, in at.%. Black curves 1 stand for the initial distribution of components, while red curves 2 and blue curves 3 are for t = 4.5 × 10 10 and 3.0 × 10 11 , respectively. Since periodic boundary conditions are used, we have periodic structure with alternating layers of EP975 and VKNA-25 alloys. It can be seen that Ni and Al move from VKNA-25 into EP975, while Cr in the opposite direction. As a result, the distribution of elements becomes homogeneous. This behavior is in qualitative agreement with the experimental results presented in Figure 2. Order parameters at the first coordination shell, calculated with the help of Equation (10), are shown as the functions of time in Figure 5. Only α (1) NiCr shows a noticeable change in time, while other two order parameters remain practically constant. This is explained by our choice of the interatomic potential parameters, see Equation (19). Parameter D NiCr is smaller than D NiAl and D AlCr . That is why the reconstruction of Ni-Cr bonds takes place with a higher probability during the diffusion process. The absolute value of the order parameter α (1) NiAl is considerably larger than that of α (1) NiCr and α (1) AlCr , since Ni and Al form an ordered phase Ni 3 Al.
There is a relation between simulation time and temperature as regards the kinetics of welding and order parameter evolution. With an increase in temperature, the probabilities in Equation (16) become closer, the ordering process becomes faster but the maximal values of order parameters become smaller. In addition, with an increase in temperature, the rate of vacancy jumps increases.

Conclusions
A theoretical basis and its numerical implementation were developed for modeling diffusion processes in refractory nickel-base alloys. As an example, diffusion bonding of VKNA-25 and EP975 alloys is considered. All alloy components were classified into three groups as explained in Section 2.4, which made it possible to consider multi-component alloys as ternary.
In Section 2.2, expression for the alloy energy, Equation (6), was given under the assumption of non-deformable lattice and pairwise interatomic interactions. The difference between alloy energies in any state and in the disordered state, ∆E 1 , is expressed in terms of the order parameters α KL [see Equation (9)]. On the other hand, the difference between alloy energies in any state and in the decomposed into pure components state, ∆E 2 , is expressed in terms of the parameters β (i) KL and ordering energies ω (i) KL (see Equation (12)). This means that the parameters β (i) KL determine the closeness of the alloy to decomposition into pure components.
In Section 2.3, the vacancy diffusion model was described. The model takes into account the effect of temperature and the energy gain/loss due to the vacancy jump. A vacancy jump leading to an increase in the potential energy of the alloy is possible, but its probability is less than the probability of a jump leading to a decrease in the potential energy. As the temperature rises, the probability that atoms of different sorts will jump into a vacancy is getting closer. This leads to a transition to a disordered state at a sufficiently high temperature.
It is important to note that all expressions presented in Sections 2.2 and 2.3 are valid for ternary alloys based on any lattice of any dimension. Information on the type and dimension of the lattice is taken into account by the coordination numbers N i . Overall, a model of ternary alloys was developed that allows to simulate diffusion processes in heat-resistant Ni-base alloys.
In this work, a toy model of refractory alloys based on a two-dimensional square lattice was used. In an upcoming study, the developed model will be applied to the simulation of diffusion bonding of refractory alloys using a more realistic three-dimensional fcc lattice.