THE JOURNAL OF CHEMICAL PHYSICS 139, 154102 (2013)

Hybrid approach combining dissipative particle dynamics and finite-difference diffusion model: Simulation of reactive polymer coupling and interfacial polymerization Anatoly V. Berezkin1,a) and Yaroslav V. Kudryavtsev2 1

Max-Planck Institut für Eisenforschung GmbH, Max-Planck str. 1, 40237 Düsseldorf, Germany Topchiev Institute of Petrochemical Synthesis, Russian Academy of Sciences, Leninsky Prosp. 29, 119991 Moscow, Russia 2

(Received 13 July 2013; accepted 26 September 2013; published online 16 October 2013) A novel hybrid approach combining dissipative particle dynamics (DPD) and finite difference (FD) solution of partial differential equations is proposed to simulate complex reaction-diffusion phenomena in heterogeneous systems. DPD is used for the detailed molecular modeling of mass transfer, chemical reactions, and phase separation near the liquid/liquid interface, while FD approach is applied to describe the large-scale diffusion of reactants outside the reaction zone. A smooth, selfconsistent procedure of matching the solute concentration is performed in the buffer region between the DPD and FD domains. The new model is tested on a simple model system admitting an analytical solution for the diffusion controlled regime and then applied to simulate practically important heterogeneous processes of (i) reactive coupling between immiscible end-functionalized polymers and (ii) interfacial polymerization of two monomers dissolved in immiscible solvents. The results obtained due to extending the space and time scales accessible to modeling provide new insights into the kinetics and mechanism of those processes and demonstrate high robustness and accuracy of the novel technique. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4824768] I. INTRODUCTION

Many practically important heterogeneous processes including chemical and biochemical reactions, phase transitions, self-assembly, and surface reconstruction can hardly be modeled at large space and time scales with a single method that combines both accuracy and high performance. This problem stimulated development of a number of hybrid simulation techniques, which make use of various methods in different regions of a model system and involve matching solutions at the boundaries between those regions. Implementation of hybrid methods enables reasonable allocation of computational resources so that only the most complex and interesting phenomena/regions are simulated with highly accurate but resource-consuming techniques, while the rest of the system is modeled with less detailed but robust approaches. For instance, coupling of quantum and molecular mechanics (QM/MM models) is widely used to study charge transport,1, 2 the structure of active centers in proteins and enzymatic reactions.3–6 Molecular simulations combined with finite element approaches are applied to investigate mechanical properties of materials (deformation and failure during nanoindentation, etc.).7–10 Linking molecular dynamics to finite difference (FD) mesoscale solvers, such as lattice Boltzmann method, allows consideration of hydrodynamic effects in suspensions or polymer solutions.11–13 Hybrids of Brownian dynamics (BD) with FD methods14–16 are successfully used to describe biochemical processes.17 a) Author to whom correspondence should be addressed. Electronic mail:

[email protected]

0021-9606/2013/139(15)/154102/12/$30.00

In this work we propose a novel hybrid technique combining dissipative particle dynamics (DPD) and FD diffusion model for efficient investigation of complex heterogeneous phenomena including diffusion and chemical reactions. There are two key differences between our model and similar BDFD hybrids of Refs. 14–17. First, DPD preserves the total number of particles and treats solvent explicitly, whereas BD can simulate variable number of particles in an implicit solvent. Thus, addition or removal of a solute in DPD models is less straightforward than in BD ones, but hydrodynamic interactions and entropy-driven effects, which are especially important in polymer systems, can be accounted for with DPD easier and more accurately. Another important difference between the current and preceding approaches is related to the role of a buffer region between FD and molecular model domains. In BD-FD hybrids, the FD domain is used as a reservoir maintaining Dirichlet,14–16 Neumann, or Robin17 boundary conditions for the molecular simulation box. We prefer a more flexible strategy, where boundary conditions are imposed at the outer boundary of the FD domain, far off the DPD box. Simultaneously, we use a self-consistent “sutureless” linkage of the molecular and FD domains providing free diffusion of species with minimum artifacts caused by near-wall effects and nonFickian diffusion of particles at small space and time scales. Simulations for model systems prove universality, accuracy, and numerical efficiency of the elaborated approach. The paper is organized as follows. Section II starts with describing the general structure of the DPD-FD model and also includes a specification of the constituting techniques in the domains of their use and an algorithm of their

139, 154102-1

© 2013 AIP Publishing LLC

154102-2

A. V. Berezkin and Y. V. Kudryavtsev

J. Chem. Phys. 139, 154102 (2013)

of end-coupling at a flat interface is replaced with the diffusion controlled one with decreasing the fraction of reactive chains. Such regularity was predicted in theory28 but up to date it was not clearly corroborated in laboratory or numerical experiments. In linear interfacial polymerization (IP), chains of alternating AB copolymer are obtained from monomers A and B that are dissolved in two immiscible solvents. The polymer product usually precipitates forming a thin film. This phenomenon underlies a wide range of IP applications, such as in situ preparation of membranes,29–35 microcapsules,36–39 nanocapsules,40, 41 and nanoparticles,42–44 especially interesting for biomedical applications and drug delivery. However, IP kinetics remains poorly understood being affected by a complex of diverse phenomena. The proposed approach makes it possible to observe the structure of a reacting system with molecular resolution by explicitly taking into account large-scale diffusion of reagents, solvent effects, amphiphilic nature of reaction products, their precipitation, and self-assembly. Hybrid simulations can help to understand how to control the film thickness, structure, and permeability through monomer architecture and polymerization conditions. II. HYBRID DPD-FD METHOD A. General scheme

FIG. 1. Model heterogeneous processes: (I) seminfinite diffusion of particles B transforming into solvent S upon contact with the A/B interface; (II) reactive coupling of two immiscible end-functionalized polymers A and B that results in a diblock copolymer AB formation; (III) interfacial polymerization of monomers A and B dissolved in immiscible solvents. Interfacial region is treated with DPD, while semi-infinite diffusion in the bulk is simulated within FD model.

coupling. Section III reports on the simulations of three interfacial reaction-diffusion processes that are schematically shown in Figure 1. The reaction A + B → A + S (I) is used as a test system because it admits an analytical solution for the concentration profiles of reacting species. Then, it is demonstrated that such complex processes as the reactive coupling of immiscible polymer melts (II) and interfacial polymerization (III) can be effectively simulated by the DPD-FD technique at low computational cost. Reactive coupling of immiscible polymers through in situ block copolymer formation at a polymer/polymer interface is widely used for manufacturing composite materials.18–23 Our recent simulations24–27 by DPD made it possible to delineate several kinetic regimes and determine the conditions of the interfacial instability, which leads to the blend emulsification. Though molecular modeling undoubtedly helped a lot in understanding details of the heterogeneous process, practically important partially functionalized systems could not be investigated because of the small size of the simulation box. In the present paper, the novel DPD-FD hybrid technique is applied to demonstrate that the reaction controlled kinetics

Consider a simple bilayer system I consisting of a melt of A particles and a solution of B particles, where solvent molecules (S particles) are treated explicitly. The melt and solution are immiscible, whereas the components B and S are fully miscible and physically indistinguishable. Any B particle is instantaneously transformed into S upon its contact with A, while S particles are chemically inert. This reaction induces semi-infinite diffusion of B particles from the solution to the A/B interface. The hybrid model structure is shown in Figure 2. Domain in the vicinity of the A/B interface is simulated with DPD covering the range from 0 to lz in the z-direction. DPD technique provides a detailed description of all reaction-induced phenomena and possible changes in the structure and permeability of the interfacial region, as well as in the mechanisms of diffusion and reactions therein. Diffusion of B particles from the bulk toward the interface is described with the FD model defined in the range zα ≤ z ≤ zmax , as is shown in Figure 2. The DPD and FD models overlap in the “buffer” region zα ≤ z ≤ lz , where their concentration fields are self-consistently related (see Subsection II D) to provide unhampered diffusion of the active component according to Fick’s diffusion law. B. Dissipative particle dynamics

DPD is the coarse-grained molecular dynamic technique, nowadays very popular in simulations of simple fluids, suspensions, and polymer solutions/melts/mixtures for its ability to reproduce large-scale Navier-Stokes hydrodynamics and relative simplicity of mapping onto the classical

154102-3

A. V. Berezkin and Y. V. Kudryavtsev

J. Chem. Phys. 139, 154102 (2013)

FIG. 2. Scheme of the hybrid simulation: (a) general overview and (b) enlarged buffer region. Yellow, red, and blue particles represent A, B, and S components, respectively. The broken black and blue lines correspond to the DPD number density of B particles, ρ B (z), and to the FD concentration of B component, cB (z), respectively. See text for details.

Flory-Huggins theory. The method was proposed by Hoogerbrugge and Koelman45, 46 for the simulation of liquid suspensions and extended to polymer systems by Espanol, Groot, and Warren.47, 48 DPD particles, each representing a group of small molecules or extensive molecular fragments, interact by conservative, dissipative, and random which are pair forces, R (FCij + FD wise additive. The net force fi = ij + Fij ) exerted

The equations of particle motion dri /dt = vi , dvi /dt = Fi are solved using the so called DPD-VV integration scheme49 (modified velocity-Verlet algorithm) with the time step δt = 0.04. Chemical reactions are modeled with the probabilistic approach24–26 as follows.

on a given ith particle is calculated by summation over all other particles within a certain cut-off radius rc . Let rc , particle mass m, and kB T be the unit distance, mass, √ and energy, respectively, thus defining the time unit τ = rc m/kB T . The conservative force represents the excluded volume interactions as well as the interactions of chemically bonded particles  i and j in the dimensionless form FCij = aij (1 − rij )rij − ks rij ,

1. Two active particles can react upon a contact, i.e., when the distance between their centers is smaller than the reaction radius rR (here rR = rc = 1 is taken). 2. Reaction in every pair of contacting active particles takes place every 10th time step with a certain predefined probability p0 yielding the reaction probability per unit time pR = p0 /(10δt). 3. Particle pairs are examined for the possible reaction in a random order.

where rij = ri – rj , rij = |rij |, rij = rij /rij , aij is a maximum repulsion between those particles attained at ri = rj , and ks is a spring constant taken to be equal to 4 for neighboring particles in a polymer chain and to zero for non-bonded particles. The dissipative and random forces, FD √ ij    = −γ ω(rij )2 (rij · vij )rij and FRij = σ ω(rij )rij ζ / δt, respectively, constitute together the Groot-Warren thermostat, where γ is a friction coefficient related to a thermal noise amplitude σ via the fluctuation-dissipation theorem, σ 2 = 2γ , ω(r) is a weight function, ζ is a normally distributed random variable with zero mean and unit variance, which is uncorrelated for different particle pairs, δt is the time step of an integration scheme, and vij = vi – vj is the relative velocity of ith and jth particles. Following Ref. 48 we choose σ = 3.0, ω(r) = 1 – r, the average density of particles ρ 0 = 3, and the repulsion parameter between similar particles aii = 25 (i = A, B, S). The same value used for aBS makes B and S components fully miscible, whereas A becomes highly immiscible with both of them by setting aAB = aAS = 80. The only difference between B and S is the chemical activity of B particles in the reaction with A ones. Note that the repulsion parameter for dissimilar particles aij is related to the interaction parameter χ ij of the Flory-Huggins theory as48 aij = χ ij /0.306 + 25.

In the system I, the instantaneous interfacial reaction A + B → A + S proceeds at p0 = 1, which means that it takes place at every collision of active particles. For the systems II and III, smaller reaction probabilities are used because of the experimental notion21, 22 that most of the polymer reactions involving formation and cleavage of covalent bonds are reaction controlled, at least on short time scales. Simulation box for the system I has the size lx × ly × lz = 25 × 25 × 50 rc 3 including A layer 10rc thick with 18 750 particles and B/S layer 40rc thick with 15 000 B and 60 000 S particles so that the initial number density of the reagent B is ρ B0 = 0.6. In order to minimize composition fluctuations, the initial distribution of B particles across the z-axis is leveled out. In the direction perpendicular to the zaxis and parallel to the A/B interface, the box is bounded by impenetrable walls made of particles similar to the particles of an adjacent phase. Wall-forming particles have infinitely large mass, so that they cannot move, but formally attributed non-zero velocities obeying the Maxwell-Boltzmann distribution, in order to avoid temperature drop near walls and the wall slippage effect. In the x and y directions, the periodic boundary conditions are imposed. For the other systems, the box size is different but the boundary conditions are the same.

j =i



154102-4

A. V. Berezkin and Y. V. Kudryavtsev

J. Chem. Phys. 139, 154102 (2013)

C. Finite-difference model of diffusion

In the FD domain, the concentration of B component, cB , is governed by one-dimensional diffusion along the z-axis: ∂cB (t, z) ∂ 2 cB (t, z) . =D ∂t ∂z2

(1)

