Nonequilibrium molecular dynamics simulations are used to investigate the influence of hydrodynamic interactions on vertical segregation (stratification) in drying mixtures of long and short polymer chains. In agreement with previous computer simulations and theoretical modeling, the short polymers stratify above the long polymers at the top of the drying film when hydrodynamic interactions between polymers are neglected. However, no stratification occurs under the same drying conditions when hydrodynamic interactions are incorporated through an explicit solvent model. Our analysis demonstrates that models lacking hydrodynamic interactions do not faithfully represent stratification in drying mixtures, in agreement with the recent analysis of an idealized model for diffusiophoresis. Hydrodynamic interactions must be incorporated into such models for drying mixtures in future.
Coatings or films formed by drying1,2 are relevant to many technologies, including latex paints,3,4 inkjet printing,5 and polymer nanocomposites.6 Such films are often comprised of multiple components including nanoparticles or colloids,7 polymers,8 and surfactants.4,9 It is desirable to control the distribution of components within the dried film, e.g., to ensure the dispersion of nanoparticles6,10 or to impart antimicrobial properties to a surface.11 Recent experiments showed that quickly dried mixtures of small and large spherical colloids can vertically stratify into layers.12–14 Unexpectedly, the smaller colloids accumulated at the solvent-air interface and pushed the larger colloids down into the film. Computer simulations qualitatively agreed with the experiments12–15 and formed the same small-on-top stratification in ternary and polydisperse colloid mixtures16 and in polymer–polymer and colloid–polymer mixtures.17 The stratification was more pronounced for larger size ratios and faster drying speeds, although it was found to persist even at moderate drying rates.13,15,17
Stratification has been modeled theoretically as a diffusion process7,18 characterized by the film Péclet number for each component, Pei = /Di, where H is the initial film height, is the speed of the solvent–air interface, and Di is the diffusion coefficient of component i. When Pei ≪ 1, diffusion dominates over advection from drying, and a nearly uniform distribution of component i is expected. When Pei ≫ 1, component i accumulates toward the drying interface. Trueman et al. developed a model for drying colloid mixtures19 that qualitatively predicted a maximum in the number of large colloids on top of the film when Pe1 < 1 and Pe2 > 1 (here 1 is the smaller colloid and 2 is the larger colloid).20,21 However, their model does not predict the small-on-top stratification that occurs when Pe1 > 1 and Pe2 > 1.19,20 Fortini et al. hypothesized that such stratification resulted from an osmotic pressure imbalance,12,16 while Howard et al.15,17 and Zhou et al.22 proposed models using the chemical potential as the driving force. Computer simulations of polymer mixtures and the Howard et al. model agreed quantitatively,17 and qualitative agreement was found between the Zhou et al. model and experiments for colloid mixtures in certain parameter regimes.14,21 The Zhou et al. model is equivalent to the low-density limit of the Howard et al. model.15
For computational and theoretical convenience, many previous studies have employed a free-draining approximation for hydrodynamic interactions to model stratification. This approach incorporates the Stokes drag on each particle but neglects other interactions. The computer simulations12–17 implicitly modeled the solvent using Langevin or Brownian dynamics methods that treated the solvent as a quiescent viscous background, while the theoretical models12,15–17,22 assumed that the particle mobilities were given by the Einstein relation. However, Sear and Warren pointed out in a recent analysis23 that neglecting hydrodynamic interactions overpredicts the phoretic velocity of a single large colloid in an ideal polymer solution, which must be proportional to the extent of stratification in this simple model. They argued that this discrepancy is caused by the omission of solvent backflow,23 leading to an incorrect hydrodynamic mobility.24 Their argument is qualitatively supported by differences between the experiments and implicit-solvent simulations of Fortini et al.12 To our knowledge, though, there has been no direct test of how hydrodynamic interactions influence stratification at finite concentrations or size ratios between components, which are more challenging to analyze theoretically.23,24
In this article, we systematically demonstrate the influence of hydrodynamic interactions on stratification using a drying polymer mixture as a model system. Computer simulations are ideal tools for such a study because, unlike in experiments, it is possible to artificially remove hydrodynamic interactions between particles in the simulations. Molecular dynamics simulations of drying were performed using both explicit and implicit solvent models. We show that although the implicit-solvent model stratifies for the drying conditions investigated, the explicit-solvent model does not due to the presence of hydrodynamic interactions. Our simulations show that it is critical to incorporate hydrodynamic interactions into models and simulations in order to reliably predict stratification in drying mixtures.
II. MODELS AND METHODS
We performed nonequilibrium molecular dynamics simulations of drying polymer mixtures using two bead–spring models: one with an explicitly resolved solvent (Sec. II A) and a corresponding one with an implicit solvent (Sec. II B). While the explicit model included full hydrodynamic interactions, the implicit model only incorporated the Stokes drag on each polymer, thus omitting hydrodynamic interactions between polymers. We have chosen to investigate the influence of hydrodynamic interactions on stratification using a polymer mixture rather than the colloid mixture originally studied by Fortini et al.12 because the polymer mixtures are more computationally amenable to the analysis in Sec. III. However, the results and conclusions presented here extend to other types of drying mixtures. Details of both polymer models are presented in the following.
A. Explicit solvent model
Polymers were represented as linear chains of monomers with diameter σ and mass m immersed in an explicit solvent of beads of equal size and mass. All monomers and solvent particles interacted with the Lennard-Jones potential,
where ε is the interaction energy and r is the distance between particles. The monomer–monomer interactions were made purely repulsive by truncating the potential at its minimum and shifting it to zero.25 The solvent–solvent and monomer–solvent interactions included the attractive part of Unb and were truncated at 3.0 σ with the energy and forces smoothed to zero starting from 2.5 σ. Bonds between monomers were represented by the finitely extensible nonlinear elastic potential,26
with the standard parameters r0 = 1.5 σ and κ = 30 ε/σ2. This choice prevents unphysical bond crossing because the equilibrium bond length is approximately 0.97 σ at a temperature of T = 1.0 ε/kB, where kB is Boltzmann’s constant. This model corresponds to good solvent conditions for the polymers.
The drying process was simulated using the method developed by Cheng and co-workers.27–29 The polymer solution was supported by a smooth structureless substrate modeled by a Lennard-Jones 9–3 potential at the bottom of the simulation box,
where z is the distance between the particle and the substrate and εw = 2 ε is the strength of the interaction. Interactions were truncated for z > 3.0 σ. The solvent vapor above the liquid film was confined by an additional potential of the same form as Uw at the top of the simulation box, which was made purely repulsive by truncating it for z > (2/5)1/6σ. The solvent was evaporated by deleting a small number of randomly chosen solvent particles from the top 20 σ of the vapor. To maintain temperature control, monomers and solvent in a slab of height 20 σ above the substrate were weakly coupled to a Langevin thermostat with T = 1.0 ε/kB and friction coefficient 0.1 m/τ,30–32 where is the unit of time.
All simulations were performed using the HOOMD-blue simulation package (version 2.2.2) on multiple graphics processing units33,34 with a time step of Δt = 0.005 τ. The solvent coexistence densities were 0.664 σ−3 and 0.044 σ−3 for the liquid and vapor phases, respectively. The simulation box was periodic in the x and y directions with a length of L = 50 σ. The height of the box was Lz = 540 σ, where H = 500 σ was the initial film height and the remaining 40 σ was filled with solvent vapor. We investigated a mixture of short polymers of M1 = 10 monomers and long polymers of M2 = 80 monomers at an initial monomer density of roughly 0.06 σ−3 each. The simulations consisted of 730 000 solvent particles, 900 long chains, and 7100 short chains, resulting in a total of 873 000 particles. By deleting one solvent particle every 1.25 τ, we obtained an evaporative flux of 3.2 × 10−4/σ2τ. The drying rate for this model solvent is estimated to be considerably faster than for solvents like water that are often used in experiments, but is close to the slowest rate that was computationally accessible to us. However, the simulated films were also much thinner than in experiments, resulting in realistic Péclet numbers.
B. Implicit solvent model
The implicit-solvent model was constructed by matching the structure of the polymer chains in the explicit solvent at infinite dilution. Interactions between the monomers and the solvent swelled the polymers beyond the size of a chain with only the monomer–monomer interactions, so additional effective interactions were required. We measured the distance distribution between monomers separated by two bonds along a chain surrounded by solvent and a chain without solvent. We fit the negative logarithm of the ratio of the explicit and implicit distribution with a spline potential,35
The fit resulted in parameters n = 3.8, εs = 0.68 ε, rs = 1.32 σ, and rc = 2.0 σ. This potential, modeling the soft repulsion between monomers due to the solvent, was added to the bare monomer–monomer interactions. We validated this effective model by measuring the polymer structure over a range of concentrations and compositions up to a total monomer number density of 0.375 σ−3, finding overall good agreement between the implicit and explicit models (see the supplementary material).
The polymer long-time dynamics in the explicit solvent were matched in the implicit solvent using Langevin dynamics simulations.30–32 This technique incorporates the effects of Brownian motion and Stokes drag from the solvent, but neglects hydrodynamic interactions between monomers. The monomer friction coefficients were adjusted for each polymer to give the same polymer center-of-mass diffusion coefficient at infinite dilution as in the explicit solvent. This approach gives the correct long-time dynamics, but distorts the internal relaxation modes of the polymers.17 We measured diffusion coefficients in the explicit solvent of D1 = 0.0342 σ2/τ for the M1 = 10 polymers and D2 = 0.0071 σ2/τ for the M2 = 80 polymers from the polymer center-of-mass mean squared displacement in a cubic box with edge length L = 50 σ, giving 2.92 m/τ and 1.76 m/τ for the friction coefficients, respectively (see the supplementary material for details).
The liquid–vapor interface was modeled by the repulsive part of a harmonic potential.15,17 (The complete form of the potential can be found in Ref. 15.) In order to closely match the explicit-solvent simulations, the spring constant and position of the interface were adjusted to obtain similar initial density profiles to the explicit-solvent model. We used spring constants of 0.1 ε/σ2 and 0.18 ε/σ2 for the monomers of the short and long chains, respectively. The position of the interface was measured throughout the explicit-solvent simulations and used directly in the implicit-solvent simulations. The minimum of the potential was offset for the long polymers by −1 σ.
III. RESULTS AND DISCUSSION
We performed 25 drying simulations for both models and measured the monomer density profiles in the film. Figure 1 shows the average profiles at three times: before evaporation, at a film height of roughly 3H0/4, and at a film height close to H0/2. We used the average interface speed and the diffusion coefficients at infinite dilution to estimate the film Péclet numbers of Pe1 ≈ 7 and Pe2 ≈ 33. (Throughout this discussion, 1 denotes the shorter polymer, while 2 is the longer polymer.) Initially, the polymers were nearly uniformly distributed in the film. When drying began, both the long and short polymers accumulated at the moving interface, as expected from their Péclet numbers. However, there was a noticeable difference in the distribution of chains within the film between the two models. Consistent with the previous simulations and the Howard et al. theoretical model,17 the implicit-solvent model stratified with the shorter polymers on top of the longer polymers [Fig. 1(b)]. In stark contrast, the explicit-solvent model showed no small-on-top stratification [Fig. 1(a)] and in fact more long polymers than short polymers accumulated immediately below the liquid–vapor interface.
The qualitatively different behavior of the two models suggests a significant difference in the relative migration speeds of the polymer chains, which was previously shown to give rise to stratification for the implicit-solvent model.12,15,17,22 At isothermal conditions, the relative migration speed Δu predicted by the Howard et al. model15 is
where u1 is the migration speed of the short polymers, μi is the chemical potential of component i computed from the equilibrium free-energy functional, and is the effective mobility for component i. We previously estimated from the equilibrium diffusion coefficient Di using the Einstein relation, , within the free-draining approximation of hydrodynamics.15,17 The Howard et al. model was previously shown to give quantitative agreement with implicit-solvent simulations of a polymer mixture.17 Stratification of the short polymers requires that Δu > 0, which occurs when the ratio of the chemical potential gradients is larger than the ratio of mobilities. Our simulations and Eq. (5) then suggest several possible explanations for the different drying behavior:
temperature gradients from evaporative cooling affect μi or the migration velocities in the explicit-solvent model, but such gradients are absent from the implicit-solvent model,
differences in the chemical potentials of the two models lead to different diffusive driving forces,
differences in the equilibrium diffusion coefficients change for the two models, and/or
hydrodynamic interactions mediated by the explicit solvent alter the migration velocities.
We probed each of these effects in turn. As will be shown, we found that only hydrodynamic interactions (D) had a significant enough effect to explain the lack of stratification in the explicit-solvent model.
Solvent evaporation can lead to cooling at fast drying rates.27 For our explicit-solvent model, the temperature is expected to be fixed at T = 1.0 ε/kB by the Langevin thermostat close to the substrate and, at pseudo-steady state, to decrease linearly in distance from the substrate with a slope proportional to the evaporative flux. The temperature profile T(z) measured in the simulations, shown in Fig. 2(a), is consistent with this expectation. Here, the local temperature T(z) is defined by
where the sum is taken over the Nk particles in a slab k centered around z with thickness 2Δz = 1 σ and is the speed of particle j. On the other hand, there were no temperature gradients in the implicit-solvent model because the solvent was treated as an isothermal viscous background.
Temperature gradients can give rise to mass flux (i.e., the Soret effect),36 and the mobility is temperature dependent through the viscosity.37 It has recently been shown that temperature gradients can influence the stratification of nanoparticles.38 In order to test how the presence of temperature gradients influenced stratification in the explicit-solvent model, we removed the Langevin thermostat and instead coupled all particles to a dissipative particle dynamics (DPD) thermostat32,39,40 to eliminate temperature gradients and performed a simulation at the same drying rate. The DPD thermostat applies pairwise random and dissipative forces to all particles and so conserves momentum and preserves hydrodynamic interactions. We set the DPD friction coefficient to be 0.5 m/τ, which has been shown to have a negligible effect on the viscosity of a fluid of nearly hard spheres.40 Figure 2(a) shows that although there was originally roughly 12% evaporative cooling in the explicit-solvent model, the DPD thermostat maintained a constant temperature throughout the film.
We then compared the distribution of polymers in the film with and without evaporative cooling at the same drying rate. The temperature gradient in the film led to a corresponding gradient in the local total density (solvent and polymer), which was higher in colder regions [Fig. 2(b)]. This in turn decreased the total film height compared to the isothermal case for an equal amount of the solvent evaporated. To account for this difference, we compared the monomer density profiles at equal film height. The monomer density profiles for the M1 = 10 and M2 = 80 polymers are shown in Figs. 2(c) and 2(d) at film height 360 σ, when the temperature gradient was most pronounced. These monomer density profiles were indistinguishable with and without evaporative cooling. Accordingly, we excluded temperature effects as a possible reason for the different stratification behavior in the implicit- and explicit-solvent models.
B. Chemical potential
In the absence of pressure and temperature gradients, diffusive flux is driven by chemical potential gradients within the regime of linear response.41 A mismatch in the density or composition dependence of the chemical potentials could result in different stratification behavior for the explicit- and implicit-solvent models. We accordingly measured the relevant chemical potentials for both models using test insertion methods. We first generated trajectories of polymer mixtures in cubic periodic simulation boxes over a wide range of mixture compositions (see the supplementary material), saving 100 independent configurations for each composition.
For the explicit-solvent mixtures, the chemical potential of the solvent μ0 was calculated by Widom’s insertion method,42
where β = 1/kBT, is the chemical potential of an ideal gas of particles at the same temperature and density as the solvent in the mixture, and ΔU is the change in the potential energy on insertion of a test particle. The ensemble average is taken over configurations and insertions. We performed 5 × 105 insertions per configuration, which we found to be sufficient to obtain a converged value for μ0.
The chemical potential of the polymers was estimated by the chain increment method43 for both the explicit- and implicit-solvent models. The incremental chemical potential μ+(M) to grow a chain of length M + 1 from a chain of length M is
where the ensemble average is now taken over configurations and test insertions onto the ends of chains of length M. The incremental chemical potential μ+(M) converges to a constant value μ+(∞) for sufficiently large M.43 The greatest deviations from μ+(∞) occur for a single monomer, μ+(0), which can instead be measured by Widom insertion.42 The chemical potential of the chain was then approximated as44
where is the chemical potential of an ideal gas of chains of length Mi at the same chain density and temperature as in the mixture.
In order to improve convergence of the ensemble averages, test particles were inserted onto the chain ends at positions consistent with the bond length distribution.44,45 The test particle position was drawn with random orientation relative to the chain end at distance r distributed according to r2 exp(−βϕ(r)) in the range 0.7 σ ≤ r ≤ 1.3 σ, where ϕ(r) is the interaction potential between bonded monomers given by Eqs. (1) and (2). An additional analytical contribution to the chemical potential was then included to account for the weighted sampling.44 To measure μ+(10) and μ+(80), we performed 100 trial insertions per chain end in each configuration and, finding that μ+(10) ≈ μ+(80) within statistical uncertainty, took μ+(∞) ≈ μ+(10). The monomer excess chemical potential, μ+(0), was measured using 5 × 105 random insertions per configuration, as for the solvent.
In our implicit-solvent theoretical model [Eq. (5)],15,17 the diffusive driving force was given by the gradient of μi. However, this expression must be modified to account for the presence of the solvent in the explicit model. If the mixture is approximately incompressible, the solvent and polymer concentrations are not independent. Assuming that the volume occupied by one polymer chain of length Mi is roughly equal to Mi solvent particles, the chemical potential μi is replaced in Eq. (5) by the exchange chemical potential Δμi = μi − Miμ0.46,47 For the implicit model, the solvent is effectively incorporated into the free energy of the polymers, giving Δμi = μi and recovering Eq. (5), as expected. Batchelor showed that the diffusion process resulting from gradients of Δμi is equivalent to sedimentation of the solute (polymers) in a force-free solvent.47,48
Although Δμi was determined across a wide range of mixture compositions (Fig. 3), we found that it was essentially only a function of the total monomer density ρm. As expected, Δμi increased for larger ρm. This is in qualitative agreement with our analysis of stratification in mixtures of hard-sphere chains, where the excess chemical potential per monomer depended primarily on the total monomer packing fraction and more weakly on the total chain number density.17 Most importantly, Δμi is in excellent agreement between the two models for both the short and long polymers. We then expect both the explicit- and implicit-solvent models to give equivalent driving forces for diffusion and stratification.
C. Equilibrium diffusion coefficient
Having found good agreement in the exchange chemical potential between the two models, we considered dynamic effects in Eq. (5). The polymer mobility relates the diffusive driving force on the polymers to the migration velocity. In our theoretical model, we assumed that could be estimated from the equilibrium diffusion coefficient Di. In our simulations, Di for the implicit-solvent model was matched to the explicit-solvent model in the limit of infinite dilution, but the models may deviate at higher concentrations. Differences in the concentration-dependence of Di, and hence , may accordingly influence the relative migration speeds of the components during drying.
We measured Di for both the explicit- and implicit-solvent models across a wide range of concentrations and compositions from the mean-squared displacement of the polymer centers of mass. The values of all measured diffusion coefficients are available in the supplementary material. At low total monomer densities, the agreement between the two models was good, as expected from the model fitting, with larger discrepancies at higher monomer densities. Over the range of compositions relevant to the drying simulations, the agreement between the two models was quite good, with a deviation of at most 20% for the short polymers and 30% for the long polymers relative to the explicit-solvent model.
According to Eq. (5), a larger value of , and so approximately D2/D1, should result in a larger difference in migration velocities for the same chemical potential gradient and lead to stronger stratification. For all concentrations and compositions investigated, D2/D1 was larger for the explicit-solvent model (Fig. 4). This suggests that stronger stratification might be expected for the explicit model on the basis of the diffusion coefficients, which is the opposite of the observed behavior. Accordingly, we concluded that differences in the equilibrium diffusion coefficients were not likely the reason for the lack of stratification in the explicit-solvent model.
D. Hydrodynamic interactions
Given that temperature gradients, chemical potentials, and equilibrium diffusion coefficients do not explain the difference in stratification in the implicit- and explicit-solvent models, it remains now to test the influence of hydrodynamic interactions. Hydrodynamic interactions are inherently present in the explicit-solvent model, but are absent from the implicit-solvent Langevin dynamics simulations. The recent analysis by Sear and Warren demonstrated how solvent backflow reduces the phoretic motion of a large colloid in an ideal polymer solution.23 Brady showed from a rigorous microscopic perspective how this fundamentally results from a neglect of hydrodynamic interactions, which modify the mobility.24
We performed a simple simulation to test the effect of hydrodynamic interactions on polymer migration. We constructed a simulation box with L = 40 σ and Lz = 80 σ. Half of the box (z > 0 σ) was filled with pure solvent, while the other half (z < 0 σ) was filled with a solution of long and short polymers. A semipermeable membrane at z = 0 σ separated the two compartments, and walls were placed at z = ±40 σ. Interactions with the walls were given by Eq. (3) truncated at distances greater than 3 σ with εw = 1 ε, while the monomer interactions with the membrane were modeled by Eq. (3) truncated at (2/5)1/6 σ with εw = 1 ε. The solvent was allowed to exchange across the membrane, coming to equal chemical potential, but the polymers were confined to z < 0 σ. For the explicit-solvent model, there were 78 720 solvent particles in the box and the initial total monomer density of the polymers in the compartment was ρm = 0.2 σ−3 (12 800 monomers) with half of the monomers belonging to each type of polymer (ρ1 = ρ2 = 0.1 σ−3). The total density of the solution in z < 0 σ roughly corresponded to the measured equilibrium density (see the supplementary material). The polymers in the explicit solvent were used as starting configurations for the implicit-solvent model.
We subsequently removed the membrane and allowed the polymers to diffuse. At infinite dilution, the polymers would be expected to diffuse at roughly the same rate in both models because their center-of-mass diffusion coefficients were matched at these concentrations. At finite concentrations, however, hydrodynamic interactions between the polymers, qualitatively associated with solvent backflow, may alter the diffusive flux out of equilibrium.24 Since the initial gradient in Δμi across the membrane should be comparable for both models based on Fig. 3, the diffusive driving force is the same in the explicit and implicit solvent, and any differences in the measured flux should be due to hydrodynamic interactions out of equilibrium. Figure 5 shows the fraction of monomers which moved into the opposite side of the box (z > 0 σ) at a given time, averaged over 25 independent simulations. At the steady state, half of the monomers, summed over both short and long polymers, should have z > 0 σ.
It is clear that the dynamics of the implicit-solvent model are faster than the explicit-solvent model, in agreement with the considerations in Refs. 23 and 24. This is also consistent with the measured total monomer density profiles in the drying films, computed from Fig. 1, that show more monomers near the liquid–vapor interface for the explicit-solvent model. Moreover, the long (M2 = 80) polymers and short (M1 = 10) polymers responded differently to the same density gradient between the two models. In the implicit-solvent model, the long polymers migrated faster than the short polymers, consistent with Eq. (5) and stratification in the drying simulations. On the other hand, the long polymers migrated much slower than the short polymers in the explicit-solvent model, explaining the lack of stratification in the drying simulations.
We demonstrated the influence of hydrodynamic interactions on stratification in a drying polymer mixture using explicit-solvent and implicit-solvent computer simulations. The implicit-solvent model predicted stratification in agreement with the previous simulations and theoretical considerations.17 However, no such stratification was found for the explicit-solvent model under the drying conditions considered. Despite good mapping of the equilibrium bulk properties (chemical potential, diffusion coefficient) between the explicit- and implicit-solvent models, hydrodynamic interactions out of equilibrium were shown to alter the polymer diffusion in a way that is consistent with the lack of stratification. Our analysis directly tests and confirms that implicit-solvent simulations and theoretical models lacking hydrodynamic interactions are not capable of quantitatively predicting stratification, in agreement with the analysis of Sear and Warren23 for the special case of a large colloid in an ideal polymer solution. In future, hydrodynamic interactions must be incorporated into any simulations aiming to study stratification, e.g., through explicit-solvent molecular dynamics or with an appropriate mesoscale simulation method.
While there are currently no experiments available for stratification in polymer mixtures, colloid mixtures have been shown to stratify in experiments.12–14 The extent of colloid stratification in the experiments12 and recent explicit-solvent computer simulations38 appears weaker than predicted by models lacking hydrodynamic interactions, consistent with our analysis. Larger chemical potential gradients on the larger component would be required to increase the stratification, which could be induced by, for example, additional cross-interactions between the components.
See supplementary material for the single-chain structure for both models and equilibrium bulk properties of mixtures.
We gratefully acknowledge the use of computational resources supported by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology’s High Performance Computing Center and Visualization Laboratory at Princeton University. Financial support for this work was provided by the Princeton Center for Complex Materials, a U.S. National Science Foundation Materials Research Science and Engineering Center (Award No. DMR-1420541).