4604
J. Phys. Chem. B 2008, 112, 4604-4612
Effects of Polymer Chain Length and Stiffness on Phase Separation Dynamics of Semidilute Polymer Solution Yung-Hsu Wu,† Da-Ming Wang,*,† and Juin-Yi Lai‡ Department of Chemical Engineering, Institute of Polymer Science and Engineering, National Taiwan UniVersity, Taipei, Taiwan 10617, and Research and DeVelopment Center for Membrane Technology, Department of Chemical Engineering, Chung Yuan UniVersity, Chungli, Taiwan 320 ReceiVed: June 7, 2007; In Final Form: January 28, 2008
Three-dimensional dissipative particle dynamics (DPD) simulations were performed to investigate the phase separation dynamics of semidilute polymer solutions with different polymer chain length and stiffness. For the polymer solution composed of shorter and more flexible chains, a crossover of the domain growth exponent from 1/3 to 2/3 was observed during the course of phase separation, indicating that the growth mechanism altered from diffusion to interfacial-tension driven flow. When the chain flexibility was kept the same but the chain was lengthened to allow for the chain entanglement to occur, the growth exponent changed to 1/4 in the diffusion-dominating coarsening regime while the growth exponent remained 2/3 in the flow-dominating regime. When the chain length was kept short but the stiffness was increased, the growth exponent became 1/6 in the diffusion-dominating regime and little effect was observed in the flow-dominating coarsening regime. The slow down of the phase separation dynamics in the diffusion-dominating coarsening could be explained by that the polymer chains could only perform wormlike movement when chain entanglements occurred or when the chain motion was limited by chain stiffness during phase separation. Moreover, when both the effects of chain length and stiffness were enhanced, polymer networks composed of longer and stiffer chains appeared and imposed an energy barrier for phase separation to occur. As a result, the polymer solution with stiffer and longer chains required a larger quench depth to initiate the phase separation and caused the delay in crossover of the coarsening mechanism from diffusion to flow.
I. Introduction Phase separation dynamics has received considerable research attention1-9 due to its significant impacts on the final morphology of a phase-separated polymer solution. Before phase separation reaches its equilibrium state, the interfacial energy between the two phase-separated domains is dissipated by the coarsening process,10 which plays a vital role in determining the final morphology. To characterize the mechanism of coarsening, a power law of domain growth,1 R ) Ktn, was introduced to quantify the relationship between the domain size (R) and the elapsed time step of phase separation (t), where K is the rate constant and n the growth exponent. The coarsening process is dominated by diffusion in the early stage and by interfacial-tension driven flow in the later stage.3 Lifshitz and Slyozov1 (hereafter referred to as L-S) showed that the growth exponent was 1/3 in the diffusion-dominating stage. For the domain growth in the flow-dominating stage, Siggia3 advocated the domain growing as R ∝ t1 due to the balance between viscous stress and interfacial tension, whereas Furukawa4 proposed the growth exponent to be 2/3 in consideration of the balance between interfacial tension and inertia stress. An important research topic regarding to phase separation dynamics is how the chain entanglement influences the dynamics of phase separation.11-13 A polymer chain can only move in reptation when it is entangled with other chains.14 Using the tube model11 to elucidate the entanglement effect, de Gennes * Corresponding author. E-mail address:
[email protected]. † National Taiwan University. ‡ Chung Yuan University.
concluded that the growth exponent (n) was 1/5 when the domain size was much smaller than the end-to-end distance of the polymer chain, and the domain growth followed the L-S mechanism, namely n ) 1/3, when the domain size was much larger than the end-to-end distance. Similarly, using the tube model11 and treating entangled polymers as a highly stretched network, Muratov13 analyzed theoretically the transient gel discussed in the Tanaka’s experiment15 and derived a growth exponent of 1/ with the assumption that the domain size was smaller than 6 the distance between the entanglement points. Simulation techniques were also employed to investigate the phase separation dynamics of polymer solution16-18 and blends.19,20 By using a two fluid model,16 it was shown that the viscoelastic property of a polymer solution remarkably influenced the lifetime of the network structure in the solution, which was in good agreement with the experimental results.21,22 In addition, molecular dynamics (MD)17 and Monte Carlo (MC)18 simulations have been performed to investigate the domain growth during phase separation. The results indicated that the domain growth was suppressed by the polymer-rich belt bridging between clusters, resulting in a growth exponent of 1/4 during the early stage of phase separation. Recently, researchers have begun to adopt the mesoscale simulation to investigate phase separation dynamics. For example, the mesoscale simulation based on the timedependent Ginzburg-Landau (TDGL) equation was carried out and showed that polymer entanglements influenced the coarsening of polymer blends only in the early stage but not in the later stage.19 In the present work, a mesoscale simulation technique bridging the gap between microscopic and mesoscopic simula-
10.1021/jp074403y CCC: $40.75 © 2008 American Chemical Society Published on Web 03/27/2008
Chain Length and Stiffness Effects on Phase Seperation Dynamics tions, dissipative particle dynamics (DPD),23-25 was used to provide insight into the effect of chain length and chain stiffness on phase separation dynamics. The DPD offers conformational analysis of phase separation, and its physical mapping to the Flory-Huggins theory has been shown to be nontrivial.25 The DPD techniques have been successfully utilized to investigate the rheology of polymer solution23 and the phase separation of polymer mixtures26 and viscous fluids.27 It is applied here to study the thermally induced phase separation of polymer solution. In addition to the effect of chain entanglement on coarsening dynamics, another topic investigated in this work is how the polymer chain stiffness affects phase separation. As mentioned above, the chain entanglement affects the phase separation dynamics by confining the chains motion to only reptation. The increase in stiffness of polymer chains can also limit the chains motion via the reptation mechanism.28-31 Therefore, the chain stiffness should also play a role in the coarsening dynamics, which has not yet been clearly understood.
J. Phys. Chem. B, Vol. 112, No. 15, 2008 4605
FDij ) - γwD(rij)
( )
rij rij ‚vij rij rij
(4)
where γ is interpreted as the friction coefficient, wD(rij) is a rij-dependent weight function, and vij is the velocity difference between particles i and j, vij ) vi - vj. The random force FRij results from the thermal fluctuation and is expressed as
FRij ) σwR(rij)
rij θ rij ij
(5)
where σ is the noise amplitude, wR(rij) is a rij-dependent weight function, and θij(t) is a white noise with Gaussian statistic: 〈θij(t)〉 ) 0, and 〈θij(t),θkl(t′)〉 ) (δikδjl + δilδjk)δ(t - t′). The random force FRij can also be expressed in terms of a short interval of time ∆t
FRij ) σwR(rij)
II. Simulation methods
rij ζij rij x∆t
(6)
1. DPD Theory. The dissipative particle dynamics (DPD) method introduced by Hoogerbrugge and Koelman23,24 is a coarse-grained and particle-based simulation technique. In the DPD simulation, hundreds of molecules are grouped into a particle. Therefore, the DPD method can provide analysis in larger length and time scales in comparison with the molecular dynamics (MD) simulation. The DPD method combines the basic concepts of MD and takes the hydrodynamic behavior into account.23 In this approach, Newton’s law of motion is used to describe the movement of particles.
where ζij is a random number with zero mean and unit variance.25 In order to control the temperature via the fluctuationdissipation theorem, the weight functions wD(rij) and wR(rij), the friction coefficient γ, and the noise amplitude σ obey the following relationships:
dri dvi ) Vi, mi ) fi dt dt
σ2 ) 2γkBT
(1)
where ri and vi denote the position and velocity of particle i, and mi and fi are respectively the mass of the particle and the total force acting on it. The forces exerted on each particle by other particles in DPD include the conservative force (FCij ), the dissipative force (FDij ), and the random force (FRij ).32 Thus, fi ) Σj*i(FijC + FijD + FijR). In order to save computation time, the three forces between neighboring beads were neglected when their distance is beyond the cutoff radius rc.25 The forms for the three forces proposed in previous works25,32 are briefly described below. The conservative force FCij is given by
FCij
)
{
( )
aij 1 0
rij rij if rij < rc rc rij if rij g rc
(2)
where rij ) ri -rj, rij ) |rij|, rc is the cutoff radius, and aij denotes a maximum repulsion between particles i and j (i * j). The mapping of aij to the Flory-Huggins interaction parameter (χij) has been established,25 and when the reduced particle density F)3, aij can be expressed as
aij )
75 + 3.497χij F
(3) -3.
where aij is reduced by kBT/rc and F is reduced by rc The dissipative force FDij denoting the drag force between particles is given by
wD(rij) ) wR(rij) ) 2
{( ) 1-
0
rij rc
2
if rij < rc
(7)
if rij g rc (8)
DPD simulations are performed at reduced quantities, and rc, mi, and kBT are used to define the reduced length, mass, and energy. In the present study, the reduced particle density (F) was set to be 3 and a friction coefficient (γ) of 4.5 and a noise amplitude (σ) of 3 were chosen, so that the simulations were performed at reduced temperature equal to 1. In addition, periodical boundary conditions were used in our simulations. To describe the constraint between the bonded segments in polymer chains, an extra spring force is introduced.25 Hookean chain is usually assumed and a constant (C) is used to describe the flexibility of polymer chains.
FSij ) Crij
(9)
The chain stiffness can be varied by using different values for C. A modified velocity-Verlet algorithm25 is usually used in DPD, integrating the equations of motion to obtain the positions and velocities of all particles in each time step.
1 ri(t + ∆t) ) ri(t) + ∆tWi(t) + (∆t)2fi(t) 2 v˜ i(t + ∆t) ) Wi(t) + λ∆t fi(t) fi(t + ∆t) ) fi(r(t + ∆t), W˜ (t + ∆t)) 1 Wi(t + ∆t) ) Wi(t) + ∆t(fi(t) + fi(t + ∆t)) 2
(10)
In the present work, the simulations were carried out with λ ) 0.5 and a time step (∆t) of 0.04. With the time integration, the
4606 J. Phys. Chem. B, Vol. 112, No. 15, 2008
Wu et al.
TABLE 1: The Number of Beads Per Polymer Chain (Nb), Spring Constant (C), and Interaction Parameter (χ) Used in the Six Simulation Cases Nb C χ
case (a)
case (b)
case (c)
case (d)
case (e)
case (f)
10 2 1.43
100 2 1.22
10 10 1.43
100 10 1.22
100 10 1.33
100 2 1.33
trajectory of each particle in the simulation box was determined and the morphology evolution was thus obtained. 2. Simulation Details. Polymer solutions with different polymer chain length and stiffness were simulated to investigate their phase separation dynamics. For all the simulation cases, the initial concentration of polymer beads (φo) were set to be 0.3 and the number of beads per polymer chain (Nb) was set to either 10 or 100. The overlap threshold concentration33 (φ*) in athermal solution was estimated by φ* ) Nb-4/5, and thus φ* ) 0.16 for Nb ) 10 and φ* ) 0.03 for Nb ) 100. Since φ0 is larger than φ* for all simulation cases, the polymer solutions in the present study were all semidilute. In case (a), the polymer chain was composed of ten connected beads (Nb ) 10) with a spring constant C ) 2. In order to investigate the effect of chain length, case (b) simulated longer polymer chains (Nb ) 100) with the same spring constant (C ) 2) as in case (a). In case (c), the polymer chain length was the same as case (a) (Nb ) 10), while the spring constant was set as C ) 10 to study the effect of chain stiffness. To study the combined effects of chain stiffness and length, the spring constant in case (d) was set as C ) 10 and the chain length was Nb ) 100. The same quench depth (χ - χb) ) 0.55 was kept in cases (a)-(d), where χ was the interaction parameter between the polymer and solvents in the simulation, and χb is the smallest interaction parameter needed to initiate the phase separation of the polymer solution when the polymer concentration is 0.3, which was calculated by using the Flory-Huggins theory.34 Since χb ) 0.88 for Nb ) 10 and χb ) 0.67 for Nb ) 100, χ was set as 1.43 in cases (a) and (c) and 1.22 in cases (b) and (d). However, it was found that, for case (d), the interaction parameter of 1.22 was not enough to bring about phase separation, indicating that merely the Flory-Huggins theory cannot appropriately describe this simulation case. In order to study the phase separation dynamics for Nb ) 100 and C ) 10, an interaction parameter of 1.33 was adopted (case (e)), which was large enough to initiate the phase separation of the polymer solution with stiffer and longer chains. Besides, to study the chain stiffness effect on phase separation dynamics for the longer chain case (Nb ) 100), case (f) was simulated with the same Nb (100) and χ (1.33) as case (e), but with different spring constant (C ) 2). The simulation parameters (Nb, C, χ) for the six cases are listed in Table 1 for comparison. The polymer solution was modeled in a three-dimensional box with a size of 40 × 40 × 40 containing 134 400 solvent particles in each case, 5760 polymer chains in cases (a) and (c), and 576 polymer chains in cases (b), (d), (e), and (f). In all the cases, the athermally homogeneous polymer solution models (χ ) 0) were first prepared at the cost of 5 × 104 simulation steps, and phase separation was then generated by setting χ to the values listed in Table 1. According to eq 3, the reduced repulsion parameter between polymer (denoted as p) and solvents (denoted as s), can be estimated with aps ) 25 + 3.497χ. The phase separation dynamics was investigated with the information about how the polymer and solvent particles moved in the simulation box. DPD simulations were performed
on SGI Origin 3800 workstations with software Cerius2 4.9 (Accelrys Inc.). The mechanism of coarsening during phase separation is usually characterized by the growth rate of phase-separated domains. In the present work, the domain size of the polymerrich region was determined by analyzing the distribution of polymer concentration in the simulation box. The threedimensional simulation box contained totally 40 × 40 × 40 cells, with 40 divisions in each dimension. The polymer concentration was calculated by the number ratio of polymer particles to total particles. And the polymer-rich domain was defined as a region where the polymer concentration (φp) was higher than the initial value (φp g 0.3) during the course of phase separation, indicating that polymer tended to aggregate in the region. To characterize the size of the polymer-rich domain, the following procedures were carried out. First, an X-Y plane was sectioned from the simulation box by specifying a Z value, as shown in Figure 1A. The polymer concentration (φp) in all the cells on the sectioned plane was analyzed, as depicted in Figure 1B. Then, by choosing a certain Y value, the polymer concentration distribution in the X direction was plotted, as shown in Figure 1C. The width of the peak area where φp g 0.3, illustrated by the distance between the arrows in Figure 1C, was defined as the size of the polymer-rich domain, and it possesses a distribution. By varying Y from 1 to 40, and repeating the above procedures, we obtained the domain size distribution in the X direction. Similarly, by choosing X from 1 to 40 and plotting the polymer concentration distribution in the Y direction, we obtained the domain size distribution in the Y direction. With the above analysis, the domain size distribution on the sectioned plane was determined. Four sectioned planes (Z ) 10, 20, 30, 40) were analyzed to provide the domain size distribution in the simulation box. The distribution probability with regard to a specific domain size (Ri) can then be calculated and a distribution probability function P(Ri) obtained, and ∫40 0 P(Ri) dRi ) 1. Subsequently, the average size of the polymerrich regions, 〈R〉, can be calculated with 〈R〉 ) ∫40 0 RiP(Ri) dRi. The 〈R〉 at different elapsed time steps were calculated to analyze the time dependency of 〈R〉 in each simulation case, which helps to investigate the mechanism of the coarsening of polymerrich regions. III. Results and Discussions 1. Verification of the Definition of the Polymer-Rich Domain. The polymer-rich domain in the present DPD study was defined as the region where the polymer concentration was equal to or higher than 0.3. The definition was verified by comparing the average polymer concentration of the polymerrich domains after 5 × 104 time steps with the equilibrium concentration of the polymer-rich phase calculated by using the Flory-Huggins theory.34 The equilibrium polymer concentration of the polymer-rich phase was 0.81 for the case of shorter polymer chains (Nb ) 10, χ ) 1.43) and 0.77 for that of longer polymer chains (Nb ) 100, χ ) 1.22). The polymer concentration in each of the 40 × 40 × 40 cells in the simulation box was calculated from the DPD simulation results, as discussed in the previous section. The polymer concentration distribution in the simulation box was thus determined and the distribution probability function P(φp) was obtained. The average polymer concentration of the polymerrich phase 〈φpr〉 was then calculated with the following equation:
∫0.31 φpP(φp) dφp 〈φpr〉 ) ∫0.31 P(φp) dφp
(11)
Chain Length and Stiffness Effects on Phase Seperation Dynamics
J. Phys. Chem. B, Vol. 112, No. 15, 2008 4607
Figure 1. Determination of the size of polymer-rich domains. (A) Sectioned X-Y planes from the DPD simulation box, where the blue regions represent the polymer-poor phase (φp < 0.3) and the other regions stand for the polymer-rich phase (φp g 0.3). (B) The distribution of polymer concentration on the sectioned X-Y plane. (C) Determination of the domain size by measuring the width of the peak area where φp g 0.3.
After 5 × 104 time steps, 〈φpr〉 ) 0.79 under χ ) 1.43 for the cases with Nb ) 10 and 〈φpr〉 ) 0.74 under χ ) 1.22 for the cases with Nb ) 100. The 〈φpr〉 values agreed well with the equilibrium concentration calculated by the Flory-Huggins theory, which provides the evidence supporting that the definition of polymer-rich domain is reasonable. In the simulation, the mapping of DPD to the Flory-Huggins theory was used to calculate the repulsion parameter from the interaction parameter as denoted in eq 3. The accuracy of the above mapping has been investigated for 2 < Nb < 10,25 but not for Nb ) 100. The results presented here could offer some information about the applicability of the mapping to Nb ) 100. The agreement of the 〈φpr〉 after 5 × 104 time steps with the equilibrium concentration calculated by the Flory-Huggins theory revealed that the aij obtained from the above mapping for Nb ) 100 did not introduce much error. 2. Evolution of Concentration Distribution during Phase Separation. Figure 2A describes the evolution of the concentration distribution in the simulation box during phase separation. The concentration distribution at the beginning of phase separation was Gaussian, and the average concentration was 0.3. With an increase in the elapsed time of phase separation, the distribution became bimodal, one peak corresponding to the polymer-rich domain and the other to the polymer-poor domain. As discussed in the previous section, after 5 × 104 time steps, the average concentration of the polymer-rich domain was close to the equilibrium concentration of the polymer-rich phase calculated by the Flory-Huggins theory. Similarly, the average
concentration of the polymer-poor domain also agreed well with the equilibrium concentration of the polymer-poor phase. Figure 2B shows that the concentration distribution was bimodal for cases (a)-(c) after 5 × 104 time steps. However, for case (d), the distribution was not bimodal, indicating that the solution had not yet demixed to two phase-separated domains. No bimodal concentration distribution was observed even with longer elapsed time steps. The absence of phase separation for case (d) was further evidenced by the observation that the morphology of the solution did not change with the elapsed time steps of simulation, which indicated that phase separation of polymer solution was hindered by longer and stiffer polymer chains (Nb ) 100, C ) 10). An increase in quench depth was employed to force the phase separation of the solution composed of the longer and stiffer polymer chains (Nb ) 100, C ) 10). With an larger interaction parameter of 1.33, phase separation of the solution occurred, evidenced by the bimodal concentration distribution shown in Figure 2C (case (e)). For comparison, simulation was performed for case (f) with Nb ) 100, C ) 2, and χ ) 1.33 to investigate the effect of chain stiffness. The concentration distributions for cases (e) and (f) were quite similar. The average concentration of the polymer-rich domain (φp g 0.3) was 0.74 for case (e), similar to that for case (f), 0.76. Both agree well with the equilibrium concentration of the polymer-rich phase (0.82) calculated by the Flory-Huggins theory. With an interaction parameter of 1.22, the effect of chain stiffness can be seen clearly in Figure 2B by comparing the
4608 J. Phys. Chem. B, Vol. 112, No. 15, 2008
Figure 2. (A) Time dependence of polymer concentration distribution probability functions for case (a). (B) Polymer concentration distribution probability functions at t* ) 5 × 104 for cases (a), (b), (c), and (d). (C) Polymer concentration distribution probability functions at t* ) 5 × 104 for cases (e) and (f).
concentration distributions for case (b) (Nb ) 100, C ) 2) and case (d) (Nb ) 100, C ) 10). The enhancement of chain stiffness could hinder the phase separation of the polymer solution. Nonetheless, with an interaction parameter of 1.33, little effect of the chain stiffness was observed. As shown in Figure 2C, the effect of chain stiffness and its relationship with quench depth will be further discussed in the following. 3. Effect of Chain Stiffness on Phase Separation of Polymer Solution. The results shown in Figure 2B revealed that, with an interaction parameter of 1.22, the repulsion between polymer and liquid was not large enough to bring about phase separation of the polymer solution consisting of longer and stiffer chains (Nb ) 100, C ) 10). Such a suppression on phase separation is believed to be related to the chain stiffness, since a system with the same chain length (Nb ) 100) but a lower chain stiffness (C ) 2) could undergo phase separation with the same interaction parameter. The suppression effect of chain stiffness could not be seen for the system of shorter chains, as the phase separation of the polymer solution consisting of shorter chains (Nb ) 10) can still be observed after an increase in the chain stiffness from C ) 2 to C ) 10. The above results indicated that, when the stiffer chains were long enough, the
Wu et al. elastic effect of polymer chains could stabilize the polymer solution, and a larger interaction parameter was thus needed to initiate the phase separation. The elastic effect on phase separation has been investigated by Cahn,35 showing that the elastic stress of a material could stabilize its solution and retard the phase separation. In addition, Flory and Rehner36-38 studied how the elastic free energy resulting from the cross-linked polymers influenced the critical interaction parameter and the equilibrium concentration of the polymer-rich phase. A semidilute polymer solution could present its elasticity more explicitly when there are sufficient entanglements in the polymer solution, and the elasticity would be enhanced by an increase in the chain stiffness.39,40 Therefore, a stronger repulsion is needed to overcome the elastic barrier imposed by the stiff chains and to initiate phase separation, which accounts for the difference between cases (b) and (d). However, the polymer entanglements are not permanent crosslinkings. Once the interaction parameter is large enough to result in phase separation, disentanglement occurs so as to alleviate the elastic barrier, which explains why the polymer concentration in the polymer-rich domains agreed well with the equilibrium calculated by the Flory-Huggins theory in case (e). The above discussion revealed that the stiffness of entangled polymer chain could impose an energy barrier to retard the phase separation of a polymer solution. In the following sections, how the chain stiffness influences the morphology evolution and the growth exponent during the course of phase separation will be discussed. 4. Morphology Evolution during Coarsening. Figure 3 presents the three-dimensional morphology evolutions for the simulation cases (a)-(c) and (e)-(f). The black regions stand for the polymer-rich domain, where the polymer concentration (φp) was larger than 0.3. The results of case (a) with Nb ) 10, C ) 2 are shown in Figure 3A. It can be seen that the phase separation resulted in an interconnected structure and the phaseseparated domains grew with the elapsed time steps of phase separation. 4.1. Chain Length Effect. The morphology evolution of case (b) (Nb ) 100, C ) 2) is shown in Figure 3B. The comparison of Figure 3A with Figure 3B reveals that an increase in the chain length slowed down the domain growth. At the same elapsed time step (t*), the polymer-rich domain for the case with longer chain was smaller and the degree of interconnectivity was higher. 4.2. Chain Stiffness Effect. The morphology evolution of case (c) (Nb ) 10, C ) 10) is presented in Figure 3C. The results demonstrate that an increase in chain stiffness remarkably retarded the aggregation of the polymer-rich domains during coarsening. Polymer-rich fractures, not seen in cases (a) and (b), were found in case (c), providing evidence to show that the domain growth for case (c) was via the mechanism of local aggregation. More dramatic effects of chain stiffness were observed for the cases with longer chains. With Nb ) 100, an increase in chain stiffness, from C ) 2 (case (c)) to C ) 10 (case (d)), completely suppressed the phase separation of the polymer solution, as discussed in the previous section. Although the effect of chain stiffness on phase separation was clearly demonstrated, its effect on phase separation dynamics could not be evaluated since no phase separation was observed for case (d). The quench depth was enlarged from 0.55 to 0.66 to force the phase separation of the polymer solution for case (e). The results are shown in Figure 3E. For comparison, the phase separation of the polymer solution with flexible longer chains (Nb ) 100 and
Chain Length and Stiffness Effects on Phase Seperation Dynamics
J. Phys. Chem. B, Vol. 112, No. 15, 2008 4609
Figure 3. Morphology evolution in simulation cases (a), (b), (c), (e), and (f) during the course of phase separation. The black regions represent the polymer-rich regions, and the white regions represent the polymer-poor regions.
C ) 2) under the quench depth of 0.66 was simulated, and the morphology evolution is shown in Figure 3F. Panels E and F of Figure 3 give more evidence showing that an increase in chain stiffness would slow down the domain growth. The above discussion was based on the morphology evolution shown in Figure 3, which gives only a qualitative description of the coarsening process during phase separation. To obtain clearer insight into the effect of chain length and stiffness on phase separation dynamics, quantitative analysis of domain growth was performed and the results are discussed in the following section. 5. Domain Growth during Coarsening. The growth exponent (n) of the domain growth curve, indicative of the coarsening mechanism, was determined as described below. The average size of the polymer-rich domains, 〈R〉, was plotted against the elapsed time step (t*) on a log-log scale as shown in Figure 4A (case (a)-(c)) and Figure 4B (case (e)-(f)). The straight lines on the log-log plots reveal that the domain growth could be described by the power law R ) Ktn, and the slopes of the lines gave the growth exponents (n). The results in Figure 4A,B show that two exponents are needed to characterize the domain growth for each case, indicating two mechanisms are involved in the coarsening process. It was reported that a thermally induced phase separation process usually begins with a diffusion-dominating
coarsening process, characterized by a growth exponent of 1/3, and then followed by the flow-dominating coarsening process driven by interfacial tension,41 characterized by a growth exponent of 2/3. For all the cases illustrated in Figure 4A,B, the growth exponents in the late stage of coarsening ranged from 0.65 to 0.73, close to 2/3, corresponding well to the flowdominating coarsening process. However, the growth exponents in the early stage were not 1/3 except for case (a), the case with shorter chain length (Nb ) 10) and lower chain stiffness (C ) 2). For the cases with longer chain length (Nb ) 100) or higher chain stiffness (C ) 10), the growth exponents in the early stage were smaller than 1/3, showing that the chain length and stiffness had influence on the diffusion coarsening process. 5.1. Effect of Chain Length. It is known that chain entanglement can affect the domain growth exponent when the size of phase-separated domains is smaller than the end-to-end distance of polymer chains.11 To check if the chain entanglement effect was important in the simulation cases, the end-to-end distance of each polymer chain in the simulation box was determined and the distribution probability was then calculated accordingly. The distribution probability of the size of polymerrich domain (R) and the end-to-end distance (R0) were presented in Figure 5. For case (a), at the beginning of coarsening, most of the domain sizes were larger than the end-to-end distance of
4610 J. Phys. Chem. B, Vol. 112, No. 15, 2008
Wu et al.
Figure 4. Growth curves of the average size of the polymer-rich domains:. (A) cases (a)-(c) with (χ - χb) ) 0.55; (B) cases (e) and (f) with (χ - χb) ) 0.66. The slopes of the line are the growth exponents (ni). The change of growth exponents for each case is indicated in the figures.
polymer chains, with an average R of 3.35 and an average R0 of 1.77, where the length was in units of cutoff radius rc. Since the domain size was not smaller than the end-to-end distance, chain entanglement effect was not important and the coarsening followed the L-S1 mechanism, with a growth exponent close to 1/ in the early stage. For the case with larger chain length ((b), 3 (e), and (f)), the end-to-end distance of polymer chains was larger than the domain size. Therefore, the polymer chains entanglement played a significant role in the coarsening dynamics. Due to the chain entanglement effect, the initial growth exponents for cases (b), (e), and (f) were close to 1/4, smaller than 1/3. The exponent of 1/4 agreed well with the simulation results of the phase separation of the solution with long polymer chains obtained by using molecular dynamics (MD) simulation17 and Monte Carlo (MC) simulation.18 Nonetheless, the growth exponent did not match the result obtained by de Gennes11 using the tube model, 1/5. It was noteworthy that in de Gennes’s analysis, the domain size was assumed much smaller than the end-to-end distance (R , R0), though the result was used up to R ) R0. Although the domain size (R) at elapsed time step t* ) 500 in our simulation was smaller than the endto-end distance (R0), R and R0 were of the same order of magnitude. Since the simulation condition (R < R0) was between that for de Gennes’s theory (R , R0) and that for the L-S mechanism (R > R0), the degree of local adjustment should also be in between, which explains why the exponent we obtained (1/4) was between 1/5 and 1/3. Although the chain entanglement affected the growth exponent in the early stage, it had little effect on the growth exponent in the later flow-dominating coarsening stage, as shown in Figure 4A,B. Such an observation can be reasoned by the fact that the flow-dominating stage occurred after the diffusiondominating stage, during which the domain size had grown large enough to neglect the entanglement effect.
Figure 5. Distributions of the domain size (R) and the end-to-end distance (R0) for the simulation cases (a)-(c) and (e)-(f) at elapsed time step t* ) 500, where P(di) is the probability for R or R0 with a value of di.
5.2. Effect of Chain Stiffness. Acomparison of the domain growth curves of cases (a) and (c) shows that an increase in chain stiffness remarkably slowed down the growth of polymerrich domains and the growth exponent in the early stage of coarsening was reduced approximately from 1/3 to 1/6 in the diffusion regime. Nevertheless, in the following stage, the growth exponent was still 2/3. The results indicated that the chain stiffness significantly influenced the coarsening kinetics in the diffusion regime but not in the flow regime. In accordance to the distributions of the domain size (R) and the end-to-end distance (R0) in units of rc at t* ) 500 shown in Figure 5C, the average domain size (〈R〉 ) 3.27) was larger than the average end-to-end distance of polymer chains (〈R〉 ) 2.94), suggesting that the chains were too short to introduce significant entanglement effect. The results showed that, even when the chain entanglement effect was not significant, the growth exponent could still become smaller due to an increase in chain stiffness. It has been shown that29-31 the stiffer chains in the semidilute solution were trapped in the cage formed by the neighboring stiffer rods and the motion of chains transverse to the rod was strictly hindered. The translational motion of the stiffer chains was relatively significant, and the stiffer chains tend to perform rodlike motion. Therefore, the severe confinement caused by the chain stiffness was analogous to the constraint by the entanglements. Consequently, the coarsening of stiffer polymers was slowed down due to the confined motion of polymer chains. The growth exponent of 1/6 agreed well with the results of the theoretical work of Muratov13 and the experiments performed by Tanaka15 who proposed that the transient gel occurred during the course of phase separation kinetics. Muratov analyzed the effect of chain entanglement on domain growth under the
Chain Length and Stiffness Effects on Phase Seperation Dynamics condition that the domain size was smaller than the distance between entangled points. It is noteworthy that in de Gennes’s analysis,11 the domain size was still much larger than the distance between entangled points, though it was much smaller than the end-to-end distance. de Gennes postulated that the chain segments can still move randomly between the entangled points in the tube, but in Muratov’s analysis, the chain can only move in a wormlike pattern between the entangled points. The good agreement between the results of case (c) with Muratov’s work indicated that the chain stiffness in case (c) resulted in such a severe limitation that the polymer chain can only move in reptation. Figure 4B shows the chain stiffness effects on domain growth for the solution with longer polymer chains (Nb ) 100) under the quench depth of (χ - χb) ) 0.66. In both cases (e) and (f), the growth exponents in the diffusion regime were close to 1/4. By taking cases (b) and (d) into account, with the chain length Nb ) 100 and the quench depth (χ - χb) ) 0.55, the absence of the phase separation in case (d) and the domain growth in case (b) demonstrated the remarkable difference due to increase in the chain stiffness. While the dramatic effect of chain stiffness was not observed in case (e), implying that the increase of quench depth could alleviate the confinement originated from the chain stiffness, a less dramatic influence was shown. More research work is needed to elucidate the interplay between the quench depth and the constraints due to chain stiffness. Nonetheless, an effect of the chain stiffness can be seen from Figure 4B; the crossover in the coarsening mechanism from diffusion to flow was delayed by the chain stiffness. The result indicated that the polymer network structure formed by the stiffer, longer chains provided a barrier for the onset of interfacial-tension driven flow. Consequently, the increase in chain stiffness restricts the chain to move in reptation and intensifies the polymer structure during coarsening so that the structure is stabilized, and it can be stressed that the chain stiffness remarkably impacts the coarsening in diffusion regime as well as the chain length. Moreover, the effects of chain stiffness and chain length on growth exponents can only be observed in the diffusion regime, indicating that the interfacialtension driven flow in the late stage could overwhelm the topological confinements caused by the chain interactions. The simulation results indicated that the coarsening dynamics was influenced by the chain entanglement and an increase in chain stiffness could further strengthen such effect. Theoretically, the chain entanglement effect cannot be properly simulated by the bead-spring models due to the problem of artificial bond crossing42 caused by the lack of spring-spring repulsions.43 Two reasons are given below to account for why the chain entanglement effects were clearly seen in our simulations. First, in all the simulation cases, the distance between the adjacent beads in a polymer chain was shorter than the cutoff radius within the interaction range; therefore, the interaction forces between beads could help to avoid the artificial bond crossing. Second, the increase in chain stiffness might also reduce the occurrence of the bond crossing. There is a possibility that in our simulations the artificial bond crossing was not completely suppressed and the chain entanglement effect was thus underestimated. Investigations including the spring-spring repulsions to completely avoid the artificial bond crossing still remain a challenge. IV. Conclusion DPD simulations were carried out to investigate the effects of polymer chain length and stiffness on the phase separation
J. Phys. Chem. B, Vol. 112, No. 15, 2008 4611 dynamics of semidilute polymer solutions. For the case of shorter and more flexible chains, a crossover was observed in the growth exponent from 1/3 to 2/3, reflecting that the coarsening process was first dominated by diffusion and then by interfacialtension driven flow. With longer chains, the morphological evolution was slowed down and resulted in a growth exponent of 1/4 in the diffusion-dominating regime, regardless of chain stiffness. The slow down was resulted from the entanglement effect of polymer chains, and the growth exponent of 1/4 obtained agreed well with that determined by the MD and MC simulations. As for the case with shorter but stiffer chains, the polymer chains diffused in a rodlike pattern due to the chain stiffness. As a result, a growth exponent of 1/6 in the diffusion regime was obtained, corresponding to those obtained from the theoretical and experimental analyses in the coarsening dynamics of transient gels. When both the polymer chain length and stiffness were increased, the stiffer polymer network stabilized the solution and it required larger repulsion between polymer and liquid to break up the network and to initiate phase separation. Moreover, the crossover of the domain growth from diffusion to flow was delayed by the network consisting of stiffer and longer chains, indicating that the enhanced chain length and stiffness provided a barrier for the onset of interfacial-tension driven flow. It was also demonstrated that the chain length and stiffness had influences on the growth exponents only in the diffusion regime but not in the flow regime, revealing that the minimization of the interfacial energy overwhelmed the chain interactions in the later stage of coarsening. More research work is needed to avoid artificial bond crossing and the underestimation in entanglement effect by modifying the bead-bead repulsion. Acknowledgment. We gratefully thank the National Center for High-Performance Computing of Taiwan for providing computing time. References and Notes (1) Lifshitz, I. M.; Slyozov, V. V. J. Phys. Chem. Solids 1961, 19, 35. (2) Binder, K.; Stauffer, D. Phys. ReV. Lett. 1974, 33, 1006. (3) Siggia, E. D. Phys. ReV. A 1979, 20, 595. (4) Furukawa, H. Phys. ReV. A 1985, 31, 1103. (5) Furukawa, H. AdV. Phys. 1985, 34, 703. (6) Bray, A. J. AdV. Phys. 1994, 43, 357. (7) Gunton, J. D.; San Miguel, M.; Sahni, P. S. Phase Transitions and Critical Phenomena; Academic Press: New York, 1983; Vol. 8. (8) Mulder, M. Basic Principles of Membrane Technology, 2nd ed.; Kluwer Academic: Dordrecht, Netherlands, 1996. (9) Kesting, R. E. Synthetic Polymeric Membranes : A Structural PerspectiVe, 2nd ed.; Wiley: New York 1985. (10) Voorhees, P. W. J. Stat. Phys. 1985, 38, 231. (11) de Gennes, P. G. J. Chem. Phys. 1980, 72, 4756. (12) Pincus, P. J. Chem. Phys. 1981, 75. (13) Muratov, C. B. Phys. ReV. Lett. 1998, 81, 3699. (14) Doi, M.; Edwards, S. F. The Theory of Polymer Dynamics; Clarendon Press: Oxford, 1986. (15) Tanaka, H. Phys. ReV. Lett. 1993, 71, 3158. (16) Araki, T.; Tanaka, H. Macromolecules 2001, 34, 1953. (17) Liu, H.; Bhattacharya, A.; Chakrabarti, A. J. Chem. Phys. 1999, 111, 11183. (18) Termonia, Y. Macromolecules 1997, 30, 5367. (19) Cao, Y.; Zhang, H. D.; Xiong, Z.; Yang, Y. L. Macromol. Theory Simul. 2001, 10, 314. (20) Luo, K. F.; Gronski, W.; Friedrich, C. Macromol. Theory Simul. 2004, 13, 365. (21) Tanaka, H. J. Chem. Phys. 1994, 100, 5323. (22) Tanaka, H. J. Phys.: Condens. Matter 2000, 12, R207. (23) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Europhys. Lett. 1992, 19, 155. (24) Koelman, J. M. V. A.; Hoogerbrugge, P. J. Europhys. Lett. 1993, 21, 363. (25) Groot, R. D.; Warren, P. B. J. Chem. Phys. 1997, 107, 4423. (26) Groot, R. D.; Madden, T. J. J. Chem. Phys. 1998, 108, 8713.
4612 J. Phys. Chem. B, Vol. 112, No. 15, 2008 (27) Novik, K. E.; Coveney, P. V. Phys. ReV. E 2000, 61, 435. (28) Kholodenko, A. L. J. Chem. Soc., Faraday Trans. 1995, 91, 2473. (29) Odijk, T. Macromolecules 1983, 16, 1340. (30) Keep, G. T.; Pecora, R. Macromolecules 1988, 21, 817. (31) Tracy, M. A.; Pecora, R. Annu. ReV. Phys. Chem. 1992, 43, 525. (32) Espanol, P.; Warren, P. Europhys. Lett. 1995, 30, 191. (33) de Gennes, P. G. Scaling Concepts in Polymer Physics; Cornell University Press: Ithaca, N.Y., 1979. (34) Flory, P. J. Principles of Polymer Chemistry; Cornell University Press: Ithaca, N.Y., 1953.
Wu et al. (35) Cahn, J. W. Acta Metall. Mater. 1961, 9, 795. (36) Flory, P. J.; Rehner, J. J. Chem. Phys. 1943, 11, 521. (37) Flory, P. J.; Rehner, J. J. Chem. Phys. 1943, 11, 512. (38) Flory, P. J. J. Chem. Phys. 1950, 18, 108. (39) Tschoegl, N. W.; Ferry, J. D. J. Am. Chem. Soc. 1964, 86, 1474. (40) Ferry, J. D. Viscoelastic Properties of Polymers; John Wiley: New York, 1985. (41) Haas, C. K.; Torkelson, J. M. Phys. ReV. E 1997, 55, 3191. (42) Pan, G.; Manke, C. W. Int. J. Mod. Phys. B 2003, 17, 231. (43) Kumar, S.; Larson, R. G. J. Chem. Phys. 2001, 114, 6937.
