Introduction

Precise determination of structural and chemical properties of metallic systems is strenuous, and sometimes, even impossible, due to the technical limitations of the experimental set-up. However, a database limited to one property or even an incomplete experimental dataset can usually be sufficient to perform molecular dynamics simulations and/or use CALPHAD-type models enabling to quickly cover broad compositional dependence. Recent literature reports gave proof of MD simulations being sufficient, and then, such simulations have been extensively used for studies of structure, dynamics and thermodynamics of liquid and solid bulk materials [1,2,3,4,5,6]. Sometimes MD simulations are insufficient to obtain information that can be straightforwardly retrieved from CALPHAD calculations. This is a result of CALPHAD-type models’ ability to give a first insight into the structure and atomic transport properties of a melt without any knowledge of semi-empirical potentials [7,8,9,10].

At least four intermetallic phases (θ-Al2Cu, η-AlCu, γ-Al4Cu9, β-b.c.c.) exist in solid state of the Al–Cu system [11], which makes it interesting according to experimental and computational contributions given in the literature [12,13,14]. In contrast, only a few reports are related to the liquid phase. Rybicki et al. [15] were the first to experimentally study the structure of the liquid Al83Cu17 alloy, observing how the atomic groups are linked to each other. Later, this melt was the subject of computational investigations, including ab initio molecular dynamics (AIMD) [16] and classical molecular dynamics studies [5], employing MEAM potentials [17]. Consequently, both studies revealed short-range order presence in the liquid state of the Al80Cu20 alloy, which is responsible for anomalous structural dynamics. Similarly to the Al80Cu20 melt investigations, analogous computational techniques have been used for liquid structure study at an another Al–Cu alloy composition, i.e. Al60Cu40, including AIMD and CMD simulations performed by Wang et al. [18] and Dziedzic et al. [19], respectively.

Thermophysical properties of liquid Al–Cu alloys also caught scientific interest [19,20,21,22], including experimental and modelling contributions, which altogether build a solid background and open up a field for discussion of the new data presented in this paper. Schick et al. [20] experimentally determined viscosity for a series of liquid Al–Cu alloy compositions, excluding only the liquid Al60Cu40 alloy, using a high-temperature oscillating-cup viscometer method and supplementing it with thermodynamics-based modelling analyses. As a result, a pronounced maximum in viscosity was observed for Al–Cu alloy composition close to intermetallic phases existing in the solid state of the Al–Cu system, which could suggest “associates” influence. This issue has not been discussed in the literature with regard to structure and diffusivity investigations. The latter ones include the determination of the self-diffusion coefficient of Cu and of the mutual diffusion coefficient using neutron scattering [21] or X-ray radiography [22] for liquid Al-rich alloys, namely Al90Cu10, Al83Cu17, Al80Cu20 and Al75Cu25.

The “associate forming tendency” has been partially studied by Trybula et al. [23], using the free volume model (FVM), but this investigation provided only fragmentary insights into the structure and thermodynamics of Al–Cu liquid alloys. Unfortunately, those results were not enough to fully describe the physical phenomena occurring in the liquid state of the Al–Cu alloys due to the FVM constraints [24]. Apart from that, no detailed study on the topology of local atom environment and its relationship with chemistry and atom kinetics has been given yet, which could play a vital role in the investigation of glass or crystal state stability. As it has been shown in a recent study on liquid Al–Zn alloys [4], a combination of MD simulations with thermodynamics-based modelling can successfully account for experimental data and describe the relationship between the structure and thermophysical properties. One of the widely used thermodynamics-based models is the quasi-chemical approximation [7], which can successfully reproduce experimental data for various liquid metallic alloys [25, 26].

In this work, we propose to combine molecular dynamics simulations with thermodynamics-based models to study the structure and chemistry of liquid Al–Cu alloys. Furthermore, the results presented herein extend the applicability of the new modified embedded atom model (MEAM) parametrization established by Trybula [5] to Al–Cu alloys with higher Cu content. In this concept, we would like to initially compute properties that are available in the literature to check the validity of the MEAM potential. We then proceed to investigate the energetic stability of atomic clusters and their influence on kinetics and thermodynamics of Al–Cu alloy liquid state. Structural characterization obtained by performing MD simulations concerns the analysis of local chemical surroundings of atoms, including the computation of a structure factor and pair distribution functions, which is then supplemented by 3D imaging of local atom topology. In addition, we search for a structure–dynamics–thermodynamics relationship by considering the universal scaling law [27] and extended Enskog collision frequency theory [28]. Finally, the main scope of the present work is building a bridge between MD-based simulations and CALPHAD-type modelling based on quasi-chemical approximation for a weakly compound-forming molten alloy [29] and tests it on liquid Al–Cu alloys that will be helpful in discussion of structure and chemistry of the Al–Cu system liquid state.

Computational methods

Quasi-chemical approximation for weak compound-forming melt

Within the compound formation model (CFM), the presence of an appropriate chemical compound of \( A_{m} B_{n} \) type (m and n are integer numbers) is considered for a binary system, A–B. \( A_{m} B_{n} \) type corresponds to one of the intermetallic phases occurring in the solid state, or it is just a hypothetical one. The complete mathematical formalism of the CFM approach can be found elsewhere [7, 30].

The principal thermodynamic quantity is Gibbs free energy of mixing, Gmix, which has the form:

$$ G_{{\text{mix}}} = k_{{\rm B}} T\left[ {\left( {\frac{W}{{k_{{\rm B}} T}}} \right) \cdot C(1 - C) + \left( {\frac{{\Delta W_{\text{AB}} }}{{k_{{\rm B}} T}}} \right) \cdot \varPhi_{AB} (C) + \left( {\frac{{\Delta W_{\text{AA}} }}{{k_{{\rm B}} T}}} \right) \cdot \varPhi_{\text{AA}} (C) + \left( {\frac{{\Delta W_{\text{BB}} }}{{k_{{\rm B}} T}}} \right) \cdot \varPhi_{\text{BB}} (C)} \right]. $$
(1)

Such form of Gmix is obtained when standard thermodynamic relations are applied, and some algebraic rearrangements are made according to a receipt previously used by Novakovic et al. [10, 26] or Trybula et al. [9]. \( \varPhi_{ij} \,(i,j\,\, = \,A,B) \) are concentration functions depending on m and n values, C is concentration of one component, k B is the Boltzmann constant, and T is absolute temperature. \( W,\Delta W_{\text{AB}} ,\Delta W_{\text{AA}} ,\Delta W_{\text{BB}} \) are energy-pair parameters optimized by utilizing the least-squares method and fitting the modelled Gmix to experimental one. Different expressions can be devised for the computation of the concentration functions in relevance to the assumed type of clusters existing in the liquid state. In present consideration, we examined two variants of Φ ij functions, the first corresponds to the Al2Cu compound, while the latter variant is the Al2Cu3 phase. Therefore, Φ ij functions are formulated as follows [31]:

  1. 1.

    m = 1 and n = 2 for Al2Cu compound:

    $$ \begin{aligned} \varPhi_{\text{AB}} & = \frac{1}{6}(1 - C) + (1 - C)^{2} - \frac{5}{3}(1 - C)^{3} + \frac{1}{2}(1 - C)^{4} ; \\ \varPhi_{\text{AA}} & = 0; \\ \varPhi_{\text{BB}} & = - \frac{1}{4}(1 - C) + \frac{1}{2}(1 - C)^{2} - \frac{1}{4}(1 - C)^{4} \\ \end{aligned} $$
    (2a)
  2. 2.

    m = 2 and n = 3 for Al2Cu3 clusters:

    $$ \begin{aligned} \varPhi_{\text{AB}} & = \frac{13}{420}(1 - C) + \frac{2}{3}(1 - C)^{3} - \frac{3}{2}(1 - C)^{4} + \frac{2}{3}(1 - C)^{6} - \frac{5}{7}(1 - C)^{7} + \frac{1}{4}(1 - C)^{8} ; \\ \varPhi_{\text{AA}} & = - \frac{53}{840}(1 - C) + \frac{2}{3}(1 - C)^{3} - \frac{5}{4}(1 - C)^{4} - (1 - C)^{6} + \frac{4}{7}(1 - C)^{7} - \frac{1}{8}(1 - C)^{8} ; \\ \varPhi_{\text{BB}} & = \frac{23}{280}(1 - C) - \frac{1}{2}(1 - C)^{4} + \frac{2}{5}(1 - C)^{6} + \frac{1}{7}(1 - C)^{7} - \frac{1}{8}(1 - C)^{8} . \\ \end{aligned} $$
    (2b)

The aforementioned functions can also be used for examining the chemical short-range order in a melt. The first structure index, that is directly calculated from the CFM formalism and linked to Bhatia–Thornton (B–Th) theory [32], is the concentration fluctuation function at the long wavelength limit, SCC(0). It is an indicative of chemical short-range order, and for a sub-regular metallic system it is formulated as follows:

$$ {\text{S}}_{\text{cc}} (0) = k_{{\rm B}} T\left[ {\left( {\frac{W}{{k_{{\rm B}} T}}} \right) \cdot C(1 - C) + \left( {\frac{{\Delta W_{AB} }}{{k_{{\rm B}} T}}} \right) \cdot \varPhi_{\text{AB}} + \left( {\frac{{\Delta W_{\text{AA}} }}{{k_{{\rm B}} T}}} \right) \cdot \varPhi_{\text{AA}} + \left( {\frac{{\Delta W_{\text{BB}} }}{{k_{{\rm B}} T}}} \right) \cdot \varPhi_{\text{BB}} } \right] , $$
(3)

while for ideal solution, Scc(0, id) has the form of:

$$ S_{\text{CC}} (0,{\text{id}}) = C(1 - C). $$
(4)

The second indicator of short-range order (SRO) in a molten alloy is the Warren–Cowley parameter αtot [33], quantifying chemical SRO within the first coordination shell and it is written:

$$ \alpha_{\text{tot}} = \frac{M - 1}{[1 - M \cdot (Z - 1)]} , $$
(5)

where Z is the coordination number and M is the Scc(0)/(Scc(0, id)) ratio. A melt has a compound formation tendency for which the M value is lower than unity [34].

Molecular dynamics simulations

A recent parametrization of the modified embedded atom model (MEAM), given by Trybula [5] for liquid Al80Cu20 alloys, has been used for the investigation of Al1−xCux (x = 10, 30, 40 and 50 at%) alloys properties by means of molecular dynamics simulations (MD). A consistent and compact description of liquid Al80Cu20 alloys’ properties was obtained in reference to the previous ab initio molecular dynamics [16] (AIMD) and experimental [20, 22, 35] results. Herein, we validate this MEAM parametrization for structure and transport properties investigations over a broad Cu concentration range. All MD simulations have been performed using the LAMMPS code [36] and carried out using the Verlet algorithm with a time step of 1.5 fs. Cubic box with a total number of 1372 atoms was built, and periodic boundary conditions in three directions were imposed for each Al–Cu composition investigated. Randomly oriented atoms in the simulation box were initially heated up to 2500 K and thermalized during 100 ps. Pressure, volume and temperature of each sample during quenching were monitored to reproduce experimental densities [37] at each temperature investigated (T = 1073, 1173, 1223, 1345 and 1523 K). Before structural and dynamic properties determination, each sample was thermalized during 20 ps and the final simulations were continued for 500 ps to determine structural and transport properties of liquid Al1−xCux alloys investigated (x = 10, 30, 40 and 50 at% Cu). For structure analysis, 50 atomic configurations regularly spaced in time were selected and extracted from the MD trajectories and then were used for three-dimensional topological analysis including the Voronoi tessellation method [38]. Each atomic structure was brought to minimum energy using the steepest descent algorithm as implemented in LAMMPS.

Self-diffusion coefficient, describing the random walk motion of tagged atomic species β (β = A, B for A and B being the particles of an AB binary system), was determined from the linear slope of the mean square displacement (MSD), equivalent to the Green–Kubo formalism, given as:

$$ D_{\beta } \, = \,\mathop {\lim }\limits_{t \to \infty } \frac{1}{{N_{\beta } }}\sum\limits_{i = 1}^{{N_{\beta } }} {\frac{{\left\langle {\left| {{\mathbf{r}}_{i,\beta } (t)\, - \,{\mathbf{r}}_{i,\beta } (t = 0)} \right|^{2} } \right\rangle }}{6t}} , $$
(6)

where ri,β(t) is the position of particle i of species β at time t.

While for A–B binary systems various methods for computing the mutual diffusion coefficient exist, the first one is the simple Darken’s relationship [39] linking the interdiffusion coefficient to self-diffusion coefficients, which is written as:

$$ D_{\text{AB}} = C_{B} D_{A} + C_{A} D_{B} . $$
(7)

In practice, Eq. (7) is not valid for concentrated melts with a strong associate forming tendency; hence, a modified Darken’s formula is used instead. It incorporates an additional thermodynamic factor, ϴ [39].

The second expression that can be used for mutual diffusion coefficient computation is the universal scaling law of Dzugutov [27] that is constituted by two important assumptions. The first includes the Enskog collision frequency describing binary collisions of dense hard-sphere fluids, only restricted to short-range repulsive forces, and governing the energy and momentum transfer [28]. In our investigations, we have considered the extended Enskog’s theory on a binary liquid, where reduced self-diffusion coefficient written as D* i  = D i /χi is computed using relevant diffusion scale factor of ith species in a binary liquid [40]:

$$ \chi_{1} = 4\sigma_{1}^{4} g_{11} (\sigma_{11} )\rho_{1} \sqrt {\frac{{\pi k_{{\rm B}} T}}{{m_{1} }}} + 4\sigma_{12}^{4} g_{12} (\sigma_{12} )\rho_{2} \sqrt {\frac{{\pi k_{{\rm B}} T}}{{m_{1} m_{2} }}} \cdot $$
(8)

Here, g11(σ11), ρ1 and m1 are partial pair-correlation functions of the hard sphere, atomic density and mass, respectively. An analogous factor is defined for second species by interchanging index 1 with 2 in Eq. (8). The second assumption relates to the frequency of the local structural relaxations associated with diffusive atomic motions dominated by cage effect which is proportional to the number of accessible configurations. In turn, the normalized diffusion D* is proportional to exp(S), where S is the excess entropy expressed in NkB units and D* can be estimated by reformulating Eq. (8) as:

$$ D^{*} = \left( {\frac{{D_{i} }}{{\chi_{i} }}} \right)^{{X_{1} }} \cdot \left( {\frac{{D_{j} }}{{\chi_{j} }}} \right)^{{X_{j} }} . $$
(9)

If the Scc(0) function is known, the interdiffusion coefficient of a binary liquid, DAB, can also be calculated using Eq. (10a):

$$ D_{\text{AB}} = \frac{{S_{\text{CC}} (0)}}{{S_{\text{CC}} (0,{\text{id}})}} \cdot D_{\text{id}} , $$
(10a)
$$ S_{\text{CC}} (0) = k \cdot T \cdot \left( {\frac{{\partial^{2} G_{\text{mix}} }}{{\partial x^{2} }}} \right)^{ - 1} . $$
(10b)

D id is the chemical diffusion coefficient for an ideal mixture. Importantly, Eq (10a) explores the structure–dynamics relation involving only macroscopic thermodynamic information and links the short-range order to diffusion phenomenon [29]. Gmix is sum of Gid and Gxc, the latter is expressed in terms of thermodynamic parameters (L Liq i , i = 1, 2, 3) standing for Redlich–Kister polynomial coefficients [41].

The second important transport property, that is readily computed using molecular dynamics simulations, is shear viscosity η. It can be accessed through the Stokes–Einstein relation [42] [Eq. (11a)] or Green–Kubo formalism [43] [Eq. (11b)], both of them describe the collective motion of a macroscopic particle and are formulated as:

$$ \eta_{\text{SE}} = \frac{{k_{{\rm B}} T}}{{6 \cdot \pi \cdot D_{\text{AB}} \cdot R}} $$
(11a)
$$ \eta_{\text{GK}} = \frac{1}{{k_{{\rm B}} \cdot T \cdot V}}\int\limits_{0}^{\infty } {\left\langle {\sigma^{xz} (t) \cdot \sigma^{xz} (0)} \right\rangle \,{\text{d}}t} . $$
(11b)

Here σxz is an off-diagonal element of the stress tensor computed by averaging over time, and similar averages are computed for the equivalent yz and xy components of the stress tensor. R corresponds to the position of the first maximum of the radial distribution function.

Results and discussion

Thermodynamics-based modelling

Keeping in mind the brief description of the CFM method, the mixing Gibbs free energies (Gmix) were computed at 1345 K assuming the presence of two “associate” types, i.e. Al2Cu and Al2Cu3; the results are plotted in Fig. 1. The interaction energy parameters W ij , used for fitting Gibbs free energy to the one computed employing the Redlich–Kister (RK) polynomial (solid line in Fig. 1), are given in Table 1 in the form of their temperature dependence for both associate types considered. Here, the RK’s Gibbs energies are those computed using thermodynamic parameters (iLLiq, i = 0, 1, 2, 3) taken from Witusiewicz et al. [44] to Redlich–Kister polynomial [41] and they are collected in Table 1.

Figure 1
figure 1

Computed integral Gibbs energy of mixing as a function of Cu concentration at 1345 K using CFM formalism with assumption of: circle Al2Cu and rectangle Al2Cu3 “associate” type, solid line—RK polynomial [44]

Full size image
Table 1 Cluster types, interaction parameters W, WAB, W AA and W BB used for CFM model and thermodynamic parameters of the Gibbs energy of Al-Cu liquid phase [44]
Full size table

Two important features can be read from Fig. 1, with regard to the Gmix minimum position and the convergence of the experimental points with those obtained from CFM computations. First of all, for both atomic cluster assumptions used for CFM, the Gmix minimum position was located at the same Cu concentration (the Al40Cu60 alloy). Secondly, consideration of the Al2Cu atom cluster gives a slightly better agreement of CFM-based Gmix with the RK’s Gmix as opposed to the second atomic cluster considered in this work.

Diffusivity: CFM computation versus MD simulations

Below, we report the self-diffusion (SD) coefficients of each species composing Al–Cu alloys and the mutual diffusion coefficient. The temperature dependence of the SD coefficient for Al and Cu species is graphically presented in Fig. 2a, b, respectively. One can observe that the SD coefficients of both Al and Cu species exhibit significant deviation from a linear dependence for Cu-deficient content of Al–Cu alloys which becomes linear again for Cu-enriched alloys. An isotherm of SD coefficients computed at 1345 K for Al and Cu atoms in liquid Al–Cu alloys is plotted in Fig. 3, which includes also SD coefficients taken from Meyer et al. measurements for Al [45] and Cu [46] and represented by diamonds. Therefore, a significant difference in the SD coefficient values between Al and Cu atomic species can be easily seen from Fig. 3 for the investigated Al–Cu liquid alloys. Another striking feature is an appreciable decrease in SD coefficients between 10, 20 and 30 at% Cu content as demonstrated in Fig. 2a, b. Such changes could be attributed to the variety of intermetallic phases present in solid Cu-rich Al–Cu alloys, in which the Al atoms are more likely to be chemically bonded to Cu atoms compared to the Al-rich alloys. Regarding the concentration-dependent SD coefficient of Al, the present MEAM-MD simulations predict its almost linear decrease with increasing Cu content in liquid Al–Cu alloys. A substantial disagreement is found between the present studies and the AIMD data of Wang et al. [16], which particularly concerns the Al80Cu20 alloy, and this issue has been previously discussed by Trybula [5]. Interestingly, SD coefficient of Al in liquid Al–Cu alloys has attracted less scientific attention compared to that of copper, and this produces some difficulty in its discussion. In the case of SD coefficient of Cu species, investigated for selected binary Al–Cu alloys and drawn in Fig. 3, the MEAM-MD simulations also give its linear change and demonstrate a very good agreement with previous AIMD simulations performed for the liquid Al60Cu40 [18] and Al80Cu20 alloys [16]. Such observation does not apply to experimental data by Dahlborg et al. [21], whose values are substantially higher than both MEAM-MD and AIMD simulations results. To our best knowledge, no regular studies on diffusivity over a broad Cu concentration range using a single computational or experimental technique have been performed so far, which introduces a slight difficulty in their validation.

Figure 2
figure 2

a, b Temperature-dependent self-diffusion coefficients of Al and Cu atomic species of liquid Al–Cu alloys, respectively

Full size image
Figure 3
figure 3

SD coefficients isotherm as a function of Cu concentration computed at 1345 K. AIMD data—Al80Cu20 [16] and Al60Cu40 [18], diamonds—expt. data for pure elements [45, 46], and asterisk—simulation data taken from [5]

Full size image

The mutual diffusion coefficient of Al–Cu liquid alloys is drawn as a function of Cu concentration in Fig. 4. The linear change in observed composition dependence of self-diffusion coefficients is incompatible with that observed for the mutual one which is observed for the present MEAM-MD simulation results and CFM calculations, as it is illustrated in Fig. 4. Indeed, MEAM-MD and CFM-Al2Cu calculations give results quite close to each other compared to an appreciable difference recorded for the CFM-Al2Cu3 data. The most striking feature in DAB coefficient behaviour is seen for Al-rich alloys: it consists in a nonlinear decrease of the DAB coefficient with increasing Cu content which gradually changes over to a slight increase near the Al40Cu60 alloy composition. When CFM calculation results for DAB coefficient are confronted with the ones computed using Eq. (10a), the CFM-Al2Cu results appear to correlate well with Eq. (10a) compared to MEAM-MD simulation data exhibiting a worse agreement. Apart from that, both “associate” assumptions employed in the CFM formalism, MEAM-MD simulation and Eq. (10a), are consistent in describing the negative deviation from linearity represented by a long-dashed line in Fig. 4, corresponding to the chemical diffusion coefficient of ideal solution, Did. The extended Darken’s formulation has also been considered in our investigation for estimation of the DAB coefficient. However, it gave significantly overestimated DAB coefficients for Al–Cu liquid alloy compositions in the vicinity of the intermetallic phases existing in solid state and consequently, they have not been depicted in Fig. 4. According to Trybula et al. studies [9], the assumption of A2M3 “associate” type enables one to predict a nonlinear change in mutual diffusion coefficient, even if various arbitrary parameters are used. Importantly, its curvature is synchronized with the shape of Scc(0) function correlating with the thermodynamic driving force. As will be shown below, the two assumptions used in this work give two different descriptions of the fluctuation function for the liquid Al80Cu20 alloy. To prove the aforementioned statement, we need to have complete information derived from structural studies that will also be presented below.

Figure 4
figure 4

Mutual diffusion coefficient computed: squares—CFM formalism [Eq (7)], long-dashed line—Did; solid line—Eq. (10a), and red stars—MEAM-MD simulations. Lines are guide to the eye

Full size image

Since a direct method for the determination of the concentration fluctuation function using molecular dynamics simulations has not been given yet, and its experimental determination is extremely difficult, first we try to discuss the data obtained using CFM formalism in our previous study [9], plotted in Fig. 5. Two different formalisms have been used for the thermodynamics-based computations, that is FVM [4] and CFM (this work), which gave quite similar Scc(0) function behaviour. In our present CFM computations, other atomic clusters have also been considered; however, we dismissed them from further discussion due to significant misfits in mixing thermodynamic properties giving an unrealistic description of concentration fluctuation and dynamic properties. Since the Al2Cu “associate” type is considered, a flat maximum located around 30 at% of Cu is observed along with almost perfect agreement to the approximated experimental SCC(0) function (dashed line in Fig. 5). The latter has been assessed by using Eq. (10b) with Gmix = (Gxc + Gid) calculated using the RK polynomial. At the same time, consideration of the Al2Cu3 cluster allows us to observe a sharp maximum at 20 at% Cu content and a minimum spanned over 60–80 at% Cu. Importantly, the CFM formalism with assumption of the Al2Cu3 “associate” predicts the Al90Cu10 and Al80Cu20 alloy to be ideal solutions, while both FVM and CFM-Al2Cu computations as well as B–Th theory data contradict such behaviour. The latter have been computed as a second derivative of mixing Gibbs free energy with respect to concentration—Eq. (10b) [47]. Despite this fact, a substantial disagreement between FVM and CFM-Al2Cu calculations is seen for Al-rich Al–Cu alloys when they are confronted with B–Th results (dashed line in Fig. 5).

Figure 5
figure 5

Concentration fluctuation function in the long-wave limit, Scc(0): solid line—ideal solution, dashed line—B–Th model [Eq. (10b)] [47], whereas symbols—modelled data at 1345 K: open triangles—CFM formalism, and stars—MEAM-MD data

Full size image

The CFM-Al2Cu calculation gives a significantly smaller discrepancy from the approximated experimental data in comparison with the CFM-Al2Cu3 results. Such a good correlation undoubtedly gives the first evidence of Al2Cu associate importance, and it might suggest their preponderance in liquid Al–Cu alloys. Indeed, the trend observed in Scc(0) change is compatible with the one that can be read from the phase diagram [11], namely intermetallic phase absence on the Al-rich side of the phase diagram corresponds to the appearance of Scc(0) maximum in comparison with a smooth decrease observed for high-Cu-content Al–Cu liquid alloys. The present MEAM-MD simulations give slightly higher values of the Scc(0) function compared to the CFM-Al2Cu calculations and experimentally predicted data. However, they are still lower than line representing the ideal solution and, importantly, this does not change the structural description of liquid Al–Cu alloys. Because the Scc(0) function values for Al–Cu liquid alloys are lower than for ideal solution (Scc(0)/Scc(0, id) < 1), the Al–Cu alloys rather have a tendency for compound formation according to the definition given by Bhatia and Thornton [47]. In turn, the present CFM-Al2Cu and MEAM-MD simulation results are in line with this observation and FVM computation data [23] over the whole range of Cu concentrations investigated, as opposed to the selected ones found by our present CFM-Al2Cu3 calculation.

Structure and local atomic arrangements: molecular dynamics study

The structure factors S(q) for a series of liquid Al–Cu alloy compositions investigated in this work are plotted in Fig. 6 and confronted with X-ray experimental data measured by Brillo et al. [48] (black squares) and AIMD simulation results [16] (open circles). The S(q) values for each of the studied Al–Cu alloys exhibit a pronounced, narrow first peak becoming much more intense with increasing Cu content, while negligible differences in damped oscillations are seen for the second and third peaks. Such changes signify that Al–Cu alloys become more densely packed liquids with increasing Cu content and inform about the loss of randomness in atomic packing of Al–Cu alloys, that is also compatible with the observed decrease in Scc(0) function and the analysis of the Al–Cu phase diagram. Indeed, the existence of at least four different copper-rich intermetallic phases (θ-Al2Cu, η-AlCu, γ-Al4Cu9, β-Bcc) [44] can imply on observing an increase in alloy ordering. The MEAM-MD simulations data are in acceptable agreement with available experimental results for the Al67Cu33 liquid alloy measured at 1023 K [35] and AIMD simulations performed for the Al60Cu40 liquid alloy at 1323 K [18], as position of the first peak is considered. A slightly different scenario can be perceived from the trend in intensity change of the first peak of S(q) (Fig. 6). The MEAM-MD simulations show its increase with increasing Cu content, while AIMD simulations gave a less intense first peak compared to the experimental data.

Figure 6
figure 6

Structure factor of liquid Al1−xCux alloys: solid, dashed and dotted lines represent MEAM data for x equal to 30, 40 and 50 at%, respectively. Open circles and black squares—AIMD results for liquid Al60Cu40 [18] and X-ray experimental data for liquid Al67Cu33 alloy [35], respectively

Full size image

Interestingly, MEAM-MD computations allow us to observe a weak pre-peak occurrence at short q for liquid Al1−xCux alloys for x = 0.3, 0.4 and 0.5 which can also be seen for X-ray and AIMD data in comparison with its absence in the Al80Cu20 liquid alloy [16]. Analogous observations emerge from the analysis of total pair-correlation function of liquid Al1−xCux alloys for x = 0.3, 0.4 and 0.5 drawn in Fig. 7a. Some insignificant mismatch between present MEAM-MD results and available literature data is found, concerning the intensity of the first peak and shape of the second and higher-order peaks, whereas the position of the first peak is moving towards shorter interatomic distance with Cu content growth and approaching to that measured by Shtablavyi et al. [15] for pure Cu element. For clarity of Fig. 7a, MEAM and AIMD data for the liquid Al80Cu20 alloy are not appended; they are discussed in detail elsewhere [5].

Figure 7
figure 7

a Total and b-d partial pair distribution functions for selected liquid Al–Cu alloys. Solid, dashed and dotted lines—MD simulations for liquid Al70Cu30, Al60Cu40 and Al50Cu50, respectively; symbols (open circles)—AIMD data [16] and (full squares)—experiment [35]

Full size image

Partial pair-correlation functions, contributing to the total correlation function, g(r), that is gAlAl(r), gAlCu(r), gCuCu(r), are graphically presented in Fig. 7b-d, correspondingly. As it can be seen, MEAM-MD data also show high consistency with X-ray diffraction [35] and AIMD simulations data [16, 49] with exception of the gCuCu partial. According to our recent observation made for the Al80Cu20 liquid alloy [5], an adequate inconsistency in the position and shape of the high-order peaks in gCuCu(r) partial has also been found. Therefore, it might mean that the atomic pair interactions described within MEAM method potential for Cu element need to be treated more attentively, which is beyond the scope of the present work. Apart from the slight inconsistency exhibited by the MEAM-MD data, any splitting of the second peak has not been seen for each of the g(r) partials. If it had been observed, such peculiarity would have been suggestive of the icosahedral short-range ordering [16]. Interestingly, both experimental and AIMD data support the SCC(0) function discussed above, which provides information that will be essential in further discussion concerning the topological local atom surrounding study.

By integration of each partial g ij (r) up to the first minimum, average partial and total coordination numbers (CN) were computed; they are collected in Table 2. The average total CN for the investigated liquid Al–Cu alloys increases with increasing Cu content, and importantly, and such trend is compatible to that displayed by AIMD simulations for Al80Cu20 and Al60Cu40 liquid alloys [16, 18]; however, it is less pronounced. A significant rise in total coordination number obtained by performing MEAM-MD simulations can signify the occurrence of structural changes in atomic packing of liquid Al–Cu alloys with Cu content growth. As it can also be seen, MEAM potential overestimates both experimental and AIMD data resulting from highly coordinated Al central atoms. It needs to be stressed that the MEAM potential usually gives higher coordination numbers with regard to experiment or AIMD data due to its sensitivity to the first minimum position of g ij (r) [5]. Importantly, recent experimental reports delivered a quite different look into the relationship between CN and Cu composition. The CNs increase substantially with increasing Cu content between 17 and 34% of Cu–Al83Cu17 and Al66Cu34 alloys [35]. A similar relationship to that observed for our present MEAM-MD simulations has found confirmation in a previous study of Al-Ni liquid alloys structure for which total and partial CNs increase with increasing Ni content [1]. Due to the similarity in thermodynamic and/or physicochemical properties’ concentration and temperature dependences [50], similar trend in CNs change is also expected to be seen for the Al–Cu liquid alloys.

Table 2 Symbol, method and coordination number
Full size table

Taking into account the information presented above, the atomic short-range order can be quantified by computing the Warren–Cowley parameter (α). Therefore, we compared all thermodynamics-based estimations—Eq. (5) and B–Th model [Eq. (10b)] with MEAM-MD simulations. For the latter case, we used random-packed model [51] defining the parameter αtot as:

$$ \alpha_{\text{tot}} = \frac{{c_{\text{Cu}} (Z_{{{\text{Al}} - {\text{Al}}}} - Z_{{{\text{Cu}} - {\text{Al}}}} ) + c_{\text{Al}} (Z_{{{\text{Cu}} - {\text{Cu}}}} - Z_{{{\text{Al}} - {\text{Cu}}}} )}}{{Z_{\text{tot}} }}. $$
(12)

The quantified chemical short-range order (CSRO) is gathered in Table 3 for MEAM-MD simulations data, CFM computations and available literature data [30] computed at 1345 K.

Table 3 Methods used and chemical short-range order parameter (CSRO)
Full size table

CFM computations performed assuming Al2Cu and Al2Cu3 give negative CSRO parameter, and therefore, its sign confirms the presence of chemical ordering in the liquid Al–Cu alloys which does increase with growing Cu content. As it can be seen, CFM-Al2Cu calculations predict slightly higher CSRO parameter compared to the one estimated using Eq. (10b) and experimental thermodynamic data [30]. For further comparison, a recent study of Das et al. [52]. provided some interesting observations related to chemical description of liquid Al–Ni alloys, involving diffusion and CSRO. Basically, the nonlinear compositional dependence of diffusion coefficients corresponds to increasing CSRO, moving towards the Ni-rich side. Note that a similar relationship is exhibited by the analysed liquid Al–Cu alloys, namely the growth in CSRO with increasing Cu content reflects in decreasing self- and mutual diffusion coefficients. By definition, the CSRO parameter describes chemical local atom environment; therefore, its increase informs that liquid structure becomes densely packed reflecting in a significant deviation from random atom ordering. Such behaviour is a good example of the hard-sphere-like liquids behaviour, for which mass transport is strongly correlated with packing fraction of atoms [53, 54]. To have a complete view onto the mass transport property, we decide to extend our consideration on viscosity investigation, which will be discussed in the following section.

Three-dimensional analysis based on the Voronoi tessellation was performed, in which a geometrical construction divides the investigated space into polyhedra, containing a single atom, called Voronoi cells (VC). The local atom topology as a function of composition and temperature is inspected. Within Voronoi analysis, the total number of the VC faces corresponds to the coordination number assigned to the central atom. Voronoi index \( n_{3} , n_{4} .n_{5} , n_{6} , \ldots . \) is used to differentiate between the VC types, where \( n_{i} \) is the number of i-edged sides in the polyhedron. This structural tool allows to distinguish body- and face-centred cubic (b.c.c. and f.c.c.), as well as hexagonal close-packed (h.c.p.) structures. Further details on the Voronoi tessellation method, they can be found elsewhere [55, 56]. Ten types of Voronoi cells (VC) are found to be in reasonable abundance among the 100 different types of cells around Al and Cu atomic species in liquid Al1−xCux (x = 10, 30, 40 and 50) alloys at 1345 K. VCs abundance and their geometrical construction are drawn in Fig. 8.

Figure 8
figure 8

The abundance of Voronoi cells (top) found at 1345 K around Al and Cu species in liquid Al–Cu alloys and their geometrical construction (bottom)

Full size image

High coordination numbers around Al and Cu species correspond to high-indexed Voronoi cells (VCs). As compared to the AIMD data at 1345 K for the Al60Cu40 liquid alloy [49], some similarities have been found concerning abundance and cluster types distinguished in the Al–Cu liquid alloys examined. At first, icosahedral clusters of < 0 0 12 0 > index—ideal icosahedron are found to be not the most abundant index compared to distorted icosahedra of {< 0 1 10 2 > , < 0 3 6 4 >}. It is argued that local icosahedral structures have a positive influence on stabilization of undercooled liquid state [57] resulting from their compact structure; therefore, we can expect that this issue is also applied to liquid Al–Cu alloys. On the other hand, their preponderant abundance has a beneficiary influence on slowing down relaxation processes distinguished during glass transition. Compact structures, like icosahedra derivatives, can be responsible for slow structural dynamics compared to loosely packed clusters undergoing fast structural dynamics [58, 59] resulting in significantly displacing local atom regions. Recognizable discrimination in VCs as it has been found between Al and Cu local surrounding could have indicated that liquid Al–Cu alloys undergo two relaxation processes. Invoking at this moment the previously discussed functions, representing chemistry of local atom surrounding: structure factor and distribution function, any evidence of middle-range ordered clusters has not been found. Therefore, it can be stated that Al–Cu liquid alloys undergo processes which enable liquid materials to crystallize rather than have an amorphous structure, excepting only the Al-rich alloys, due to the existence of a number of interatomic phases stabilizing the solid state of the Al–Cu system [11]. Quite interestingly, such a conclusion was confirmed by Wang et al. [16]; in particular, the liquid Al80Cu20 alloy can be classified as a strong liquid according to Angel classification [60] with fragility value lower than 10 that provides a straightforward suggestion of weak glass forming ability.

Viscosity

For viscosity determination, at first, we confine our investigation to test the validity of Stokes–Einstein (S–E) [42] relation for Al–Cu liquid alloys using only MEAM-MD data. Its compositional dependence is shown in Fig. 9. As compared to the experimental data collected by Brillo [48], the inability of S–E relation to accurately reproduce the experimental data (red stars) can be easily seen which corresponds to a substantial underestimation of experiment. Such behaviour eliminates the S-E equation consideration from further studies of liquid Al–Cu alloys in contrast to the simple liquids, for which this relation is obeyed. Consequently, we employ the Green–Kubo (G–K) formalism [43] [Eq. (11b)] as previously Asta et al. [61] used for liquid Al–Ni alloys investigation. Importantly, the G–K approach (circles) allows to get an overall 5–15% agreement with the experiment over a broad Cu concentration range.

Figure 9
figure 9

MEAM MD computed (full and empty squares) and measured (stars) viscosity at 1345 K as a function of Cu concentration. Experimental data were taken from Brillo et al. [48]. Lines are guide to the eye

Full size image

The nonlinear behaviour of viscosity is strongly correlated with nonlinear diffusion coefficients change as a result of the CSRO degree increasing with Cu content growth. Quite interestingly, the experiment also exhibits the same behaviour [21], compared to an almost complete failure given by the S–E relation. On the Al-rich side, where no intermetallic phases exist, viscosity slightly decreases towards the Al80Cu20 alloy and suddenly increases as the Al2Cu compound starts to form [11]. In turn, an estimation of the fragility parameter [60, 62, 63] is not necessary due to a negligible evidence of icosahedral motifs or Wang et al. remarks [16] as well as significant Al2Cu compound influence stimulating the solidification of liquid Al–Cu alloys [14, 64,65,66].

Excess entropy

To define the structure–dynamics relationship in liquid Al–Cu alloys, we considered the universal scaling law relating diffusion coefficient to compute the excess entropy as formulated in the original work by Dzugutov [27]. In the frame of this formulation, the excess entropy (expressed in k B units) is restricted to a two-body contribution:

$$ S_{2} \, = \, - 2\pi \rho \int\limits_{0}^{\infty } {\{ g(r)\ln [g(r)]\, - \,[g(r) - 1]\} \,r^{2} dr} , $$
(13)

for ρ being the number density. In other words, Dzugutov scaling law is an important contribution combining the kinetics of a liquid to a static thermodynamic quantity, S. In our case, it is S2 function coupled to an accessible structural property of a liquid—partial pair-correlation function, whereas the total excess entropy is given by S2 = XAlS Al2  + XCuS Cu2 , where S Al2 and S Cu2 are partial molar entropy of Al and Cu atoms, respectively. The partial and total S2 function computed at 1345 K are plotted in Fig. 10a.

Figure 10
figure 10

a Computed total and partial entropies (S2) as a function of Cu composition at 1345 K and b scaled diffusion versus entropy for liquid Al–Cu alloys and pure elements: Al [27] and Cu [61]

Full size image

For the investigated Al–Cu alloys, we observe a decrease in total excess entropy (S2) with increasing Cu content as well as substantial discrimination in partial entropy between Al and Cu elements on the Al-rich side which decreases with increasing Cu content. Displayed compositional dependences of S2 and partial entropies correlate with the chemical short-range order (CSRO) changes discussed before. Basically, low CSRO degree on the Al-rich side corresponds to a low negative S2 value compared to middle Cu composition for which significant degree of CSRO and abundance of icosahedral motifs have been recorded. This is a good demonstration of S2 function ability to appropriately reflect the local atom arrangements within the nearest neighbour environment, which can be characterized by CSRO or ISRO (icosahedral short-range ordering), in spite of negligence of many-body correlations. Similarly to the impact of CSRO and ISRO on the S2 function compositional behaviour of Al–Cu liquid alloys, we also notice such confirmation for self-diffusion coefficients’ changes. The decrease in both SD coefficients corresponds to respective decreases in S Al2 and S Cu2 entropies. In relation to this, the significant differentiation in SD values and average coordination number between Al and Cu species has a big influence on the trend in excess and partial entropies.

To clarify the effect of diffusion on excess entropy, we inspect the scaled diffusion coefficient (D*) using Enskog theory which is presented in Fig. 10b as a function of excess entropy S2 at 1345 K. It has been computed for the all Al–Cu alloys investigated and confronted with the literature data for pure liquid Al [27] and Cu [61] metals. Normalized diffusion D* computed, using Eq. (9), scales adequately to those displayed by pure liquid Al and Cu atoms, excepting only the liquid binary Al50Cu50 system exhibiting a significant coefficient shift. As it can be seen, D* decreases with increasing (-S2) function and its decrease is compatible with that exhibited in Fig. 10a.

Conclusions

In conclusion, the structure and chemistry of liquid Al–Cu alloys have been studied by means of molecular dynamics, using modified embedded atom model potential and CALPHAD-type modelling based on a quasi-chemical approximation for weakly compound-forming melts (CFM formalism).

The CFM computations provided important knowledge on the interplay between thermodynamics, kinetics and structure over a broad Cu concentration range which allowed to supplement the database of Al–Cu liquid alloy properties. The assumption of Al2Cu compound dominance within the CFM formalism allowed to accurately reproduce the available experimental data, and consequently explain the experimental structural characterization of two liquid Al–Cu alloys performed by Brillo et al. [35].

Compositional dependence of concentration fluctuation function, Scc(0), displayed significant deviation from regular solution comparable with predictions of the CFM formalism and experimental data. In turn, the nonlinear change of diffusion coefficients (both self and mutual) and viscosity is associated with pronounced chemical short-range ordering (CSRO) increasing with Cu growth confirmed by present molecular dynamics simulations. A significant effect of increasing dense-packing liquid Al–Cu alloys on excess entropy corresponding to an appreciable decrease in partial and total entropies has been found.

To summarize, using the hybrid approach presented herein, we gave a compact and detailed description of the structure, thermodynamics and kinetics of liquid phase, valuable for manufacturing of solid materials based on Al–Cu alloys. The studies presented in this work are the first ones concerning the Al–Cu alloys and the used methodology, which allowed us to find another hard-sphere-like metallic system whose transport properties are strongly affected by packing effects. Thermodynamics-based modelling allowed to quickly cover a broad concentration range compared to molecular dynamics simulations employing semi-empirical interatomic forces fitted to ab initio calculation results. Although the former calculations do not explicitly consider the atomic level, they need to be supported by atomistic computations to gain deeper insight into the physical phenomenology of any studied melt.