196 O.H.Zeng et al.Prog.Polym.Sci.33 (2008)191-269 There are many algorithms ation of mo don using f difference interface and also toth an commonly used in MDsi tions [17.All algorithms assume that the atomic 2.1.2.Monte Carlo position,velocities and accelerations a can be MC te also called me olis r approximated by a Taylor series expansion: don +80)=0+0t+5a0r2+…, numbers to generate a sample population of the (3) system from which one can calculate the properties of interest.A MC simulation usually consists of iu+0)=i(0+a0t+5b0)6r+, thre (4) e typical steps.In the first step.the physical problem unde mo In th at+80)=a)+bu)8t+. emodel is so Generally speaking.a good integration algorithm ical should conserve the total energy and momentum ined dat In the d h and be time- statistical methods.MCp ovides only the informa eversible.It should lso be easy to tion on equilibrium pronerties (e g free enerov step.The vene phase equilibrium).different from Md which gives mit a bly出 gorithn and s the positions-8t)from the pre ious sten t-t to calculate the new positions r(t+8r)at t+8t. at lly mo ing one ator According to Taylor series expansion of Eq.(3). we have 7u+60=70+t081+5a0ǒr2+…, △H=HU-H), (11) where H()and the Hamiltonian associated original an gurat 7u-60)=70-(01+5a062- (7) s.If△H Adding Eqs.(6)and(7)together gives cording to the ould g rule t Tu+80)=270-7t-6)+au)ǒ2 of lower energy Hence the move ent is immedi (8) ately accented and the displaced atom remains in its The velocities at time t and t+t can be new position.If AH0,the move is accepted only respectively estimated with a certain probability pwhich is given by 元(0=[T(u+8)-7t-6/28, (9) ep) (+58)=+0-0/8 (10) where is the Bot al between 0 and I and determine the new configuration according to the following rule: (uVT),microcanonical (NVE),canonical (NVT) and isothermal-isobaric (NPT).The constant tem m(←) the move is accepted perature and pressure can be controlled by adding over an and b the move is not acepied. respe tively. MD int allows us to in configuration is rejected.one the o and peats the polymer structure and dynamics in the vicinity of
There are many algorithms for integrating the equation of motion using finite difference methods. The algorithms of verlet, velocity verlet, leap-frog and Beeman, are commonly used in MD simulations [17]. All algorithms assume that the atomic position ~r, velocities v * and accelerations a * can be approximated by a Taylor series expansion: r *ðt þ dtÞ ¼ r *ðtÞ þ v *ðtÞdt þ 1 2 a *ðtÞdt 2 þ , (3) v *ðt þ dtÞ ¼ v *ðtÞ þ a *ðtÞdt þ 1 2 bðtÞdt 2 þ , (4) a *ðt þ dtÞ ¼ a *ðtÞ þ bðtÞdt þ . (5) Generally speaking, a good integration algorithm should conserve the total energy and momentum and be time-reversible. It should also be easy to implement and computationally efficient, and permit a relatively long time-step. The verlet algorithm is probably the most widely used method. It uses the positions ~rðtÞ and accelerations a *ðtÞ at time t, and the positions~rðt dtÞ from the previous step tdt to calculate the new positions ~rðt þ dtÞ at t+dt. According to Taylor series expansion of Eq. (3), we have r *ðt þ dtÞ ¼ r *ðtÞ þ v *ðtÞdt þ 1 2 a *ðtÞdt 2 þ , (6) r *ðt dtÞ ¼ r *ðtÞ v *ðtÞdt þ 1 2 a *ðtÞdt 2 . (7) Adding Eqs. (6) and (7) together gives r *ðt þ dtÞ ¼ 2 r *ðtÞ r *ðt dtÞ þ a *ðtÞdt 2 . (8) The velocities at time t and t þ 1 2dt can be respectively estimated v *ðtÞ¼½r *ðt þ dtÞ r *ðt dtÞ=2dt, (9) v * t þ 1 2 dt ¼ ½r *ðt þ dtÞ r *ðtÞ=dt. (10) MD simulations can be performed in many different ensembles, such as grand canonical (mVT), microcanonical (NVE), canonical (NVT) and isothermal–isobaric (NPT). The constant temperature and pressure can be controlled by adding an appropriate thermostat (e.g., Berendsen, Nose, Nose–Hoover and Nose–Poincare) and barostat (e.g., Andersen, Hoover and Berendsen), respectively. Applying MD into polymer composites allows us to investigate into the effects of fillers on polymer structure and dynamics in the vicinity of polymer–filler interface and also to probe the effects of polymer–filler interactions on the materials properties. This will be discussed in Section 3.3. 2.1.2. Monte Carlo MC technique, also called Metropolis method [19], is a stochastic method that uses random numbers to generate a sample population of the system from which one can calculate the properties of interest. A MC simulation usually consists of three typical steps. In the first step, the physical problem under investigation is translated into an analogous probabilistic or statistical model. In the second step, the probabilistic model is solved by a numerical stochastic sampling experiment. In the third step, the obtained data are analyzed by using statistical methods. MC provides only the information on equilibrium properties (e.g., free energy, phase equilibrium), different from MD which gives non-equilibrium as well as equilibrium properties. In a NVT ensemble with N atoms, one hypothesizes a new configuration by arbitrarily or systematically moving one atom from position i-j. Due to such atomic movement, one can compute the change in the system Hamiltonian DH: DH ¼ HðjÞ HðiÞ, (11) where H(i) and H(j) are the Hamiltonian associated with the original and new configuration, respectively. This new configuration is then evaluated according to the following rules. If DHo0, then the atomic movement would bring the system to a state of lower energy. Hence, the movement is immediately accepted and the displaced atom remains in its new position. If DHX0, the move is accepted only with a certain probability pi-j which is given by pi!j / exp DH kBT , (12) where kB is the Boltzmann constant. According to Metropolis et al. [19] one can generate a random number x between 0 and 1 and determine the new configuration according to the following rule: xp exp DH kBT ; the move is accepted; x4exp DH kBT ; the move is not accepted: If the new configuration is rejected, one counts the original position as a new one and repeats the process by using other arbitrarily chosen atoms. ARTICLE IN PRESS 196 Q.H. Zeng et al. / Prog. Polym. Sci. 33 (2008) 191–269
Q.H.Zeng et al Prog.Polym.SeL 33(2008)191-269 In a uVT ensemble,one hypothesizes a new and time scales currently inaccessible by the classical configuration/by arbitrarily choosing one atom MD methods and proposing that it can be exchanged by an atom of a different kind. This procedure affects the 2.2.1.Brownian dynamics chemical composition of the system.Also,the move BD simulation is similar to MD simulations [20] is accepte with a cer pro ity.How However,it introduces a few new approximations the hange e energy c that allow one to pertorm simulations on the to the follo AU2o the move of sitional own up to a e e MD accepted.However.if AU>0.the move is acce pted d with an de with a certain probability which is given by △U are typically ignored,allowing a much larger time- Pi exp(- (13 step than that of MD.Therefore,BD is particularly useful for systems where there is a large gap of time where AU is the change in the sum of the mixing energy and the chemical potential of the mixture.If example,in polym the new configuration is rejected one counts the st m original configuration as a new of th ndes of the using some othe step.However.if the detailed motion of the solvent mer na ompos olecrhots M molecules is concerned they may he removed from the simulation and their effects on the polymer are rious factors represented by dissipative(-yp)and random ((t)) orce terms. the forces in the governing Eq.(1)is replaced by a Langevin equation 2.2.Microscale methods F0=F听-m+aa0 (14) The modeling and simulation at microscale aim to bridge molecular methods and continuum where F is the conservative force of particle j methods and avoid their shortcomings.Specifically acting on particle i,y and o are constants depending in nanoparticle-polymer systems. the study of on the system,P;the momentum of particle 1,and structural evolution (i.e.,dynamics of phase separa- (1)a Gaus con the the ine on st degree th interact that the macrosconic hehavior of the system will no be hydrodynamic.In addition.the effect of one solute molecule on another through the flow of methods.In contrast.the interactions betweer solvent molecules is neglected.Thus,BD can only components can be examined at an atomistic level reproduce the diffusion properties but not the but are usually not straightforward to incorporate hydrodynamic fle the simulation at the continuum level.Therefore,various s does not obey the Navier okes equations aluate and structure and particle by H BD.DPD.LB.time-depe I anda and Koel 11t n si hath and non-Newtonian fluids.including polymer melts methods.a polymer system is usually treated with and blends,on microscopic length and time scales. a field description or microscopic particles that Like MD and BD,DPD is a particle-based method. incorporate molecular details implicitly.Therefore However.its basic unit is not a single atom or they are able to simulate the phenomena on length molecule but a molecular assembly (i.e..a particle)
In a mVT ensemble, one hypothesizes a new configuration j by arbitrarily choosing one atom and proposing that it can be exchanged by an atom of a different kind. This procedure affects the chemical composition of the system. Also, the move is accepted with a certain probability. However, one computes the energy change DU associated with the change in composition. The new configuration is examined according to the following rules. If DUo0, the move of compositional change is accepted. However, if DUX0, the move is accepted with a certain probability which is given by pi!j / exp DU kBT , (13) where DU is the change in the sum of the mixing energy and the chemical potential of the mixture. If the new configuration is rejected one counts the original configuration as a new one and repeats the process by using some other arbitrarily or systematically chosen atoms. In polymer nanocomposites, MC methods have been used to investigate the molecular structure at nanoparticle surface and evaluate the effects of various factors. 2.2. Microscale methods The modeling and simulation at microscale aim to bridge molecular methods and continuum methods and avoid their shortcomings. Specifically, in nanoparticle–polymer systems, the study of structural evolution (i.e., dynamics of phase separation) involves the description of bulk flow (i.e., hydrodynamic behavior) and the interactions between nanoparticle and polymer components. Note that hydrodynamic behavior is relatively straightforward to handle by continuum methods but is very difficult and expensive to treat by atomistic methods. In contrast, the interactions between components can be examined at an atomistic level but are usually not straightforward to incorporate at the continuum level. Therefore, various simulation methods have been evaluated and extended to study the microscopic structure and phase separation of these polymer nanocomposites, including BD, DPD, LB, time-dependent Ginsburg–Landau (TDGL) theory, and dynamic DFT. In these methods, a polymer system is usually treated with a field description or microscopic particles that incorporate molecular details implicitly. Therefore, they are able to simulate the phenomena on length and time scales currently inaccessible by the classical MD methods. 2.2.1. Brownian dynamics BD simulation is similar to MD simulations [20]. However, it introduces a few new approximations that allow one to perform simulations on the microsecond timescale whereas MD simulation is known up to a few nanoseconds. In BD the explicit description of solvent molecules used in MD is replaced with an implicit continuum solvent description. Besides, the internal motions of molecules are typically ignored, allowing a much larger timestep than that of MD. Therefore, BD is particularly useful for systems where there is a large gap of time scale governing the motion of different components. For example, in polymer–solvent mixture, a short time-step is required to resolve the fast motion of the solvent molecules, whereas the evolution of the slower modes of the system requires a larger timestep. However, if the detailed motion of the solvent molecules is concerned, they may be removed from the simulation and their effects on the polymer are represented by dissipative (gp) and random (sz(t)) force terms. Thus, the forces in the governing Eq. (1) is replaced by a Langevin equation, FiðtÞ ¼ X jai F C ij gpi þ sziðtÞ, (14) where F C ij is the conservative force of particle j acting on particle i, g and s are constants depending on the system, pi the momentum of particle i, and z(t) a Gaussian random noise term. One consequence of this approximation of the fast degrees of freedom by fluctuating forces is that the energy and momentum are no longer conserved, which implies that the macroscopic behavior of the system will not be hydrodynamic. In addition, the effect of one solute molecule on another through the flow of solvent molecules is neglected. Thus, BD can only reproduce the diffusion properties but not the hydrodynamic flow properties since the simulation does not obey the Navier–Stokes equations. 2.2.2. Dissipative particle dynamics DPD was originally developed by Hoogerbrugge and Koelman [21]. It can simulate both Newtonian and non-Newtonian fluids, including polymer melts and blends, on microscopic length and time scales. Like MD and BD, DPD is a particle-based method. However, its basic unit is not a single atom or molecule but a molecular assembly (i.e., a particle). ARTICLE IN PRESS Q.H. Zeng et al. / Prog. Polym. Sci. 33 (2008) 191–269 197
198 O.H.Zeng et al.Prog.Polym.Sci.33 (2008)191-269 ir mas of solid particles.The LB method is originated from which ,dissipative F号 and randon Fa=FC+F+FR (15) lattice with particles residing on the nodes.A set of Boolean variable ncx.t)(i=1.. .M)describing FG=Ⅱoocv)e, (16 the particle occupation is defined,where M is the number of directions of the particle FB=-yonU)(ei·Plei (17) n nod automaton is as follow FR=(rg) (18) n(x+e,1+1)=n(,)+2,》 is 1=0.1..... (20) coefficient.a noise am nplitude anda randon where e,are the local particle velocity.During the noise term with zero mean fie -0)and evolution.each time-step can be divided into two -steps particl unit variance.oco and o are the weight functions for each interaction force [221.While the mo t wit h two particles,1 and then change thei r velocity directions according to seattering rules The main feature of the LB method is to replace to decr linearly the particle occupation varables n,(Boolean var with ble unctions (rea increasing r Revond a certain cut-off se =( a r.the weight functions and thus the forces are all dual partick zero. a and -pa where Therefore,the total force Ft)acting on particle i at time is ation fro either the disc del o (19) the boltzmann kinetic equation and to derive th macroscopic Navier-Stokes equations from the LB Because the forces are pairwise and momentum is equation.For example,based on the discrete kinetic the LB equation for a multicomponen conserved, the macroscopic behavior directly in- corporates Navier-Stokes flow can be expres ed as hydrodynamics. ever, 6 fc+e△x,t+△)-fx,)=2Wfx). h dis ve and rce ten (i=0,1...,0=1.2,..,S (21) the effects of Brow ger lenat where the left-hand side of the equation describes scales.DPD has MD o changes in the distribution function due to free particle motion and the right-hand side contains the with far fewer particles than required in a MD collision operation describing change simulation because of its larger particle size.Be- ays s sides,its force forms allow larger time steps to be Bhatnagar- rook Iorm taken than those in MD 2x)=Ux,0-g9x,), B25 (22 cale method that (i=0,1,,MG=1,2,, suited for whe Ar is the dimensionless collision s It has ently relaxation time of ent of the the phase separation of binary fluids in the presence fluid.is the lattice velocity vector.the subscript
DPD particles are defined by their mass Mi, position ri and momentum pi. The interaction force between two DPD particles i and j can be described by a sum of conservative F C ij , dissipative FD ij and random forces F R ij [22–24]: Fij ¼ F C ij þ FD ij þ F R ij , (15) F C ij ¼ P0oCðrijÞe^ij, (16) FD ij ¼ goDðrijÞðe^ij pijÞe^ij, (17) F R ij ¼ szijoRðrijÞe^ij, (18) where rij ¼ jri rjj; e^ij ¼ r^ij=rij; P0 is a constant related to the fluid compressibility, g is a friction coefficient, s a noise amplitude and zij a random noise term with zero mean (i.e., ozij4 ¼ 0) and unit variance. oC, oD, and oR are the weight functions for each interaction force [22]. While the interaction potentials in MD are high-order polynomials of the distance rij between two particles, in DPD the potentials are softened so as to approximate the effective potential at microscopic length scales. The form of the conservative force in particular is chosen to decrease linearly with increasing rij. Beyond a certain cut-off separation rc, the weight functions and thus the forces are all zero. Therefore, the total force Fi(t) acting on particle i at time t is FiðtÞ ¼ X jai F C ij þX jai FD ij þX jai F R ij . (19) Because the forces are pairwise and momentum is conserved, the macroscopic behavior directly incorporates Navier–Stokes hydrodynamics. However, energy is not conserved because of the presence of the dissipative and random force terms which are similar to those of BD, but incorporate the effects of Brownian motion on larger length scales. DPD has several advantages over MD, for example, the hydrodynamic behavior is observed with far fewer particles than required in a MD simulation because of its larger particle size. Besides, its force forms allow larger time steps to be taken than those in MD. 2.2.3. Lattice Boltzmann LB [25] is another microscale method that is suited for the efficient treatment of polymer solution dynamics. It has recently been used to investigate the phase separation of binary fluids in the presence of solid particles. The LB method is originated from lattice gas automaton which is constructed as a simplified, fictitious molecular dynamic in which space, time and particle velocities are all discrete. A typical lattice gas automaton consists of a regular lattice with particles residing on the nodes. A set of Boolean variable ni(x, t) (i ¼ 1,y, M) describing the particle occupation is defined, where M is the number of directions of the particle velocities at each node. The evolution equation of the lattice gas automaton is as follows: niðx þ ei;t þ 1Þ ¼ niðx;tÞ þ Oiðnðx;tÞÞ, ði ¼ 0; 1; ... ; MÞ, ð20Þ where ei are the local particle velocity. During the evolution, each time-step can be divided into two sub-steps: (i) streaming, in which each particle moves to the nearest node in the direction of its velocity and (ii) collision, which occurs when particles arriving at a node interact with each other and then change their velocity directions according to scattering rules. The main feature of the LB method is to replace the particle occupation variables ni (Boolean variables), by single-particle distribution functions (real variables) fi ¼ /niS and neglect individual particle motion and particle–particle correlations in the kinetic equation, where / S indicates an ensemble average. There are several ways to obtain the LB equation from either the discrete velocity model or the Boltzmann kinetic equation, and to derive the macroscopic Navier–Stokes equations from the LB equation. For example, based on the discrete kinetic equation, the LB equation for a multicomponent flow can be expressed as f s i ðx þ eiDx;t þ DtÞ f s i ðx;tÞ ¼ Oiðf s i ðx;tÞÞ, ði ¼ 0; 1; ... ; M; s ¼ 1; 2; ... ; SÞ, ð21Þ where the left-hand side of the equation describes changes in the distribution function due to free particle motion and the right-hand side contains the collision operation O, describing changes due to pairwise collisions which are always simplified to the linear Bhatnagar–Gross–Krook form: Oiðf s i ðx;tÞÞ ¼ 1 ts ðf s i ðx;tÞ f s;eq i ðx;tÞÞ, ði ¼ 0; 1; ... ; M; s ¼ 1; 2; ... ; SÞ, ð22Þ where ts ¼ ls =Dt is the dimensionless collision– relaxation time constant of component s of the fluid, ei is the lattice velocity vector, the subscript i ARTICLE IN PRESS 198 Q.H. Zeng et al. / Prog. Polym. Sci. 33 (2008) 191–269
Q.H.Zeng et al Prog.Polym.SeL 33(2008)191-269 199 represents the lattice velocity direction.fct)is the distribution function of the particles of component 器a8 mh with velocity e along direction iat latticex and is minimized to simulate a temperature quench from time t,and;(x,/is the corresponding equili- the miscible region of the phase diagram to the brium distribution function. odynamic behavior immiscible region.Thus,the resulting time-dependent In the LB methoc the blend md equatio a d (23) this method to where UGL()is the Ginzburg-Landu free energy, For an incompressible binary polymer blend,the equation of motion can be written as follows: and U)is a potential that describes the coupling interaction between the ith particle and the surrounding fluid, D=Mi】+. (26 oo(r.t) UGL()=dr( (24) Here()sa thermal represen For a binary k>0 are phe ara ()describes a single homo neous phase while al=「drU】+F 27 k。T k。T <0,it yields two-or multiphases coexistence. where UFH()is the Flory-Huggins free energy of U.()=d山rVoe-r-s(p)-s》}, (25) mixing.given by UFH( KRI 2=忌n+“-+x1-9 the of th %>0 (28) ents its coup ginlCnaci8nand Here x is the enthalpic interaction parameter An important advantage of the LB method is that are the the two polymer componen V4 and N microscopic physical interactions of the fluid numbe n polymer Aa particles can be conveniently incorporated into the ely Th 8 with /2 numerical model. Compared the er Stokes LB me handle the system is ce the critical co s am centration whereas if and re it is immiscible. To simplify TDGL,Oono and co-workers pro- and posed cell dynamic method (CDM)which wa level.However.its main disadvanta is that it i the TDG typically not guaranteed to be numerically stable and aplacian term is rep by it may lead to physically unreasonable results,for nd.the free instance,in the case of high forcing rate or high iven by interparticle interaction strength UCDM[(r)] 2.2.4.Time-dependent Ginzhu TDGL is a microscale method for simulating the -4 In(cosh structural evolution of phase-separation in polymer 29) blends and block copolymers.It is based on the Cahn-Hilliard-Cook (CHC) nonlinear diffusion This model reproduces the growth kinetics of the equation for a binary blend and falls under the more TDGL model,demonstrating that such quantities
represents the lattice velocity direction, f s i ðx; tÞ is the distribution function of the particles of component s with velocity ei along direction i at lattice x and time t, and f s;eq i ðx; tÞ is the corresponding equilibrium distribution function. In the LB method, the thermodynamic behavior can be described by the following free-energy functional: UðfÞ ¼ UGLðfÞ þX N i¼1 UiðfÞ, (23) where UGL(f) is the Ginzburg–Landu free energy, and Ui(f) is a potential that describes the coupling interaction between the ith particle and the surrounding fluid, UGLðfÞ ¼ Z dr t 2 f2 þ u 4 f4 þ k 2 jrfj 2 , (24) where the term k 2jrfj 2 represents the free-energy cost of forming fluid–fluid interface, and u40 and k40 are phenomenological parameters. If t40, UGL(f) describes a single homogeneous phase, while to0, it yields two- or multiphases coexistence. UiðfÞ ¼ Z dr V0eðjrsij=r0Þ ðfðrÞ fðsiÞÞ2 , (25) where si is the position of ith particle and N the total number of particles. The constant V040 characterizes the strength of the coupling interaction and r0 represents its range. An important advantage of the LB method is that microscopic physical interactions of the fluid particles can be conveniently incorporated into the numerical model. Compared with the Navier– Stokes equations, the LB method can handle the interactions among fluid particles and reproduce the microscale mechanism of hydrodynamic behavior. Therefore it belongs to the MD in nature and bridges the gap between the molecular level and macroscopic level. However, its main disadvantage is that it is typically not guaranteed to be numerically stable and may lead to physically unreasonable results, for instance, in the case of high forcing rate or high interparticle interaction strength. 2.2.4. Time-dependent Ginzburg– Landau method TDGL is a microscale method for simulating the structural evolution of phase-separation in polymer blends and block copolymers. It is based on the Cahn–Hilliard–Cook (CHC) nonlinear diffusion equation for a binary blend and falls under the more general phase-field and reaction-diffusion models [26–28]. In the TDGL method, a free-energy function is minimized to simulate a temperature quench from the miscible region of the phase diagram to the immiscible region. Thus, the resulting time-dependent structural evolution of the polymer blend can be investigated by solving the TDGL/CHC equation for the time dependence of the local blend concentration. Glotzer and co-workers have discussed and applied this method to polymer blends and particle-filled polymer systems [29]. For an incompressible binary polymer blend, the equation of motion can be written as follows: qfðr;tÞ qt ¼r Mr dU½fðr;tÞ dfðr;tÞ þ zðr;tÞ. (26) Here z(r, t) is a thermal noise term with zero mean, U is the free energy, and the mobility M may depend on the order parameter f. For a binary polymer blend, one typically takes U½fðrÞ kBT ¼ Z dr UFH ½fðrÞ kBT þ kðfÞjrfðrÞj2 , (27) where UFH(f) is the Flory–Huggins free energy of mixing, given by UFH ðfÞ kBT ¼ f NA ln f þ ð1 fÞ NB lnð1 fÞ þ wfð1 fÞ. (28) Here w is the enthalpic interaction parameter between the two polymer components. NA and NB are the number of segments in polymer A and polymer B, respectively. The Flory–Huggins model has a critical point at fc ¼ N1=2 B =ðN1=2 A þ N1=2 B Þ and wc ¼ ðN1=2 A þ N1=2 B Þ 2 =2NANB. If wowc the system is miscible at the critical concentration whereas if w4wc it is immiscible. To simplify TDGL, Oono and co-workers proposed cell dynamic method (CDM) which was derived from the discretized TDGL equation [30,31]. Here, the Laplacian term is replaced by its isotropic discretized counterpart. And, the freeenergy functional in Eq. (27) is given by UCDM½fðrÞ kBT ¼ Z dr A lnðcosh fÞ þ 1 2 f2 þ D 2 jrfðrÞj2 . ð29Þ This model reproduces the growth kinetics of the TDGL model, demonstrating that such quantities ARTICLE IN PRESS Q.H. Zeng et al. / Prog. Polym. Sci. 33 (2008) 191–269 199
200 O.H.Zeng et al.Prog.Polym.Sci.33 (2008)191-269 are in is the calar order parameter. the di PA-PB ee-energy term. d p investigate the phase-separation of er nand M is the mobility for the order p Variable describes the asymmetry of the diblock. nanoparticles [32-36]. is the noise field.and parameter I determines the thickness of the lamellar domains which is related to 22 5 Dynamic DET method the length of the diblock.The free energy U()is Dynamic DFT method is usually used to model given by the local free energy U and polymer- the dynamic behavior of polymer systems and has particle coupling interaction term Up: U(6)=U:++Ucp (31 m th with TDGL mo el U=d/(,》+号v,2 (32) d However.in contrast to traditional phenomenolo gical free-e nergy expansion methods employed in Up=Cdr∑Vv-R(,)-中，)， (33) the TDGL approach,the free energy is not truncated at a certain level,and instead retains the =1 is the value of the full polymer path integral numerically. At the enging comput un expen The motion of particle is described by the me mply Hug Langevin equation addition.viscoelasticity.which is no =M,,-a0)+n (34 OR included in TDGL approaches,is included at the level of the Gaussian chains.A similar DFT th approach has been developed by Doi and co workers [38,39]and forms the basis represe ts a pa e no their nev software tool Simulation Utilities for Soft and Hard Inter 2.3.Mesoscale and macroscale methods at s40 Despite the importance of understanding the The of dyn nic DeT method is that the molecular structure and nature of materials,their instantaneous unique conformation distribution ca behavior can be homogenized with respect to be obtained from the off-equilibrium density profile different ects which can be at different scales by coupling a fictitious external potential to the the observed macrosc behavior Is Hamiltonian.Once such distribution is knowr the the energy is ther calculated by standar The orce me.The continuur m rial is th ned to gra ent of t fur an a subiected to the pody forces such as a ity and surface forces equations for both nolymer and particle in the Generally speaking.the macroscale methods(or diblock polymer-particle composites 1351 alled cont uum methods hereafter)obey the The dynamics of microphase separation in a melt undamental ws of:(1)continuity. derived from of diblocks is described through the following onservatio of ma equation [341: 8r,0 chang (30) an arbitrary oint is equal to the resultant moment
are insensitive to the precise form of the double-well potential of the bulk free-energy term. The TDGL and CDM methods have recently been used to investigate the phase-separation of polymer nanocomposites and polymer blends in the presence of nanoparticles [32–36]. 2.2.5. Dynamic DFT method Dynamic DFT method is usually used to model the dynamic behavior of polymer systems and has been implemented in the software package MesodynTM from Accelrys [37]. The DFT models the behavior of polymer fluids by combining Gaussian mean-field statistics with a TDGL model for the time evolution of conserved order parameters. However, in contrast to traditional phenomenological free-energy expansion methods employed in the TDGL approach, the free energy is not truncated at a certain level, and instead retains the full polymer path integral numerically. At the expense of a more challenging computation, this allows detailed information about a specific polymer system beyond simply the Flory–Huggins parameter and mobilities to be included in the simulation. In addition, viscoelasticity, which is not included in TDGL approaches, is included at the level of the Gaussian chains. A similar DFT approach has been developed by Doi and coworkers [38,39] and forms the basis for their new software tool Simulation Utilities for Soft and Hard Interfaces (SUSHI), one of a suite of molecular and mesoscale modeling tools (called OCTA) developed for the simulation of polymer materials [40]. The essence of dynamic DFT method is that the instantaneous unique conformation distribution can be obtained from the off-equilibrium density profile by coupling a fictitious external potential to the Hamiltonian. Once such distribution is known, the free energy is then calculated by standard statistical thermodynamics. The driving force for diffusion is obtained from the spatial gradient of the first functional derivative of the free energy with respect to the density. Here, we describe briefly the equations for both polymer and particle in the diblock polymer–particle composites [35]. The dynamics of microphase separation in a melt of diblocks is described through the following equation [34]: qfðr;tÞ qt ¼ Mr2 dU½fðr;tÞ dfðr;tÞ G½fðr;tÞ l þ xðr;tÞ, (30) where f is the scalar order parameter, f ¼ rArB, i.e., the difference of the respective local densities of the A and B segments of AB diblocks. The constant M is the mobility for the order parameter field. Variable l describes the asymmetry of the diblock. x is the noise field, and parameter G determines the thickness of the lamellar domains which is related to the length of the diblock. The free energy U(f) is given by the local free energy Ul and polymer– particle coupling interaction term Ucpl: UðfÞ ¼ Ul þ þUcpl, (31) Ul ¼ Z dr f 1ðfðr; tÞÞ þ D 2 ðrfðr; tÞÞ2 , (32) Ucpl ¼ C Z dr X i Vðr RiÞðfðr;tÞ fsÞ 2 , (33) where fs ¼ 1 is the value of the order parameter at the particle surface, and V(r) is the potential functional. The motion of particle is described by the Langevin equation R_i ¼ Mp f i qUðfÞ qRi þ Zi, (34) where Mp is the particle mobility, fi is the force acting on the particle from other particles, and Z represents a Gaussian white noise. 2.3. Mesoscale and macroscale methods Despite the importance of understanding the molecular structure and nature of materials, their behavior can be homogenized with respect to different aspects which can be at different scales. Typically, the observed macroscopic behavior is usually explained by ignoring the discrete atomic and molecular structure and assuming that the material is continuously distributed throughout its volume. The continuum material is thus assumed to have an average density and can be subjected to body forces such as gravity and surface forces. Generally speaking, the macroscale methods (or called continuum methods hereafter) obey the fundamental laws of: (i) continuity, derived from the conservation of mass; (ii) equilibrium, derived from momentum considerations and Newton’s second law; (iii) the moment of momentum principle, based on the model that the time rate of change of angular momentum with respect to an arbitrary point is equal to the resultant moment; ARTICLE IN PRESS 200 Q.H. Zeng et al. / Prog. Polym. Sci. 33 (2008) 191–269