Multiphase Phase-Field Lattice Boltzmann Method for Simulation of Soluble Surfactants

: This paper proposes a phase-ﬁeld model for the lattice Boltzmann method which has discretized symmetrical directions of velocities in a cartesian grid, to simulate the soluble surfactant in a Multicomponent multiphase system. Despite other existing phase-ﬁeld models following Langmuir relation, the interfacial tension can be calculated analytically in this proposed model. Parameters playing roles in the models and controlling the surfactant’s strength and interaction with other phases are obtained directly from a given initial interfacial tension and bulk surfactant. Consequently, there is no further need for trial-and-error simulations, and a real system, e.g., oil-water-surfactant, can be simulated with given initial parameters. The model is validated with the analytical result for a planar oil–water-surfactant system. Furthermore, the method for reobtaining numerical interfacial tension for ﬁve different cases is tested and compared with the given initial values for an oil droplet surrounded by water and surfactant. The results show that the obtained interfacial tension from the method is in good agreement with the given initial interfacial tension. Furthermore, the spurious velocity of the model is calculated and seen that the magnitude of spurious velocities is proportional to interfacial tension.


Introduction
Surfactants are amphiphilic compounds containing hydrophobic and hydrophilic groups that cause an imbalance at the interface and reduce surface tension [1]. Surfactants are interfacially active agents which selectively adhere to fluid interfaces creating a buffer zone to reduce the system energy. It plays an essential role in different industries such as process industry, pharmaceutical industry, food processing [2] or can be used in microfluidic applications [3]. In microfluidic systems, surfactants are often used to thermodynamically stable systems such as water-in-oil and oil-in-water emulsions [4]. Moreover, surfactants are expected to alter droplet dynamical behavior in the microfluidic devices significantly.
There are studies in binary mixture related to the surfactant theoretically and experimentally. The effect of surfactant on the droplet dynamics is experimentally investigated in [5,6]. By experimentally investigating the equilibrium oil-water mixture, Mulqueen and Blankschtein [5] proposed a theoretical model for the absorbed surfactant at the oil-water interface. The proposed model is able to predict the interfacial properties of the solution in the presence of surfactant without performing any further experiments on the mixture. Zhang et al. [6] showed that the surfactant concentration in the interface of a dielectric oil droplet can enhance the drop deformation under electrical force. Besides experimental and theoretical studies, the behavior of surfactants on the droplet can be numerically investigated. Two approaches of sharp and diffuse interfaces have been used in immiscible binary mixtures [7,8]. In the sharp interface, the interface is located with accurate calcula-tion of curvature and in the diffuse interface methods, the interface has a finite thickness which is small in comparison to the smallest length scale in the domain. Jesus et al. [9] applied a hybrid method in which the sharp interface method is coupled with an Eulerian method to investigate soluble surfactant in the mixture. Similar methods have been used by other researchers [10,11]. The phase-field method uses a free energy function to capture the interface in the multiphase system as a diffusive interface model [12]. The method is successfully applied for the simulation of immiscible binary mixtures [13][14][15][16][17]. A lattice Boltzmann method is proposed by Van der Sman and van der Graaf to simulate multiphase system with surfactant [15]. Later, the method has been improved in [18] by a better approximation of the delta function.
Although several researchers proposed models for the simulation of the surfactant using the phase-field method, these models still suffer from problems that have not been solved so far. One problem that is addressed in this paper is that there is no full control on the systems with surfactants. Each model has some parameters determining interfacial tension and bulk surfactant concertation, but it is impossible to set interfacial tension as an input value for the system and calculate the rest of the required parameters. The only way to find these parameters is using trial-and-error simulations, which leads to run the program several times to find appropriate values. The reason that causes us to lose full control is that the models have complicated mathematical forms, have to follow Langmuir relation, and consequently, an analytical solution cannot be obtained. In this paper, we propose a model that allows us to fully control the system and determine the required parameters directly by giving an initial interfacial tension and bulk surfactant.
In the proposed model, the surfactant concentration of the system is shifted to the bulk surfactant concentration, and the interfacial surfactant concentration is controlled with only one parameter. Since there is an analytical solution for interfacial tension of the proposed method, the parameter playing a role in determining interfacial tension in the model can be calculated directly without any trial-and-error simulation. The required bulk surfactant concentration of the system is set in the model directly. A simple thermodynamic analysis of the system in the equilibrium can prove that the model agrees well with the Langmuir relation. The numerical instability, which is observed in the logarithmic form (at low surfactant concentration), does not occur in the present model. In this paper, the model is described, and mathematical derivation is demonstrated. Furthermore, validations are given to test the accuracy of the model. We would leave the study of dynamics of surfactant layer for further publications.
In this paper, in Section 2, the numerical method is explained, and the mathematical derivation of the proposed method is described. The thermodynamics of the system in the equilibrium based on the model is discussed in Section 2.2. The LBM as a fluid solver is explained in Section 2.3. The model is compared with an analytical solution in Section 3. Finally, the conclusion of the paper is given in Section 4.

Numerical Method
A phase-field model capturing thermodynamic and hydrodynamics effects associated with surfactants in realistic systems is proposed and presented here. The mathematical derivation is given in the following.

Proposed Phase-Field Model for Immiscible Fluids Including Surfactants
The Landau-Ginzburg free energy functional is used to determine the thermodynamics of a binary mixture including surfactant. The interfacial tension of binary mixture is reduced by adding surfactant in which it adheres to the interface. Extra terms are defined and added to the original Landau-Ginzburg free energy functional to account for the surfactant effect. The proposed free energy functional is introduced as: is the order parameter representing the relative concentration of the local compositions. The bulk phase behavior for pure mixture binary is described by the first two terms with a minimum φ b = » A B for component one and two. In this paper, we set φ b = 1, so that parameters A and B are equal. The third term in the equation determines the interfacial of clean system. The remaining terms specify the behavior of surfactant in which the system is shifted to the bulk surfactant concentration ψ b . The forth term is the entropic part of free energy of mixing of the surfactant with the bulk phase. The pronominal form is substituted for common logarithmic form. This simplification allows us to get around the complexity of logarithmical form although we still see the Langmuir relation in the model. The surfactant order parameter ψ is in fact normalized and it is equal to one for fully saturated interface with surfactant. The parameters K B and T are Boltzmann constant and tempreture of the system, respectively. The last term controls the loading of surfactant at the interface. Parameter d determines the interfacial tension of a mixture binary with surfactant.
The interfacial tension of the whole domain can be calculated by integrating the excess free energy per unit interface area which can be defined by W: which satisfies the Euler-Lagrange minimization equations: After substituting Equations (1)-(4), the three basic equations are obtained: The surfactant concentration can be calculated by the use of Equation (7) as: With the use of Equations (6) and (7), the interfacial thickness and a solution for the order parameter is obtained as: By integrating the excess free energy per unit interface area (Equation (5)) in the whole domain, the interfacial tension is obtained: The chemical potentials µ φ and µ ψ can then be obtained via the variational derivatives of the free energy functional with respect to φ and ψ: The divergence of pressure tensor can be calculated by Excess chemical potential gradients from the Gibbs-Duhem equality: with the pressure tensor P given by where I and P 0 are the second-order unit tensor and the scalar part of the pressure tensor, respectively. It can be calculated by the thermodynamic relation as:

Thermodynamic Equilibrium
The chemical potentials in the equilibrium are constant throughout the entire system. In this section, the prediction of the equilibrium properties of the surfactant adsorption as described by the Langmuir isotherms of the proposed model is analyzed. There is no need to assume that the bulk surfactant concentration is much smaller than unity for further simplification despite of many proposed model. From Equation (13), the relation between the chemical potentials at bulk µ b and at interface µ o is obtained: This equation can be also obtained by Equation (8) in which φ = 0 at the interface.

Lattice Boltzmann Method
In order to fully characterize and demonstrate the hydrodynamics of an immiscible binary fluid, four main equations, continuity equation, the Navier-Stokes equation and the Cahn-Hilliard equations describing the order parameter and surfactant concentration are used. They can be summarized as follows: in the above equations, ρ, u, P and µ are the fluid density, velocity, pressure tensor and dynamic viscosity, respectively. The mobilities for the order parameter, φ, and the surfactant concenteration, ψ, are M φ and M ψ , respectively.
The symmetrical lattice Boltzmann algorithm is used to solve Equations (18)- (21) to determine the macroscopic properties. The method has second order accuracy in space and time which can recover Navier-stokes accurately. Three particle distribution functions of f i (x, t), g i (x, t) and h i (x, t) are used to describe multiphase flow with surfactant. The time evolution equations for the LB equation for these three distribution functions can be written as: where f eq , g eq and h eq are the equilibrium distribution functions and τ f , τ g and τ h are the relaxation times. The macroscopic variables are related to the particle distribution functions as follows: The relaxation times in the lattice Boltzmann algorithm are related to the physical variables as follows: The equilibrium distribution functions are defined below so as to obtain the continuum equations related to an immiscible binary fluid: where the subscripts α and β represent the components along x and y directions, respectively, and δ αβ is the Kronecker delta.

Model Validation
The proposed method for predicting the profile of surfactant concentration and the order parameter at a planar oil-water interface is tested and the numerically obtained results are compared with analytical solution. A symetrical 2D computational domain with 200 × 10 lattice cell is considered with an oil phase initially located at 50 < x < 150 and periodic boundary conditions are used for all the boundaries. The remaining parameters are σ 0 = 0.04, ξ 0 = 2∆x, K B T = 0.01, M φ = M ψ = 1, τ φ = τ ψ = 1. Two different cases are tested where case one is for σ = 0.03, ψ b = 0.312 and case two is for σ = 0.025, ψ b = 0.03. An excellent agreement is observed between our numerical results and the analytical solution, as it is demonstrated in Figure 1 for both the surfactant concentration and the order parameter. The profiles of the surfactant concentration are compatible with the analytical solution even for the relatively large value of ψ. This figure confirms that the order parameter is independent of the surfactant loading. It is assumed that surfactant concentration is low in most models. They use this assumption to simplify their models when they mathematically derive their methods. Consequently, their results are largely deviated from analytical solution at the high surfactant concentration [46]. However, an analytical solution is obtained in the proposed model without such an assumption or simplification.
The influence of interfacial thickness on the simulation results is studied. In addition to the result presented in Figure 1, in which interfacial thickness is equal two, we simulate the same setup with similar interfacial tensions and bulk surfactant concentrations but with ξ = 1.14∆x and ξ = 3∆x. The results are demonstrated in Figure 2. A small deviation is observed for ξ = 1.14∆x which is similar for a pure mixture binary system while the excellent result is obtined for ξ = 3∆x. The sharp profile of surfactant concentration can be correctly captured across the interface with a reasonably thick interface ξ = 2∆x demonstrated in Figure 1.
The presence of surfactants changes the interfacial tension. We can control the interfacial tension of droplets with a specific surfactant concentration in the emulation with this proposed model. To exam the accuracy of the model, we consider five samples in which their surfactant concentration and interfacial tension are given in Table 1. The samples are selected so that the Langmuir relation is satisfied. These samples can be extracted from experimental results.    An oil droplet resting at the water in a periodic boundary condition is considered for further validation of the model. The computational domain, the interfacial tension and thickness interface without surfactant are 100 × 100 lattice cell, σ 0 = 0.4 and ξ = 2∆x, respectively. The radius of droplet R is 30 lattice cell centered in the middle of the computational domain. The resting parameters are K B T = 0.01, M φ = M ψ = 1, τ φ = τ ψ = 1. It must be noted that parameter d is calculated by using Equations (10) and (11) and rearranging them as follows: where σ is the interfacial tension with the presence of surfactant concentration. When the droplet reaches its equilibrium, the interfacial tension is calculated by the Laplace law (for two dimension): where ∆p is the pressure difference across the droplet interface. For the five samples, the interfacial tension is obtained by the proposed method and compared with the given interfacial tension at the given surfactant concentration. The deviation of the numerical solution from the given interfacial tension in percentage is calculated and called error in the last column of Table 1. It is observed for all samples, the maximum error is less than 3%. The errors in the table show that the model can predict interfacial tension in equilibrium state very good and can be useful when we want to reobtain experimental data in which the availble information are only interfacial tension and surfactant concentration at the bulk. Surfactant loading at the interface can also be computed numerically and compared with the theory (Equation (17)) for the similar setup and samples when the system reaches the equilibrium. The comparison has been demonstrated in Figure 3 where the good agreement is observed. As it is expected from the analytical expression, a linear relation between the surfactant concentration at bulk and interface is shown. Time evolution of the surfactant concentration for sample three with a randaom distribution is demonstrated in Figure 4. The surfactant concentration is increased at the interface till it reaches the euqlibrium. In fact, the divergence of chemical potential creates a force and pushs surfactant to attach to the interface. while the surfactant concentration at the interface is increasing, the concentration at the bulk is decreased.

Time step=1
Time step=5 Time step=20 Time step=100 Time step=1000 Time step=20000 Spurious velocities at the phase interface in many numerical methods in multiphase flows cause artifacts in the simulation. The magnitude of the spurious velocities is proportional to the interfacial tension for a clean droplet in quiescent fluid [47]. The contour of spurious velocity for sample three is shown in Figure 5. Spurious velocities can be de-creased by using other collision operators of LBM in which multiple relaxation parameters are used [48]. The larger surfactant concentration in the system, the lower interfacial tension is expected and consequently it leads to bigger reduction of spurious velocities. This behavior is displayed in Figure 6 in which it is also proportional to initial interfacial tension. In this figure, the influence of the interfacial tension on the spurious velocities has been shown. The maximal spurious velocities for three different initial interfacial tension σ 0 = 0.002, σ 0 = 0.02 and σ 0 = 0.1 at bulk surfactant concentration ψ b = 0.5 are demonstrated. The surfactant concentration itself does not affect spurious velocities as long as the interfacial tension is constant. The model has a great potential to simulate multiphase flows with the presence of the surfactant. However, the model needs to be more investigated in various fluid flows in which the dynamics of the droplet is important in order to extract the limitations of the model. In this paper, we focus on the validation of the model and the ability of the model to reobtaining numerical interfacial tension. In the future works, the authors intend to study the limitation of the model in various fluid flow.

Conclusions
A phase-field model based on free energy is proposed to simulate the surfactant concentration on a binary mixture system. The proposed model allows us to have an analytical solution for the interfacial tension of the system in the presence of the surfactant concentration. The mathematical derivation has been demonstrated. The parameters determining the loading surfactant at the interface can be given as input values in this model. The full control of the system can be obtained using the proposed method in which the only required parameters are the interfacial tension and the bulk surfactant concentration of the system. This model can be used to reobtain the experimental results where there are only these two measured parameters. In fact, this model has promising features that can bridge between numerical and experimental data, which is not available explicitly. The LBM has been used to solve transport equations (the Cahn-Hilliard equation) for the order parameter and the surfactant concentration. The models can be used to simulate both symmetric and asymmetric fluid flow.
The order parameter and the surfactant concentration of numerical results are compared with the analytical solution at a planar oil-water interface. A good agreement is observed in this comparison. Furthermore, the effect of interfacial thickness is investigated, and it is shown the numerical results for the interfacial thickness larger than 2 are in good agreement in comparison with the analytical solution, although for ξ = 1.14∆x the numerical results are acceptable.
Furthermore, the ability of the method for reobtaining numerical interfacial tension for five different samples are examined and compared with the given initial values for an oil droplet surrounded by water and surfactant. The results show that the obtained interfacial tension from the proposed method agrees well with the given initial interfacial tension. The chemical protentional in the equilibrium obtained by the analytical anumerical solution are compared, and a good agreement is observed. In order to study the ability of the model to predict the time evolution of surfactant, a simulation has been carried out with a random distribution of surfactant concentration. The surfactant moves into the interface. The process has been demonstrated time-dependently.
The maximal value of spurious velocities as a function of initial interfacial tension for contaminated droplets has been shown. It is seen that the spurious velocity of the model is proportional to the interfacial tension. The study of the dynamics of a droplet with the presence of the surfactant is left for future works.