Impact of Values of Diffusion Coefficient on Results of Diffusion Modelling Driven by Chemical Potential Gradient

were


Introduction
Finding relations for various phenomena and describing them in terms of mathematical equations in order to use them later to predict what could happen in a similar situation is a wellestablished approach in science.As an example can serve the problem of forecasting weather or modelling the spread of disease, both with their own journals i.e.Weather and Forecasting published by American Meteorological Society and Infectious Disease Modelling published by Elsevier B.V. where such models are shared.But these are just two of many very complex problems that need some assumptions, drawn boundaries or simplifications to be solvable.Every numerical solution bears the risk of an error and every simplification can make that error even bigger or narrow the use of the model [1].But using such shortcuts is not a bad thing.It is often a matter of being able to acquire some results or none at all by the means of modelling.
Simulation of diffusion process makes a very good example in the case: some results or none.It could be and was done with the use of empirical Fick's laws, where chemical composition gradient was the driving force of mass transport.For example, Fick's First Law is exact for steady state diffusion in isothermal conditions.If some process conditions differ only a little from those restrictions then this approach still can be used to obtain modelling results that agree with the experiment within some predefined margins.However, the more conditions vary from those assumed in Fick's law, the bigger deviation of modelling results from the experimental ones may become.To ignore the limits of Fick's law, a non-simplified approach would have to be used, where mass transport is related to the gradient of chemical potential, as it is the realcorelated with the definition of equilibriumdriving force of diffusion.The ability to calculate chemical potentials, not just by obtaining the limited number of results but also by performing many calculations for everchanging conditions in the system, was the limiting factor that prevented the use of this method.Therefore, it was preferable to obtain results with the use of methods based on chemical composition gradient, which was more or less accurate than no results at all.Nowadays the Calphad approach can be used to obtain such a big number of calculated values of chemical potentials and enables applying of methods based on chemical potential gradient.However, there is still research being done, where composition gradients are employed even in ternary systems, for example [2].With the occurrence of the uphill diffusion in the system, considered in the mentioned paper, a mathematical trick that would change the sign of diffusion coefficient was used.This shows the supremacy of methods based on chemical potential gradients, where diffusion against composition gradient can be acquired in a straightforward way.
The presented model of diffusion, in the aspect of the driving force, is not a simplified one.It calculates fluxes using chemical potential-based equations coupled with the Calphad method to obtain up-to-date values of chemical potentials as the composition changes throughout the system during modelling time.The first results obtained using the model, that have been presented earlier [3] employed other simplification: the values of mobility and diffusion coefficient were considered constant.In this paper, the comparison of results for the constant value of diffusion coefficients versus those dependent on temperature, chemical composition or both will be assessed.

Model description
The model is based on a phenomenological flux equation that relates mass transfer to the gradient of chemical potential [4,5].
Starting with the divergence theorem and applying the Finite Difference Method (FDM) one can derive a difference equation for the change of mass for each element in the system.Such derivation was presented in previous work [3].The final form of the equation for finite difference calculation scheme, defining change of mass for i-th element in considered cell (Δm0,i), for one timestep, is given by: where Bdir,i is the mean mobility, of i-th element, in the area between two neighbouring domains (dir, 0), Δx, Δτstep in space and time, ptotal number of neighbours (p = 2 for 1D calculations).Indexes used in equation ( 2) have the form: (direction, element).Figure 1 presents exemplary values of indexes.
Fig. 1. Outline of FDM cell currently considered for calculations (red) together with two neighbours in 1D model.Cells contain possible values of indexes -direction "dir", element "i"and exemplary indications for composition (within the dashed frame).
Integer values of dir are introduced to sum in equation ( 2) Knowing the Nernst-Einstein relation that connects diffusion coefficient (Di) and mobility [6][7][8], together with Tabsolute temperature and Rgas constant: allows calculating Bi with the use of literature data for Di.The only variable that still needs to be determined, to calculate the change in mass, is µi.
As stated in [3], a very convenient way of acquiring such data is using the Calphad method.If the relation for Gibbs free energy (since in most cases this fundamental equation is available [9]) of a phase for varying values of composition and temperature is known, µ can be calculated from the minimization of the total Gibbs energy [10].Utilising Calphad-based software allows to calculate equilibrium and determine values of chemical potentials whether or not one has access to Gibbs free energy relations and is even more straightforward.In the presented work Thermo-Calc software (version 2019a, thermodynamic database: TCFE7) with TQ-Interface was used to obtain values of chemical potentials of each element for given conditionstemperature, composition (and constant pressure 101325 Pa).

