A Symmetric Free Energy Based Multi-Component Lattice Boltzmann Method

# A Symmetric Free Energy Based Multi-Component Lattice Boltzmann Method

## Abstract

We present a lattice Boltzmann algorithm based on an underlying free energy that allows the simulation of the dynamics of a multicomponent system with an arbitrary number of components. The thermodynamic properties, such as the chemical potential of each component and the pressure of the overall system, are incorporated in the model. We derived a symmetrical convection diffusion equation for each component as well as the Navier Stokes equation and continuity equation for the overall system. The algorithm was verified through simulations of binary and ternary systems. The equilibrium concentrations of components of binary and ternary systems simulated with our algorithm agree well with theoretical expectations.

## I Introduction

Multicomponent systems are of great theoretical and practical importance. An example of an important ternary system, that inspired the current paper, is the formation of polymer membranes through immersion precipitation Akthakul et al. (2005). In this process a polymer-solvent mixture is brought in contact with a non-solvent. As the non-solvent diffuses into the mixture, the mixture phase-separates, leaving behind a complex polymer morphology which depends strongly on the processing conditions. The dependence of the morphology on the parameters of the system is as yet poorly understood. Preliminary lattice Boltzmann simulations of this system exist Akthakul et al. (2005). However, this work did not recover the correct drift diffusion equation. A general fully consistent lattice Boltzmann algorithm with an underlying free energy to simulate multicomponent systems is still lacking. This paper strives to bring us a step nearer to achieving this goal.

There are several previous lattice Boltzmann methods for the simulation of multi-component systems. There are three main roots for these approaches. There are those derived from the Rothmann-Keller approach (2); (3) that attempt to maximally phase-separate the different components. A second approach by Shan and Chen is based on mimicking the microscopic interactions (4); (5); Shan and Doolen (1996) and a third approach after Swift, Orlandini and Yeomans Osborn et al. (1995); (8) is based on an underlying free energy. All of these have different challenges. Since we are interested in the thermodynamics of phase-separation we find it convenient to work with a method based on a free energy. This allows us to easily identify the chemical potentials of the components. This is convenient since the gradients of the chemical potentials drive the phase separation as well as the subsequent phase-ordering.

The challenge for the LB simulation of a multicomponent system lies in the fact that momentum conservation is only valid for the overall system but not for each component separately, and diffusion occurs in the components. For a binary system of components and with densities and , the simulation usually traces the evolution of the total density and the density difference (8). Although this scheme is successful in the simulation of a binary system (8); Wagner and Cates (2001), its generalization for the LB simulations of systems with an arbitrary number of components is asymmetric. For instance, to simulate a ternary system of components , , and with densities , and , the total density of the system, , should be traced, and the other two densities to be traced may be chosen as, e.g., and Lamura et al. (1999). This approach is likely to be asymmetric because the three components are treated differently as is the case of Lamura’s model Lamura et al. (1999). If an LB method is not symmetric, it will lose generality an will only be adequate for special applications. In this paper, we established a multicomponent lattice Boltzmann method based on an underlying free energy that is manifestly symmetric.

## Ii Macroscopic Equations for Multicomponent System

The equation of motion for a multicomponent system are given by the continuity and Navier-Stokes equations for the overall system and a drift diffusion equation for each component separately. The continuity equation is given by

 ∂tρ+∇⋅J=0, (1)

where is the mass density of the fluid, is the mass flux which is given by , and is the macroscopic velocity of the fluid. The Navier-Stokes equation describes the conservation of momentum:

 ∂t(ρuα)+∂β(ρuαuβ)=−∂βPαβ+∂βσαβ+ρFα, (2)

where and are the pressure and viscous stress tensors respectively, is the component of an external force on a unit mass in a unit volume, and the Einstein summation convention is used. For Newtonian fluids, the viscous stress tensor is given by

 σαβ=η(∂βuα+∂αuβ−d2δαβ∇⋅u)+μBδαβ∇⋅u, (3)

where is the shear viscosity, and is bulk viscosity, and is the spacial dimension of the system.

Free energy, chemical potential, and pressure are key thermodynamic concepts to understand the phase behavior of a system. The chemical potential of each component can be obtained by a functional derivative as

 μσ=δFδnσ, (4)

where is the chemical potential of component ; is the number density of component ; and is the total free energy of the system.

The pressure in a bulk phase in equilibrium is given by

 p=∑σnσμσ−ψ. (5)

The pressure tensor is determined by two constraints: in the bulk and everywhere.

In multicomponent systems, there are two mechanisms for mass transport: convection and diffusion. Convection is the flow of the overall fluid, while diffusion occurs where the average velocities of components are different. The velocity of the overall fluid is a macroscopic quantity because it is conserved, but the average velocities of the components are not. The macroscopic velocity of the fluid can be expressed in terms of the density and velocity of each component in the form of

 u≡∑σρσuσ∑σρσ. (6)

With the notation

 Δuσ≡uσ−u, (7)

the flux of each component can be divided into a convection part and a diffusion part :

 Jσ≡ρσuσ=ρσ(u+Δuσ)=Jσc+Jσd. (8)

Because mass conservation still holds for each component, the continuity equation for each component is valid:

 ∂tρσ+∇⋅Jσ=0. (9)

Substituting Eq. (8) into Eq. (9), the convection diffusion equation for a component can be obtained.

 ∂tρσ+∇⋅Jσc=−∇⋅Jσd. (10)

From Eqs. (6) and (7), we see that

 ∑σJσd=0, (11)

which ensures the recovery of the continuity equation for the overall system. The diffusion process between two components is related to the difference of the chemical potential of the two components, which is also called the exchange chemical potential Jones (2002). Recognizing that the gradient of the exchange chemical potential determines the diffusion processes, we obtain a first order approximation for the diffusion flux of one component into all other components as

 Jσd=−∑σ′Mσσ′∇(μσ−μσ′), (12)

where and enumerate the components; and are the chemical potentials of components and ; and is a symmetric positive definite mobility tensor.

A simple model for the diffusion process assumes that a diffusion flux between two components is proportional to the overall density and the concentration of each component. Then mobility tensor can be expressed as

 Mσσ′=kσσ′ρσρσ′ρ, (13)

where is the constant diffusion coefficient between components and . It depends on components but is independent of the total densities and concentration of each component. Substituting Eq. (13) into Eq. (12), we have

 Jσd=−∑σ′kσσ′ρσρσ′ρ∇(μσ−μσ′). (14)

Substituting Eq. (14) into Eq. (10), the general form of a convection diffusion equation is obtained as

 ∂tρσ+∇(ρσu)=∇∑σ′kσσ′ρσρσ′ρ∇(μσ−μσ′). (15)

## Iii Lattice Boltzmann for Multicomponent System

To simulate a multicomponent fluid using LB we set up a LB equation for each component. The LBE for a component of a multicomponent system is given by

 fσi(r+viΔt,t+Δt)−fσi(r,t) (16) = Δt[1τ(fσei(r,t)−fσi(r,t))+Fσi],

where is the particle distribution function with velocity for component , is its equilibrium distribution and is the forcing term of component due to the mean potential field generated by the interaction of the component with the other components. The main task in setting up this lattice Boltzmann method is to determine the correct form of the forcing term which will recover the convection diffusion equation (15).

The density of each component and the total density are given by

 ρσ = ∑ifσi, (17) ρ = ∑σρσ. (18)

The average velocity of one component and the overall fluid can be defined as

 ρσuσα ≡ ∑ifσviα, (19) ρuα ≡ ∑σρσuσα, (20)

where is the average velocity of the component , and is the average velocity of the overall fluid.

The moments of equilibrium distributions for one component are chosen to be

 ∑ifσei = ρσ, ∑ifσeiviα = ρσuα, ∑ifσeviαviβ = ρσ3δαβ+ρσuαuβ, ∑ifσeviαviβviγ = ρσ3(uαδαβ+uβδαγ+uγδαβ) (21)

The moments for the forcing terms of one component are

 ∑iFσi = 0 (22) ∑iFσiviα = ρσaσα, (23) ∑iFσiviαviβ = ρσ(aσαuσβ+aσβuσα), (24) ∑iFσiviαviβviγ = 13ρσ(aσαδβγ+aσβδαβ+aσγδαβ). (25)

To utilize the analysis of the one component system we can establish a LB equation for the total density by defining

 ∑σfσi = fi, ∑σFσi = Fi, ∑σρσaσα = ρaα. (26)

Similar to the counterparts of the one-component system, the moments for the overall equilibrium distribution function are given by

 ∑ifei = ρ, ∑ifeiviα = ρuα, ∑ifeiviαviβ = 13ρδαβ+ρuαuβ, ∑ifeiviαviβviγ = 13ρ(uαδβγ+uβδαγ+uγδαβ) (27) +ρuαuβuγ+Qαβγ.

The moments for the overall force terms are then given by

 ∑iFi = 0, ∑iFiviα = ρaα,

Using Eq. (24), we obtain

 ∑iFiviαviβ (28) = ∑σ(aσαuσβ+aβuσα) = ρ(aαuβ+aβuα)+∑σ(aσαΔuσβ+aσβΔuσα),

where the second term of Eq. (28) is of a higher order smallness than the first terms, and therefore does not enter the hydrodynamic equations to second order. For the third moment we have

 ∑iFiviαviβviγ=13ρ(aαδβγ+aβδαβ+aγδαβ). (29)

By summing Eq. (16) over , an effective LB equation for the total density is

 fi(r+viΔt,t+Δt)−fi(r,t) (30) = Δt[1τ(fei(r,t)−fi(r,t))+Fi],

This is identical to the LB equation for a system of one component. Therefore, the continuity equation and the Navier Stokes equation of the overall fluid of a multicomponent system are recovered as

 ∂tρ+∂α(ρUα)=0+O(ϵ3), (31)

where is the macroscopic velocity of the fluid. The Navier Stokes equation for the overall fluid is:

 ∂t(ρUβ)+∂α(ρUαUβ) (32) = −∂α(13ρδαβ) +∂α(w3ρ(∂αUβ+∂βUα))+ρaβ −w∂γ∂αρuαuβuγ.

To recover the convection diffusion equation of each component, we performed a Taylor expansion on the left of Eq. (16) to second order:

 Δt(∂t+viα∂α)fσi+(Δt)22(∂t+viα∂α)2fσi+O(ϵ3) (33) = Δt(1τ(fσei−fσi)+Fσi).

Because of the recursive nature of Eq. (33), can be expressed by and derivatives of as

 fσi = fσei+τFσi (34) −τ(∂t+viα∂α)(fσei+τFσi)+O(ϵ2).

Substituting Eq. (34) into the left side of Eq. (33) we obtain

 (∂t+viα∂α)(fσei+τFσi) (35) − w(∂t+viα∂α)2(fσei+τFσi)+O(ϵ3). = 1τ(fσei+τFσi−fσi).

Summing Eq. (35) over gives,

 ∂tρσ+∂α(ρσuα)+τ∂α(ρσaσα) (36) − w∑i(∂t+viα∂α)2(fσei+τFσi)=O(ϵ3)

The first moment of and are not identical, and the continuity equation cannot be obtained. Eq. (36) shows that is of order , and is of order . Therefore is of order , and we get

 w∑i(∂t+viα∂α)2(fσei+τFσi) = w∂β[∂t(ρσuβ)+∂α(ρσ3+ρσuαuβ)]+O(ϵ3).

So Eq. (36) can be simplified to

 ∂tρσ+∂α(ρσuα)+τ∂αρσaσα+O(ϵ3) −w∂β[∂t(ρσuβ)+∂α(ρσ3+ρσuαuβ)]=0. (37)

Eqs. (32) yields

 ∂tUβ=−Uα∂αUβ−1ρ∂α(ρ3δαβ)+aβ+O(ϵ2). (38)

From Eq. (36) it follows that

 ∂tρσ=−∂α(ρσUα)+O(ϵ2). (39)

Inserting Eqs. (38) and (39) into Eq. (37) we get

 ∂t(ρσuβ)= − ∂α(ρσUαUβ)−ρσρ∂α(ρ3δαβ) (40) + ρσaβ+O(ϵ2).

Substituting Eq. (40) into Eq. (37) results in

 ∂tρσ+∂α(ρσUα) = ∂α[τρσaα−τρσaσα+w∂α(ρσ3) (41) −wρσρ∂α(ρ3)],

From this we deduce that the correct form of the forcing term is

 ρσaσα ≡ −wτρσ∂α(μσ−13lnρσ) (42) = −wτ(ρσ∂αμσ−13∂αρσ),

where the coefficient is . This coefficient approaches 1 as approaches 0, as one would expect from the continuum limit. This constitutes the main result of this paper. Plugging Eq. (42) into Eq. (41), we then obtain the convection diffusion equation

 ∂tρσ+∂α(ρσUα) = ∂α[w∑σ′ρσρσ′ρ∂α(μσ−μσ′)]. (43)

The diffusion flux of component is

 Jσd=−w∑σ′ρσρσ′ρ∂α(μσ−μσ′). (44)

So that the in Eq. (44) is equivalent to in Eq. (14).

## Iv Numerical validation

We examined the equilibrium behavior of phase separated binary and ternary systems. We used the Flory-Huggins free which is a very popular model to study polymer solutions. It is given by

 F = ∫(−∑σθnσmσ+∑σnσθlnϕσ (45) +∑σ∑σ′12χσσ′θmσnσϕσ′)dV,

where is the polymerization of the component , is its number density, and is its volume fraction. It is defined as

 ϕσ=mσnσ∑σ(mσnσ)=ρσρ, (46)

where is the mer density of component and is the mer density of the system, which is a constant in the Flory-Huggins model. To validate our algorithm we compared the binodal lines obtained by our algorithm to the theoretical ones obtained by minimizing the free energy. We used the interfacial tension parameter in all our LB simulations of binary and ternary systems, because there is an intrinsic surface tension in the LB simulation due to higher order terms Wagner (2006), which did not appear explicitly in the second order Taylor expansion presented in this paper. Since we are only evaluating the phase-behavior here we use a one dimensional model known as D1Q3. This model has the velocity set . This is an important test since all other frequently used higher dimensional model have this D1Q3 model as a projection.

We consider two binary systems: a monomer system with and , and a polymer system with and . For both systems, the total density was . Throughout this paper we choose the self interaction parameters to vanish: . The critical volume fractions for the monomer system are and and for the polymer systems are and . To induce phase separation a small sinusoidal perturbation was added in the initial conditions. The amplitude of the perturbation is 0.1 and its wavelength is the lattice size. The initial volume fraction of component A is given by . The initial volume fraction of component B is given by .

The monomer system was simulated with different inverse relaxation times. In Figure 1 we show results for , and . We see that the equilibrium densities have only a very slight dependence on the relaxation time, although the range of stability depends noticeably on the relaxation time. The polymer system was simulated with only one inverse relaxation time of . Starting from the critical point and increasing the value for each initial condition until the simulations were numerically unstable, we obtained a pair of binodal points for each initial condition. The system reached a stable state after about 5000 time steps. The measurement were taken after 50000 time steps to be sure that an equilibrium state had been reached.

For the polymer system, Figure 2 shows the comparison of the total density, and the volume fractions and chemical potentials of each component to the corresponding theoretical values. The total density of a system in equilibrium by LB is essentially constant with a variation of . The volume fractions of each component in the LB simulation agree well with the theoretical values. The chemical potential of each component by the LB simulation was very close to the theoretical value. The chemical potential , corresponding to the polymer component, varied slightly with a difference for the bulk values of about and a variation in the interface of about . This is the underlying reason for the small deviation from the theoretically predicted concentration. The potential was nearly constant with a variation of less than in the bulk and a variation of about at the interface. For large values of this discrepancy increases leading to the noticeable variation of the equilibrium densities of the polymer system as shown in Figure 1.

We also performed LB simulations with two ternary systems: a monomer system with , , and and a polymer system with , , and . The parameters for both systems were , , and . The other parameters were zero. The inverse relaxation time constant for both simulations was . The critical point for the monomer system was , , and . The critical point for the polymer system was , , and . The initial state of each simulation were set from the critical points towards the end point (, , ). Initially a small sinusoidal wave perturbation of an amplitude of 0.1 and wavelength of the lattice size was superimposed on the initial volume fraction of the A component. This perturbation was subtracted from the B component and the C component was constant. I performed a LB simulation for each set of initial volume fractions and obtained the volume fractions of the two phases in the equilibrium state, resulting in two binodal points. The simulation reached a stable state after about 20,000 time steps. The measurements were taken after 200,000 time steps to make sure the equilibrium state was reached.

Figure 3 shows the comparison of the binodal points by LB simulation to the theoretical binodal lines of both systems. The binodal points obtained by the LB simulation agree fairly well with the theoretical binodal lines for the monomer and polymer systems. The simulation becomes unstable when is close to zero, i.e. when one component is nearly depleted. In this region the simulation results also deviate noticeably from the theoretical binodal lines. Immediately near the critical point, the evolution of the system becomes extremely slow so the slight deviation between the binodal points obtained through the LB simulation and the theoretical ones probably indicates that the LB simulation was not yet fully equilibrated.

For the polymer system Figure 4 shows a comparison of the volume fractions and chemical potentials of each component. The total density of the system is again nearly constant with variation of less than in the bulk. At the interface there is a small variation of . The volume fractions of each phase in the simulation were very close to their theoretical values. The chemical potential of component A was slightly different in two phases with a variation of about , while the chemical potentials of components B and C were much closer in the two phases with a variation of less than .

## V Outlook

We have presented a general lattice Boltzmann algorithm for systems with an arbitrary number of components which is based on an underlying free energy. In this algorithm the key thermodynamic quantities such as the chemical potentials of the components are immediately accessible. It is also manifestly symmetric for all components. We tested the equilibrium behavior of the new algorithm for two and three component systems in each case examining both the case of monomer and polymer mixtures with an underlying Flory-Huggins free energy. We obtained to expected phase-diagrams to good accuracy and the chemical potentials were constant to good accuracy for the monomer systems. Polymer systems were more challenging to simulate but still obtained acceptable results for . Higher polymerizations, however, become increasingly difficult to realize with the current algorithm.

There are three directions in which we hope to extend this algorithm in the future. The current algorithm does not allow for component dependent mobility. We are working on developing an algorithm that can recover an arbitrary mobility tensor . The chemical potential is only approximately constant. Recent progress for liquid-gas systems Wagner (2006) makes us hopeful that we will be able to ensure that the chemical potential is constant up to machine accuracy. And lastly we hope to extend to model so that it can simulate polymer systems with significantly larger polymerization.

### References

1. A. Akthakul, C. Scott, A. Mayes, and A.J. Wagner, J. Membrane Sci. 249, 213 (2005).
2. D.H. Rothman and J.M. Keller, J. Stat. Phys. 52, 1119 (1988).
3. T. Reis and T.N. Phillips J. Phys. A 40, 4033 (2007).
4. X. Shan and H. Chen, Phys. Rev. E 47, 1815 (1993).
5. X. Shan and H. Chen, Phys. Rev. E 49, 2941 (1994).
6. X. Shan and G. Doolen, Phys. Rev. E 54, 3614 (1996).
7. W. Osborn, E. Orlandini, M.R. Swift, J.M. Yeomans, and J.R. Banavar, Phys. Rev. Lett. 75, 4031 (1995).
8. M.R. Swift, E. Orlandini, W.R. Osborn, J.M. Yeomans, Phys. Rev. E 54, 5041(1996).
9. A.J. Wagner and M. Cates, Europhys. Lett. 56, 556 (2001).
10. A. Lamura, G. Gonnella, and J.M. Yeomans, Europhys. Lett. 45, 314 (1999).
11. R. Jones, Soft Condensed Matter (Oxford University Press, 2002).
12. A.J. Wagner, Phys. Rev. E 74, 056703.1 (2006).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters