Mesoscale modelling of miscible and immiscible multicomponent fluids.

A mesoscopic simulation method based on the integration of dissipative particle dynamics (DPD), smoothed particle hydrodynamics (SPH) and computational thermodynamics (CT) has been developed. The kinetic behaviours of miscible and immiscible fluids were investigated. The interaction force between multicomponent mesoscopic particles is derived from the system free energy. The diffusivity of the components in non-ideal solution is determined by the chemical potential. The proposed method provides convincing predictions to the effects of convection, diffusion and microscopic interaction on the non-equilibrium evolution of engineering fluids, and demonstrates a potential to simulate more complicated phenomena in materials processing.

Introduction

Multicomponent fluids contain miscible and/or immiscible elements. The mixing and demixing phenomena are important in chemical engineering and material fabrication[1,2]. Many advanced materials are fabricated by precise control of kinetic evolution. For example, a core-shell nanostructured metallic glass with supreme plasticity was fabricated by well-controlled demixing processing[3], while an aluminium-bismuth alloy with uniformly distributed Bi particle relies on well-controlled mixing procedure to improve its bearing ability[4]. On the contrary, a poorly controlled process degrades the materials properties[5].

The kinetic evolution in mixing and demixing is affected by convection, diffusion and thermodynamics of the materials. The diffusion-convection equation illustrates the coherent interaction among the factors[6]where D is the diffusivity, the convection velocity and R is the sources or sinks of the composition. R is dependent on the reactive thermodynamics. Materials microscopic interaction affects the viscosity and diffusivity[7,8]. Convection influences distribution of compositions, which subsequently affects the phase transformation. The redistribution of composition affects the force field between compositions and hence to affect flow pattern. The kinetic evolution of multicomponent material should be described by an integrated model compromising at least the diffusion, convection and thermodynamics.

Many numerical methods have been developed to tackle this problem under various approximations. Phase field model has been implemented to study the mixing and demixing processes either without convection[9] or with a largely simplified hydrodynamic description of flow behaviour[10]. Finite element method has been integrated with phase-field model to study the micro-segregation in the solidification of binary alloys[11]. The approach faces difficult in description of realistic cases such as that in blast furnace, where the flow is laminar at one location and tubular in another with some complex geometries of the phase boundaries. Lattice Boltzmann equation has demonstrated the capability in description of the mixing and demixing process[12,13], which is applicable to multiphase and multicomponent systems[14,15] but facing challenges in dealing with non-isothermal, high density-ratio and high Reynolds number systems. The dissipative particle dynamics (DPD) and smooth particle hydrodynamics (SPH) are off-lattice hydrodynamic methods[16,17] capable of simulating multiphase fluid with no need of explicit interface tracking[18,19]. However, their capability to the incorporation of engineering thermodynamics needs to be developed further. There is no generic method reported in literature to allow the integration of the mesoscopic simulation method to commercial thermodynamic database that has been implemented widely in engineering.

This work aims to develop an integrated method based on the frameworks of DPD, SPH and computational thermodynamics, and use the method to investigate the effect of convection, diffusion and thermodynamics on the materials kinetic behaviour in mixing and demixing processes. The purpose of the research was to simulate materials processing in engineering cases.

Results

The material is represented approximately by a pack of mesoscopic interactive particles. Each particle contains many atoms and is a subsystem that obeys the statistical thermodynamics. Inheriting the basic idea from DPD method, the interaction between a pair of mesoscopic particles contains the conservative, dissipative and random forces[20]where the sub-index represents the particle i and j respectively. The dissipative force F represents the viscosity. The random force F represents the thermal fluctuation. The fluctuation-dissipation theorem can be satisfied if F and F are correlated adequately[20]. The kinetic behaviour of an ideal solution can be reproduced by Eq. (2) using . represents the interaction between particles. Given the fact that many engineering materials are processed in an open environment under approximately constant pressure, the state thermodynamic quantity in such an environment is Gibbs free energy. The equivalent mesoscopic driving force in kinetic evolution can be derived from the thermodynamic quantity as[21]where G = G−G is the Gibbs free energy difference between two points in the space. Equation (3) means that the mesoscopic driving force intends to drive the particle to move between two different points. This can be considered as two nearby mesoscopic particles i and j. , where is the spatial vector at the mass centre of particle i. , where g(c, T) is the Gibbs free energy density and c = {c1, c2, … c} for the material containing n components. For a material without ferromagnetic element one has[22]where g(c, T) = g(c, T) + g(c, T) is the free energy density of ideal material, is the free energy of bulk pure elements, is the free energy due to ideal mixing and g(c, T) is called the excess free energy density. R is the gas constant. g(T) is the free energy density of pure element α at temperature T, which is available for all the elements[22]. Substituting Eq. (4) into (3) and comparing with (2), one haswhere . Equation (5) has the same format as that of derived by Pagonabarraga and Frenkel[23]. In the computational thermodynamics, the excess free energy density is expressed as Redlich-Kister expansion[24]where is the m-order interactive coefficient between element k and l and is dependent on the temperature. The value of are available for many materials in commercial thermodynamic databases[25]. For a binary fluid, Eq. (6) reduces into . For binary regular solution, it further reduces into Porter equation . The convection of an engineering fluid can be calculated according to the derived mesoscopic interactive force.

The governing equation for solute diffusion was derived by Cahn as  [26], where M is a positive coefficient called mobility and the diffusivity is defined as in comparison with Fick’s law. D ≥ 0 when , which corresponds to solute mixing in miscible fluid. D < 0 when , which corresponds to the case of spinodal decomposition. The more generic case was derived from the irreversible law of thermodynamics where the effect of inter-diffusion on the volume change was taken into account[7]. In the present approximation, materials are represented by moving mesoscopic particles. The chemical composition and temperature in different mesoscopic particles might be different to reflect the heterogeneous non-equilibrium materials. The free energy density and hence the diffusivity in each mesoscopic particle can be in different values. The diffusion between two adjacent particles i and j requires to determine its diffusivity D according to the respective diffusivity D and D. To this purpose one considers a fictitious mid-point position r between i and j positions, the diffusivity from r to r is D and from r to r is D. The net solute flux from both i and j particles to r should be zero according to mass conservative, i.e.Equation (7) gives . The flux from i to j is . This gives

The diffusivity can be in positive or negative values that are dependent on the chemical constitution, composition and temperature. Therefore, it is possible that D + D = 0 in some cases. Equation (8) goes infinite when D = −D. In such case, The mass conservation in Eq. (7) cannot be satisfied if c ≠ c with D ≠ 0 and D ≠ 0. To guarantee the mass conservation, it requires D = 0 and D = 0. This will enable Eq. (8) to be not infinite because DD is a higher order infinite small quantity than D + D. In summary, the diffusivity between adjacent particles i and j isWith Eq. (9) in mind and also to represent the diffusion equation into following format

Numerical solution of Eq. (10) can be achieved using SPH methodology. The latter uses interpolation method and introduces a kernel to calculate the quantity and the associate derivatives. Equation (10) in SPH is represented as[26]where m and ρ are the mass and density of particle i. h is the smooth length and W is the kernel function. There are several possible formats of W available to choose from[17].

In the numerical calculation, one chooses and  [20], where , . ζ is the white noise between −0.5 and 0.5 and with 0 mean. γ and σ are coefficients satisfying σ2 = 2γkT, k and T are the Boltzmann constant and temperature, respectively. ω and ω are weight functions satisfying

The theoretical derivations presented in the earlier paragraphs have demonstrated that the method proposed in the present work is based on an integration of DPD, SPH and CT. The kinetic evolution of the system was described following a DPD strategy, where conservative, dissipative and random forces are implemented to describe the driving force. The relationship between dissipative force and random force are defined according to classical DPD theory. However, the chemical constitution and its evolution are described according to the SPH framework, where interpolation method and kernel function are implemented. The conservative force and the diffusivity are derived and obtained according to the CT theory and the associate format of database. The method developed in the present work is based on integration of DPD, SPH and CT.

If the chemical diffusion is neglected by letting D = 0 and conservative force is defined as with a a constant and ω(r) the same format as that defined in Eq. (12), the newly developed method in the present work reduces to the classical DPD method that was discussed by Groot and Warren[20]. Their results can be reproduced if the same coefficients are chosen because the governing equations for both models in such conditions are identical[20]. The numerical results for particle demixing at various ratios of particle numbers have been demonstrated in Supplementary Figures. On the other side when the coefficients of dissipative and random forces are defined as γ = 0 and σ = 0, the newly developed method will reduce to conventional SPH where kernel interpolation is applied to calculate all the quantities (e.g. mass of particle) and their spatial derivations (e.g. gradient of particle density). The thermodynamic quantity represented in Eq. (6) is widely applied in description of condensed matters such as solid and liquid metallic materials. For gaseous phases, activities and partial pressure can be used. All of those are well described in computational thermodynamics.

To implement the integrated model developed in the present work, one considers a regular binary solution with solute A and solvent B. For a mesoscopic particle containing N atoms of regular binary elements with solute mass composition c, the excess free energy is and diffusivity is . Its mass is obtained by , where m and m are the atom mass of element A and B respectively, N is the Avogadro’s constant. Its density, without considering of the mixture-induced volume change, can be obtained by , where ρand ρ are mass density of pure element A and B respectively.

For total N number of mesoscopic particles in a cubic box, one defines the dimensionless units for the convenience of numerical calculation. The dimensionless mass reduces into , where c is the average composition. The temperature reduces into where T is the melting temperature of a binary mixture with the average composition. The free energy density is reduced to . The smooth length reduces to h* = h/h, where and L, L and L are number of unit cells in each coordinates. Based on the above definition, one has t* = t/τ, where . It gets and .

For the parameters of N = 4000, L = L = L = 5 and c = 0.45, N = 10 gives average mass of mesoscopic particle of m = 1.01 × 10−24 kg and h = 15.35 Å for Cu-Ni liquid miscible alloy and m = 1.945 × 10−24 kg and h = 21.06 Å for Al-Pb liquid immiscible alloy. Cu-Ni alloy is well-known for its unlimited miscibility of one element in another in liquid or solid states. The alloy has been developed for many decades and used as coinage, electrical components and other occasions[27]. Al-Pb is an immiscible binary alloy[28]. The mechanical mixture of Al-Pb is desirable for making wear parts in automotive industry. Both alloys have been studied intensively by experiments. The melting temperature for Cu-Ni is 1573 K and for Al-Pb about 1643 K at the given solute composition of c = 0.45. One chooses σ = 3 and γ = σ2/2 for random force and dissipative force according to the recommendation in literature[20], and selects for miscible alloy and for immiscible alloy. The selection of in the present simulation didn’t follow the actual value of specific materials for the purpose to prove general applicability of the proposed method. This is because is not only materials-related but also temperature-dependent. More importantly, several interactive coefficients are required to determine the non-regular and non-ideal materials. However, the values are in the range of real conditions. The computational thermodynamic methodology is followed strictly in the following studies. The calculation of phase diagram (CALPHAD) based on computational thermodynamics has been proved to agree with experimental characterization of phase diagram in many engineering cases. Giving the fact that the bulk free energy doesn’t usually affect the phase equilibrium, one plotted g(c, T) + g(c, T) in Fig. 1. It shows clearly that corresponds to a miscible fluid at any solute composition but represents an immiscible fluid. Spinodal decomposition happens when the average solute composition is between c and c. To selection of diffusivity is after the reference of true values of realistic system, e.g. the diffusivity of dilute Cu-Ni alloy of 3.27 × 10−9 m2s−1 and that of dilute Al-Pb of 1.0 × 10−6 m2s−1 respectively[9,29], one has M* = 0.0144 and M* = 5.84 for the two systems. Periodical boundary condition was implemented in the simulations.

Figure 1

The change of gid(c, T) + gex(c, T) as a function of solute composition shows that is miscible and is immiscible. ∂2g/∂c2 = 0 reveals the spinodal decomposition area between c and c. For the immiscible case, the cross points between common tangent line and Gibbs free energy curve gives the equilibrium compositions and . The cross point between the comment tangent line and the average chemical composition determines the minimum Gibbs free energy.

The histogram of particle concentration distribution for miscible system is demonstrated in Fig. 2. At the beginning, an average solute composition of 0.45 with a range of variation between 0.4 and 0.5 was assigned to the mesoscopic particles randomly. The solute composition axis from 0 to 1 was divided into 100 parts and the mesoscopic particles were counted within each solute part by a statistical code program. For example, those particles with solute composition between 0.4 and 0.41 were considered as the same solute composition of 0.405 during the plotting of Fig. 2. As the time evolves, more and more particles have their solute composition to approach to the average value. The variation range decreases and the distribution peak becomes sharp, which indicates a gradual mixing of solute in the solvent. After 105 time steps, a single peak in the histogram was nearly formed. There are 2970 particles with composition between 0.44 and 0.45, and 1030 particles with composition between 0.45 and 0.46. This suggests that the mixing was almost reached and solute were distributed almost homogeneously.

Figure 2

Distribution of particle numbers with different solute compositions at 0, 1 × 104 and 1 × 105 time steps for miscible fluid.

The 3D configurations of particles at 0 and 1 × 104 time steps were plotted respectively using MatVisual, as are presented in Fig. 3(a,b). The colour of the particles is determined by the solute composition in a range between 0.4 and 0.5, with 0.4 in blue and 0.5 in red at the colour spectrum. The colours were distributed randomly in Fig. 3(a), which demonstrates the random initialization of the solute composition to mesoscopic particles. The colours are converged toward that in the middle of spectrum after 1 × 104 time steps, as shown in Fig. 3(b). The other time steps demonstrated in Fig. 2 were also plotted but not presented here because the Fig. 3(b) has already shows the tendency of solute homogenization clearly.

Figure 3

Configurations of mesoscopic particles at (a) 0 and (b) 1 × 105 time steps. The colour spectrum corresponds to the solute concentration from 0.4 to 0.5.

The histogram for the distribution of particle numbers with various solute concentrations for immiscible fluids at 0, 100, 1000 and 1 × 105 time steps was presented in Fig. 4(a). The initial solute distribution was assigned following the same initialization procedure as that of in Fig. 2, i.e. the solute was randomly distributed around the average value of 0.45 with a range of slight fluctuation altitude of 0.05. In the time evolution, the range of solute variation becomes wider firstly and then two distinct peaks appear after 1000 steps. The two distinct peaks represent a solute-rich and another solute-poor liquid phases respectively. The solute compositions at the peaks are in agreement with the minimum Gibbs free energy positions indicated in Fig. 1, which is the cross point between the comment tangent line and the average chemical composition. It can be seen clearly from Fig. 4(a) that the two peaks at 105 time steps are not in the same height. This is due to the initial average solute composition is 0.45 rather than 0.5 and the equilibrium solute composition positions is in symmetrical position regarding to 0.5 rather than 0.45. To prove that the level rule is valid in the method developed in the present work, another case with identical simulating conditions as that in Fig. 4(a) apart from the average composition value of 0.5 has been studied. The solute concentrations at 0, 100, 1000 and 1 × 105 time steps was plotted in Fig. 4(b). The heights of two peaks at 1 × 105 time steps are the same within fluctuation. The distribution of the particle numbers obeys the level rule in equilibrium phase diagram. The presence of two peaks in concentration distribution indicates the happening of solute demixing. This is in agreement with the thermodynamic prediction indicates in Fig. 1.

Figure 4

Distribution of particle numbers with various solute compositions at 0, 100, 1000 and 1 × 105 time steps for immiscible fluids. (a) The average solute composition is 0.45. (b) The average solute composition is 0.5.

The statistical calculation indicates some particles are with the solute composition between that of two equilibrium values, or called intermediate solute composition. The number of particles with intermediate solute concentration reduces drastically during the time evolution. Figure 5 presents the fraction of such particles at various time steps. In the beginning, the number of particles with intermediate solute composition reduces in a scaling law with a scale index of −0.47 according to the trend line fixing. It can be seen that that a large reduction of intermediate phase occurs at early stage and the reduction rate decreases drastically after the two distinct peaks appear. Some particles may be trapped at the interface between the solute-rich and solute-poor phases, and take longer time to evolve toward the peak area along with the phase coarsening and domain coalescence. At 105 time steps the ratio of intermediate phase is less than 1%. The intermediate particle can be seen clearly in the 3D particle distribution figures shown in Fig. 6, where the solute compositions were indicated by particles’ colours from blue for 0.25 to red for 0.75. It can be seen clearly that intermediate particles are located in the interface area.

Figure 5

Evolution of the fraction of particles with intermediate solute composition. Trendline-fitting reveals scaling factor around −0.47.

Figure 6

Particle solute distributions at 0, 103, 104, and 105 time steps. The colour spectrums are from blue (0.25) to red (0.75).

Discussion

In order to understand how the initial speed of mesoscopic particles and substance diffusivity affect fluid kinetic behaviours, several pseudo simulations are performed for the immiscible fluids. For the low speed case, the initialization of particle speed components in 3 Cartesian coordinates were assigned randomly between −0.5 and 0.5 speed unit with average macroscopic speed of zero. In the high initial speed simulations, the variation was 1000 times larger as that of the low speed case. The solute configurations among particles at 100, 200 and 1000 time steps are presented in Fig. 7. Larger speed causes longer evolution time toward the equilibrium. This is understandable because larger fluctuation requests longer time to dissipate its kinetic energy[30]. However, it has been noticed that after 1000 time steps, the solute distributions in particles are not distinguishable for the two speed initialization cases. This is due to the definition of dissipative force in DPD method. Larger fluctuation of speed causes larger dissipative force, as is shown in for larger relative speed . The initial speed condition will be smeared out. To keep the high speed, driving force will be required to the system, which is beyond the scope of the present work.

Figure 7

Effect of different initial random speed on the solute demixing at 100, 200 and 1000 time steps.

The effect of diffusivity on the system evolution has also been investigated. Figure 8(a) demonstrated the solute distribution for a low experimental diffusivity and another high speculate diffusivity (10 times of the true diffusivity) at 100 and 1000 time steps. It can be seen that diffusivity makes huge differences. Large diffusivity leads to a faster evolution toward thermodynamic equilibrium. Figure 8(b,c) demonstrate the comparison of effects of initial speed and diffusivity on the solute decomposition at 100 and 1000 time steps.

Figure 8

(a) Effect of diffusivity on solute demixing at 100 and 1000 steps. Comparison of the effects of initial speed and diffusivity on solute decomposition at (b) 100 and (c) 1000 time steps.

Figure 9 plotted the particles solute distribution at 20000 time steps at four different conditions, namely the high speed with low diffusivity, high speed with high diffusivity, low speed with low diffusivity and low speed with high diffusivity. In the stage after many time steps of evolution, the microstructure evolution is not mainly due to diffusion but due to coalescence. The effect of diffusivity is not demonstrated clearly. Figure 9(a,b) are with high particle random speed and (c) and (d) with low speed. (a) and (c) have low diffusivity but (b) and (d) have high diffusivity. The high random speed helps to destroy the structures. This is similar to that of the effect of temperature on the fluid[31].

Figure 9

(a) High speed low diffusivity; (b) high speed high diffusivity; (c) low speed low diffusivity and (d) low speed high diffusivity. The spectrum is between 0.25 to 0.75.

In summary, a generic model that integrated the hydrodynamics, diffusion and thermodynamics has been developed. The method implemented the spirit of DPD and SPH and linked the theoretical frame into thermodynamics. Using the model, many engineering materials fabrication can be simulated. The parameters are determined by the real system. The results are in agreement with the phase equilibrium prediction made by computational thermodynamics. Their respective structural evolutions can also be shown in 3D space particle plots and some fine structure can be captured during the structure evolution. The method developed in the present work has some advantage and disadvantage aspects in comparison with other methods. It is better than lattice Boltzmann method in simulation of adiabatic condition and high Reynolds number fluids but worse than that in simulation of complex geometric boundaries[8,13]. It provides less accuracy than that of the phase field method in simulation of spinodal decomposition by diffusion but can deal with mixing and demixing in convection, in which phase-field method needs to integrate with another hydrodynamic method to achieve the task[10,32]. In comparison with conventional computational fluid dynamics, the method developed in the present work can provide information on microstructure evolution.

An understanding of how fluid speed and substance diffusivity affect phase separation/mixing as a case research has been given as well. The flow speed does not influence the concentration exchanges very much and then only have a subtle influence on the speed of phase separation/mixing while the substance diffusivity influences concentration exchanges a lot and have a big influence on the speed of phase separation/mixing. In terms of structural evolution, the substance diffusivity does not influence the structure evolution while the high speed flow violates forming fine structure at early stage and contributes to forming loose structure at late stage.

Several further works need to be done. An investigation of the influence of phase behaviours on fluid flows is still left out. A spatial temperature variation can be added to the model to accommodate temperature effects on fluid flow with phase behaviours. Moreover, as the framework of this model is generic, by replacing the thermodynamic law of phase behaviours with chemical reaction laws, a promising model for reactive fluids can be potentially developed, which will be very useful and meaningful in engineering applications.

Methods

Numerical algorithm

In combination of the modified velocity-Verlet numerical algorithm that was used in most DPD simulation[20] and the leap-frog numerical algorithm that was applied in SPH simulation[33], following numerical algorithm was implemented in the present work, where the iterations of velocity and solute concentration were divided into two half steps separated by the acceleration iteration.where is the acceleration of particle i and Δt is the time step.

Software

The three-dimensional configurations of mesoscopic particles with various solute concentration were plotted using MatVisual software.

Supplementary Figure