Relations for diffusion coefficients and the choice of calculation parameters
It was pre-defined that the paper will present results of single dimension modelling for Fe-based binary systems.It was also decided to simulate diffusion in a single FCC phase.Literature data for diffusion coefficients was the limiting factor for the choose of systems and calculation parameters.By referring to phase diagrams (calculated using Thermo-Calc version 2019a, thermodynamic database: TCBIN) and analysing data for Di, the following systems were selected: Fe-Si, Fe-Mn and Fe-C.
From the available relations of diffusion coefficients as a function of T, ci or both the final calculation parameters, defining modelling tasks, for each system were chosen.If Di(T) relation was available then composition across modelled specimen was uniform and there was a linear temperature profile with Tmin on one edge of the specimen and Tmax on the other.If Di(ci) was found then the composition profile was linear with cmin and cmax on edges whereas temperature across the specimen was uniform.For the case where Di(ci,T) relation was available linear profiles of ci and T with reversed slopes were used: cmin, Tmax on one edge and cmax, Tmin on another edge.
The number of space grid points was equal for each calculated system.The aim was to obtain a similar relative change of composition, in the same modelled time.For that reason, the total length of the specimen varied to reflect different values of mobilitythe lower mobility the shorter specimen's length was simulated.However, for Fe-Si system it would require to use of extreme temperature gradients and though modelling allows to choose such input, it was decided to use more realistic conditions.

Fe-Si system
Figure 2 presents a fragment of Fe-Si phase diagram that covers the whole FCC area.It shows that data for diffusion coefficients should be searched roughly for temperature from the range 1000÷1300°C and cSi for up to 1.5 wt.%.For this system parameters for the Arrhenius relation [7], given by equation ( 4), were found.

Fig. 2. Fragment of Fe-Si phase diagram covering FCC phase range
Used values of frequency factor D0 and activation energy Q in this system, within considered T and ci limits: For DSi: D0 = 0.07 cm 2 •s -1 , Q = 243.0kJ•mol -1 [12].
The above allowed to examine how adding temperature dependency to diffusion coefficients influences modelled diffusion.The following initial conditions for modelling tasks were chosen: cSi = 0.8 wt.%, the temperature range for D(T) 1150÷1250°C, the single temperature used to calculate D const.T = 1200°C.The total modelled length of the specimen equals 10 mm.

Fe-Mn system
The large region of single FCC phase in Fe-Mn system, for up to 50 wt.% Mn, is presented in figure 3. It spreads throughout a huge range of temperatures and compositions.
Relations for individual diffusion coefficients, for this system, as a function of cMn (at.%) were found in [13].Dependencies were digitalized, recalculated with respect to weight percent and approximated by polynomials.Both curves were divided into two segments for better fitting.Figure 4 presents digitalized points with polynomial approximation curves.Table 1 contains parameters of polynomials, which were used to describe relations for DFe(cMn), DMn(cMn) utilised in one of the calculation tasks.Another calculation task, for single, constant values of diffusion coefficients used the data from another literature source [11].Relation for chemical diffusion coefficient, D , in Fe-Mn system (for 5 at.% = 4.923 wt.% Mn, T = 1170°C) was utilised.Arrhenius equation parameters for this case are the following: D in Fe-Mn (5 at.%): The chemical diffusion coefficient is related to the individual diffusion coefficient of species A and B in binary alloy [8] as: where NA, NB are molar fractions of two components.DFe was calculated using Arrhenius equation parameters for self diffusion of Fe, given earlier for Fe-Si system, for temperature 1170°C.
Rewriting equation ( 5) gives a relation for DMn: ( ) Based on the presented data, the influence of adding composition dependency to the diffusion coefficient on the modelled diffusion was tested.The following initial conditions for modelling tasks were chosen: composition range 0÷10 wt.% Mn, temperature 1170°C, single Mn composition used to calculate D const: cMn = 5.0 at.%.Total modelled length of specimen: 1.0 mm.

Fe-C system
Figure 5 presents a fragment of Fe-C phase diagram containing the whole FCC single phase range.Concentration and temperature dependencies of DC in FCC iron for temperature range 800-1100°C and carbon concentration 0.0÷1.4wt.% was found in [14].The general form of this relation is given by equation (7) where A(T) and B(T) are polynomials of third degree, given by equations (8)(9).Table 2 contains values of polynomials' parameters.For DFe the same temperature dependency as for Fe-Si system was used.Concentration dependency on DFe was not taken into account.

Fe-Si system
Table 3 contains input values defining two calculation tasks: for constant diffusion coefficients (A) and temperature-dependent one (B).Figure 6 presents the initial distribution of Si and a comparison between the modelling results after 25, 50 and 100 days for tasks A & B. Change in concentration is very low and can be observed almost only at the edges of the specimen.For that reason figure 6 presents enlarged, 1.5 mm long, segments on the margins.For task A the cSi curves are close to being symmetricalafter 100 days cSi loss at x = 0 equals 128.9E-4 wt.% whereas at x = 10 mm cSi gain equals to 130.3e-4 wt.%.For task B diffusion at colder edge (x = 0) is slower than for hotter one (x = 10 mm).The shape of the curves suggests that after 100 days system is still far from equilibrium and only the initial stage of diffusion has been calculated.
Within modelled 100 days absolute difference in cSi at the edges between tasks A & B, calculated using equation (10), continuously increases as modelling progresses in time.Relative change at edges calculated using equation (11), with task A being reference tents to 100%.for xmin or xmax (11) It is noteworthy that when considering the time needed for the system to reach a certain point (for example some fixed composition on the edge or at a distance from it) the differences in presented simulations are considerable.Figure 6 gives a good example showing that curve (b) for task A almost overlaps curve (c) for task B, meaning that the system needs nearly twice much time to reach the same state when employing a different approach to diffusion coefficients

Fe-Mn system
Table 5 contains input values defining two calculation tasks: for constant diffusion coefficient (A) and concentration-dependent one (B).Figure 7 presents the initial distribution of Mn and a comparison between modelling results after 25, 50 and 100 days for tasks A & B.
The86iffererence between tasks A and B at x = 0 mm is close to zero whereas on the other edge (x = 1.0 mm) this difference is much bigger (factor of 10 for t = 100 days) with faster diffusion for task A. Close to the middle of specimen we observe higher values of cMn for task A.
Table 7 presents the initial values of diffusion coefficients and the absolute and relative difference in manganese concentration for t = 25, 50, 100 days.In the presented system diffusion coefficients for task B change in timethey are being recalculated for updated cMn values during modelling.At x = 0 initially cMn is higher for task A, between t = 50 and 100 days however due to increase of Di, manganese concentration for task B becomes greater: ΔcMn, rel crosses 100%, probably to reach a maximum but at t = 100 days it still slowly increases.At x = 1.0 mm there is also an increase of Di for task B but DMn stays lower than the coefficient for task A and ΔcMn, rel tends to 100%.7. Initial values of diffusion coefficients at the edges (task A: 5 at.% Mn, data from ref. [11], task B: extreme Mn compositions, data from ref. [13]) and comparison of modelled Mn composition after 25, 50 and 100 days t = 25 days t = 50 days t = 100 days  8. Initial values of diffusion coefficients at the edges (mean composition for task A and extreme cC for task B, data from ref. [11,14]) and comparison of modelled C composition after 25, 50 and 100 days t = 25 days t = 50 days t = 100 days The difference between tasks A and B for all three presented modelling times, throughout the entire specimen length, is very low.At the edges, there is 0.25 wt.% difference between cC, used to calculate DC, for two tasks and this causes 7% difference in the carbon diffusion coefficient.Since the coefficient for task B is being recalculated for updated values of carbon composition DC(cC) tends to values of task A. For t = 25 days there is still 1.4% ΔcC, rel at x = 0. Table 8 presents the initial values of diffusion coefficients and the absolute and relative difference in carbon concentration, according to equations (10) and (11), for t = 25, 50, 100 days.

Constant Di versus Di(T)
Table 9 contains input values defining two calculation tasks: for constant diffusion coefficients (A) and temperature-dependent Di (B). Figure 9 presents the initial distribution of C and a comparison between modelling results after 25, 50 and 100 days for tasks A & B. Similarly to Fe-Si system cC is lower for task B in the middle of the specimen and higher on the edges.Diffusion for temperature-dependent coefficients is again faster for higher temperatures and slower for lower.The difference between tasks A and B, both absolute and relative, is the biggest (for presented data) for t = 25 days.Table 10 contains initial values of diffusion coefficients and the absolute and relative difference in carbon concentration, according to equations ( 10) and (11), for t = 25, 50, 100 days.At the edges, there is 50°C difference in temperature used to calculate Di for two tasks.it causes a much bigger difference in Di then in the previous simulation and consequently greater differences in modelled cC between two tasks.Since, as presented in 0, carbon diffuses to the hotter edge, having the initial composition and temperature gradient in the same direction would mean that system is already quite close to equilibrium.To obtain bigger values of change in composition, and simultaneously have composition impact on diffusion coefficient (rejection of initial uniform composition) the composition gradient has been reversedthe higher value of cC is located at edge x = 0 whereas lower at x = 50 mm.
Table 11 contains input values defining two calculation tasks: for constant diffusion coefficients (A) and DC(cC,T) with DFe(T) (B). Figure 10 presents the initial distribution of C and a comparison between modelling results after 25, 50 and 100 days for tasks A & B. Table 10.Values of Diffusion coefficients at the edges (mean temperature for task A and extreme temperatures for task B, data from ref [11,14]) and comparison of modelled C composition after 25, 50 and 100 days t = 25 days t = 50 days t = 100 days The fastest diffusion closer to the edges causes the composition profile to have two extremes for significant simulation time (curves b and c in fig.10).Initially, a bigger relative difference between tasks A and B can be observed for a hotter edge, x = 50 mm: for t = 25 days deviation from 100% is still higher there.In later stages, however, composition equalises faster on that edge due to being closer to equilibrium and the value of DC being even higher than initially, because of its composition dependency.Table 12 presents initial values of diffusion coefficients and the absolute and relative difference in carbon concentration, according to equations ( 10) and ( 11), for t = 25, 50, and 100 days.

Conclusions
At the beginning of each simulation, the system is the furthest from equilibrium and the driving force of diffusion (gradient of chemical potentials) is at its maximum.For that reason, the biggest difference in modelled composition, for the varying approach to the calculation of diffusion coefficients, can be found in these early stages.Since each system tends to equilibrium at its own pace, the duration of that period should be identified with the amount of mass already transferred rather than elapsed time.
For simulations that are concluded during the early diffusion (example of Fe-Si system), taking into account temperature or composition dependency on diffusion coefficient is strongly recommended.In later stages, when the flux is getting smaller as the driving force tends to zero, the difference between the two modelling tasks slowly vanishes.For simulations that would be performed well beyond the initial stage, when the system approaches equilibrium (simulated 100 days of Fe-C diffusion), the use of the constant diffusivity simplification in the model is desirable.For cases in between (examples of Fe-Mn and modelled 25 or 50 days of Fe-C systems), the decision about the way of determining Di should be made according to the current requirements: simplicity or accuracy of the model.The results show that the smaller the difference between constant and varying diffusion coefficient is the faster the system reaches the point when the distinction in results between two approaches becomes

6 .
Concentration profile cSi(x) for the whole length of specimen (I), for 1.5 mm on the left edge (II) and for 1.5 mm on the right edge (III).Initial silicon concentration (a), and modelled distribution after: 25 days (b), 50 days (c), 100 days (d).Dashed lines: constant values of Ditask A, solid lines: DSi(T), DFe(T)task B

Table 9 .
Parameters defining tasks A and B for Fe-C calculations (Di(T

Table 1 .
Parameters of polynomials used to calculate values of DFe(cMn) and DMn(cMn)

Table 3 .
Parameters defining tasks A and B for Fe-Si calculations

Table 4 .
[11,12]of Diffusion coefficients at the edges (mean temperature for task A and extreme temperatures for task B, data from ref.[11,12]) and comparison of modelled Si composition after 25, 50 and 100 days

Table 11 .
Parameters defining tasks A and B for Fe-C calculations (DC(cC, T))

Table 12 .
[11,14]of diffusion coefficients at the edges (mean temperature and composition for task A and extreme temperature and composition for task B, data from ref[11,14]) and comparison of modelled cC after 25, 50 and 100 days