In general, the diffusion equations can be written in terms of the chemical potentials and Onsager coefficients, or even the Nernst-Planck equation can be numerically solved, but here we limit ourselves to Eq. (1) for simplicity. Since the diffusivity D should be the same in the DPD and FD models, we adopt the value D = 0.3029 ± 0.0002 estimated for the DPD liquid of identical non-bonded particles. Equation (1) can be numerically integrated with standard methods. We apply here the unconditionally stable implicit second-order Crank-Nicolson discretization:50

tD cB (t + t, z) = cB (t, z) + [cB (t, z − z) 2( z)2 − 2cB (t, z) + cB (t, z + z) + cB (t + t, z − z) − 2cB (t + t, z) + cB (t + t, z + z)],

(2)

on a two-dimensional grid with the space step z = rc and time step t = δt using the sweep algorithm. The initial and boundary conditions read cB (0, z ≥ 10) = cB0 , cB (t, zα ) = ρB (t, zα ),

cB (t, zmax ) = cB0 = ρB0 ,

(3)

where cB0 is the initial reagent B concentration in the solution, which should be equivalent to the initial average number density of B particles, ρ B0 , in the DPD model. D. Linking the DPD and FD models

The DPD and FD models overlap in the buffer region zα ≤ z ≤ lz (Figure 2(b)), where the local number density of DPD particles, ρ B (t, z), is matched to the concentration cB (t, z) of the FD model in a self-consistent manner. To this end, the value ρ B (t, zα ), which is calculated from the DPD model by counting all B particles in the “sampling interval” [zα − δz/2; zα + δz/2], is substituted for the left boundary condition cB (t, zα ) in Eq. (3) before each FD integration step. In its turn, the distribution of B particles is periodically redefined within the “back-coupling interval” zβ ≤ z ≤ lz using the concentration field cB (t, z) taken from the FD model. The general algorithm of self-consistent matching the DPD and FD solutions is as follows. 1. (DPD step) Integrate the equations of motion from time t1 to t2 = t1 + t. 2. Calculate the local density of B particles, ρ B (t2 , zα ), within the sampling interval [zα − δz/2; zα + δz/2]. 3. (FD step) Integrate the diffusion Eq. (2) from time t1 to t2 with the boundary condition given by Eq. (3), where cB (t, zα ) = ρ B (t, zα ). 4. Use the function cB (t2 , z) from step 3 to adjust the distribution of B particles so that ρ B (t2 , z) ≈ cB (t2 , z) in

the back-coupling interval zβ ≤ z ≤ lz . This is done by considering all particles from the mentioned z-range in an arbitrary order and redefining their type according to the Metropolis rule. At given cB (t2 , z), the local fraction of B particles equals pB (t2 , z) = cB (t2 , z)/ρ 0 (do not confuse pB with ρ B ). Therefore, if a random number ζ , which is uniformly distributed from zero to unity, obeys the inequality ζ < pB (t, z), then a chosen particle with the coordinate z adopts type B, and type S otherwise. 5. Collect the data and go to step 1. In our model, it is assumed that only diffusive mass transfer takes place, while convection is negligible. Thus, only concentrations, not fluxes, should be matched in the buffer region. In a more general situation, our method requires extension in the spirit of the recent hybrid models,51, 52 where non-equilibrium molecular dynamics was coupled with continuum fluid dynamics. The described means of linking DPD and FD domains requires defining the widths of the buffer region lz − zα , sampling interval δz, and back-coupling interval lz – zβ . The last one flanks the wall of the DPD simulation box, whereupon the particle density there noticeably fluctuates (see the broken curve in Figure 2(a)). Therefore, the back-coupling interval, where the particle distribution is redefined, should be thick enough. On the other hand, the distance zα – zβ between the regions, where particles are inserted and where concentrations are measured, should be large compared to the characteristic distance, at which the diffusive regime of particle motion is established. This is especially important for large reactive molecules such as polymers. In Sec. III we discuss the choice of hybrid model parameters in more detail based on the simulation results. III. RESULTS AND DISCUSSION A. Simple heterogeneous reaction

Let us start with considering the system I, in which B particles are instantaneously transformed into S ones upon contact with A particles. The problem can be treated analytically in terms of the semi-infinite diffusion model. It means that the concentration of B particles obeys Eq. (1) with the following initial and boundary conditions: cB (0, z ≥ z0 ) = cB0 ,

t = 0,

cB (t, z0 ) = 0, cB (t, ∞) = cB0 ,

t > 0,

(4)

where z0 is the coordinate of the A/B interface. The solution reads 

z − z0 cB (t, z) = cB0 erf √ 2 Dt

 ,

2 erf(x) = π

x exp(−k 2 )dk. 0

(5) It is straightforward to found a time-dependent consumption rate of B component at the interface    dc D J =D = cB0 , (6) dz z=z0 πt

154102-5

A. V. Berezkin and Y. V. Kudryavtsev

J. Chem. Phys. 139, 154102 (2013)

FIG. 3. Concentration of B component along the z-axis at different moments of time given by Eq. (5) (the black curves) and calculated by hybrid simulations (the DPD and FD profiles are shown with red and blue, respectively). The model parameters are zβ – zα = 5, δz = 1, and lz – zβ = 10.

FIG. 5. Influence of varying one of the model parameters (shown) on the deviation between the simulated ns (t) and theoretical n(t) conversions for the system I.

which can be integrated over time to get the consumption of B particles per unit interfacial area by the time t as a measure of the reaction conversion  t Dt . (7) n(t) = J dt = 2cB0 π

and it well reproduces the analytical dependence n(t) given by Eq. (7), whereas the pure DPD simulation in the small box with lz = 50 predicts saturation caused by the exhaustion of B particles. Using the theoretical value of n(t) given by Eq. (7), a fractional standard deviation of the simulated conversion ns (t) can be expressed in terms of the integral function   t    ns (τ ) − n(τ ) 2 1 dτ , (8) σn (t) = t n(τ )

0

Simulations with the hybrid model were carried out within the time span t = 35 000. The constant density of B particles cB (zmax ) = ρ B0 = 0.06 was kept at zmax = 2000. In test runs, a standard set of the model parameters zβ – zα = 5, δz = rc = 1, and lz – zβ = 5 was used. As is seen from Figure 3, the obtained concentration profiles at different moments of time are close to the analytical solution given by Eq. (5) despite high fluctuations of the local concentration in the DPD model. Figure 4 demonstrates the dependence of the conversion n on square root time at different zβ – zα values. It is seen that the DPD-FD model is not very sensitive to this parameter

FIG. 4. Reaction conversion vs. square root time according to Eq. (7) (the black dashed curve) and simulated for the pure DPD model (gray dashed) and the DPD-FD model for different zβ – zα values (color). δz = 1, lz – zβ = 5.

0

where the integrand [(ns (τ ) − n(τ ))/n(τ )]2 is an instantaneous squared error of the conversion and the integration gives the value averaged over the time interval [0, t]. It was calculated for different values of the hybrid model parameters zβ – zα , δz, and lz – zβ (two of them were kept fixed at their standard values, while the remaining one was designated as x and varied) and plotted in Figure 5 for t = 35 000. It is seen that in all the cases the deviation varies from 1.3% to 2.5%. More detailed analysis reveals that it mostly comes from an initial stage of the simulation, when the conversion is low compared to fluctuations of the particle density and surface area in the DPD model. Figure 5 also demonstrates that a contact between the regions, where the concentration of particles is measured (the sampling interval) and then redefined (the back-coupling interval), i.e., zβ – zα = 0, produces unacceptably high deviation of 45%. In could be explained by short-time non-Fickian dynamics of the reactive particles. This raises the question, how large should be the distance zβ – zα for polymeric reactants? In an unentangled polymer melt, the diffusivity of chain ends coincides with that of the center of mass at times t τ R , where the Rouse time scales as τ R ∼ N2 . The time it takes for a chain to diffuse a distance of zβ – zα is tαβ ∼ (zβ – zα )2 /DN ∼ (zβ – zα )2 N, where DN is the chain diffusivity, which scales in the Rouse model as DN ∼ N−1 . If N is increased, one should keep the constant ratio tαβ /τ R to guarantee that the chain diffusion over the distance zβ – zα will follow Fick’s law with the diffusivity

154102-6

A. V. Berezkin and Y. V. Kudryavtsev

DN . To meet this condition, zβ – zα should scale as ∼N1/2 , i.e., similarly to the chain end-to-end distance, and always exceed it. The value zβ – zα = 5 satisfies those requirements in all our simulations for the systems II and III, which contain polymers. There are two possible origins of errors in calculating local concentrations within the sampling interval. The first one comes from stochastic distribution of B particles in space. If the sampling volume Vs = lx ly δz contains n = Vs ρ 0 particles in total, then the solute density obeys the binomial distribution with the mean ρ B and standard deviation σ1

= ρB (ρ0 − ρB )/(ρ0 lx ly δz), which decreases as δz grows. The second type of errors is produced by averaging the nonlinear concentration fields within the sampling interval. Here the deviation is proportional to the second derivative of concentration, i.e., σ 2 ∼ = δz2 (∂ 2 ρ B /∂z2 ) and it increases with δz. Thus, there should be an optimal δz value that minimizes the sum σ12 + σ22 . It follows from Figure 5 that the highest accuracy is attained when δz = 1, which coincides with the cutoff radius for all non-bonded DPD interactions, rc , and (incidentally) with the space step of the FD scheme, z. An example of the near-wall effect is visible in Figure 3, where the density of B particles ρ B drops down in the vicinity of the DPD box boundary at z ≈ lz . Since the total particle density does the same (a known DPD artifact), the local fraction of B particles, pB , which is used for redefining particle types within the back-coupling interval, is much less affected. That is why the system I kinetics weakly depends on the backcoupling interval width lz – zβ . In order to avoid polymeric near-wall effects, such as the entropy-driven accumulation of end groups at interfaces,26 it is worth taking the value of lz – zβ larger than the coil size, as is done here. Further we use the values zβ – zα = 5, δz = z, and lz – zβ = 5, which are called the standard set, which seems to be optimal for all the processes simulated in this study.

J. Chem. Phys. 139, 154102 (2013)

FIG. 6. Dependence of the copolymer coverage nc on the scaled time for different densities of the reactive end groups B. Full and dashed color lines (coinciding at ρ B0 = 0.60) describe the results for the DPD-FD and pure DPD models, respectively. Black dashed lines show the linear (reaction controlled) and square root (diffusion controlled) trends.

adjusting their distribution when linking the DPD and FD domains, equals pB (t, z) = cB (t, z)/ρ B0 . All the data were averaged over 8 independent runs. Results of the simulations are presented in Figures 6 and 7, where the interfacial coverage nc (here defined as the number of copolymer molecules per unit area of the initial interface) is plotted vs time. At the earliest stage the linear kinetics nc ∼ t is observed in Figure 6 indicating about the reaction controlled process. Then, the copolymer formation rate gets slower and, provided only a small fraction (10%) of B chains are functionalized, a diffusion controlled regime nc ∼ t1/2 sets

B. Reactive coupling

Here we consider the system II that contains two adjacent layers of melted, strongly immiscible (the repulsion parameter aAB = 50) polymers A and B. All chains consist of five DPD particles (NA = NB = 5) connected with Hookean springs of zero equilibrium length and the spring constant ks = 4. Each A chain has a functional group at one of its ends (ρ A0 = 0.6), while only part of B chains bear the same (ρ B0 is varied from 0.6 to 0.06). Irreversible coupling between functionalized A and B end groups at the interface takes place with the probability pR = 0.25 and results in the formation of AB diblock copolymer. Simulation boxes of two different sizes 30 × 30 × 30 rc 3 (with the A/B boundary at z0 = 10rc ) and 20 × 20 × 80 rc 3 (z0 = 20rc ) are used in the intervals from 0 to 104 and from 200 to 2 × 105 units of dimensionless time, respectively. Other parameters are as for the system I. Diffusivity of end groups in the FD model is assumed to be equal to that of the whole chain DB ≈ DN = 5 = 0.06942 ± 0.00003, which was found by the direct DPD simulation. The local fraction of reactive B particles, which is used for

FIG. 7. The same curves as in Figure 6 plotted in other coordinates to illustrate the saturation regime of end-coupling. Black dashed line shows the logarithmic (reaction controlled) trend.

154102-7

A. V. Berezkin and Y. V. Kudryavtsev

in. By formally inserting into Eq. (7) the numerical values of cB0 and D = DN = 5 for the current model, we can find the nc (t) dependence, which is plotted as a dashed straight line in the double logarithmic coordinates of Figure 6. Though Eq. (7) does not account for the copolymer brush at the interface, it perfectly reproduces the simulation trend and only slightly (by 8%) overestimates the simulated polymer coverage nc (t) plotted with the red line. On the opposite, the system with fully functionalized B chains demonstrates a transition from the linear to logarithmic kinetics nc (t) ∼ (ln t)1/2 , which can be explained by saturating the interface with the copolymer product that keeps the process under reaction control.53, 54 Thus, switching between diffusion and reaction controlled regimes in the polymer end-coupling is possible by varying the concentration of reactive end groups. One finds from Figures 6 and 7 that fast exhaustion of the reactive chains leads to kinetic artifacts (false plateaus) in the pure DPD model. This becomes more evident for the systems with low concentration of reactants, where largescale diffusion is simulated. Since the new hybrid model is free from those problems, its use for such situations is highly recommended. The blend, where all A and B chains are end-functionalized, reveals the same coupling kinetics for the both DPD-FD and DPD models thus admitting no doubt in our previous results obtained by the latter model.24–27

C. Interfacial polymerization

Modeling of IP has a long history reviewed in Ref. 55. It was shown that the polymerization kinetics is governed by a complex interplay of reaction and diffusion near the interface. Analytical solution of the corresponding infinite set of partial differential reaction-diffusion equations is impossible. Theories that are based on the linearization of concentrations in the reaction zone fail to correctly predict the monomer conversion, not to speak about the molecular mass distribution, without empirical fitting. As a consequence, the most accurate method is nowadays related to the numerical solution of partial differential reaction-diffusion equations. Its realization is most evident if polymer molecules are soluble and do not influence diffusivities and reactivities of each other. This so-called “phantom polymer approximation” was first introduced by Sharikov et al.,56 who studied a frequent situation, when IP takes place in an organic phase, while an adjacent aqueous phase only provides convective diffusion of a second monomer. Later one of us used the same assumptions in the detailed investigation of IP mechanism and macrokinetics.55 It was shown that long polymer chains can be formed only if the polymerization is very fast and the process is controlled by the diffusion of monomers in the organic phase. If it is limited by diffusion in the aqueous phase or the polymerization is slow, then IP yields oligomers only. The theory55 qualitatively reproduces numerous experimentally known57 peculiarities of IP: weak dependence of the average degree of polymerization (DP) on the molar ratio of monomers, high polydispersity index, which grows with DP, non-linear dependence of DP on time and monomer concentration, etc. Besides, it was shown

J. Chem. Phys. 139, 154102 (2013)

that monomer concentrations strongly fluctuate near the interface. There is a narrow region (reaction zone), where they are nearly equal and DP attains its maximum, whereas in the most part of the interfacial region there is a high excess of one of the monomers. The reaction zone expands and simultaneously drifts with time due to monomer diffusion and depletion. That is why the kinetics of IP cannot be correctly described with simple models assuming stationary mass transfer rate, linear or uniform distribution of reagents inside the reaction zone or polymer film. Since phantom polymer approximation becomes inaccurate if the polymer product precipitates (adsorbs) onto the interface, the existing FD techniques were recently extended by Freger and Srebnik58 and by Yashin and Balazs59 in order to take into account the mentioned effect. It was concluded that the polymer film at the interface creates a non-uniform barrier to monomer diffusion, the thickness of which is proportional to the total amount of polymer and its miscibility with monomers and solvents. Both the polymer density and DP decrease from the film center to its periphery. Despite a lot of valuable information provided by the pure FD models,55, 56, 58, 59 they have serious limitations that hinder realistic simulations of IP. One-dimensional FD models cannot correctly represent the structure of a dispersion formed during phase separation, while 2D and 3D models of IP are computationally too expensive. FD modeling also involves complicated problems of phase separation and interdiffusion in polydisperse multicomponent systems that are insufficiently studied in theory. In face of these obstacles, the natural next step was taken by Srebnik et al.,60, 61 who performed the first molecular simulation of IP, as molecular models naturally incorporate all mentioned phenomena thus avoiding the corresponding theoretical problems. However, the effects originating from a small size of the simulation box were not ruled out. Indeed, in Ref. 61 the concentration of monomers at the simulation box boundaries was kept constant according to the simple, but crude two-film Lewis and Whitman’s theory62 of mass transfer. The correctness of this assumption in the course of IP is questionable because a thickness of the diffusion zone is expected to grow due to mechanical stabilization of the interface by a polymer film, whereas in simulations61 it can only shrink as the film thickens in a finite-size box. Thus, it is hardly possible to expect realistic long-term polymerization kinetics from the model of Ref. 61. Our DPD-FD technique is free from such problems, since non-stationary semi-infinite diffusion of monomers is considered in the spirit of more comprehensive Higbie’s penetration theory.62 The hybrid model takes into account polymer precipitation and swelling in solvents, monomer sorption by the polymer film, and, with some modifications, can encompass glass transition and crystallization effects. Consider the system III (Figure 1), in which an alternating copolymer is formed from two immiscible bifunctional monomers A and B dissolved in two immiscible solvents SA and SB , respectively. Each monomer is represented by a single bifunctional particle, and solvent molecules are inert particles of the same size. The simplest reaction scheme can be written

154102-8

A. V. Berezkin and Y. V. Kudryavtsev

J. Chem. Phys. 139, 154102 (2013)

in the form 4pR

Am + Bm → Ae -Be , 2pR

R-Ae + Bm → R-A-Be , 2pR

Am + Be -R → Ae -B-R, pR

R1 -Ae + Be -R2 → R1 -A-B-R2 . Here indexes “m” and “e” denote monomers and polymer end groups, and R, R1 , R2 stay for arbitrarily long AB sequences. The reaction between end groups proceeds at their contact with a probability pR resulting in the formation of a new bond with the spring constant ks = 4 (or ks = 0, if a phantom polymer is simulated). Other reactions have larger probabilities because of twice higher functionality of monomers in comparison to end groups. In contrast to our previous FD models,55 we assume that the solubility of A particles in SB solvent is the same as that of B particles in SA solvent, therefore, there is no phase, where the reaction proceeds preferentially. This assumption enables comparing the results with those reported in Ref. 59, which is up to date the most accurate and relevant FD study of linear IP accompanied by polymer precipitation. More practical situation, when IP preferably proceeds in one of the contacting phases (usually it is an organic phase), can easily be simulated with our approach by a corresponding choice of the DPD repulsion parameters. Let particles A (B) be fully miscible with SA (SB ) and A (SA ) be highly immiscible with B (SB ), which is provided by the following values of the repulsion parameters: aAA = aBB = aASA = aBSB = 25, aAB = aSA SB = aASB = aBSA = 45, 55, or 65 that correspond to the Flory-Huggins parameter χ AB = 6.12, 9.18, or 12.24, respectively.48 The DPD periodic box of sizes 20 × 20 × 80 containing 96 000 particles at the initial densities of A and B monomers ρ A0 = ρ B0 = 0.6 in the corresponding phases is used. The reaction probability is pR = 0.025. Since the system is symmetric relatively to the interfacial plane, the origin of the z-axis coincides with the center of the DPD simulation box (z0 = 0) and the box boundaries are situated at z = ±lz /2 = ±40rc (positive sign corresponds to the phase B). Two FD domains flank the box at zα = ±30rc , where the sampling is performed, and the Dirichlet boundary conditions cA (t, −zmax ) = ρ A0 , cB (t, zmax ) = ρ B0 are imposed at the outer boundaries z = ±zmax = ±2000rc . IP is a multistage process evolving through several kinetic regimes that depend on monomer functionalities, concentrations, diffusivities, reaction rates, polymer properties, and other factors.55 At the very first stage of IP (at t ≤ 200), when monomer concentrations are nearly constant and equal to their initial values, the process is reaction controlled and qualitatively similar to a homogeneous ideal binary polymerization. For example, the conversion n, which is defined as the number of monomer units linked into polymer chains per unit interfacial area, linearly grows with time n ∼ pR ρA0 ρB0 t, as is seen in Figure 8.

FIG. 8. Time dependence of the monomer conversion n during IP at the DPD repulsion parameter aAB = 45 and 65 in the cases of (solid lines) precipitating and (color dots) phantom polymers.

The number-averaged (Pn ) and weight-averaged (Pw ) degrees of polymerization during the homogeneous ideal process read63 Pn = 2(1 + kcA0 t),

Pw /Pn = 2.

(10)

Figure 9 demonstrates that at t < 200 IP produces similar trends Pn − 2 ∼ Pw − 2 ∼ t, while the polydispersity index Pw /Pn does not exceed the value of 2. Thus, an assumption about the Flory molecular mass distribution of the IP product made by Yashin and Balazs59 stays reliable at the linear stage. However, it is worth mentioning that the linear stage of IP is rather short at the scale of the whole process and it yields only oligomers, as follows from Figure 9. Further kinetics strongly depends on the polymer architecture and solubility.55, 61 In the simplest case of phantom (non-precipitating) polymer, the reaction slows down

(9) FIG. 9. Time dependences of the number average, Pn , and mass average, Pw , DP of the copolymer AB formed during IP at different aAB values.

154102-9

A. V. Berezkin and Y. V. Kudryavtsev

because of the exhaustion of monomers near the interface. A crossover to the diffusion controlled regime should be√ gin when the diffusive flux at the interface ρB0 DB /(π t) (Eq. (6)) becomes comparable to the homogeneous reaction rate pR ρ A0 ρ B0 (Eq. (9)), whence one finds the char2 ) ∼ 400. As is seen from acteristic time tD ∼ DB /(πpR2 ρA0 Figure 8, the monomer conversion rate is actually decreased at that time. The color dotted n(t) curves asymptotically approach the square root trend, which can be found similarly to Eq. (7):



(11) n = 2ρA0 DA t/π + 2ρB0 DB t/π , where DA and DB are A and B monomer diffusivities in the corresponding solutions. In the case of precipitating polymer, monomer diffusion within the polymer film becomes a rate-controlling step and the IP kinetics gets even slower. The time when the precipitation becomes important, strongly depends on the incompatibility between comonomers: at aAB = 65 the solid and dotted blue curves diverge at t ≈ 200, whereas at aAB = 45

J. Chem. Phys. 139, 154102 (2013)

phantom and precipitating polymers lead to the same kinetics until t ≈ 2000. Since real polymer/solvent segregation is usually less than in our DPD model, the phantom polymer approximation described in Ref. 55 is expected to satisfactorily describe IP kinetics much longer than the linear regime lasts. In general, after the linear stage (at t > 200) the dependence of the number-average DP on time takes the form Pn ∼ tx , where x is less than unity being determined by the value of aAB (1/3 < x < 1/2 for 45 ≤ aAB ≤ 65 in Figure 9). Surprisingly, the linear dependence of Pw ∼ t for any aAB survives until t ∼ 7000, when monomer diffusion through the polymer film surely dominates the IP process. The polydispersity index considerably exceeds its homogeneous value, which is consistent with a number of classical experiments57 but in no way with the assumption of Yashin and Balazs.59 The detailed analysis reveals that a high polydispersity of the resulting polymer stems from the large fraction of oligomers, which was predicted in Ref. 55. Polymer accumulation at the interface is illustrated in Figure 10 for aAB = 45. It is seen that the retardation of

FIG. 10. Polymer and monomers concentration profiles across the initial interface at different moments of time. aAB = 45. Dotted lines correspond to the concentrations of polymer end groups.

154102-10

A. V. Berezkin and Y. V. Kudryavtsev

FIG. 11. Dependence of the local number average DP, Pn , on the zcoordinate at different moments of time (aAB = 45).

the kinetics at t ≈ 2000 (Figure 8) corresponds to the formation of a continuous copolymer film, which later thickens at nearly constant density. Polymer concentration in the film grows with aAB but its distribution within the film remains uniform. We suppose that leveling the film density becomes possible due to the high content of mobile oligomers. The local number average DP Pn (z) is calculated as an inverse sum of inverse chain lengths for all monomer units within the layer [z, z+1]. As is seen from Figure 11, the dependence of Pn on z is highly non-monotonous, it has a global maximum at the film center and rapidly decreases toward film boundaries. This fact reflects the difference in the apparent reaction rate across the film and qualitatively agrees with the results of Refs. 58–61. The surrounding solutions contain only small amounts of the most soluble short chains with Pn ≈ 3, which is seen from Figures 10 and 11. For applications, it is important to know the location of end groups within the polymeric membrane, prepared by IP,64 in order to predict its permeability. Corresponding distributions are shown in Figure 10 with the dashed lines. Generally, the concentration of end groups A (B) attains a maximum at the film boundary facing SA (SB ) solvent and gradually decreases to zero at the film center, similarly to the concentration of the corresponding monomer. Figure 12 presents the time dependences of the film thickness h, which is measured as a distance between the points where the polymer concentration reaches half of its maximum value. It is seen that after the film formation the power law h ∼ tα is established, with the exponent α varying from 0.505 to 0.571 for different values of the repulsion parameter aAB . When the smaller aAB is applied, the trend h ∼ t1/2 predicted in Ref. 59 is better reproduced. On the other hand, since the film density is nearly constant, its thickness has to be proportional to the amount of polymer. Indeed, for n > 20, h ≈ βn with β = 0.49, 0.44, and 0.42 at aAB = 45, 55, and 65, respectively, (for a dry polymer film, β = 1/ρ 0 = 0.333 is expected). Thus, a film consisting of linear poly-

J. Chem. Phys. 139, 154102 (2013)

FIG. 12. Dependence of the film thickness h on the simulation time. Values of the repulsion parameter aAB and power-law exponents are shown.

mers formed by IP retains the ability to grow with time, oppositely to a rather limited thickening of the crosslinked polymer film.58, 60, 61 Let us roughly evaluate the rate of IP at its late stage, when the reaction rate is controlled by monomer diffusion within the polymer film. As we already know, the film is nearly uniform so that its thickness is proportional to the amount of polymer, i.e., h = βn. Moreover, as is seen from Figure 10, at t ≥ 5000 the monomer concentrations everywhere outside the film become close to their bulk values ρ A0 , ρ B0 , while inside the film they decrease to zero values at the film center. Let us denote the monomer concentrations inside the film near its boundaries as ρ Af , ρ Bf , the monomer diffusivities in the film as DAf , DBf and introduce the distribution coefficients KA = ρ Af /ρ A0 , KB = ρ Bf /ρ B0 relating the monomer concentrations in the solution and film. The rate of polymer formation in the diffusion controlled regime reads ∂ρA ∂ρB ∂n = DAf + DBf . ∂t ∂z ∂z

(12)

Since A and B monomers react in an equimolar ratio, their diffusion fluxes are equal DAf ∂ρA /∂z = DBf ∂ρB /∂z.

(13)

Assuming a linear decrease of the concentrations within the film (this is a very crude approximation, which is used here to evaluate the mass transfer rate, but is inapplicable for the calculations of DP) one can rewrite Eq. (12) in the form ρAf ρBf ∂n = DAf + DBf , ∂t hA hB

(14)

where hA + hB = h and hA /hB = DAf ρ Af /(DBf ρ Bf ), since DAf ρ Af /hA = DBf ρ Bf /hB according to Eq. (13).

154102-11

A. V. Berezkin and Y. V. Kudryavtsev

After simple transformations we get ∂n 2 DAf KA ρA0 + DBf KB ρB0 = ∂t h 2  DAf KA ρA0 + DBf KB ρB0 . = βn

J. Chem. Phys. 139, 154102 (2013)

(15)

Integration of Eq. (15) yields the dependence of the conversion on time

(16) n = 2 (DAf KA ρA0 + DBf KB ρB0 )t/β. Equation (16) describes the diffusion controlled regime, which is revealed at t ≥ 104 in Figures 8 and 9. Possible glass transitions, polymer crosslinking, or crystallization can additionally decrease permeability of the film and promote earlier establishing of this regime in laboratory experiments. Summarizing, it is worth noting that the diffusion control over the kinetics in the considered model system can stem either from low permeability of the polymer film or from depletion in the monomer concentrations in solutions outside it. As is seen from Figure 10, in our simulations these effects are competing because the depletion, which is easily seen at t = 500, becomes much less pronounced with the polymer film thickening at t = 5000. As a result, no clear n ∼ t1/2 trend is found for the solid kinetic curves in Figure 8 until the very late stage, when the mature polymer film dominates the process. In practice one can also consider the dependencies of n on the initial monomer densities ρ A0 , ρ B0 , which are square root in Eq. (16) and linear in Eq. (11), in order to find whether the rate of the whole IP process is controlled by diffusion in the film or solutions. For example, the experimental data on microcapsule preparation by interfacial polymerization of diamines were successfully interpreted65 in terms of the simple theoretical model that included a scaling relation similar to Eq. (16), which means that the late stage of the process was fully controlled by diffusion in the polymer film. IV. CONCLUSIONS

We proposed the novel hybrid of dissipative particle dynamics and finite-difference diffusion model (DPD-FD) that describes heterogeneous reaction-diffusion phenomena with molecular resolution and simultaneously treats the large-scale reagent diffusion via computationally cheap numerical solution of partial differential equations. Despite its simplicity, the new approach is well applicable to studying complex and practically important processes of polymer synthesis and modification. For the macromolecular reaction of coupling between two partially end-functionalized immiscible polymers, it was shown that at the intermediate stage of the reaction either diffusion or reaction controlled kinetics can be observed depending on the fraction of reactive end groups, in agreement with earlier theoretical predictions.28 For the linear interfacial polymerization of two incompatible monomers dissolved in immiscible solvents, the molecular simulation, which accounts for semi-infinite diffusion of monomers toward the interface and precipitation of the polymer product onto it, was carried out for the first time. After

a very short initial stage, when the process proceeds like in a homogeneous system, the IP kinetics can pass through two diffusion controlled regimes, where the limiting processes are monomer diffusion through the solution bulk and through the polymer film near the interface, respectively. The polymer forms a film of nearly uniform density at the interface, which is thickened proportionally to the conversion. The new technique provides an opportunity to study practically important details of the local film structure. The mass average polymerization degree grows linearly with time much longer that the number average DP does. As a result, the molecular mass distribution reveals a high polydispersity, in accordance with the published experimental data.57 Finally, we would like to emphasize that the proposed hybrid approach is in no way limited to polymeric systems only. Moreover, it could be extended to couple the finite element method with more detailed coarse-grained molecular or even atomistic models.

ACKNOWLEDGMENTS

The work was supported by the Russian Foundation for Basic Research (project 12-03-00817a). We are also grateful for the opportunity to use computational cluster Chebyshev at the Moscow State University. 1 L.

Wang, D. Beljonne, L. Chen, and Q. Shi, J. Chem. Phys. 134, 244116 (2011). 2 S. Y. Kim and S. Hammes-Schiffer, J. Chem. Phys. 124, 244102 (2006). 3 J. Aaqvist and A. Warshel, Chem. Rev. 93, 2523 (1993). 4 F. Ahmadi, N. Jamali, S. Jahangard-Yekta, B. Jafari, S. Nouri, F. Najafi, and M. Rahimi-Nasrabadi, Spectrochim. Acta, Part A 79, 1004 (2011). 5 P. M. Zimmerman, M. Head-Gordon, and A. T. Bell, J. Chem. Theory Comput. 7, 1695 (2011). 6 J. Kang, Y. Hagiwara, and M. Tateno, J. Biomed. Biotechnol. 2012, 236157 (2012). 7 S. Ogata, E. Lidorikis, F. Shimojo, A. Nakano, P. Vashishta, and R. K. Kalia, Comput. Phys. Commun. 138, 143 (2001). 8 R. E. Miller and E. B. Tadmor, MRS Bull. 32, 920 (2007). 9 R. E. Miller and E. B. Tadmor, Modell. Simul. Mater. Sci. Eng. 17, 053001 (2009). 10 N. Bernstein, J. R. Kermode, and G. Csanyi, Rep. Prog. Phys. 72, 026501 (2009). 11 G. Giupponi, G. De Fabritiis, and P. V. Coveney, J. Chem. Phys. 126, 154903 (2007). 12 A. Ladd, J. Fluid Mech. 271, 285 (1994). 13 P. Ahlrichs and B. Dünweg, J. Chem. Phys. 111, 8225 (1999). 14 T. Geyer, BMC Biophys. 4, 7 (2011). 15 T. Geyer, C. Gorba, and V. Helms, J. Chem. Phys. 120, 4573 (2004). 16 C. Gorba, T. Geyer, and V. Helms, J. Chem. Phys. 121, 457 (2004). 17 P. Bauler, G. A. Huber, and A. J. McCammon, J. Chem. Phys. 136, 164107 (2012). 18 Reactive Extrusion, edited by M. Xanthos (Hanser Publ., Munich, 1992). 19 Polymer Blends, edited by D. R. Paul and C. B. Bucknall (Wiley, New York, 2000). 20 C. Koning, M. Van Duin, C. Pagnoulle, and R. Jerome, Prog. Polym. Sci. 23, 707 (1998). 21 A. D. Litmanovich, N. A. Platé, and Y. V. Kudryavtsev, Prog. Polym. Sci. 27, 915 (2002). 22 C. W. Macosko, H. K. Jeon, and T. R. Hoye, Prog. Polym. Sci. 30, 939 (2005). 23 G. J. Jiang, H. Wu, and S. Y. Guo, Polym. Eng. Sci. 50, 2273 (2010). 24 A. V. Berezkin and Y. V. Kudryavtsev, Macromolecules 44, 112 (2011). 25 D. V. Guseva, Y. V. Kudryavtsev, and A. V. Berezkin, J. Chem. Phys. 135, 204904 (2011). 26 A. V. Berezkin, D. V. Guseva, and Y. V. Kudryavtsev, Macromolecules 45, 8910 (2012).

154102-12 27 A.

A. V. Berezkin and Y. V. Kudryavtsev

V. Berezkin and Y. V. Kudryavtsev, Macromolecules 46, 5080 (2013). O’Shaughnessy and U. Sawhney, Phys. Rev. Lett. 76, 3444 (1996). 29 V. Enkelmann and G. Wegner, Appl. Polym. Symp. 26, 365 (1975). 30 J. E. Cadotte, R. J. Petersen, R. E. Larson, and E. E. Erickson, Desalination 32, 25 (1980). 31 B. H. Jeong, E. M. V. Hoek, Y. S. Yan, A. Subramani, X. F. Huang, G. Hurwitz, A. K. Ghosh, and A. Jawor, J. Membr. Sci. 294, 1 (2007). 32 N. Y. Yip, A. Tiraferri, W. A. Phillip, J. D. Schiffman, and M. Elimelech, Environ. Sci. Technol. 44, 3812 (2010). 33 W. J. Lau, A. F. Ismail, N. Misdan, and M. A. Kassim, Desalination 287, 190 (2012). 34 K. P. Lee, T. C. Arnot, and D. Mattia, J. Membr. Sci. 370, 1 (2011). 35 C. Feng, K. C. Khulbe, and T. Matsuura, J. Appl. Polym. Sci. 115, 756 (2010). 36 P. W. Morgan, J. Macromol. Sci., Chem. A 15, 683 (1981). 37 R. Arshady, J. Microencapsul. 6, 13 (1989). 38 T. L. Whateley, in Microencapsulation: Methods and Industrial Applications, edited by S. Benita (Marcel Dekker, New York, 1996). 39 H. Uludag, P. De Vos, and P. A. Tresco, Adv. Drug Delivery Rev. 42, 29 (2000). 40 C. Mayer, Int. J. Artif. Organs 28, 1163 (2005). 41 P. Couvreur, G. Barratt, E. Fattal, P. Legrand, and C. Vauthier, Crit. Rev. Ther. Drug Carrier Syst. 19, 99 (2002). 42 M. Barrere and K. Landfester, Polymer 44, 2833 (2003). 43 L. Torini, J. F. Argillier, and N. Zydowicz, Macromolecules 38, 3225 (2005). 44 J. P. Rao and K. E. Geckeler, Prog. Polym. Sci. 36, 887 (2011). 45 P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett. 19, 155 (1992). 28 B.

J. Chem. Phys. 139, 154102 (2013) 46 J.

M. V. A. Koelman and P. J. Hoogerbrugge, Europhys. Lett. 21, 363 (1993). 47 P. Espanol and P. B. Warren, Europhys. Lett. 30, 191 (1995). 48 R. D. Groot and P. B. Warren, J. Chem. Phys. 107, 4423 (1997). 49 G. Besold, I. Vattulainen, M. Karttunen, and J. M. Polson, Phys. Rev. E 62, R7611 (2000). 50 J. Crank and P. Nicolson, Proc. Cambridge Philos. Soc. 43, 50 (1947). 51 T. Werder, J. H. Walther, and P. Koumoutsakos, J. Comput. Phys. 205, 373 (2005). 52 R. Delgado-Buscalioni, K. Kremer, and M. Praprotnik, J. Chem. Phys. 131, 244107 (2009). 53 B. O’Shaughnessy and U. Sawhney, Macromolecules 29, 7230 (1996). 54 G. H. Fredrickson and S. T. Milner, Macromolecules 29, 7386 (1996). 55 A. V. Berezkin and A. R. Khokhlov, J. Polym. Sci., Part B: Polym. Phys. 44, 2698 (2006). 56 Yu. V. Sharikov, M. I. Fedotova, and L. B. Sokolov, Vysokomol. Soedin., Ser. A 15, 982 (1973). 57 P. W. Morgan and S. L. Kwolek, J. Polym. Sci. 40, 299 (1959). 58 V. Freger and S. Srebnik, J. Appl. Polym. Sci. 88, 1162 (2003). 59 V. V. Yashin and A. C. Balazs, J. Chem. Phys. 121, 11440 (2004). 60 R. Nadler and S. Srebnik, J. Membr. Sci. 315, 100 (2008). 61 R. Oizerovich-Honig, V. Raim, and S. Srebnik, Langmuir 26, 299 (2010). 62 R. B. Bird, W. E. Stewart, and E. N. Lightfoot, Transport Phenomena, 2nd ed. (John Wiley & Sons, New York, 2007). 63 P. J. Flory, Principles of Polymer Chemistry (Cornell University Press, Ithaca, NY, 1953). 64 A. A. Hussain, S. K. Nataraj, M. E. E. Abashar, I. S. Al-Mutaz, and T. M. Aminabhavi, J. Membr. Sci. 310, 321 (2008). 65 D. Poncelet, B. Poncelet, and R. J. Neifeld, J. Membr. Sci. 50, 249 (1990).

The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp