ScienceDirect ELsevier Prog.Polym.Sai.33(2008)191-269 www.cl vier.com/locate/ppolys Multiscale modeling and simulation of polymer nanocomposites Q.H.Zeng,A.B.Yu.*G.Q.Lub ARC Centre of Excellence for Functional Nanomaterials.The University of Queensland.Brisbane.OLD 4072.Australia Polymer nanocor the reinforcement nd ti ng and simt that have been applied to polymer nanocomposites,covring from cale (e.g molecular dynamics,Monte Carlo) (e.g., rownian dynam dissipa tive particle dynamic lattice Boltzmann,time-dependent Ginzburg- polymer nanocomposites,including the thermodynamics and kinetics of formation,molecular structure and dynamics morphology.proc The present review aims to summariz the recent advance in the fundamental understandin of polymer nanocomposites 207ti,nanorube.day patelets)and stimate fuher ihsarea. Keywords:Polymer nanocomposite:Multiscale:Modeling:Simulation:Nanostructure:Nanoparticles Contents 18 2.1.Molecular scale methods 194 olecular dynamics Monte Carlo.......................................................... 196 Abbreviations:BD:Br nod:CHC:Cah onal theory D T629385442gRn60295956 E-mail address a.yu@unsw.eduau (A.B.Yu). 0079-6700/S-see front m er2007 Elsevier Ltd.All rights reserved. doi:10.1016/j-progpolymsci.2007.09.002
Prog. Polym. Sci. 33 (2008) 191–269 Multiscale modeling and simulation of polymer nanocomposites Q.H. Zenga , A.B. Yua,, G.Q. Lub a Centre for Simulation and Modeling of Particulate Systems and School of Materials Science and Engineering, The University of New South Wales, Sydney, NSW 2052, Australia b ARC Centre of Excellence for Functional Nanomaterials, The University of Queensland, Brisbane, QLD 4072, Australia Received 13 December 2006; received in revised form 5 September 2007; accepted 5 September 2007 Available online 3 December 2007 Abstract Polymer nanocomposites offer a wide range of promising applications because of their much enhanced properties arising from the reinforcement of nanoparticles. However, further development of such nanomaterials depends on the fundamental understanding of their hierarchical structures and behaviors which requires multiscale modeling and simulation strategies to provide seamless coupling among various length and time scales. In this review, we first introduce some computational methods that have been applied to polymer nanocomposites, covering from molecular scale (e.g., molecular dynamics, Monte Carlo), microscale (e.g., Brownian dynamics, dissipative particle dynamics, lattice Boltzmann, time-dependent Ginzburg–Landau method, dynamic density functional theory method) to mesoscale and macroscale (e.g., micromechanics, equivalent-continuum and self-similar approaches, finite element method). Then, we discuss in some detail their applications to various aspects of polymer nanocomposites, including the thermodynamics and kinetics of formation, molecular structure and dynamics, morphology, processing behaviors, and mechanical properties. Finally, we address the importance of multiscale simulation strategies in the understanding and predictive capabilities of polymer nanocomposites in which few studies have been reported. The present review aims to summarize the recent advances in the fundamental understanding of polymer nanocomposites reinforced by nanofillers (e.g., spherical nanoparticles, nanotubes, clay platelets) and stimulate further research in this area. r 2007 Elsevier Ltd. All rights reserved. Keywords: Polymer nanocomposite; Multiscale; Modeling; Simulation; Nanostructure; Nanoparticles Contents 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192 2. Modeling and simulation techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194 2.1. Molecular scale methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194 2.1.1. Molecular dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194 2.1.2. Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196 ARTICLE IN PRESS www.elsevier.com/locate/ppolysci 0079-6700/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.progpolymsci.2007.09.002 Abbreviations: BD: Brownian dynamics; BEM: boundary element method; CDM: cell dynamic method; CHC: Cahn–Hilliard–Cook; CNT: carbon nanotube; DFT: density functional theory; DPD: dissipative particle dynamics; FEM: finite element method; LB: lattice Boltzmann; MC: Monte Carlo; MD: molecular dynamics; MM: molecular mechanics; NMR: nuclear magnetic resonance; ODT: order–disorder phase transition; PCL: poly (e-caprolactone); PMMA: polymethylmethacrylate; RVE: representative volume element; SCF: self-consistent field; SWCN: single-walled carbon nanotube; TDGL: time-dependent Ginsburg–Laudau Corresponding author. Tel.: +61 2 93854429; fax: +61 2 93855956. E-mail address: a.yu@unsw.edu.au (A.B. Yu).
192 O.H.Zeng et al.Prog.Polym.Sci.33 (2008)191-269 2.2 oscale method Brownian dynamics 187 Dissipative particle dynamics.. attice B ndent Cinburg-Landau method. 225 Dynamie DFT method 2.3. e and macroscale methods 132 self-similar ppr 2.3.3 Finite element method. 204 3.Modeling and simulation of polymer nanocomposites.................................... 32 Nanocomposite thermodynamics........................ Nanocomposite molecular structure and dynamic properties. 213 osite morphology 342 3 Nanocomposite rheological and processing behaviors. 222 6 al properties .......... 。。。，。。。。。。 36 Continuum models 3.6.3 Equivalent-continuum and self-similar models 42 Sequential and concurrent approaches 4.3 Current research status......................... 155 Ks. References 261 1.Introduction 2 the hierarchical characteristics of the structure and dynamics of polymer nanocomposites ran- Polymer materials reinforced with nanoparticles ging from molecular scale,microscale to mesos- (e.g..nanosphere,nanotube,clay platelet)have recently cale and macroscale,in particular,the molecular received tremendous attention in both scientific and structures and dynamics at the interface between industrial communit s due to their extraordinary nanoparticles and polymer matrix: 0J.However,fror point on the ion of nanopart fabrication of polym sites.Tpu mechan ment of such materials is still largely empirical and a isms of nanoparticles in polymer nanocomposites. finer degree of control of their properties cannot be achieved so far.Therefore.computer modeling and simulation will play an ever-increasing role in predict The purpose of this review is to discuss the application of modeling and simulation techniques to ing and designing material properties,and guiding suck polymer nanocomposites.This includes a broad subject covering methodologies at various length and time ntal is We les and many aspects o for th 1.the thermodynamics and kinetics of the forma an he ughly divided int tion of polymer nanocomposites:
2.2. Microscale methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197 2.2.1. Brownian dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197 2.2.2. Dissipative particle dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197 2.2.3. Lattice Boltzmann . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 198 2.2.4. Time-dependent Ginzburg–Landau method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199 2.2.5. Dynamic DFT method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 200 2.3. Mesoscale and macroscale methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 200 2.3.1. Micromechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201 2.3.2. Equivalent-continuum and self-similar approaches. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202 2.3.3. Finite element method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 204 3. Modeling and simulation of polymer nanocomposites. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205 3.1. Nanocomposite thermodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205 3.2. Nanocomposite kinetics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 209 3.3. Nanocomposite molecular structure and dynamic properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213 3.4. Nanocomposite morphology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222 3.4.1. Homopolymer nanocomposites . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222 3.4.2. Block copolymer nanocomposites. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 227 3.5. Nanocomposite rheological and processing behaviors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233 3.6. Nanocomposite mechanical properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 236 3.6.1. Molecular models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 237 3.6.2. Continuum models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 244 3.6.3. Equivalent-continuum and self-similar models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 251 4. Multiscale strategies for modeling polymer nanocomposites. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253 4.1. Challenges . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253 4.2. Sequential and concurrent approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 254 4.3. Current research status . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 255 5. Concluding remarks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 260 Acknowledgment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 261 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 261 1. Introduction Polymer materials reinforced with nanoparticles (e.g., nanosphere, nanotube, clay platelet) have recently received tremendous attention in both scientific and industrial communities due to their extraordinary enhanced properties [1–10]. However, from the experimental point of view, it is a great challenge to characterize the structure and to manipulate the fabrication of polymer nanocomposites. The development of such materials is still largely empirical and a finer degree of control of their properties cannot be achieved so far. Therefore, computer modeling and simulation will play an ever-increasing role in predicting and designing material properties, and guiding such experimental work as synthesis and characterization. For polymer nanocomposites, computer modeling and simulation are especially useful in addressing the following fundamental issues: 1. the thermodynamics and kinetics of the formation of polymer nanocomposites; 2. the hierarchical characteristics of the structure and dynamics of polymer nanocomposites ranging from molecular scale, microscale to mesoscale and macroscale, in particular, the molecular structures and dynamics at the interface between nanoparticles and polymer matrix; 3. the dependence of polymer rheological behavior on the addition of nanoparticles, which is useful in optimizing processing conditions; and 4. the molecular origins of the reinforcement mechanisms of nanoparticles in polymer nanocomposites. The purpose of this review is to discuss the application of modeling and simulation techniques to polymer nanocomposites. This includes a broad subject covering methodologies at various length and time scales and many aspects of polymer nanocomposites. We organize the review as follows. In Section 2, we introduce briefly the computational methods used so far for the systems of polymer nanocomposites which can be roughly divided into three types: molecular scale methods (e.g., molecular dynamics (MD), Monte Carlo ARTICLE IN PRESS 192 Q.H. Zeng et al. / Prog. Polym. Sci. 33 (2008) 191–269
Q.H.Zeng et al.Prog.Polym.Sci.33 (2008)191-269 193 Nomenclature S average compliance of polymer compo- site strain-concentration tensor of polymer T temperature composite U interaction potential A area of interface interaction △U change of potential energy acceleration of atom UGL()Ginzburg-Landu free energy B stress-concentration tensor of polymer UCDM free energy calculated in cell dynamical composite method O stiffness tensor of polymer composite UK kinetic energy in finite element model Cr elastic stiffness tensor of the effective Uv Hookian potential energy in finite ele- fiber in polymer composite ment model E total modulus of polymer composite Um total potential energy of molecule model E11 longitudinal modulus of polymer com- U total potential energy of truss model posite U Udh energies associated with covalent E22 transverse modulus of polymer compo- bonding stretching,bond-angle bending, site and van der Waals interactions,respec- Er modulus of filler in polymer composite tively Em modulus of matrix in polymer composite UU energies associated with truss ele- ei local particle velocity in lattice gas ments that represent covalent bonding automaton method,or lattice velocity stretching,bond-angle bending,and van vector in lattice Boltzmann method der Waals interactions,respectively force acting on the ith atom or particle displacement field in finite element model Fyc conservative force of particle j acting on V(r) potential functional particle i Vt volume of filler in polymer composite dissipative force of particle j acting on volume of matrix in polymer composite particle i friction coefficient FyR random forces of particle j acting on interface binding energy particle i E volume-average strain of polymer com- f(x,t)distribution function of the particles of posite component a with velocity e;along 可 volume-average strain of filler in poly- direction i at lattice x and time t mer composite (x,)equilibrium distribution function of the 高m volume-average strain of matrix in poly- particles of component a mer composite G chain gyration tensor infinitesimal strain tensor of filler in H Hamiltonian associated with a configura- polymer composite tion infinitesimal strain tensors of matrix in △H change of system Hamiltonian or free polymer composite energy random noise term with zero mean h,ho final and initial clay platelet separation (r,t) thermal noise term with zero mean kB Boltzmann constant Gaussian white noise M mobility of the order parameter field enthalpic interaction parameter Mp particle mobility 0 volume fraction of filler in polymer 171 atomic mass composite nx,t)particle occupation variable with direc- Um volume fraction of matrix in polymer tion i at distance x and time t composite P average segment orientation velocity of atom Pi momentum of particle 5 random number between 0 and I which is Pi-→j probability of accepting a new config- to determine the acceptance or rejection uration ( of a new configuration position of ith atom or particle 5 shape parameter depending on filler geometry and loading direction
ARTICLE IN PRESS Nomenclature A strain–concentration tensor of polymer composite Ai area of interface interaction a * acceleration of atom B stress–concentration tensor of polymer composite C stiffness tensor of polymer composite Cf elastic stiffness tensor of the effective fiber in polymer composite E total modulus of polymer composite E11 longitudinal modulus of polymer composite E22 transverse modulus of polymer composite Ef modulus of filler in polymer composite Em modulus of matrix in polymer composite ei local particle velocity in lattice gas automaton method, or lattice velocity vector in lattice Boltzmann method F ~i force acting on the ith atom or particle Fij C conservative force of particle j acting on particle i Fij D dissipative force of particle j acting on particle i Fij R random forces of particle j acting on particle i f s i ðx;tÞ distribution function of the particles of component s with velocity ei along direction i at lattice x and time t f s;eq i ðx; tÞ equilibrium distribution function of the particles of component s Gij chain gyration tensor H Hamiltonian associated with a configuration DH change of system Hamiltonian or free energy h, h0 final and initial clay platelet separation kB Boltzmann constant M mobility of the order parameter field Mp particle mobility mi atomic mass ni(x, t) particle occupation variable with direction i at distance x and time t P2 s average segment orientation pi momentum of particle pi-j probability of accepting a new configuration (j) ~ri position of ith atom or particle S average compliance of polymer composite T temperature U interaction potential DU change of potential energy UGL(f) Ginzburg–Landu free energy UCDM free energy calculated in cell dynamical method UK kinetic energy in finite element model UV Hookian potential energy in finite element model Um total potential energy of molecule model Ut total potential energy of truss model Ur , Uy , Uvdw energies associated with covalent bonding stretching, bond-angle bending, and van der Waals interactions, respectively Ua , Ub , Uc energies associated with truss elements that represent covalent bonding stretching, bond-angle bending, and van der Waals interactions, respectively u_ displacement field in finite element model V(r) potential functional Vf volume of filler in polymer composite Vm volume of matrix in polymer composite g friction coefficient gi interface binding energy ¯ volume–average strain of polymer composite ¯f volume–average strain of filler in polymer composite ¯m volume–average strain of matrix in polymer composite e f infinitesimal strain tensor of filler in polymer composite e m infinitesimal strain tensors of matrix in polymer composite zij random noise term with zero mean z(r, t) thermal noise term with zero mean Z Gaussian white noise l enthalpic interaction parameter vf volume fraction of filler in polymer composite vm volume fraction of matrix in polymer composite v * velocity of atom x random number between 0 and 1 which is to determine the acceptance or rejection of a new configuration xf shape parameter depending on filler geometry and loading direction Q.H. Zeng et al. / Prog. Polym. Sci. 33 (2008) 191–269 193
194 Q.H.Zeng et al.Prog.Polym.Sci.33 (2008)191-269 Pmax maximum stress along the tube axis dimensionless collision-relaxation time noise amplitude constant of component o of the fluid G volume-average stress of polymer com- x°=10/△t posite order parameter Gf volume-average stress of filler in polymer 中 value of the order parameter at particle composite surface Gm volume-average stress of matrix in poly- asymmetry of diblock copolymer mer composite 2 collision operation total stress of filler in polymer composite weight function for conservative force total stress of matrix in polymer compo- WD weight function for dissipative force site OR weight function for random force (MC)),microscale methods (e.g.,Brownian dynamics clusters as the basic units considered.The most (BD),dissipative particle dynamics (DPD),lattice popular methods include molecular mechanics Boltzmann (LB),time-dependent Ginzburg-Lanau (MM),MD and MC simulation.Modeling of method,dynamic density functional theory (DFT) polymer nanocomposites at this scale is predomi- method),and mesoscale and macroscale methods (e.g., nantly directed toward the thermodynamics and micromechanics,equivalent-continuum and self-similar kinetics of the formation,molecular structure and approaches,finite element method(FEM)).We do not interactions.The diagram in Fig.I describes the aim to provide a detailed description of each method equation of motion for each method and the typical but its basic principles,strengths and weaknesses,and properties predicted from each of them [11-16].We potential applications.Interesting readers can refer to introduce here the two widely used molecular scale relevant books,reviews and research articles for details. methods:MD and MC. In Section 3,we discuss in detail the applications of these modeling and simulation methods in some specific aspects of polymer nanocomposites,including 2.1.1.Molecular dynamics the thermodynamics and kinetics of the formation, MD is a computer simulation technique that allows molecular structure and dynamics,morphologies (i.e., one to predict the time evolution of a system of phase behaviors),rheological and processing behaviors, interacting particles (e.g.,atoms,molecules,granules, and mechanical properties.We pay more attention to etc.)and estimate the relevant physical properties clay-based polymer nanocomposites because of their [17,18].Specifically,it generates such information as importance in polymer nanocomposites and our own atomic positions,velocities and forces from which the research interests.Of course,we also refer to research macroscopic properties (e.g.,pressure,energy,heat activities in polymer systems reinforced by nanospheres capacities)can be derived by means of statistical and nanotubes,an area developed rapidly in the past mechanics.MD simulation usually consists of three few years.In Section 4,we highlight the importance of constituents:(i)a set of initial conditions (e.g.,initial multiscale strategies of modeling and simulation in positions and velocities of all particles in the system); understanding and predicting the hierarchical structure (ii)the interaction potentials to represent the forces and behaviors arising from the polymer and nanopar- among all the particles;(iii)the evolution of the ticle together with the properties observed at various system in time by solving a set of classical Newtonian scales.We discuss the current applications of some equations of motion for all particles in the system. multiscale methods in polymer nanocomposites.Final- The equation of motion is generally given by ly,we conclude the review by emphasizing the current d27 challenges and future research directions. 下(0=mdr (1) 2.Modeling and simulation techniques where F;is the force acting on the ith atom or particle at time t which is obtained as the negative gradient of 2.1.Molecular scale methods the interaction potential U,m;is the atomic mass and Fi the atomic position.A physical simulation involves The modeling and simulation methods at mole- the proper selection of interaction potentials,numer- cular level usually employ atoms,molecules or their ical integration,periodic boundary conditions,and
(MC)), microscale methods (e.g., Brownian dynamics (BD), dissipative particle dynamics (DPD), lattice Boltzmann (LB), time-dependent Ginzburg–Lanau method, dynamic density functional theory (DFT) method), and mesoscale and macroscale methods (e.g., micromechanics, equivalent-continuum and self-similar approaches, finite element method (FEM)). We do not aim to provide a detailed description of each method but its basic principles, strengths and weaknesses, and potential applications. Interesting readers can refer to relevant books, reviews and research articles for details. In Section 3, we discuss in detail the applications of these modeling and simulation methods in some specific aspects of polymer nanocomposites, including the thermodynamics and kinetics of the formation, molecular structure and dynamics, morphologies (i.e., phase behaviors), rheological and processing behaviors, and mechanical properties. We pay more attention to clay-based polymer nanocomposites because of their importance in polymer nanocomposites and our own research interests. Of course, we also refer to research activities in polymer systems reinforced by nanospheres and nanotubes, an area developed rapidly in the past few years. In Section 4, we highlight the importance of multiscale strategies of modeling and simulation in understanding and predicting the hierarchical structure and behaviors arising from the polymer and nanoparticle together with the properties observed at various scales. We discuss the current applications of some multiscale methods in polymer nanocomposites. Finally, we conclude the review by emphasizing the current challenges and future research directions. 2. Modeling and simulation techniques 2.1. Molecular scale methods The modeling and simulation methods at molecular level usually employ atoms, molecules or their clusters as the basic units considered. The most popular methods include molecular mechanics (MM), MD and MC simulation. Modeling of polymer nanocomposites at this scale is predominantly directed toward the thermodynamics and kinetics of the formation, molecular structure and interactions. The diagram in Fig. 1 describes the equation of motion for each method and the typical properties predicted from each of them [11–16]. We introduce here the two widely used molecular scale methods: MD and MC. 2.1.1. Molecular dynamics MD is a computer simulation technique that allows one to predict the time evolution of a system of interacting particles (e.g., atoms, molecules, granules, etc.) and estimate the relevant physical properties [17,18]. Specifically, it generates such information as atomic positions, velocities and forces from which the macroscopic properties (e.g., pressure, energy, heat capacities) can be derived by means of statistical mechanics. MD simulation usually consists of three constituents: (i) a set of initial conditions (e.g., initial positions and velocities of all particles in the system); (ii) the interaction potentials to represent the forces among all the particles; (iii) the evolution of the system in time by solving a set of classical Newtonian equations of motion for all particles in the system. The equation of motion is generally given by F ~iðtÞ ¼ mi d2 ~ri dt2 , (1) where F ~i is the force acting on the ith atom or particle at time t which is obtained as the negative gradient of the interaction potential U, mi is the atomic mass and ~ri the atomic position. A physical simulation involves the proper selection of interaction potentials, numerical integration, periodic boundary conditions, and ARTICLE IN PRESS rmax maximum stress along the tube axis s noise amplitude s¯ volume–average stress of polymer composite s¯f volume–average stress of filler in polymer composite s¯ m volume–average stress of matrix in polymer composite sf total stress of filler in polymer composite sm total stress of matrix in polymer composite t s dimensionless collision–relaxation time constant of component s of the fluid, ts ¼ ls =Dt f order parameter fs value of the order parameter at particle surface w asymmetry of diblock copolymer O collision operation oC weight function for conservative force oD weight function for dissipative force oR weight function for random force 194 Q.H. Zeng et al. / Prog. Polym. Sci. 33 (2008) 191–269
Q.H.Zeng et al Prog.Polym.SeL 33(2008)191-269 9 QM MD MC otion of electron,nuclei motion of atoms motion of atom chroainger equat 70 ●3 ctronis distribution anoparticle surface Bonding and unbonding micdifasioncocfec hterealaioakineties Mechanical properties ethods and lo n of a We-VCYn The interaction potentials together with their parameters,i.e.,the so-called force field,describe in detail how the particles in a system interact with each other,i.e.,how the potential energy of a system depends on the particle coordinates.Such a force may be obt me (e.g., (2) me0ho order potential).The criteria for selec The first four terms represent bonded interac- field include the accuracy.transferability and tions,i.e.,bond stretching ond,bond-angle bend computational speed.A typical interaction potential fhreioanierqio and dihedral angle torsion and U may consist of a number of bonded and non- whil the last two bonded interaction terms interactions 1.e.,van de energy and el energy a U,2,,fN)= involved in a given interaction:NN and N. stand for the total numbers of these respective interactions in the simulated system; and linrersion uniquely specily an indi idual interaction of each type;i and j in the
the controls of pressure and temperature to mimic physically meaningful thermodynamic ensembles. The interaction potentials together with their parameters, i.e., the so-called force field, describe in detail how the particles in a system interact with each other, i.e., how the potential energy of a system depends on the particle coordinates. Such a force field may be obtained by quantum method (e.g., ab initio), empirical method (e.g., Lennard–Jones, Mores, Born-Mayer) or quantum-empirical method (e.g., embedded atom model, glue model, bondorder potential). The criteria for selecting a force field include the accuracy, transferability and computational speed. A typical interaction potential U may consist of a number of bonded and nonbonded interaction terms: Uð~r1;~r2; ... ;~rNÞ ¼ N Xbond ibond Ubond ðibond ;~ra;~rbÞ þ N X angle iangle Uangleðiangle;~ra;~rb;~rcÞ þ N Xtorsion itorsion Utorsionðitorsion;~ra;~rb;~rc;~rd Þ þ NXinversion iinversion Uinversionðiinversion;~ra;~rb;~rc;~rd Þ þ N X1 i¼1 X N j4i Uvdwði; j;~ra;~rbÞ þ N X1 i¼1 X N j4i Uelectrostaticði; j;~ra;~rbÞ. ð2Þ The first four terms represent bonded interactions, i.e., bond stretching Ubond, bond-angle bend Uangle and dihedral angle torsion Utorsion and inversion interaction Uinversion, while the last two terms are non-bonded interactions, i.e., van der Waals energy Uvdw and electrostatic energy Uelectrostatic. In the equation, ~ra;~rb;~rc, and ~rd are the positions of the atoms or particles specifically involved in a given interaction; Nbond, Nangle, Ntorsion and Ninversion stand for the total numbers of these respective interactions in the simulated system; ibond, iangle, itorsion and iinversion uniquely specify an individual interaction of each type; i and j in the van der Waals and electrostatic terms indicate the atoms involved in the interaction. ARTICLE IN PRESS Fig. 1. Molecular modeling and simulation methods commonly used for polymer nanocomposites [11–16]. Reproduced from (i) Lee, Baljon and Loring by permission of American Institute of Physics; (ii) Smith, Bedrov, Li and Byutner by permission of American Institute of Physics; (iii) Smith, Bedrov and Smith by permission of Elsevier Science Ltd.; (iv) Zeng, Yu, Lu and Standish by permission of American Chemical Society; (v) Vacatello by permission of Wiley-VCH; and (vi) Zeng, Yu and Lu by permission of Institute of Physics. Q.H. Zeng et al. / Prog. Polym. Sci. 33 (2008) 191–269 195
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 (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
