A study of chemical reactions in coarse-grained simulations

Abstract

We introduce a reaction model for use in coarse-grained simulations to study the chemical reactions in polymer systems at mesoscopic level. In this model, we employ an idea of reaction probability in control of the whole process of chemical reactions. This model has been successfully applied to the studies of surface initiated polymerization process and the network structure formation of typical epoxy resin systems. It can be further modified to study different kinds of chemical reactions at mesoscopic scale.

Keywords

coarse-grained simulation / reaction probability / surface initiated polymerization / curing reaction

Introduction

As a feasible route for validating experimental results and predicting new phenomena, computer simulation becomes especially important on presenting physical insights in modern polymer science. Because of the broad coverage of time and length scales in most of the polymer systems, coarse-grained (CG) models from atomistic simulations are often utilized to describe the polymers on mesoscale, for example, the coarse-grained molecular dynamics (CGMD) and dissipative particle dynamics (DPD) methods, etc. In some cases, independent kinetic processes, e.g., chemical reactions, may accompany the dynamic processes of complex fluids. A lot of studies have shown that the reactions may greatly disturb the behavior of diffusion [1-4]. Therefore, how to associate hydrodynamic process of the fluids with chemical reactions becomes an important problem in the study of complex fluids. Especially in polymeric systems, it is of great significance to simulate the reaction equilibrium at mesoscopic level, since most of the dynamical problems of polymers are generally related to the mesoscopic length and time scales. Moreover, designing the polymer chains with controllable chain length distribution is a critical issue in polymer sciences. These problems urgently demand computational models in which the chemical reactions can be taken into account.

In recent years, more and more researchers have focused on developing new reaction models for use in polymeric systems, especially on how to design suitable models so that the typical polymerization processes could be reasonably reproduced. One of the attempts on this issue was done by Lísal and the coworkers [5,6], who had proposed a reaction ensemble DPD method to study the polydisperse polymer systems and supramolecular diblock copolymers. In their model, based on DPD equations of motion, they considered the random changes in the coupling parameter to govern the time evolution of the reaction system, which would be accepted with a probability derived from grand canonical partition function for the expanded ensemble. Another idea was completely based on the Monte Carlo simulations, suggested by, for example, He et al. [7]. In their model, first they presented all kinds of possible reactions in the polymeric systems (e.g., the chain initiation, propagation, termination, as well as the chain transfer, etc.), and then normalized the reaction rate coefficients so that they could obtain a reaction probability for each kind of reaction. Thus, in each reaction step, they selected one kind of the reactions to occur according to the reaction probability. This model ignored real kinetic reaction processes, thus it was not very suitable for the study of the structural characteristic variations during the reaction. For reproducing the real propagation of polymer chains, some models are related to the design of many local processes with the active growth ends connecting to free monomers in their capture radii. Müller-Plathe and the coworkers [8] employed two parameters to control this process, i.e., the characteristic delay time separating two successive reaction steps, and the capture radii defining the geometrical conditions for chain initiation (propagation). In each step, the active end will be connected with the closest free monomer in its capture radius. They applied this reaction model to study the chain propagation of styrene to polystyrene. In fact, the capture radii in their model should be a function of reaction energy barrier, and the delay time should be inversely proportional to the reaction rate. Thus it was not easy for this model to bring intuitive definition for the parameters. In Ref. [9], the authors proposed the way of polymer chain generation using radical-like polymerization. In their model, they let the active ends randomly search one of the first monomer neighbors to connect so that the chains could grow. For typical free radical polymerization, the reaction is very fast as compared to the diffusion of particles in CG simulations. Thus the reaction turns out to be a pure probability issue of judging whether the polymerization can take place within one time step or not. It implies that we may construct a probability-based reaction model in the simulations. Genzer et al. [10,11] introduced the reaction and the motion probabilities P_x to model the key steps governing the free radical polymerization. They also admitted that it was challenging to derive exact connections between the reaction probability and the reaction rate coefficient. Several other research groups also applied the reaction probability idea to construct the reaction protocols, e.g., He et al. [12] recently proposed a diffusion-limited reaction model to study the structure of hyperbranched polymers prepared by slow monomer addition to a core molecule, in which they introduced the absolute reaction probability to control the reaction. Berezkin et al. [13], Ding et al. [14] and Mukherji et al. [15] also adopted similar reaction probability concept to establish their respective reaction models. However, notably, few of them could rationally present the derivation of the physical reaction rate from their local reaction probability parameters. On this issue, Toxvaerd and the coworkers [16] suggested their polymerization model in which they presented a derivation from the preset parameters to the reaction rate coefficients. In this paper, we employ a similar idea as that in Ref. 16, and then establish our distinctive reaction model. We also construct the relationship between the preset reaction probability and the reaction rate, which offers possible ways to obtain the reaction parameters in real experiments.

We demonstrate our model by applying it to the studies of surface initiated polymerization (SIP) processes and the network structure formation in typical epoxy resin systems. In the SIP model, by choosing different initiator densities and different polymerization probabilities, we are able to obtain the polymer brushes with different degrees of polydispersity and different chain length distributions. In the study of network structure formation in typical epoxy resin systems, we design a CG scheme to map from all-atom molecular model to the mesoscopic DPD model, from which we can successfully construct a chemically reasonable network structure with conversion comparable to experiments. Moreover, we introduce a few applications of our CG reaction model on the studies of some other typical complex fluid systems. In Section 2, we introduce two common coarse-grained simulation methods adopted in our study, the reaction model, and the derivation of reaction rate. In Section 3, we mainly show the applications of the model, and in Section 4, we present the concluding remarks.

Simulation methods

Coarse-grained simulation method

Two commonly used simulation methods are employed here, i.e., the DPD method and the CGMD method. In DPD, the motion of interacting CG particles is governed by Newton’s equations of motion [17]. Interparticle interactions are characterized by pair wise conservative, dissipative, and random forces acting on particle i by particle j. They are given by

(1)F_{ij}^C=-{\alpha }_{ij}{\omega }^C(r_{ij})e_{ij}

(2){\mathbf{F}}_{ij}^D=-\gamma {\omega }^D\left(r_{ij}\right)\left({\mathbf{v}}_{ij}·{\mathbf{e}}_{ij}\right){\mathbf{e}}_{ij}

(3){\mathbf{r}}_{ij}={\mathbf{r}}_i-{\mathbf{r}}_j

where {\mathbf{r}}_{ij}={\mathbf{r}}_i-{\mathbf{r}}_j, r_{ij}=\left|{\mathbf{r}}_{ij}\right|, {\mathbf{e}}_{ij}={\mathbf{r}}_{\iota j}/r_{ij}, and {\mathbf{v}}_{ij}={\mathbf{v}}_i-{\mathbf{v}}_j. {\xi }_{ij} is a random number with zero mean and unit variance. For easy numerical handling, the cutoff radius, the particle mass, and the temperature are often set to be the units, i.e., r_c=m=k_{B}T=\mathrm{1}. {\alpha }_{ij} is the repulsion strength which describes the maximum repulsion between interacting particles. {\omega }^C, {\omega }^D, and {\omega }^R are three weight functions for the conservative, dissipative and random forces, respectively. For the conservative force, {\omega }^C(r_{ij})=\mathrm{1}-r_{ij}/r_c for and {\omega }^C(r_{ij})=\mathrm{0} for r_{ij}\ge r_c. {\omega }^D(r_{ij}) and {\omega }^R(r_{ij}) have a relation according to the fluctuation–dissipation theorem [18],

(4){\omega }^D(r_{ij})={[{\omega }^R(r_{ij})]}^{\mathrm{2}},

(5){\sigma }^{\mathrm{2}}=\mathrm{2}\gamma k_{B}T.

Here we choose a simple form of {\omega }^D and {\omega }^R due to Groot and Warren [19],

(6)

Groot–Warren-velocity Verlet algorithm [17,19] is used here to integrate the Newton’s equations of motion,

{\mathbf{r}}_i\left(t+\mathrm{\Delta }t\right)={\mathbf{r}}_i\left(t\right)+\mathrm{\Delta }t{\mathbf{v}}_i\left(t\right)+1/2{\left(\mathrm{\Delta }t\right)}^2{\mathbf{f}}_i\left(t\right)

{\tilde{\mathbf{v}}}_i\left(t+\mathrm{\Delta }t\right)={\mathbf{v}}_i\left(t\right)+\lambda \mathrm{\Delta }t{\mathbf{f}}_i\left(t\right)

{\mathbf{f}}_i\left(t+\mathrm{\Delta }t\right)={\mathbf{f}}_i\left[\mathbf{r}\left(t+\mathrm{\Delta }t\right),\tilde{\mathbf{v}}\left(t+\mathrm{\Delta }t\right)\right]

(7){\mathbf{v}}_i\left(t+\mathrm{\Delta }t\right)={\mathbf{v}}_i\left(t\right)+1/2\mathrm{\Delta }t\left({\mathbf{f}}_i\left(t\right)+{\mathbf{f}}_i\left(t+\mathrm{\Delta }t\right)\right)

Here \lambda =\mathrm{0.65} according to Ref. [19]. In DPD, polymers are constructed by connecting the neighboring particles together via harmonic springs {\mathbf{F}}_i^S={\sum }_{j}C{\mathbf{r}}_{ij}. We choose the spring constant C=\mathrm{10} according to Ref. [19]. In general, we choose the time step \Delta t=\mathrm{0.05} in the DPD simulations.

