A predictive model for the diffusion of a highly non-ideal ternary system

Diffusion plays a central part in many unit operations. The Maxwell–Stefan model is the dominant model for both gaseous and liquid diffusion. However, it was developed from the kinetic theory of gases, raising the question of whether it can be extended to non-ideal liquid systems. The dynamic fluctuation model is an alternative model based on the Cussler theory and predicts a smaller thermodynamic influence relative to the linear influence of the Maxwell–Stefan model due to dynamic concentration fluctuations. Since the dynamic fluctuation model, which uses the scaling factor a, had improved performance relative to the Maxwell–Stefan model for a wide range of binary systems, it is postulated that this improved performance should also be observed for a ternary system. In this work, the dynamic molecular fluctuation model was extended to a highly non-ideal ternary system, using the same scaling factor a, through matrix manipulation. Using self-diffusion data measured by NMR, mutual diffusion predictions of the developed model and the Maxwell–Stefan model were compared to experimental mutual diffusion data of the partially miscible system ethanol/toluene/n-decane. It is demonstrated that the dynamic fluctuation model gives improved predictions relative to the Maxwell–Stefan approach, consistent with previous observations on binary systems, showing that the reduced thermodynamic influence of the dynamic fluctuation model is an improvement. In addition, we show that the use of local mole fractions, to account for molecular association, in both the dynamic fluctuation and Maxwell–Stefan models, results in improved diffusion predictions for the ternary system. The results confirm that the dynamic fluctuation model improves predictions of mutual diffusion in liquid mixtures, suggesting a non-linear correction to the thermodynamic correction factor. The results also suggest that that the key assumptions in the Maxwell–Stefan model and its derivation, rooted in the kinetic theory of gases, are not entirely accurate for highly non-ideal liquid systems. The optimum a for the ternary system studied here is approximately 0.45, similarly to the optimum a of 0.40 to 0.80 for a range of binary systems previously studied, suggesting that the use of the a scaling factor, which is grounded in scaling laws theory, is of general validity.


Introduction
Diffusion is an important natural phenomenon in chemical processes. Nearly all industrial separations are influenced by liquid component diffusion, and in certain cases, diffusion may be the rate limiting process. Therefore, the ability to accurately estimate the diffusion coefficients is necessary for optimum unit design.
The theories governing liquid diffusion are incomplete and limited experimental diffusion data are available. This problem is exacerbated for multicomponent liquid diffusion. [1][2][3][4] Considering the industrial relevance of liquid diffusion, an accurate diffusion model that requires easy to measure parameters must be developed.
The Maxwell-Stefan (MS) model is the dominant diffusion model for gaseous systems and it is also applied to liquid systems due to its simple, linear dependence on the chemical potential gradient. This model was developed using the kinetic theory of gases 5,6 and it has been extended unaltered to highly non-ideal liquid systems without rigorous experimental validation. Due to both a lack of experimental data, and the MS diffusion coefficients being impossible to measure directly, exhaustive molecular dynamic simulations have been pursued by some researchers. 7,8 These exhaustive calculations are impractical for general engineering purposes. An alternative diffusion model is the dynamic fluctuation model, based on the Cussler theory of diffusion. [9][10][11][12] This model accounts for dynamic concentration fluctuations, which results in a smaller thermodynamic influence on diffusion relative to the MS model. The dynamic fluctuation model has been proven to be a significant improvement over the MS model for binary systems. [9][10][11] The purpose of this paper is to extend this dynamic fluctuation model to a ternary system and compare it to the MS model. A partially miscible ternary system is used to evaluate the two models, as this highly non-ideal system facilitates model discrimination.

Theory
The understanding of diffusion is dominated by two theories, the Maxwell-Stefan model and Fick's law. The key principles and requirements for the two theories are outlined in the sections that follow, with a particular focus on the potential inadequacies of the Maxwell-Stefan model for highly non-ideal systems.

Maxwell-Stefan model
The Maxwell-Stefan (MS) model is developed from the kinetic theory of gases, by assuming that the force acting on an individual particle is in equilibrium with the resistance to fluid motion due to collisions with other species. Fig. 1 shows the MS model of diffusion for species 1 through species 2.
There is extensive experimental evidence showing that the MS model is accurate for multicomponent, dilute gases. 5,13 For a gas, the pressure gradient can be related to the chemical potential gradient, rm. The MS model with a rm driving force is: 14 where x i (or x j ) and J i (or J j ) are the mole fraction and molar flux, respectively, of component i (or j), Ð ij is the MS diffusion coefficient between components i and j, T is the system temperature, R is the universal gas constant and C t is the total molar concentration. Considering the fluid resistance analogy, Ð ij can be viewed as an inverse drag coefficient. Eqn (1) shows how rm i is the driving force for diffusion, as required by thermodynamics. 5 It is important to stress that thermodynamics only requires rm i = 0 at equilibrium and that species move from a higher to a lower chemical potential; it does not place a requirement on the form of the chemical potential driving force for diffusion. The linear driving force in eqn (1) is a convenient result and not a thermodynamic requirement. For an n component system, rm i is related to the activity coefficient, g i , by: where d ij is the Kronecker delta and G ij are the elements of the matrix G, the latter being the thermodynamic correction factor for the system. 5,6 For an n component system the MS model can be expressed in terms of the n À 1 component molar flows as: where rx is a column vector of the mole fraction gradients for the n À 1 components.

Fick's law
Fick's law proposes a linear dependence between the diffusion flux and the concentration gradient based upon experimental observations; it is not a theoretical development like the MS model. The generalised Fick's law for an n component system is written as: 5 where [D] is an n À 1 Â n À 1 matrix of mutual diffusion coefficients. In an n component system there are only n À 1 independent fluxes and so [D] implicitly contains information on the final component. Thus, the multicomponent diffusion coefficients are not directly related to the binary diffusivity coefficients, making it difficult to develop a diffusion model based on the binary coefficients alone.

Estimating ideal mobility coefficients
Since both Fick's law and the MS model describe the same physical process, they are related. Comparing eqn (3) and (5), the [D] predicted by MS model is: This shows that the mutual diffusivities obtained from the MS model are composed of two distinct parts, namely hydrodynamic ([B] À1 ) and thermodynamic ([G]) influences. [D] can be obtained experimentally, as it is related to the concentration gradient. Therefore the MS model's accuracy can be determined by comparing its predicted [D] to that experimentally measured. However, Ð ij is not accessible by experiments and unlike gases, there is no universal kinetic theory for liquids that allows Ð ij to be predicted. 5 Hence eqn (6) has been assumed to be correct and Ð ij have been calculated using [G] and the measured [D].
Due to both experimental and simulation difficulties, numerous semi-empirical relationships for predicting Ð ij have been developed. However, all these models are based on the assumption that the MS model is correct for liquid systems, which need not be the case. The two foremost means to estimate Ð ij are the Darken and Vignes models.
The Vignes equation is a logarithmic interpolation that only requires the infinite dilution values, 15 whilst the Darken model is a linear interpolation. 16,17 In previous work using pulsed gradient stimulated echo nuclear magnetic resonance (PGSTE-NMR) of the nearly ideal ternary system of methanol/butan-1-ol/propan-1-ol, 18 we have shown that the multicomponent Darken model of Liu et al. 8 is more accurate than that of Krishna and van Baten. 7 However, the average error for the system was 26%. This large error is primarily a result of large errors at low compositions, where the system was modelled as ideal, and erroneous mutual diffusion data, the latter being prone to high errors due to the challenges associated with mutual diffusion in ternary systems. The multicomponent Darken mobility model is: 8 where D mix is the effective diffusivity of the mixture. The Vignes model is an empirical model. 19 Since the objective of this paper is to compare the MS model and the dynamic fluctuation model for a ternary system, and due to there being a physical basis for the Darken model, eqn (7) will be used to estimate the ideal mobility. It is stressed that this work does not set out to assess the accuracy and the validity of the Darken model relative to the Vignes model.

Dynamic concentration fluctuation model
For binary mixtures, the MS model predicts a linear relationship between D and G: However, this is not experimentally observed for hexane/ nitrobenzene as the consulate point is approached and for numerous other non-ideal binary systems. [9][10][11][12] Previous work on binary systems has demonstrated that a model based on critical point scaling laws 20 successfully predicts mutual diffusion in non-ideal binary mixtures, which is given by: were D ij is the ideal mobility between species i-j and a is a parameter from dynamic scaling theory. 21 The parameter a is a system independent constant and is taken to be around twothirds. 5,22,23 The parameter a accounts for dynamic concentration fluctuations within the system. These occur because the driving force for the removal of these fluctuations is rm i , which tends to zero at the critical point. Hence, the processes removing such fluctuations become slow and the fluctuations become important in determining the physical properties of the mixture. At any point in a mixture with significant concentration fluctuations a diffusing species may experience a concentration either a little higher or a little lower than the mean concentration; the result will be that locally the thermodynamic correction factor is higher than that which would be expected from the mean concentration. 24 This provides a theoretical explanation as to why a o 1, hence the non-linear dependence. Further work showed that the dynamic fluctuation model, with a = 0.64, was accurate in predicting the diffusivities of 14 non-ideal binary liquid mixtures, systems which exhibited both positive and negative deviations from Raoult's law and were far away from the consulate point. 10,11 This provides further evidence that the MS model inaccurately accounts for the rm i driving force.
An explanation for the MS model's failure to accurately describe diffusion in highly non-ideal liquid mixtures is that the intermolecular forces become significant in such systems. These forces may result in clusters forming, with clusters and individual particles transported in the diffusion process, 12 rather than the single spherical particles undergoing random elastic collisions proposed by the kinetic theory of gases, as is modelled in the MS model. 5 Therefore, the extension of the MS model from dilute gases to highly non-ideal liquid systems does not have a strong theoretical foundation. The difference between the MS model and the underlying diffusion of highly non-ideal systems is shown in Fig. 2.
The success of the binary dynamic fluctuation model suggests that for a highly non-ideal ternary system a scaling power adjustment is required on G.
[B]'s dependence on Ð ij , eqn (4) and (7), is expected to remain unchanged, as the dynamic fluctuation model only adjusts G due to concentration fluctuations; hence, the hydrodynamic influence is unaffected. However, for the dynamic fluctuation model Ð is replaced by D, the ideal mobility, with the only difference between these two coefficients being their physical interpretation in that D is an ideal mobility term whilst Ð is an inverse drag coefficient.
Modelling the mutual diffusion of a ternary system faces an additional complexity, as the entire system is a matrix and raising the matrix G to a fractional power cannot be done elementwise. Raising a matrix to a fractional power requires diagonalisation; a derivation is provided in the Appendix for the interested reader. Using matrix diagonalisation, the ternary dynamic fluctuation model becomes:

Local composition dynamic fluctuation model
An implicit assumption in the Darken model, eqn (7), is that D can be determined from the measured self-diffusivities, D*, of single molecules. If, however, species' motion is highly correlated due to strong association, this assumption is physically inconsistent. This is because the measured motion of a labelled molecule is accompanied by the motion of one or more unlabelled molecules. This problem can be overcome by using the local composition in the Darken model rather than the bulk composition. 23 The local composition form of the dynamic fluctuation model accurately predicts the mutual diffusion coefficients of both associating and non-associating species for binary systems. 23 The applicability of the model to a wide variety of thermodynamic systems suggests that the dynamic fluctuation model has a strong predictive capacity. Indeed, such a model has also been recently tested and assessed by Guevara-Carrion et al. 25 in a recently published work on mutual diffusion in binary mixtures, who concluded that this model was the most accurate predictive tool for mutual diffusion in binary systems.
Based on the success of using local mole fractions to account for association in binary systems, 23 the local mole fraction model is extended to ternary systems. For the general multicomponent system, the NRTL (Non Random Two Liquid) model relates the local mole fraction of molecule j and k in the immediate neighbourhood of molecule i. 26 To account for molecular association, the mole fraction of a component in its own immediate vicinity, x ii , replaces its bulk mole fraction x i in the Darken mobility model, eqn (7), to yield:D This impact on the ideal mobility only changes the hydrodynamic influence. The hydrodynamic influence for local composition, B*, is: The dynamic concentration fluctuation model using local compositions to account for molecular association is thus:

Experimental
Approaching the two-phase region of a partially miscible system will result in maximally non-ideal diffusion, and is hence the best means to identify the appropriate form of the thermodynamic correction factor. To test the ternary diffusion models across a wide-range of non-idealities, the coexistence curve was approached from the most ideal binary subsystem, toluene/n-decane. 27 This pathway results in the ideality of the system continuously decreasing. Eight sample points were evaluated, as seen in Fig. 3.

Thermodynamic data
The ternary system ethanol (1)/toluene (2)/n-decane (3) is a partially miscible system at 25 1C, and has phase equilibrium data available. 28 The thermodynamically consistent NRTL parameters for the system are presented in Table 1. 29 These values were used to determine G, the thermodynamic correction factor, in the diffusion model.

Self-diffusion measurements
To undertake the ideal mobility analysis, D* measurements are required. The pulsed gradient stimulated echo nuclear  magnetic resonance (PGSTE-NMR) method allows D* to be measured as a function of a species' NMR signal decay. 30 The details of the method have been reported elsewhere. 9,24,31 All NMR experiments were performed at a temperature of 298 K on a Bruker DMX 300 operating at a 1 H frequency of 300.13 MHz. The NMR diffusion experiments were carried out using a diffusion probe capable of producing magnetic field gradient pulses up to 11.76 T m À1 in the z-direction. The duration of the 901 pulse of the 1 H nuclei was 4.8 ms. A Bruker Variable Temperature unit, BVT 3000, was used to set the required temperature.
Self-diffusion measurements of the ethanol (1)/toluene (2)/n-decane (3) mixtures were conducted using the PGSTE pulse sequence. 32 Measurements were performed by holding the gradient pulse duration (d) constant and varying the gradient strength (g). The observation time (D) was set to 50 ms. All NMR spectra were referenced to bulk liquid TMS.

Mutual diffusion measurements
Mutual diffusion measurements were performed in a lab-scale equipment, constituted by a peristaltic pump, by which a solution of ethanol (1)/toluene (2)/n-decane (3) is flowing at a fixed flow-rate in a 20 m length capillary bore, characterised by an internal radius of r = 3.945 Â 10 À4 m. The pipe is set in a thermostatic bath working at 25 1C. The outlet flow is sent to a RI Detector (Knauer RI Detector K-2301) with a sensitivity of 3 Â 10 À8 RIU and a noise of AE1.5 Â 10 À8 RIU. The collected data were interpreted by applying the analytical solution of the laminar flow model reported in eqn (15) and (16).
D 1 and D 2 represent the eigenvalues of the function, while W 1 is the weight of the first peak function. a 1 is the fraction of the initial refractive index contributed by solute 1, defined as in eqn (17).
where R i (i = 1 or 2) is the derivative of the molar refractivity of solute i at constant concentration of the second solute.
The eigenvalues are weighted through a linear correlation; a and b are the related parameters. D 1 , D 2 , a and b, were obtained by parameter estimation analysis on the collected experimental data. In particular, the objective function written as in eqn (18) was minimised by using the ''particleswarm'' algorithm implemented in MATLAB libraries. 33 The mutual diffusion coefficients were finally calculated as in eqn (19)- (22).

Self-diffusion measurements
The self-diffusion measurements from the PGSTE-NMR experiments for the three components are shown below in Fig. 4. The values are reported for different ethanol mole fractions whilst the toluene/n-decane molar ratio is kept equal to one. The toluene and n-decane D* values only vary slightly over the composition range. However, D ethanol * at ethanol mole fractions less than 0.05 is substantially higher than that at ethanol mole fractions greater than 0.10. This suggests that ethanol undergoes strong molecular association at ethanol mole fractions above about 0.10. Above this mole fraction, D ethanol * decreases approximately linearly with increasing mole fraction.  4 Ethanol (1), toluene (2), and n-decane (3) self-diffusion measurements at 25 1C using PGSTE-NMR at the ethanol concentrations investigated. The toluene/n-decane molar ratio is kept equal to one.
The self-diffusion measurements are not at precisely the same mole fractions as the mutual diffusion mole fractions due to experimental difficulties involved in the measurements. To overcome this, a linear interpolation was used for D* of all the components. For ethanol, the linear interpolation was conducted in two separate sections, namely for x ethanol o 0.10 and x ethanol 4 0.10.

Ternary mutual diffusion results
A comparison of the experimental mutual diffusion data and the results from the Maxwell-Stefan model and dynamic fluctuation model, eqn (11), with a = 0.64 is provided in Fig. 5. Both models are a poor approximation for D 22 and D 21 at all the experimental points investigated. However, the most notable difference between the models is for D 22 .
Considering the small variation in D* for toluene and n-decane over the composition range of interest and D ethanol * being almost linearly related to the ethanol mole fraction system for x ethanol 4 0.10, Fig. 4, there is only limited source of uncertainty from the hydrodynamic term in both the MS, eqn (6), and dynamic fluctuation models, eqn (11). Therefore, the source of the discrepancy in the predicted D 22 and D 21 from both models must lie either in the accuracy of the thermodynamic data or in the form of the models.
It is reasonable to expect the difference between the experimental data and the two models to be smallest at the lowest x ethanol as this corresponds to a more thermodynamically ideal system. However, the ternary thermodynamic model was obtained from regressing the LLE (Liquid Liquid Equilibrium) coexistence curve. Therefore, as one moves further away from the coexistence curve, the ternary thermodynamic parameters become an extrapolation, although this influence is offset by the system becoming more ideal. Furthermore, in general the regressed ternary system parameters are unable to accurately model the binary subsystems VLE (Vapour Liquid Equilibrium). 34 However, as one proceeds to the coexistence curve with increasing x ethanol , as shown in Fig. 3, the ternary thermodynamic parameters are a more accurate representation of the system, in particular, at x ethanol = 0.309 when the system is close to the coexistence curve.
Nevertheless, eqn (11) is an improvement over the traditional MS model for the entire composition range when considering D 22 in Fig. 5. The absolute error of the models compared to the experimental data is provided in Table 2. The table clearly shows that the dynamic fluctuation model, eqn (11) with a = 0.64, is an improvement over the conventional MS model, eqn (6), except for D 21 for which the absolute improvement of the MS formulation is marginal. Consequently, it is reasonable to claim that the dynamic fluctuation model with a = 0.64 is an improvement over the conventional MS model.
The results when using the local mole fraction for both the MS and dynamic concentration fluctuation models are shown  in Fig. 6. As expected from its success in the binary systems, the local composition model provides improved performance. Comparing Tables 2 and 3 one sees that the use of the local mole fraction model gives a noticeable improvement in the accuracy of both models. The reason for this is that strong molecular self-association of a polar species is expected when it is mixed with a non-polar one. This requires modification of the average molecular mobility; for example; if a species is dimerised, the self-diffusivity of the dimerised species should be doubled in the form of doubling the tracer diffusivity of the dimerised species. 10 In this system, ethanol is expected to undergo molecular association, and this is consistent with the step change in D ethanol * at x ethanol B 0.10, Fig. 4. The use of a local mole fraction provides an a priori tool to account for the influence of cluster formation on the molecular mobility. 23 Fig. 7 shows the change in both the absolute diagonal and total error between the experimental mutual diffusion data and that of the dynamic concentration fluctuation model for varying values of a. Clearly, as one decreases the value a from 1 through to B0.40, both the total and the diagonal error decrease. However, decreasing a from 0.50 to 0.40 sees a negligible change in the model prediction, as the error terms remain roughly constant. An a below 0.40 implies a reduction of the diffusion prediction accuracy. Therefore, the optimum a is in the range of 0.40 to 0.60. This result is similar to that of binary systems, which also have an optimum a of 0.40 to 0.80. 10,11 This suggests that the scaling factor a is thus consistent with the theoretical scaling laws theory scaling factor of B0.64. Fig. 8 shows the result of the application of the dynamic concentration fluctuation model using local mole fraction with the optimal a of 0.45. There is a distinct improvement in the accuracy of the model, with this being most prominent for D 22 . Table 4 shows that the local mole fraction model is approximately an order of magnitude more accurate than the traditional MS model.  Table 3 The average absolute error for the mutual diffusivity diagonal and cross-diagonal elements using the dynamic fluctuation model, eqn (14) with a = 0.64, and the Maxwell-Stefan model with local mole fractions for ethanol (1)/toluene (2)/n-decane (3)

Conclusion and recommendations
The dynamic fluctuation model, which uses the scaling factor a to reduce the non-ideal thermodynamic impact based on scaling theory, was extended to ternary systems. The developed model and the most common model currently in use, the Maxwell-Stefan model, were compared to experimental mutual diffusion data for the partially miscible system ethanol/toluene/ n-decane. The analysis showed that the Maxwell-Stefan formulation inaccurately accounts for the thermodynamic influence on diffusion, the same conclusion as that for binary systems. [9][10][11] The inaccuracy of the Maxwell-Stefan model for both binary and ternary systems suggests that the key assumptions in its derivation, which are rooted in the kinetic theory of gases, are invalid for highly non-ideal liquid systems.
Considering the shortcomings of the Maxwell-Stefan formulation, molecular dynamic simulations are not an entirely appropriate means for further diffusion model validation. The developed model improves prediction of mutual diffusion; in particular, using such model with local mole fractions, to account for self-association of ethanol, substantially improved the match between predicted and measured diffusion coefficients. Thus, the model which best predicted experimental results was one incorporating local compositions into the hydrodynamic term, as well as a modulation of the thermodynamic correction factor due to dynamic concentration fluctuations, eqn (14). An interesting result is that optimum a for the ternary system investigated is approximately 0.45. This is similar to the optimum a of 0.40 to 0.80 for a range of binary systems previously investigated. 10,11 This suggests that the a scaling factor is uniform.
We note that further work is required to test the accuracy of the ternary dynamic fluctuation model, and determine whether or not the scaling factor a is uniform, as was the case for binary systems. 10,11 To reduce the uncertainties in thermodynamic modelling and thus improve model discrimination, ternary VLE data must be obtained for the chosen sample points at the temperature of interest. Ternary VLE data is preferred over LLE data, as VLE data can be obtained for partially miscible systems at sample points outside the two-phase coexistence region. Thus, more VLE data would allow the thermodynamics to be accurately modelled within the region of interest; it would no longer be an extrapolation from the two-phase coexistence region, as is the case with using thermodynamic parameters regressed from LLE data obtaining a high discrepancy between the experimental data and the model extrapolation.
It is important to note that ternary systems are complex to analyse, both in terms of experimental measurements of mutual diffusion and determination of accurate thermodynamic data and this is perhaps the reason why few systems have so far been investigated. However, as shown in this work, the use of the dynamic concentration fluctuation diffusion model with local compositions for binary systems 23 and its  Table 4 The average absolute error for the mutual diffusivity diagonal and cross-diagonal elements using the dynamic fluctuation model, eqn (14) with the optimum a = 0.45, and the Maxwell-Stefan model with local mole fractions for ethanol (1)/toluene (2)/n-decane (3)

Conflicts of interest
There are no conflicts to declare.

A. Appendix
Fractional power of matrix derivation Let us assume that the square matrix A can be diagonalised via eigenvalue decomposition and written as: where L is a diagonal matrix featuring eigenvalues of A along its diagonal, and V is the eigenfunction matrix of A. Then: and, in general, Since L is a diagonal matrix, it follows that: For a general real-valued square matrix, the eigenvalues can be negative or complex. In either of these cases, fractional powers of the matrix might feature complex entries. In the case of the diffusion matrix G raised to a: However, all entries are real since diffusion matrices are symmetric positive semidefinite matrices and possess real, nonnegative eigenvalues. Thus, fractional powers of diffusion matrices make physical sense.
The ternary dynamic fluctuation model becomes: The ternary dynamic fluctuation model just requires the eigen decomposition of the thermodynamic matrix, G.

Ternary local mole fraction model derivation
For the general multicomponent system, the NRTL (Non Random Two Liquid) model relates the local mole fraction of molecule j and k in the immediate neighbourhood of molecule i. 26 To account for molecular association, the mole fraction of a component in its own immediate vicinity, x ii , replaces its bulk mole fraction x i in the Darken mobility model, eqn (7), to yield: where e ij and e ik are NRTL non-randomness parameters, and g ji and g ki the energies of interactions. For each bulk composition for a ternary system, there are three local mole fraction relationships, eqn (A.7), and one local mole fraction balance, eqn (A.8). This provides a set of four equations that must be implicitly solved to yield the local mole fractions for each ternary composition. Greek letters a Dynamic fluctuation parameter (-) a 1 Fraction of the initial refractive index (-) DC i Concentration difference between the stream and the pulse of component i (mol m À3 ) e ij NRTL non-randomness parameter between i-j (-) m i Chemical potential energy of i (J mol À1 K À1 ) G Thermodynamic correction factor (-) g i Activity coefficient of i (-) t UNIQUAC i-j interaction parameter (-)