Surface Interaction of Nanoscale Water Film with SDS from Computational Simulation and Film Thermodynamics

Abstract

Foam systems have been attracting extensive attention due to their importance in a variety of applications, e.g., in the cleaning industry, and in bubble flotation. In the context of flotation chemistry, flotation performance is strongly affected by bubble coalescence, which in turn relies significantly on the surface forces upon the liquid film between bubbles. Conventionally, unusual short-range strongly repulsive surface interactions for Newton black films (NBF) between two interfaces with thickness of less than 5 nm were not able to be incorporated into the available classical Derjaguin, Landau, Verwey, and Overbeek (DLVO) theory. The non-DLVO interaction would increase exponentially with the decrease of film thickness, as it plays a crucial role in determining liquid film stability. However, its mechanism and origin are still unclear. In the present work, we investigate the surface interaction of free-standing sodium dodecyl-sulfate (SDS) nanoscale black films in terms of disjoining pressure using the molecular simulation method. The aqueous nanoscale film, consisting of a water coating with SDS surfactants, and with disjoining pressure and film tension of SDS-NBF as a function of film thickness, were quantitatively determined by a post-processing technique derived from film thermodynamics.

1. Introduction

The stability of foam depends a great deal on thin liquid films between the bubbles. A better understanding of the function and behavior of bubble coalescence requires detailed knowledge of the surface interaction of the liquid films. Surface force is one of the most fundamental thermodynamic properties that characterize the stability of thin films. The film thickness has great influence on the surface force. As liquid films become thinner (Figure 1a), these films begin to exhibit a black color, as there is nearly no light been reflected. These black films can be divided into two types based on their thickness ranges: common black films, with thickness 10~100 nm; and Newton black films, which are usually are much thinner, with thicknesses below 5 nm.

Experiments [1,2] have indicated that Newton black films (NBF) can be fairly stable in the absence of water evaporation and mechanical disturbances. Classical DLVO theory can successfully predict and explain the stability for CBF based on the contributions of electrostatic double layer (EDL) and attractive van der Waals (VDW) interactions, and CBF thickness with salt concentration could be well quantified [3] by the Poisson-Boltzmann theory. However, NBF with <5 nm thickness is sometimes not possible to interpret on the basis of EDL + VDW as suggested by conventional DLVO theory, because a remarkable repulsive interaction is observed within short range for NBFs that cannot be classified.

As an important quantity governing film stability, it is possible to determine disjoining pressure (Π) based on the sum of the so-called attractive component, e.g., VDW or hydrophobic interaction, and the repulsive part, e.g., EDL, according to extended-DLVO theory [5]. On the other hand, disjoining pressure Π (shown in Figure 2) can be regarded as the difference between the normal pressure of the film and pressure of the homogeneous bulk liquid. In fact, disjoining pressure could be regarded as the molecular interaction within the liquid film. From the perspective of molecular modeling, the molecular interaction using molecular simulation can represent the sum of the surface force components, as in classical DLVO theory.

To be specific, the net (or total) disjoining pressure of a liquid film, according to classical DLVO theory, can be divided into several components, depending on its molecular types,

$$\Pi ={\Pi }_{vw}+{\Pi }_{el}+{\Pi }_{hb}$$

(1)

where Π is the disjoining pressure, and Π_vw, Π_el, and Π_hb are the VDW, electrostatic, and hydrophobic components, respectively. The VDW can be written as:

$${\Pi }_{vw}=-\frac{A_{232}}{6\pi H^3}$$

(2)

where A₂₃₂ is the Hamaker constant between two gas phases and the liquid phase, H is thickness. The hydrophobic force has a similar form:

$${\Pi }_{hb}=-\frac{K_{232}}{6\pi H^3}$$

(3)

where K₂₃₂ is hydrophobic interaction parameter between the phases. For the component of electrostatic double-layer (EDL), the weak-overlap expression can be used:

$${\Pi }_{el}=2RT{C}_s\left(\cosh \left(\frac{e{\psi }_{mid}}{kT}\right)-1\right)$$

(4)

where R is the gas constant, T is temperature, C_s is electrolyte concentration, and ψ_mid is the potential at the mid-plane of the film. At equilibration [6], P_c is the capillary pressure.

$$2RT{C}_s\left(\cosh \left(\frac{e{\psi }_{mid}}{kT}\right)-1\right)-\frac{K_{232}}{6\pi H^3}=P_c$$

(5)

Classical DLVO theory does not consider the factors that arise from the microstructure of water molecules [7], and regards water as a uniform medium, which does not accurately identify the role of this short-ranged and strongly repulsive interaction. It has been proposed previously that such “non-conventional” interactions [8] are rooted in the density profiles of the molecules that form the double-layer NBF, and exhibit entropy characteristics. Such an interpretation provides a reasonable benchmark account for the strong stability of NBFs. With the fast advancement by computational resources, a wide range of applications related to complicated subjects has become examinable with simulation methods, e.g., Molecular Dynamics (MD) techniques, which can be applied to investigating a range of systems [9,10,11] such as carbon nanotubes or diamond nano-thread, and under various conditions (e.g., compounds onto water droplets) [12,13,14]. It has been reported [15] that, by using MD simulations, electrostatic energy could be obtained based on electrostatic potential profile, by taking into account the uniform characteristic of water. The electrostatic contribution of disjoining pressure [16,17] can thus be derived, and its value is in similar order of magnitude to that of the VDW component from experimental evaluation.

Certain quantities and liquid film properties [18,19,20,21] can be readily observed or calculated using computational simulations. Determination of disjoining pressure would be of great interest in liquid film research. However, it is one of the most challenging tasks in molecular simulations. Winter et al. [22] and Bhatt [23,24,25] developed a methodology for determining disjoining pressure, and studied the disjoining pressure isotherm for thin films with LJ fluid, water and surfactant. Their results showed that the disjoining pressure of LJ fluid and water were purely negative. However, disjoining pressure was positive for surfactant, and decreases with an increase in film thickness.

Classical DLVO theory is usually applied to interpret the stability of foams and films, and it has been commonly regarded as the sum of the van der Waals (VDW) component and electrostatic double layer (EDL) component. However, discrepancies between DLVO theory and experimental observations (e.g., the existence of Newton black films not predicted by DLVO theory, or even some common black films) have been found, indicating that DLVO forces alone are not able to explain the results. Thus, non-DLVO forces, such as steric force, hydration force and hydrophobic force have been proposed. This insufficiency is significant, especially for very thin thicknesses (e.g., 1–2 nm), as the classical Hamaker theory neglects fluid structure and entropy completely, and assumes a slab density [24]. VDW interaction, as in DLVO, is based on Hamaker theory (or, the classical Hamaker constant was estimated from Lifshitz’s formulation). Therefore, rigorous calculation of disjoining pressures requires detailed molecular simulations with physical molecular profiles.

Previously, Π-h of liquid film with salt [26] has been calculated, and it was found that VDW + EDL from DLVO theory were not able to explain the simulation results, with Π from DLVO theory being 1–2 × 10⁶ Pa lower than the results obtained from the MD simulations for film thickness at 1~3 nm. Furthermore, it was suggested that hydration force could be contributing to the difference between the simulation results and the classical DLVO theory using SPC/E-ions model. In this study, SDS is applied, as this molecule serves as a popular model for computational studies. The difference between wetting film and free film is shown in Figure 1b; in this study free-standing liquid film with SDS was investigated, to examine the disjoining pressure and film tension using a convenient alternative method to MD simulations. This provides a basis for validation of such systems and the further interpretation for experimental measurements.

2. Methodology

2.1. Theoretical Approach from Film Thermodynamics

For an aqueous thin film [27], the Gibbs free energy differential of thin film can be written as:

$$d\left(F^f+P^{\beta }{V}^f\right)=-S^{f}dT+Ahd{P}^{\beta }+2{\sigma }_{f}dA-\Pi Adh+{\displaystyle \sum _i{\mu }_i}d{N}_i^f$$

(6)

Here, h is film thickness, $S^f$ is the excess entropy of the thin film, ${\sigma }_f$ is the surface tension of the thin film, ${\mu }_i$ is chemical potential of the ith component and $N_i^f$ is mole number of component i. From the point of energy change, disjoining pressure could be thermodynamically regarded as the change of Gibbs free energy [28,29,30] within the thin film, that is:

$$-{\left(\frac{\partial \left(G/A\right)}{\partial h}\right)}_{T,P,A,N}=\Pi$$

(7)

For the surface tension of the film σ_f, it could be written as:

$${\left(\frac{\partial G}{\partial A}\right)}_{T,P,H,N}={\sigma }_f$$

(8)

By integrating Equation (6), there is the Gibbs-Duhem type Equation (9),

$$-2d{\sigma }_f=\left(\frac{S^f}{A}\right)dT-hdP+\Pi dh+{\displaystyle \sum _i{\Gamma }_i}d{\mu }_i$$

(9)

It can be seen from Equation (9) that disjoining pressure has following relationship:

$$-2{\left(\frac{\partial {\sigma }_f}{\partial h}\right)}_{T,P,\mu }=\Pi$$

(10)

Integrating Equation (10), the corresponding film tension can be given in Equation (11).

A key equation based on Equation (6) (also seen from Figure 2, pressure balance) can be obtained:

$${\gamma }_{\mathrm{f}}=2{\sigma }_{\infty }+{\displaystyle {\int }_h^{\infty }\Pi dh}+\Pi h$$

(11)

Here γ_f is film tension of the liquid film; σ_∞ is surface tension of bulk liquid. There is also another correlation:

$${\gamma }_{\mathrm{f}}=2{\sigma }_f+\Pi h$$

(12)

Combining Equations (11) and (12), disjoining pressure can be obtained, and the key process is as described below. It can be seen that the two types (energy, pressure balance) would lead to the same inner relationship for disjoining pressure and surface tension.

$${\displaystyle {\int }_h^{\infty }\Pi }dh=2{\sigma }_f-2{\sigma }_{\infty }$$

(13)

$$\Pi =-\frac{{\displaystyle {\int }_h^{\infty }\Pi }dh}{dh}=-\left(\frac{2{\sigma }_f-2{\sigma }_{\infty }}{dh}\right)=-{\left(2{\sigma }_f-2{\sigma }_{\infty }\right)}_h^{\prime }=-{\left(2{\sigma }_f\right)}_h^{\prime }$$

(14)

By performing derivation directly on both sides of Equation (13), then Equation (14) could be obtained. Please note that σ_∞ is the surface tension of bulk liquid, and it is a constant; thus, its derivation is zero. From Equation (14), disjoining pressure has a relationship with the surface tension of thin film (σ_f) and film thickness. To be more specific, the 2σ_f (double surface tension) at varied thicknesses needs to be obtained first; then, there was 2σ_f vs. h. After that, a derivative procedure was conducted on the curve of 2σ_f vs. h; subtraction of the result would result in the final disjoining pressure isotherm (Π-h). Surface tension was determined by the usual pressure tensor way [31] within GROMACS, i.e., given by

$${\sigma }_f=\frac{1}{2}{\displaystyle {\int }_0^{L_z}\left(P_N-P_T\right)}dt$$

(15)

where P_N and P_T are pressure tensors, at normal and tangential directions, respectively.

2.2. Computational Modeling System

The schematic simulation setup is shown in Figure 3, the thin film contains SDS and water molecules. The lateral dimension of the simulation box is L_x = L_y (~13 nm), and L_z = 3 L_x, two vapor areas were placed in the Z direction to form the film structure. The vapor size on the Z-axis is big enough to avoid possible interaction between the up and bottom surface.

Water model TIP4P/2005 [32] was applied using GROMACS 4.5.3 [33,34,35] under periodic boundary conditions (PBC) in all directions. Long-range corrections to the Coulomb interactions adopted the particle mesh Ewald (PME) [36] method; both Lennard-Jones (LJ) and Coulomb terms had a potential cut-off 1.4 nm. The Nose-Hoover method was applied to keep the system temperature at T = 298 K, with each running step set as 2 fs. Surface tension was obtained by running in the NVT (number of molecules, volume and temperature were kept constant). SDS was modeled by united atom force fields [16], the parameters of the SDS head group were from the AMBER force-field [37], Na⁺ parameters were from Dang [38]. The TIP4P/2005 water model was suggested to be better for describing water surface tension. Film thickness, h, was determined as the thickness (from configuration) within the NVT running process, in order to keep consistent.

To well equilibrate the simulation system and minimize the potential, 10 ns were first run, and another 10 ns trajectory for calculating average properties. Each result presented is the average of at least 3 simulations, each with a different initial configuration. Surface concentration of SDS was fixed at ~33 Ų, as previous verified by X-ray experimental data [39] and suggested that film usually contains between two and four H₂O for one SDS molecule. The value is in agreement with observation in the spontaneous formation of the SDS-NBF. Water molecules were placed at the middle of the SDS surfactants, forming a sandwiched slab geometry (Figure 3a), and Na⁺ could move around the SDS (Figure 3b), and can also interact with water molecules. The film thickness changes with the change of ratio of N_water/N_SDS, e.g., ratio of N_water/N_SDS = 3, for each monolayer 256 SDS (N_SDS = 256 for each surface, and 512 for two surfaces), N_water = 1536 is total. With the changing of N_water/N_SDS from 2 to 3.5, film thickness and the corresponding disjoining pressure and film tension will change.

3. Results and Discussion

3.1. Disjoining Pressure Isotherms

Surface tension and disjoining pressure are the key parameters that can be applied to obtain the film tension of a liquid thin film. In molecular modeling, film tension cannot be directly obtained; the surface tension and disjoining pressure need to be determined firstly. It is noted that the disjoining pressure is crucial to understanding the behavior of drainage and rupture processes of surfactant and non-surfactant thin films. In addition, the understanding of the effects of specific ions on film- or bubble-related problems [40,41,42,43] also needs a better, more concise description for the tensions within these thin liquid films. By molecular simulation, it is possible to calculate the surface tension of bulk liquids, and also possible to calculate the surface tension of thin films (i.e., of smaller thickness). However, it is not possible to calculate film tension directly by simulation.

One methodology for Π was developed by Bhatt [23,24,25], under the same chemical potential, by comparing the difference between the normal pressure and the surrounding bulk liquid pressure. This study applied film thermodynamics, and derived disjoining pressure using 2σ_f vs. h; the obtained 2σ_f vs. h curve and corresponding fitting line are shown in Figure 4.

The curve of 2σ_f vs. h is empirically fitted (shown in Figure 4a) using a single exponential decay model with three parameters, and the obtained parameters for SDS-stabilized films are given by Equation (16), below, with R = 0.99.

$$2{\sigma }_f=57.15+25130.38\times e^{-2.87h}$$

(16)

For very thin films, e.g., nanoscale black films, the surface interaction becomes stronger with a decrease in film thickness, and the corresponding surface force (in terms of disjoining pressure) will not be zero. Surface tension for SDS-stabilized nano-film was observed to decrease with increasing film thickness; this is opposite to the case of pure water film, as σ_f of water films would increase with increasing film thickness. Please note that, for liquid films with very large thickness, their σ_f is very similar to the value of the bulk liquid (i.e., very thick film with film thickness ∞). However, with the decrease in film thickness, its surface tension would change accordingly, and this change is more pronounced for thinner films, as can be seen in Figure 4a.

Equation (14) is used for the determination of Π; the corresponding Π-h below 5 nm thickness is obtained and shown in Figure 4b. In turn, the disjoining pressure isotherm derived by 2σ_f vs. h can be expressed as:

$$\Pi =72042.76\times e^{-2.87h}$$

(17)

the units of σ_f, h, and Π are mN/m, nm, and in MPa (i.e., ×10⁶ Pa), respectively. The disjoining pressure is positive across the entire given thickness range. For films with thickness range 30–100 nm, the disjoining pressure is below 10⁵ Pa, i.e., ≤1 bar, at such film thicknesses, their Π is easy to detect and measure, and cannot be negligible under experimental conditions. However, at film thicknesses of 30–100 nm, due to the small Π (≤1 bar), it is not possible for MD to capture its accurate Π value. Molecular simulations can only simulate very thin films (usually ~10 nm, below 100 nm).

The disjoining pressure can be studied from different perspectives, as shown in Figure 5. In this study, we applied a non-polarizable model, i.e., Equation (18) was applied:

$$U_{pair}={\displaystyle \sum _i{\displaystyle \sum _j\left(4\epsilon \left[{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^{12}-{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^6\right]+\frac{q_i{q}_j}{r_{ij}}\right)}}$$

(18)

LJ fluid thin film has only attractive components, that is to say, there are no other components. For LJ fluid, there is no Coulombic interaction in Equation (18). However, for non-LJ fluid, water molecules coated with SDS have LJ and Coulombic interactions (plus polarizable interactions, when the polarizable model is applied). The LJ term and Coulombic term cannot provide further insight into the surface force components. It has been reported [44] that, by fitting Electrostatic Helmholtz free energy exponentially, it is possible to derive the electrostatic contribution of Π. Disjoining pressure in Figure 4b contains all of the contributions. There is no solid theory to divide the so-called various components, and these interactions cannot be distinguished in MD simulations.

Classical theory can well interpret the observed results for most cases in usual ranges. However, that is not always satisfactory, especially under some extreme conditions, where it may, for example, underestimate VDW interactions by 1–2 orders of magnitude compared to the molecular simulation results for water films with 1–2 nm thickness.

3.2. Film Tension

Usually, there is very small (and sometimes negligible, in the case of very thick films) difference between film tension and 2σ_∞ for common free liquid films, or for CBF. Hence, experimental measurements of film tension can be within the margin of uncertainty for these relatively thick liquid films. The important properties of the liquid bulk solution from which the aqueous film is produced dictate the final state, and the corresponding tension values (surface tension, film tension, disjoining pressure) change with variations in temperature and concentration. For these CBFs, indirect measurement can be preferable, as the relationship (cosθ = γ/2σ_∞) defines the ratio of film tension with that of the surface tension of bulk solution. For NBFs, direct measurement for γ is possible. The surface tension, Π, and film tension are mainly dependent on film thickness.

The obtained disjoining pressure profile was subsequently employed to investigate the film tension for the thin films; Equation (12) is applied to determine film tension, and the results are shown in Figure 6. It can be seen that film tension decreased gradually as the film thickness increased, and the decreasing trend was more pronounced in the small thickness range. Moreover, the film tensions thus determined were positive values across the entire given thickness range. The film tension is great in the small thickness range, and becomes much smaller in the large thickness region, indicating that thick film tends to gradually fall into the bulk solution zone.

4. Conclusions

The disjoining pressure isotherm is an important thermodynamic property for aqueous thin films, and surface interaction is not zero when the liquid film is sufficiently thin, especially when it is nanoscale. Disjoining pressure was determined in this study for SDS-stabilized nanoscale liquid film by deriving the curve of 2σ_f vs. h based on film thermodynamics. The σ_f was different for varying film thicknesses, which indicates that surface interaction plays an important role at the nanoscale level. The determined Π for nanoscale liquid films with SDS was positive (the total Π is repulsive) at all film thicknesses studied, and of comparable magnitude to previous simulation results.

This paper is the first determination of NBF film tension and disjoining pressure using a convenient method, without complicated process, by molecular simulations using a post-processing method. However, it is not possible to separate the total interaction into various components at this stage. We have no solid theory to separate these components, and sometimes these individual components may not be additive. The present study has potential implications for the interpretation and validation of similar experimental measurements.