In CGMD, we consider the truncated and shifted Lennard-Jones (LJ) potential (i.e., the Weeks-Chandler-Andersen (WCA) potential) [20] between particles, in which the potential vanishes at the cutoff radius (r_c={\mathrm{2}}^{\mathrm{1}/\mathrm{6}}) so that only its repulsive part is left,

(8)U(r)=U_{LJ}(r)-U_{LJ}(r_c)

with

(9)U_{LJ}(r)=\mathrm{4}ϵ\left[{\left({\displaystyle \frac{\sigma }{r}}\right)}^{\mathrm{12}}-{\left({\displaystyle \frac{\sigma }{r}}\right)}^{\mathrm{6}}\right].

For simplicity, we define the units in the simulations by setting m=\mathrm{1},\sigma =\mathrm{1} and ϵ=\mathrm{1}. With WCA potential, the polymer chains are constructed by connecting the neighboring monomers together via finite extensible nonlinear elastic (FENE) spring potential [21,22],

(10)

where the spring constant is k=\mathrm{10}ϵ/{\sigma }^{\mathrm{2}}. Here r=\left|{\mathbf{r}}_i-{\mathbf{r}}_j\right| is the distance between neighboring monomers. The longest stretch of the bonds l_{\text{max}} is chosen according to the condition of bond uncrossibility [23,24] that \sqrt{\mathrm{2}}{r}_{\min }>l_{\max }. We can obtain the minimal distance between two particles, r_{\min }, from radial distribution function [24], then the choice of l_{\max } is clear. The DPD thermostat is employed here to control the constant temperature as that introduced above, since it was evaluated as a powerful technique for good thermostating [25,26]. In general, we choose the time step \Delta t=\mathrm{0.005} in CGMD simulations.

Combining reaction model in CG simulations

In the construction of the reaction model combined with CG simulations, we will introduce the idea of reaction probability in control of each time interval in the reaction process. In each reaction time interval \tau, if an active end meets several free monomers in the reaction radius, first it randomly chooses one of the monomers as a reacting object. Subsequently, another random number P is generated, then by checking if it is smaller than the preset reaction probability P_r, we decide whether the monomer will connect with the active end or not. This judging process in one reaction step is schematically illustrated in Fig. 1. If the bond can be formed between the active end and the reacting object, we record the connection information and update the spring forces between them.

This idea of reaction is especially suitable for the design of polymerization-type reactions. During the polymerization, the newly connected monomers then turn out to be the growth centers in the next propagation step of the same chain to connect other free monomers, so that the active end is transferred forward. We have employed this idea to construct several living polymerization models [27,28]. In our living polymerization models, the chain termination step is ignored as that in the real living polymerizations.

Assuming that we are now simulating a living polymerization-type reaction, there are basically three types of reactions included:

(11)I^{*}+M\to I-P^{*}

(12)I-P^{*}+M\to I-P-P^{*}

(13)I-P_{n-\mathrm{1}}-P^{*}+M\to I-P_n-P^{*}(n=\mathrm{2},\mathrm{3},\mathrm{4},\cdots )

in which I* and I denote the active and the dead initiator, respectively. While P* and P are the active growth chain end and the inner chain monomer, respectively, and M the free monomer. If we focus on the living polymerization, i.e., there is no chain termination reaction, the reaction rate r_p can be expressed as,

(14)r_p=-k[P^{*}][M]

where [P*] and [M] are the concentrations of active chain ends and free monomers, respectively, and k is the reaction rate coefficient.

As we know, the CG model inherently defines the time and length scales of the simulation, whereas notably the reaction may run with another independent clock. If we assume that we have also obtained the reaction time scale (in practice it is possible to get the approximate value by fitting the concentration change of one component or one kind of chemical group in the simulations to the real data in experiments), and occasionally in one time unit of the simulation, the chain propagation reaction takes place for \tau times. Based on the predefined reaction probability P_r and our uniformly generated random number, it is easy to conclude that the average decrease number of free monomers in one propagation step is N_{P^{*}}{P}_r, where N_{P^{*}} is the number of the active ends in the system. Since we are considering the living polymerization, any radical termination and bi-radical termination are omitted, thus the value of N_{P^{*}} is always constant, which is equivalent to the number of free radicals at the beginning. Therefore, we obtain the average decrease number of free monomers in one time unit as N_{P^{*}}{P}_r\tau. Consequently, the concentration change of the free monomers in one time unit is given by

(15)\text{d}[M]={\displaystyle \frac{N_{P^{*}}·P_r·\tau }{V·Na}}={\displaystyle \frac{[P^{*}]·P_r·\tau }{Na}}

where [M] is the concentration of free monomers, V is the volume of the system and Na is the Avogadro’s number. As a result, we can obtain the reaction rate as

(16)r_p=-\frac{d\left[M\right]}{\text{d}t}=\frac{\left[P^{*}\right]·P_r·\tau }{Na·t_0}

where t_{\mathrm{0}} is the real time unit in the simulation. Combining Eqs. (14) and (16), it is easy to obtain the reaction rate coefficient as,

(17)k={\displaystyle \frac{P_r·\tau }{Na·t_{\mathrm{0}}·[M]}}.

Assuming that the reaction is a classical Arrhenius-type, i.e., k=A\cdot \exp (-E_A/T) (here E_A is the activation energy), we can finally get the expression of P_r as

(18)P_r={\displaystyle \frac{A·\exp (-E_A/T)·Na·t_{\mathrm{0}}·[M]}{\tau }}.

This expression gives us two hints. First, the value of P_r is [M] dependent. In the consideration of our simulation length scale and time scale, it is reasonable to neglect the change of [M] when we are mapping the simulation onto an experiment. Thus the value of P_r should be nearly constant in the simulations. Second, there are two possible ways to get the value of P_r from experiments. For example, we may measure the activation energy E_A from experiments or use quantum chemistry calculations to estimate reaction rate coefficient k. Anyway, in this study, we assume that different reaction probabilities correspond to different reaction rates of the systems determined by different levels of reaction activation energies.

Applications and results

Surface initiated polymerization model

The first application relates to a simple CG model to simulate the surface initiated polymerization. Polymer brushes are generally defined as layers of polymer chains end-grafted to a surface. They give rise to a wide range of important technological and industrial applications, such as lubrications, oil recovery, colloid stability, modification of surfaces, adhesion, reversible tuning of wetting, biotechnology, and contact formation [25,29-32]. The structures and the polydispersity of polymer brushes are greatly dependent on some details of the surface initiated polymerization. In experiments, it is important to understand the effects of polymerization rate on the properties of polymer brushes in SIP reactions. As regard to the effects of polymerization rate (and most importantly, its interplay with the initiator density) on the properties of the polymer brushes, we present a simulation study on the SIP process by utilizing the reaction model and taking into consideration of the variation of polymerization rates in the simulations [28]. In this study, to highlight the excluded volume effect in the systems, WCA potential between interacting particles is adopted, while DPD thermostat is retained. We construct the model with two layers of planar walls at the top and the bottom of the simulation box by fixing wall particles. Some anchored active initiator sites are regularly distributed on the surface of the bottom wall layer. We then apply the polymerization probability to control the chain growth process. In our canonical ensemble simulations, the initial configuration is set to be the randomly distributed free monomers in the system, then \mathrm{1}\times {\mathrm{10}}^{\mathrm{5}} steps integration is carried out to relax the system. Consequently, a process of polymerization with \mathrm{1}\times {\mathrm{10}}^{\mathrm{6}} time steps is modeled. After a period of simulation, we can obtain the polydisperse polymer brushes with different morphologies and grafting densities, see Fig. 2.

In this study, we analyze the influence of both polymerization rate and initiator density on the grafting density, the chain conformation, the fraction of unsaturated initiators, and the polydispersity index (PDI) of the brush chains. First, we investigate the influence of the initiator density on the grafting density in this model. Fig. 3 shows the grafted monomer fraction versus the initiator density. At different polymerization rates, we find that the relationship between the initiator density and the grafting density is very similar to that reported in experiments [33-35]. An initiator density window can be identified in the range of 0.353-0.508, after which each curve arrives at a plateau, indicating a decrease of initiation efficiency. We also find that faster reaction yields higher plateau value. This is ascribed to that in the same period of polymerization, faster reaction is able to graft more free monomers.

We present the time evolution of the fraction of unsaturated initiators f_{N=\mathrm{1}} in Fig. 4. As shown in Fig. 4(a), in the early stage, for the systems with the same P_r, the curves almost coincide with each other, indicating the free polymerization stage without competition of ambient free monomers between the growing short chains. After this early stage, for the systems with higher initiator density ({\sigma }_i=\mathrm{0.508}), the decrease of f_{N=\mathrm{1}} slows down, showing that grown chains screen the unsaturated initiators from reacting to the monomers. It is ascribed to that in such systems, the monomers near the grafting surface are soon consumed out, while the grown chains form a compact layer to block the monomers far from the surface diffusing nearer. Thus, the unsaturated initiators have little opportunity to find free monomers, resulting in the retardness of the decrease of f_{N=\mathrm{1}}, as shown in Fig. 4. We then re-plot Fig. 4(a) in semi-log coordinates, as shown in Fig. 4(b). We find that at the same P_r, the curves for {\sigma }_i=\mathrm{0.127} and 0.508 can both be generally fitted by the same linear decay in the early stage of polymerization. Theoretically, if there is no screening effect of the grown chains on the unsaturated initiators, the number of unsaturated initiators after time t should be

(19)N_i^{\text{unsaturate}}=N_i{\left(1-P_r\right)}^k

where k=t/(\Delta t·\tau ) is the frequency of the polymerization since the beginning of the reaction and N_i is the total number of initiators. Therefore

(20)\log \left(f_{N=1}\right)=\log \left(N_i^{\text{unsaturate}}/N_i\right)=\left[\frac{\log \left(1-P_r\right)}{\Delta t·\tau }\right]t

We then can estimate the slopes for \log (f_{N=\mathrm{1}})~t from Eq. (20). We find in Fig. 4(b) that the slopes in the early stage are very close to the calculated values. We can imagine that in the early stage of the polymerization, there is no screening effect of grown chains on unsaturated initiators. At later times, the screening effect becomes obvious, especially for the system with higher initiator density, therefore we can find that the decay curves deviate from the theoretical slopes in Fig. 4(b).

We also study the PDI of the grafted chains after polymerization. Fig. 5 shows the PDI for each system with different initiator densities and polymerization rates. At the same initiator density, the systems with faster polymerization rates possess larger PDI, which supports the fact that strong screening effect of the grown chains occurs in these systems. For a specific polymerization rate, PDI is higher in the system with higher initiator density. Thus, if one needs to keep the polydispersity index of high density brushes as low as possible, the surface polymerization rate should be kept as slow as possible. A relatively homogeneous growth condition of the chains will benefit low polydispersity.

In this study, we find that the properties of the polymer brushes are greatly dependent on the coupling between the initiator density and the polymerization rate. By tuning the initiator density and modifying the polymerization rate, we can obtain the polymer brushes with different degrees of polydispersity. This study partially emphasizes the importance of considering the effects of polymerization rate in further investigations [28].

Crosslinking reaction in epoxy resin systems

In the first example, we construct a CG model system to study the controlling factors in SIP reaction. Here, we will start from an experiment relevant system and use our reaction model to study the crosslinking reaction in the epoxy resin systems.

In common epoxy resin molecules, there are two or more epoxy groups which can react with the active hydrogens in amines of the curing agents [36]. The products with crosslinked network structures exhibit excellent thermal stability, high static modulus, low creep, and high-temperature performance, etc. As a result, these materials are widely applied as electronic packaging, bonding agents, and in aerospace industry as coatings, composites, and adhesives. In practice, what is the network formation mechanism at microscale and how to control the network structure are the main concerns in the applications. Here we propose the CG simulation strategy to reasonably generate chemical network structures by incorporating both particle diffusion and curing reactions in the models. We use this model to solve a typical problem in polymer composite industry, i.e., the formation of epoxy resin with carbon fibers. Normally, commercial carbon fibers are covered by a layer of sizing agent, in most cases composed of epoxy molecules with high molecular weight. The existence of the sizing agent introduces new interfaces between the epoxy resin and carbon fibers, which consequently influences the properties of the produced composites. Thus the mechanisms of the epoxy network formation and the corresponding network structure with the existence of sizing agent are the most important issues to be clarified for manufacturing satisfactory composites.

Suppose the curing reaction is fast as compared to the sluggish diffusion of epoxy molecules. We then can adopt the reaction model together with DPD to study the influences of controlling factors on the epoxy resin network structures. In correspondence with the experiments [37], the curing agent is chosen as 4,4'-diaminodiphenyl sulfone (DDS), and the epoxy resin is diglycidyl ether of bisphenol A (DGEBA). Another epoxy resin is also included in this system as the sizing agent, which has the same repeat unit like DGEBA but a higher repeat number n = 9 (n is the number of repeating units). We first design the mapping project between the atomistic structures and the DPD coarse-grained particles. According to the structural characteristics of the three kinds of molecules, the curing agent DDS is represented by one DPD particle; the DGEBA molecule (abbreviated as RA in the following) is coarse-grained into N = 2 oligomer (N is the number of DPD particles of a molecule); whereas the sizing agent (abbreviated as SA in the following) is coarse-grained into N = 10 chain. As a result, as illustrated in Fig. 6, there are three types of particles, A, B, and C in our coarse-grained model, in which A denotes DDS, B denotes RA, and C denotes SA. Some of the particles are taken as reactive ones which contain amine or epoxy groups, such as the DDS particles and the terminal particles in RA or SA chains. According to the corresponding experimental heat treatment, we determine 423 K as our simulation temperature. In a general correspondence with experiments, we need to construct a 25 nm³ SA layer on top of the unmovable planar substrate to represent the carbon fiber surface. Then another 125 nm³ RA/DDS mixture system on top of the SA layer is constructed in our model. Furthermore, based on the densities of the components and their volume fractions, we obtain the average density of the mixture as 1.08 g/cm³. Since the CG scheme is clear, as shown in Fig. 6, combining the volume fractions of each kind of particle, we get the mean molar mass of the particles, 266.28 g/mol. Then the average volume of the CG particle is v = 409.42 ų. In classical DPD simulations, the reduced particle number density is chosen as {\rho }_n=\mathrm{3}, and we can define the length scale of the simulation as L={(v\cdot {\rho }_n)}^{\mathrm{1}/\mathrm{3}}=\mathrm{1.07}nm. This length scale denotes that there are 3 DPD particles in a small cell of a size L^{\mathrm{3}}. As conventionally used in DPD simulations, the interaction radius of CG beads, r_c, is set equal to L. As that in the surface initiated polymerization model, we construct the carbon fiber surface as a planar substrate by regularly arranging the unmovable particles at the bottom of the simulation box, and it is designed as an inactive carbon fiber without any functional groups on the surface. According to the length scale, we construct a carbon fiber surface with a size \mathrm{24}L\times \mathrm{24}L. The initial configuration is constructed with two separated layers on the surface, in which the upper layer is the mixture of RA and DDS (with the thickness of \mathrm{120}L), and the lower layer is the pure sizing agent (with the thickness of \mathrm{24}L).

We then estimate the parameters in DPD simulations via the calculation of solubility parameters of each component. First, for the system containing k molecules of the same species, we can obtain the non-bonded energy of the model E_{\text{nb}}^k by all atomistic MD simulations. Then each molecule is extracted to vacuum to calculate the non-bonded energy for the individual molecule. Thus the cohesive energy {E_{coh}}_{} can be calculated via [38]

(21)E_{\text{coh}}=\left({\displaystyle \sum _{i=1}^k{E}_{\text{nb}}^{\text{isolated}}\left(i\right)-E_{\text{nb}}^k}\right)/k

where E_{\text{nb}}^{\text{isolated}}\left(i\right) is the non-bonded energy for the i-th isolated molecule in vacuum. The cohesive energy densities of the components at the experimental temperature can therefore be calculated, i.e., e_{\text{coh}}=E_{\text{coh}}/V_m, where V_m is the molar volume of the component. Subsequently, the solubility parameters of the components can be estimated with \delta =\sqrt{e_{\text{coh}}}. Thus, the Flory-Huggins interaction parameters \chi at 423 K can be calculated by [39,40]

(22){\chi }_{ij}=\frac{{\left({\delta }_i-{\delta }_j\right)}^2{v}_{\text{ref}}}{k_{B}T}

where {v_{ref}}_{} is the reference volume, which can be taken as the average size of the CG particles (i.e., 409.42 ų as calculated above). The empirical relation between Flory-Huggins \chi parameters and DPD interaction parameters can be used to evaluate the DPD interaction parameters between different species:

(23)\chi =(\mathrm{0.286}\pm \mathrm{0.002})({\alpha }_{ij}-{\alpha }_{ii})\text{}(\rho =\mathrm{3})

where {\alpha }_{ii}(=\mathrm{25}) and {\alpha }_{ij} are the DPD interaction parameters between the same and different types of species, respectively. The calculated DPD interaction parameters between different components are shown in Table 1. Obviously, the compatibility between the three components is high, and hence we can expect an apparent mixing between them due to mutual diffusion.

With this CG scheme combined with the reaction model, we can study the coupling influence of diffusion and crosslinking reaction on the network structures. We also restrict the reactions between the particles representing the epoxy and the amine groups in order that it describes a real consumption process of the active hydrogen atoms in DDS. The simulations without curing reactions are conducted for a period of time, so that the molecules in these two layers can diffuse into the other layer and the functional groups can mix together. Then the curing reactions are switched on and the formation of the network structure begins. Fig. 7 shows a typical snapshot of the network structure obtained in our simulations, in which it is clear that the crosslinking points are knotting in the network structures, implying an enhanced mechanical property of the resulted material. Due to different mobilities of the molecules at varying temperatures, the composition of the local network structure is different (based on different gradient distribution of the components). It clearly manifests the competition between reaction and diffusion.

After the curing reaction is switched on in the DPD simulations, the epoxy groups consume the reactive amine hydrogen atoms and gradually form the network structure. We use the percentage of ring-opened epoxy groups in the total original epoxy groups to present the conversion of the curing reaction, as shown by the blue symbols in Fig. 8. We find that the percentage of the ring-opened epoxy groups increases quite slowly after the early stage of the curing reaction. The reason might be that, as the curing reaction proceeds, the reactants become insufficient in local regions and most of the components become part of the large network structure, with which the mobilities of the reactive particles are largely reduced and the curing reactions are effectively slowed down. We also use the percentage change of unreacted reactive groups in the total original reactive groups to characterize the curing conversion, as shown by the red symbols in Fig. 8. Its opposite trend generally proves the above observations.

Figure 9 shows the gradient distribution of the normalized density of crosslinking points and the left epoxy groups along the direction perpendicular to the surface after the accomplishment of the curing reactions, with a 4.5 × 10⁵ time steps diffusion period without reactions at the beginning of the simulations. In the regions with higher crosslink densities (z/L_z = 0.5-1), we find better mechanical properties; while in the regions with higher epoxy group densities, we can observe the diffusion frontier of the primary amine groups in DDS (z/L_z≈ 0.1).

This CG model for specific epoxy system with typical reactions shows the possibility of observing microscopic details of this complex system and clarifying the controlling factors for designing high-performance materials.

Other applications

This CG simulation coupled with a chemical reaction model can be extended to the study of other complex fluid systems. To refine the surface initiated polymerization model, we take into account the activation/deactivation process by defining both the relative activation probability and deactivation probability for the dormant and active chain ends. The influence of activation/deactivation process on the dynamics of SIP with high grafting density can then be investigated by this model. We find that the ratio of activation/deactivation probabilities also greatly influences the properties of polymer brushes. By simultaneously modifying the polymerization rate and tuning activation/deactivation process, we can obtain the polymer brushes with low polydispersity index in a relatively short time.

This simulation model can also be used to study how to control molecular weight distribution (MWD) in emulsion polymerization. We apply DPD method coupled with reaction, with both the chain propagation and bi-radical termination kinetics taken into account, to study the emulsion polymerization. We find that the monomer concentration is the key factor of controlling MWD and plays a decisive role in the formation of high molecular weight polymers. Moreover, increasing polymerization rate can slightly narrow the MWD, and it seems to be possible to obtain most appropriate MWD in emulsion polymerization through suitably tuning the surfactant chain length.

The CG simulation coupled with reaction is not limited to polymerization. Actually, it can be extended to the studies of other typical reaction systems. For example, we investigate the phase-separating system coupled with a simple reversible reaction \text{A}\leftrightharpoons \text{B} in a binary immiscible mixture due to critical quench. In this study, we again apply the reaction probability to the control of the species transformation from A(B) to B(A). We find that the system viscosity strongly influences the asymptotic relationship between the excess energy and the reaction rate. The competition between different dynamic factors results in the steady states with characteristic domain sizes.

Conclusion

Here we introduce a reaction model for use in coarse-grained simulations to study the reactions in complex fluid systems. In this model, we employ a key idea of reaction probability to control the whole process of chemical reactions. For a clear illustration, we apply our model to the studies of surface initiated polymerization and the network structure formation process in typical epoxy resin systems. In the SIP model, by choosing suitable initiator density and polymerization probabilities, we can obtain the polymer brushes with different degrees of polydispersity and different chain length distributions. In the study of the network structure formation in epoxy resin system, we first coarse-grain the all-atom structure into a DPD-suitable system, then use our reaction model to study the effects due to coupling between diffusion and crosslinking reaction on the network structures. By this model we construct a chemically reasonable network structure with a high fractional conversion comparable to experiments in the epoxy resin system. The obtained crosslinking structure implies an enhanced mechanical property of the resulted material. Moreover, this model can be extended to the studies of other types of complex fluid systems. We find that this reaction probability idea should be applicable for most polymer systems with reactions when we need to design the reaction processes.

Also we should admit that there are still some key problems to be solved in the application of this reaction model. Although we have strictly derived the reaction rate from the preset reaction probability concept and further presented possible ways to obtain the reaction probability values from experiments, it is in fact not easy work as it refers to an accurate measurement of activation energy and the reaction rate coefficient in experiments. Besides, we believe the relative values of reaction probability in a series of comparable systems should be more meaningful than their absolute values. So how the relative values and accuracy of the reaction probability will affect the whole performance of the reaction model needs further discussion. Based on the reaction probability idea, we expect further reaction models that closely correlate to the experiments and truly describe the reaction processes will be developed in the future.

References

[1] Glotzer, S. C.; Stauffer, D.; Jan, N., Phys. Rev. Lett. 1994, 72, 4109-4112

[2] Glotzer, S.; Di Marzio, E.; Muthukumar, M., Phys. Rev. Lett. 1995, 74, 2034-2037

[3] Kyu, T.; Lee, J. H., Phys. Rev. Lett. 1996, 76, 3746-3749

[4] Smith, G. W.Mol. Cryst. Liq. Cryst. Sci. Technol. Sect. A1994, 241, 77-89

[5] Lísal, M.; Brennan, J. K.; Smith, W. R., J. Chem. Phys. 2006, 125, 164905

[6] Lísal, M.; Brennan, J. K.; Smith, W. R., J. Chem. Phys. 2009, 130, 104902

[7] He, J. P.; Zhang, H. D.; Chen, J. M.; Yang, Y. L., Macromolecules1997, 30, 8010-8018

[8] Farah, K.; Karimi-Varzaneh, H. A.; Müller-Plathe, F.; Böhm, M. C., J. Phys. Chem. B2010, 114, 13656-13666

[9] Perez, M.; Lame, O.; Leonforte, F.; Barrat, J. L., J. Chem. Phys. 2008, 128, 234904

[10] Genzer, J., Macromolecules2006, 39, 7157-7169

[11] Turgman-Cohen, S.; Genzer, J., Macromolecules2010, 43, 9567-9577

[12] Wang, L.; He, X. H.; Chen, Y., J. Chem. Phys.2011, 134, 104901

[13] Berezkin, A. V.; Kudryavtsev, Y. V., Macromolecules2011, 44, 112-121

[14] Lu, W. Q.; Ding, J. D., Macromolecules2006, 39, 7433-7440

[15] Mukherji, D.; Abrams, C. F., Phys. Rev. E Stat. Nonlin. Soft Matter Phys.2009, 79, 061802

[16] Akkermans, R. L. C.; Toxvaerd, S.; Briels, W. J., J. Chem. Phys.1998, 109, 2929

[17] Groot, R. D.; Warren, P. B., J. Chem. Phys.1997, 107, 4423

[18] Español, P.; Warren, P. B., Europhys. Lett.1995., 30, 191-196

[19] Groot, R. D.; Madden, T. J., J. Chem. Phys.1998, 108, 8713

[20] Weeks, J. D.; Chandler, D.; Andersen, H. C., J. Chem. Phys.1971, 54, 5237

[21] Bird, R. B.; Armstrong, R. C.; Hassager, O., Dynamics of Polymeric Liquids; Wiley: New York, 1977; Vols. 1 and 2.

[22] Kremer, K.; Grest, G. S., J. Chem. Phys.1990, 92, 5057

[23] Nikunen, P.; Vattulainen, I.; Karttunen, M., Phys. Rev. E Stat. Nonlin. Soft Matter Phys.2007, 75, 036713

[24] Liu, H.; Xue, Y. H.; Qian, H. J.; Lu, Z. Y.; Sun, C. C., J. Chem. Phys.2008, 129, 024902

[25] Pastorino, C.; Binder, K.; Kreer, T.; Müller, M., J. Chem. Phys.2006, 124, 064902

[26] Soddemann, T.; Dünweg, B.; Kremer, K., Phys. Rev. E Stat. Nonlin. Soft Matter Phys.2003, 63, 046702

[27] Liu, H.; Qian, H. J.; Zhao, Y.; Lu, Z. Y., J. Chem. Phys.2007, 127, 144903

[28] Liu, H.; Li, M.; Lu, Z. Y.; Zhang, Z. G.; Sun, C. C., Macromolecules2009, 42, 2863-2872

[29] Pal, S.; Seidel, C., Macromol. Theory Simul.2006., 15, 668-673

[30] Chen, C. M.; Fwu, Y. A., Phys. Rev. E Stat. Nonlin. Soft Matter Phys.2003, 63, 011506

[31] Goujon, F.; Malfreyt, P.; Tildesley, D. J., ChemPhysChem2004, 5, 457-464

[32] Malfreyt, P.; Tildesley, D. J., Langmuir2000, 16, 4732-4740

[33] Yamamoto, S.; Tsujii, Y.; Fukuda, T., Macromolecules2000, 33, 5995-5998

[34] Yamamoto, S.; Ejaz, M.; Tsujii, T.; Fukuda, T., Macromolecules2000, 33, 5608-5612

[35] Tsujii, Y.; Ohno, K.; Yamamoto, S.; Goto, A.; Fukuda, T., Adv. Polym. Sci.2006., 197, 1-45

[36] Girard-Reydet, E.; Riccardi, C. C.; Sautereau, H.; Pascault, J. P., Macromolecues1995, 28, 7599-7607

[37] Dai, Z. S.; Zhang, B. Y.; Shi, F. H.; Li, M.; Zhang, Z. G.; Gu, Y. Z., Appl. Surf. Sci.2011, 257, 8457-8461

[38] Li, Z. W.; Lu, Z. Y.; Sun, Z. Y.; Li, Z. S.; An, L. J., J. Phys. Chem. B2007, 111, 5934-5940

[39] Fermeglia, M.; Cosoli, P.; Ferrone, M.; Piccarolo, S.; Mensitieri, G.; Pricl, S., Polymer (Guildf.)2006, 47, 5979-5989

[40] Fan, Z. J.; Williams, M. C.; Choi, P., Polymer (Guildf.)2002, 43, 1497-1502

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg