Computational study of imperfect networks using a coarse-grained model Yelena R. Sliozberg1,2,a) and Tanya L. Chantawansri2 1 2

Bowhead Science and Technology, 15163 Dahlgren Rd, King George, Virginia 22485, USA US Army Research Laboratory, Aberdeen Proving Ground, Maryland 21005, USA

(Received 8 October 2013; accepted 4 November 2013; published online 19 November 2013) The structural and mechanical properties of imperfect entangled polymer networks with various fractions of elastically active chains are studied using a generic coarse-grained model. Network topology is analyzed at various degrees of cross-linking and correlated with the mechanical response under uniaxial deformation at various strain rates. We found excellent agreement between results obtained from the structural analysis and from fitting to stress relaxation data. The relaxation tensile modulus at various engineering strains was also calculated as a function of the fraction of active strands. Results indicate that the mechanical and viscoelastic properties of entangled polymer networks are susceptible to variation in the network structure, where defects can affect the mechanical response especially at low strain rates and the relaxation behavior at long times. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4832140] I. INTRODUCTION

Chemically cross-linked polymer networks represent a remarkable class of soft materials, whose mechanical properties can be tailored for a wide range of applications in areas such as the military,1 industry,2–4 and biomedicine.5, 6 The outstanding mechanical properties of polymer networks under macroscopic deformation are due to the entropic elasticity of long chain molecules linked together by permanent junctions called cross-links. In addition to cross-links, high molecular weight polymers can also form topological constraints or entanglements due to topological restrictions that constrain the motion of the polymer chains. Entanglements and cross-links impose qualitatively similar effects, where they both suppress long-range conformational rearrangements of the strands of the polymer network. However, entanglements can slip when the network is strained. To determine the effect of each contribution to the stress, they are often assumed to be independent and additive.7, 8 The time-dependent stress relaxation modulus, G(t), is typically used to describe the viscoelastic properties of entangled polymers and networks. For entangled polymers, G(t) flattens at the melt plateau modulus, G0N , which is linked with the average monomer spacing between entanglements. At longer time G(t) decreases quickly and ultimately vanishes due to complete relaxation of the polymer chains. However, chemically cross-linked polymer networks behave as viscoelastic solids, where after a fast initial decay, the timedependent relaxation modulus G(t) approaches a finite value forming a equilibrium zone for a large range of time. That is, G(t) ≈ Ge , where the equilibrium modulus, Ge , is treated by the theory of rubberlike elasticity. If the polymer chains forming the network are long enough to have entanglements, they can become permanently trapped during the cross-linking a) Author to whom correspondence should be addressed. Electronic mail:

[email protected]

0021-9606/2013/139(19)/194904/8/$30.00

process and Ge will contain a contribution from both crosslinks and entanglements. In this instance, the entanglement contribution will dominate and control the modulus of the polymer network.9 The imperfections in a polymer network can have a large impact on its mechanical and viscoelastic properties. Topological defects in networks include dangling ends and loops, where the former occurs when a chain is connected at only one cross-link within a network. Occurrences of these defects reduce the concentration of elastically active strands and junctions in the network. For a multi-functional polymer network, a junction is considered to be elastically active when three or more of its sites are attached to the network. Network strands connecting the junction are elastically active when both ends are attached to elastically active junctions. Consequently, imperfections in the network structure and the concentration of elastically active strands and junctions are greatly dependent on the network architecture. For a tetrafunctional cross-linker, Gottlieb and Macosko found that only 75% of the strands are elastically active in a 90% cross-linked polydimethylsiloxane (PDMS) network.10 The formation of elastically inactive strands and junctions contributes to a lower modulus compared to a “perfect” network due to possible slippage of entanglements of these incompletely attached to the network structures9 and thus G(t) is expected to have a higher time dependency. These dangling structures vary in complexity from a linear chain to a highly branched star like structure and the entanglements of these dangling imperfections are not permanently trapped. Computational simulations are an excellent tool to study the effect of imperfections on the mechanical and viscoelastic properties of entangled polymer networks. Experimentally, it is difficult to accurately extract the number of active strands in an entangled network due to the simultaneous presence of trapped and slipping entanglements and structural effects such as cross-links. Several computational investigations have been done to elucidate the structure and relaxation mechanism

139, 194904-1

© 2013 AIP Publishing LLC

194904-2

Y. R. Sliozberg and T. L. Chantawansri

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

of end-linked and randomly cross-linked polymer networks including the effect of entanglements and cross-links on the mechanical properties of entangled networks.11–15 They have demonstrated the importance of trapped entanglements in determining the elastic and relaxational properties of networks and have shown that they dominate the relaxation modulus when the strand length between cross-link is large.11–13 It has also been shown that for both randomly cross- and end-linked networks that the entanglement contribution to the stress is in excellent agreement with the predictions of the Rubinstein-Panyukov slip-tube model.14, 15 Recently, Jensen et al. studied the viscoelastic properties of entangled polymer network with balanced stoichiometry and a network containing imperfections by using the slip-link model.16 Such mean-field models, however, involve a priori assumptions regarding the architecture, molecular weight, and concentration of network imperfections. In addition, the model disregards fluctuation of the cross-links and existence of soluble fractions.16 From this study they found that Ge is a linear function of the density of cross-links plus trapped entanglements and that Ge ∝ φs2 Z, where φ s is the fraction of elastically active strands and Z is the average number of entanglements of a polymer chain. However, real networks have more intricate dangling end defects which are not only polydisperse but can vary from linear to highly branched structures that each contribute differently to the mechanical behavior of the material. In this study we use coarse-grained molecular dynamics simulations, where the polymer is modeled using the beadspring model, to study the effect of imperfections in an entangled polymer network on the mechanical and viscoelastic properties. This method is very general because it is based only on chain connectivity and on the excluded volume interaction between monomers. Even so, the model does not require a priori information regarding network imperfections and entanglements making it predictive. We will also demonstrate how simulation results can be mapped to experimental data such as PMDS gels.

and

2 aF ENE 2 r R0 ln 1 − , UF ENE (r) = − 2 R0

(2)

where UFENE is finite extensible nonlinear elastic potential. We have used standard parameter values for R0 = 1.5a and aFENE = 30U0 /a2 . In this work, both U0 and a are set to 1 and all quantities are given in conventional LJ units. III. MAPPING TO EXPERIMENTAL POLYMER NETWORK

To compare with experimental data, it is necessary to define an energy and length scale. Here we demonstrate this mapping for PDMS. The Kuhn length lk was chosen as the unit of length, in agreement with standard convention,19 where lk is a feature of a single polymer chain. It is characterized as the length of an single step of a freely jointed chain with the same mean-square end-to-end distance, R2 = lk L, where L is the contour length of the polymer chain, L = Nb, b is the bond length, and N is the number of monomers. A value of a = 0.72 nm is estimated using the definition lk = C∞ b, where the characteristic ratio C∞ = 1.88 and b = 0.9655a for our model,19, 20 in addition to lk = 1.3 nm for PDMS chains at 413 K.7 To estimate units of pressure, U0 /a3 , U0 can be taken to be ∼400kB , giving U0 /a3 ∼15 MPa. The local structure of the polymer melts is characterized by two independent length scales: the Kuhn length and the packing length, p. The packing length, p = (ρ ch R2 )−1 , indicates the length at which the polymers start to interpenetrate, where chain number density is defined as ρ ch = ρ/N. The ratio of Kuhn and packing length lk /p has be used to compare experimental and simulation systems,19 where the model and experimental PDMS chains have a lk /p values are equal to 2.7 and 3.2, respectively.7 Using a = 0.72 nm the model monomer density of ρ = 0.85a−3 can also be converted to real units, 0.867 g/cm3 , which is similar to the experimental density of PDMS at 413 K, 0.895 g/cm3 .

II. POLYMER MODEL

IV. SIMULATION METHODS

Polymer chains are represented as a sequence of coarsegrained beads connected with springs, where we used the Kremer-Grest model.17, 18 The standard truncated LennardJones pair potential is used to describe the pair interaction between topologically nonconnected particles:

In this section we will discuss simulation details which include the initial equilibration of the entangled polymer chains, network preparation, and deformation of the entangled imperfect network. All simulations were executed using LAMMPS, which is a molecular dynamics program from Sandia National Laboratories.21, 22

ULJ (r) = 4U0

a 12 r

−

a 6 r

−

a rc

12

+

a rc

6 ,

(1)

where U0 is the depth of the potential well, rc represents the cutoff distance, and a is the distance where the interparticle force is zero. A rc = 21/6 a is chosen, producing the Weeks-Chandler-Andersen excluded volume potential, UWCA . The standard FENE/Lennard-Jones bonded potential is used to model the interaction between topologically bound monomers, UFENE/LJ . UF ENE/LJ (r) = UF ENE (r) + UW CA (r)

A. Melt equilibration

The first step of our simulation protocol was to produce a well-equilibrated polymer melt composed of the pre-polymer chains and four-functional cross-linkers (f = 4, where f is the functionality of cross-links). The reactant mixture consists of Mx = 500 cross-linkers modeled as four-arm stars each with an arm length of Na = 1 and Mp = 1000 linear polymer precursor chains with length of Np = 246. Equilibration of entangled polymers is nontrivial even for a case of

194904-3

Y. R. Sliozberg and T. L. Chantawansri

coarse-grained chains because of slow reptation dynamics exhibited by high molecular weight chains. In this paper, we have used our fast equilibration protocol, which has been shown to produce well-equilibrated polymer melts for both unentangled and entangled systems.23 This method involves generation of initial configurations whose structures are as close as possible to equilibrated structures at large-length scales. In addition, the chains are allowed to pass through each other to speed up the polymer dynamic through the use of Dissipative Particle Dynamics (DPD).24, 25 In a DPD simulation, particles interact with each other via a soft pairwise, two-body, short-ranged force, F, which is the sum of a conservative force (FC ), dissipative force (FD ), and random force (FR ): FCij + FD FRij . (3) Fi = ij + j =i

j =i

j =i

The conservative force, FC , can be derived from the excluded volume potential UWCA . FC also contains a contribution from bonded particles (FFENE ). The thermostat in the DPD method is provided by the two remaining forces, FD and FR . The dissipative force slows down the particles by decreasing kinetic energy, and this effect is balanced by the random force due to thermal fluctuations. After this initial equilibration, a standard “push-off” procedure was performed to remove overlap and introduce excluded volume prior to switching to the full Kremer-Grest model. A detailed description of this method can be found in Ref. 23.

B. Network preparation

After the precursor melt was equilibrated by our fast equilibration protocol, the end-linked polymer network was prepared by curing the linear chains and cross-linkers, where the latter is modeled as a four-arm star. End particles in the stars and linear chains react during the MD simulation, where FENE bonds are formed during the simulation when the separation distance between the two different end particles is less than 1.2a. In our simulation, the ends of the pre-polymer only react with two different cross-linkers to exclude loop formation. This was excluded due to the low probability of their formation in long entangled chains. We found that removing this restriction resulted in less than 1% of the precursor chains creating loops, and therefore their creation is negligible. Previous simulation and theoretical work have shown that the fraction of loops decreases with increasing chain length due to a reduction in intramolecular reactions with increasing chain length.26, 27 Recent computer simulations by Gilra et al.27 indicate that the fraction of loops produced for a stoichometrically balanced network consisting of strands of 50 units each is less than 3%. In addition, to increase conversion a truncated LennardJones potential (Eq. (1)) with rc = 2.5a is imposed between the end particles of the stars and linear chains. During the course of the simulation an imperfect network is formed, where the number of particles in the network strands, dangling ends and free chains are equal to Ns = 250, Nd = 248, and Nf = 246, respectively. The curing simulation

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

is performed using a time step of t = 0.01τ LJ . A Langevin thermostat with a damping constant of 1.0τ LJ was used to maintain a temperature of T = 1.0 U0 /kB . Since the topological and rheological entanglement spacing Ne is equal to ∼50 and ∼8628 for a Kremer-Grest flexible chain, the polymer chains which composed the network are entangled. C. Deformation

The relaxation modulus in simple extension is evaluated from a stress-relaxation simulation, which consists of a volume-conserving elongation of the sample followed by a long relaxation. An elongation strain e with −1 is imposed, a constant true strain rate of e˙ = 10−4 τLJ which expands the simulation box in both the positive and negative z-directions. The box dimension in the z direction L ˙ where changes uniaxially with time as L(t) = L0 exp(et), L0 is the initial box size. We assume, incompressibility or constant volume, i.e., extension in z is accompanied by lateral compression since changes in volume caused by extension of the elastomers are insignificant in comparison with shape changes.9 Deformation simulations were performed up to a strain of 2.8. At this strain, the lateral dimensions were significantly greater than the particle diameter. The elongational stress σ in the system is calculated from differences in the normal pressure: 1 (4) − σ = Pzz − (Pxx + Pyy ). 2 After applying an initial step strain (we use engineering strain ε from 0.1 to 1.5), the stress decays very slowly, where a total simulation time of t = 2.3 × 106 τ LJ was used to reach equilibrium. We have also evaluated the stress-strain relationships for our networks, where uniaxial tensile deformation simulations were performed at strain rates ranging from e˙ = 5 × 10−7 −1 − 10−4 τLJ under continuous uniform strain (every time step, t = 0.0075 τ LJ ). To maintain temperature, the deformation simulations are performed with a Langevin thermostat with a damping time constant of 1.0 τ LJ . Every time step, t, the ˙ positions of the monomers were rescaled by a factor Lz Lt where zi is the portion of the position vector in the direction of deformation, and L˙ is the deformation velocity. This rescaling of the cell dimensions, particle coordinates, and the resulting stress allows for a uniform strain to be imposed at each step. V. RESULTS AND DISCUSSION

In this work we considered an imperfect entangled polymer network at various degrees of cross-linking. The results shown are the average of five replicas, unless otherwise stated. A. Network structure

Since imperfections in the polymer network have a large impact on the mechanical properties, it is important to know the fraction of elastically active strands, φ s , for a given fractional conversion, φ c . The fractional conversion is given by φ c = kr /k, where kr and k are the number of reacted ends and

194904-4

Y. R. Sliozberg and T. L. Chantawansri

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

1

simulation model

0.8

φs

0.6 0.4 0.2 0 0.5

0.6

0.7

0.8

0.9

1

φc FIG. 2. The fraction of elastically active strands, φ s , as a function of conversion, φ c . Modeling relations are taken from Ref. 8.

FIG. 1. Schematic representation of an imperfect polymer network. The elastically active (f = 3, 4) and inactive junctions (f = 1) are shown in black and white, respectively. “Semi-active” junctions with f = 2 are shown in red. The elastically active and inactive strands are red and green, respectively. Complex elastically active strands have both ends anchored by elastically active junctions and can be composed of both blue polymer chains and red junctions, where the curvilinear distance is shown as a dashed red line.

the total number of ends of the pre-polymer chains, respectively. A schematic of the polymer network is shown in Figure 1. The connectivity of the polymer networks were analyzed based on the structure searching method proposed by Duering et al.11 For a tetrafunctional polymer network, a cross-link junction is elastically active only when three or more of its sites are independently attached to the network, f = 3 or 4. Network strands are regarded as elastically active when both ends are attached to elastically active junctions. An active strand can be composed not only of a single-polymer chain, but of multiple-polymer chains, which occurs when one or two ends of a chain is connected to a “semi-active” junction with f = 2. These “complex” active strands have both ends anchored by elastically active junctions and the curvilinear distance between these junctions composes the length of a complex active strand. Strands, such as dangling and free chains, are elastically inactive if at least one end is not attached to an active junction. Like active strains, inactive strands can also be composed of multiple-polymer chains. Results of the connectivity analysis of the polymer networks during the curing process are shown in Figure 2. Simulation results were compared with the Miller and Macosko relations,10 where we observed excellent agreement. This finding indicates that the probability of forming each bond during the simulated curing is independent of the probability of creating other bonds. Agreement between theoretical prediction and simulation results suggests that the number of initial configurations and our system size are adequate for this calculation. To study the mechanical and viscoelastic properties of imperfect networks, we have considered polymer networks with various fractional conversion values. The initial configurations were taking from the curing simulations at φ c = 0.7, 0.77, 0.83, 0.89, and 0.95, which correspond to φ s = 0.12,

0.30, 0.50, 0.70, and 0.87. For simplicity, we will round up our notation to φ s = 0.1, 0.3, 0.5, 0.7, and 0.9 for the remainder of the paper. Fractions of active junctions ϕ f with functionality f = 3 and 4 for a given φ c for these networks also shown in Table I. When φ c > pc , where pc is the critical extent of reaction, a network is known to behave as a viscoelastic solid. To determine which of our systems will behave as a viscoelastic solid pc is evaluated from10, 29 rpc2 =

1 , (f − 1) (g − 1)

(5)

where r = 1 is the stoichiometric ratio of cross-linker to prepolymer and g is the functionality of the pre-polymer. For our networks f = 4 and g = 2 which gives pc = 0.58. Since all of the networks considered satisfy φ c > pc , they are expected to behave as viscoelastic solids. 1. Networks conformation prior to deformation

To ensure that our systems are well-equilibrated, we compare structural properties of our network with a reference structure composed of a pre-polymer melt with N = 500. The reference structure was obtained by long brute-force equilibration23 carried out for significantly long time such that the mean-square displacement of the inner monomers scales linearly with time. To characterize the conformational change of the linear chains for various degrees of crosslinking we have computed their mean square internal distances, R2 (n)/n, where all of the network strands, dangling ends, and free chains were considered in this calculation. TABLE I. Fractional conversion, φ c , fraction of elastically active strands, φ s , fraction of active junctions, ϕ f , where f is equal to 3 or 4 and total fraction of active junctions, ϕ a . φc 0.7 0.77 0.83 0.89 0.95

φs

ϕ3

ϕ4

ϕa

0.12 0.30 0.50 0.70 0.87

0.14 0.29 0.40 0.43 0.32

0.01 0.08 0.20 0.38 0.61

0.15 0.37 0.60 0.81 0.93

Y. R. Sliozberg and T. L. Chantawansri

2

2

0.25

φs=0.1 φs=0.3 φs=0.5 φs=0.7 φs=0.9

1.6

φs=0.1 φs=0.3 φs=0.5 φs=0.7 φs=0.9

0.2 3

1.8

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

σ, U0 /a

194904-5

0.15 0.1

1.4

0.05

1.2

0

(a) 0

0.5

1

1.5

2

2.5

3

ε

1 0.25 100

φs=0.1 φs=0.3 φs=0.5 φs=0.7 φs=0.9

1000

n FIG. 3. Mean square internal distances R2 (n)/n for various degrees of cross-linking, φ s .R2 (n)/n for the reference melt of N = 500 equilibrated for t = 6 × 106 τ LJ shown in black dots for comparison.23

0.2 3

10

σ, U0 /a

1

0.15 0.1 0.05

(b)

0

2. Mechanical properties

We performed a series of simulations of uniaxial tensile deformation at strain rates ranging from e˙ = 5 −1 . We present our results in Figure 4 using × 10−7 − 10−4 τLJ the dimensionless strain rate or the Weissenberg number. The ˙ d , where τ d is the Weissenberg number is defined as Wi = eτ relaxation time of a linear unconstrained polymer chain. Since the connection between different time scales of the simulation and experiment cannot be obtained directly due to the faster dynamics observed in the coarse-grained systems, combination of the elongation rate and the length-scale by using Wi is useful for comparison between simulation and experimental results.31 We estimated a value of τ d from τ d = 0.39N2 (1 + N/Ne ) τ LJ ,32 where Ne = 8628 and found that τ d ≈ 105 τ LJ . Figure 4 shows the stress-strain relationships for three values of Wi : 10, 1.0, and 0.05. The stress-strain curves show the characteristic rubbery behavior with a very low elastic modulus. When the strain rate is greater than the relaxation time of unconstrained chain Wi > 1, untrapped entanglements are unable to relax, and consequently the entangled pendant and free chains significantly contribute to the stress. In addition, both active and inactive network strands will contribute equally to the stress since for a long polymer chain the entanglements contribution to the stress is much larger than the cross-links contribution. This is especially evident at Wi = 10 where our results clearly show that the stress is essentially equal for all degrees of cross-linking (Figure 4(a)). For Wi = 1, entanglements in the network “imperfections”

0 0.25

0.5

1

1.5 ε

2

2.5

3

φs=0.1 φs=0.3 φs=0.5 φs=0.7 φs=0.9

3

0.2 σ, U0 /a

R2 (n)/n is averaged over all possible combinations of segments of size n = |i − j| along the chains, where i < j ∈ [1, N] are monomer indices.30 This metric is an excellent indicator of chain configuration at all length scales. Figure 3 shows that the structural properties of the networks obtained using our curing simulation protocol are very close to the reference melt calculations.23 This means that the network strands have the same conformation as chains not constrained by cross-links and that our networks are not deformed prior to the deformation simulation.

0.15 0.1 0.05

(c)

0 0

0.5

1

1.5

2

2.5

3

FIG. 4. True stress, σ , versus engineering strain, ε, under uniaxial tension deformation for polymer networks with various fraction of elastically active chains, φ s , for Weissenberg number, Wi , equal to (a) 10, (b) 1.0, and (c) 0.05.

start to relax and as a result, the stress becomes a function of the network structure (Figure 4(b)). Note that relaxation time of these complex branched “imperfections” are longer than τ d of a linear unconstrained polymer chain and consequently some entanglements of incompletely attached structures are still not relaxed at Wi = 1. At Wi 1, these entanglements of incompletely attached network structures are relaxed, and consequently the stress strongly increases with increasing φ s (Figure 4(c)). Our results indicate that the stress response for a polymer network at a specific φ s drops with decreasing strain rate (Figure 5). A faster decrease in σ is observed at higher strain rates, while it becomes slower at lower strain rates. However, a polymer network of φ s = 0.1, which contains a large fraction of incompletely attached structures, exhibits a stronger rate dependent mechanical response than a near perfect polymer network with φ s = 0.9 (Figure 5). 3. The relaxation modulus

We have measured the relaxation tensile modulus, E(t), at an engineering strain, ε, equal to 0.1, 0.3, 0.5, 1.0, and 1.5. Values of ε were chosen to be below the onset of strain hardening and corresponds to the nonlinear regime with the exception of ε = 0.1, which is expected to produce a near linear relaxation response. E(t) is plotted versus time for a

Y. R. Sliozberg and T. L. Chantawansri

0.25

0.15

(a)

3

3

1⋅10-1

-4

1*10-5/τLJ 5*10 /τLJ -5 1*10-6/τLJ 5*10-6/τLJ 1*10-7/τLJ 5*10 /τLJ

0.2 σ, U0 /a

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

E (t), U0 /a

194904-6

0.1

ε=0.1 ε=0.3 ε=0.5 ε=1.0 ε=1.5

0.05

(a) 0 0 0.25

1.5 ε

2

2.5

3 1⋅10

1⋅10 1⋅10-1

3

4

1⋅10 t, τLJ

1⋅10

1⋅104

1⋅105 t, τLJ

1⋅106

1⋅10

5

6

1⋅10

7

(b)

3

0.15

-2

0.1

E (t), U0 /a

3

1

1*10-4 /τLJ 5*10-5 /τ -5 LJ 1*10-6/τLJ 5*10 /τLJ -6 1*10-7/τLJ 5*10 /τLJ

0.2 σ, U0 /a

0.5

0.05

(b)

ε=0.1 ε=0.3 ε=0.5 ε=1.0 ε=1.5

0 0

0.5

1

1.5 ε

2

2.5

3

1⋅10

FIG. 5. True stress, σ , versus engineering strain, ε, for the uniaxial tension of polymer networks at various strain rates for fraction of the elastically active chains (a) φ s = 0.1 and (b) φ s = 0.9.

1⋅103

1⋅107

FIG. 6. Relaxation tensile modulus E(t) plotted against time for networks with a fraction of active strands (a) φ s = 0.1 and (b) φ s = 0.9.

time to approach the plateau modulus. However we found that E(t) for this network is almost independent of time (E(t) ∼ t−0.038 ) during the prescribed time range (2.25 × 106 − 4.5 × 106 τ LJ ) and thus we can assume that E(t) reaches its equilibrium value. Similar results are observed for all φ s considered, confirming that our networks are viscoelastic solids. In this study, Ee is taken as the asymptotic value of E(t) at long time (t > 2.25 × 106 τ LJ ). Since the relaxation time of the dangling ends and unreacted pre-polymer chains, also called sol, are less than the -1

3

1⋅10

E (t), U0 /a

polymer network of φ s equal to 0.1 and 0.9 in Figure 6, where due to the finite system size of our simulations, the lowest strain produces the largest noise in E(t). E(t) is evaluated from a stress-relaxation experiment (step-strain) described previously in Sec. IV C. Results indicate that E(t) of a network composed of 10% of elastic strands flattens somewhat and then continues to decrease (Figure 6(a)). This behavior can be attributed to the relaxation of entanglements on branched structures which are not completely attached to the network. However, the network composed of 90% of elastic strands behaves more like a solid where after a fast initial decay, E(t) approaches a finite value called the long time or equilibrium modulus, Ee , forming a equilibrium zone for a large range of time (Figure 6(b)). The Ee in polymer networks is very sensitive to changes in the network structure, and Figure 7 shows that the presence of pendant structures and unreacted chains in the network lead to a decrease of the equilibrium modulus. However, at short times, values of E(t) are similar for all φ s considered. For networks with a higher amount of the active chains such as φ s = 0.7 and φ s = 0.9, we observe faster convergence of E(t) to Ee . Lower values of φ s contain a larger quantity of soluble and complex pendant structures, which require longer times to completely relax. Although all of the networks considered satisfy the relation φ c > 0.58 (see Sec. V A) and thus are expected to behave as viscoelastic solids, a more accurate method to determine the critical gel concentration is the time dependency of the relaxation modulus. The network with the lowest value of active chains, φ s = 0.1, requires the longest

-2

φs=0.1 φs=0.3 φs=0.5 φs=0.7 φs=0.9

1⋅10-2 1⋅103

1⋅104

1⋅105 t, τLJ

1⋅106

1⋅107

FIG. 7. Relaxation tensile modulus, E(t), plotted against time for network at strain ε = 1.0 for various fractions of active strands φ s.

Y. R. Sliozberg and T. L. Chantawansri

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

simulation time, we assumed that we have reached equilibrium and could evaluate the equilibrium tensile stress Ee and σ . Since strands of the simulated networks are long enough to have entanglements, the network shear modulus Ge contains both cross-links and entanglement contribution and thus Ge = Gc + G0N , where Gc and G0N are contributions from cross-links and entanglements, respectively. The dangling ends and sol do not contribute to the equilibrium modulus because they cannot store elastic energy at equilibrium. Thus the cross-linked contribution, Gc , is given by the phantom model,7 Gc = (ν − μ)kB T ,

where μf is number density of the junctions with functionality f = 3 and 4. The entanglement contribution is given by7 G0N = Te ρe kB T ,

(8)

where Te is trapping factor defined as a fraction of the entanglements permanently trapped in the network and ρ e is the entanglement density of the perfect network. ρ e is comM N puted from ρe = V pNe , where V is the volume of the system and Ne is taken as the topological entanglement length which is equal to 50.28 We have fitted our simulated stress relaxation data of the network with φ s = 0.9 to two models: classical stress-elongation and slip-tube model.7, 8 The classical stresselongation model predict variation of the true stress with uniaxial elongation as 1 , (9) σ = Ge λ2 − λ where Ge = Gc + G0N and λ is the elongation ratio. We evaluated Gc from Eq. (7) using μ3 and μ4 evaluated from the structural analysis of the networks and G0N from Eq. (8) using

Slip-Tube Model Classical Model

0.14

3

0.12 σ U0/a

φs Te

0.8 0.6 0.4 0.2 0 0.5

0.1 0.08 0.06

0.6

0.7

0.8

0.9

1

φc

(6)

where ν and μ are the number density of the elastically active strands and cross-links, respectively. Since there are f/2 networks strands per cross-link, μ = 2ν/f and Gc could be rewritten as f Gc = − 1 μf kB T , (7) 2 f =3,4

0.16

1

φs, Te

194904-7

FIG. 9. The fraction of elastically active strands, φ s , and trapping factor, Te , as a function of conversion φ c . The symbols represent the results obtained by fitting the stress relaxation data to Eqs. (8) and (10). The line represent φ s found from the connectivity analysis.

Te = ϕ s = 0.9. Our results indicate that the classical theory is applicable only to very small deformations and since entanglements slip during deformation, the model clearly deviates from the simulation results (Figure 8). To describe larger deformation we have used non-classical stress-strain relation proposed by Rubinstein and Panyukov in their non-affine slip-tube model:8 G0N 1 2 λ . (10) σ = Gc + − 0.74λ + 0.61λ−1/2 − 0.35 λ Our results demonstrate that the slip-tube model accurately predicts the simulation results and can be used to capture linear viscoelastic properties at low strain. This is advantageous since the viscoelastic properties obtained directly from a particle-based simulation at low strain is generally not very accurate due to large noise. We have also used the slip-tube model to compute the trapping factor, Te for φ s = 0.1, 0.3, 0.5, 0.7, and 0.9 at strain ε = 0.1. For this purpose, we computed Gc from Eq. (7) and found G0N by fitting our relaxation data to Eq. (10). Then the trapping factor Te was extracted from Eq. (8). Our results demonstrate that Te ≈ φ s for networks with φ c > 0.77 (Figure 9). Networks with φ c ≤ 0.77 have a slightly higher amount of effectively trapped entanglements than entanglements made up of elastically active strands, Te > φ s . This discrepancy stems from the fact that complex branched pendant structures exhibited by these networks have a very long relaxation time and their entanglements are still not fully relaxed during the long simulation runs (t = 4.5 × 106 τ LJ ). Even so, we observe excellent agreement between results obtained from structural analysis and from fitting to stress relaxation data.

0.04 0.02

VI. CONCLUSIONS

0 0

0.5

1

1.5

2

ε

FIG. 8. Elongational true stress, σ , as a function of strain for network with φ s = 0.9. The symbols represent results from our simulations and the lines represent fits to the classical stress-elongation (blue) and slip-tube model (green).

Our results demonstrate that viscoelastic and mechanical properties of the entangled polymer networks are very sensitive to changes in the network structure. While at short times the relaxation modulus is almost independent of the defects in the network connectivity, the modulus at long times increases with increasing fraction of active chains. The

194904-8

Y. R. Sliozberg and T. L. Chantawansri

stress-strain relationships follow a similar trend. When the strain rate is faster than the relaxation time of unconstrained chains the elongational stress of the entangled networks is essentially the same for all degrees of cross-linking. However, at lower strain rates, the stress strongly increases with the degree of cross-linking due to relaxation of the entanglements in the incompletely attached network structures. Our results demonstrate that the slip-tube model accurately predicts the stress relaxation for imperfect networks. Results of this study can be further used to study the stress optical behavior of polymer networks under uniaxial elongational flow. ACKNOWLEDGMENTS

The authors will like to thank Dr. Joseph Lenhart, Dr. Randy Mrozek, Dr. Jan Andzelm, and Dr. Timothy Sirk for useful discussion. Calculations were performed on Department of Defense High Performance Computing site at the Air Force Research Laboratory through the use of the supercomputer allocation “Challenge Project C5M.” 1 J.

L. Lenhart, P. J. Cole, B. Unal, and R. C. Hedden, Appl. Phys. Lett. 91, 061929 (2007). 2 A. Moradi-Araghi, J. Pet. Sci. Eng. 26, 1 (2000). 3 S. Ozkan, T. W. Gillece, L. Senak, and D. J. Moore, Int. J. Cosmet. Sci. 34, 193 (2012). 4 M. Rutkevicius, S. K. Munusami, Z. Watson, A. D. Field, M. Salt, S. D. Stoyanov, J. Petkov, G. H. Mehl, and V. N. Paunov, Mater. Res. Bull. 47, 980 (2012). 5 H. Xu, J. Wu, C. C. Chu, and M. L. Shuler, Microdevices 14, 409 (2012). 6 I. E. Palama, S. D’Amone, A. M. L. Coluccia, M. Biasiuccia, and G. Gigli, Integr. Biol. 4, 228 (2012). 7 M. Rubinstein and R. H. Colby, Polymer Physics (Oxford University Press, Oxford, 2003). 8 M. Rubinstein and S. Panyukov, Macromolecules 35, 6670 (2002).

J. Chem. Phys. 139, 194904 (2013) 9 J.

D. Ferry, Viscoelastic Properties of Polymers (Wiley, New York, 1980). Gottlieb, C. W. Macosko, G. S. Benjamin, K. O. Meyers, and E. W. Merill, Macromolecules 14, 1039 (1981). 11 E. R. Duering, K. Kremer, and G. S. Grest, Prog. Colloid Polym. Sci. 90, 13 (1992). 12 E. R. Duering, K. Kremer, and G. S. Grest, Macromolecules 26, 3241 (1993). 13 E. R. Duering, K. Kremer, and G. S. Grest, J. Chem. Phys. 101, 8169 (1994). 14 G. S. Grest, M. Pütz, R. Everaers, and K. Kremer, J. Non-cryst. Solids 274, 139 (2000). 15 C. Svaneborg, R. Everaers, G. S. Grest, and J. G. Curro, Macromolecules 41, 4920 (2008). 16 M. K. Jensen, R. N. Khaliullin, and J. D. Schieber, Rheol. Acta 51, 21 (2012). 17 K. Kremer and G. S. Grest, J. Phys. Condens. Matter 2, SA295 (1990). 18 K. Kremer and G. S. Grest, J. Chem. Phys. 92, 5057 (1990). 19 K. Kremer, S. K. Sukumaran, R. Everaers, and G. S. Grest, Comput. Phys. Commun. 169, 75 (2005). 20 M. Kröger, W. Loose, and S. Hess, J. Rheol. 37, 1057 (1993). 21 S. Plimpton, J. Comp. Phys. 117, 1 (1995). 22 Sandia Laboratories Home Page. LAMMPS Molecular Dynamics Simulator. http://lammps.sandia.gov (accessed on 11/13/2013). 23 Y. R. Sliozberg and J. W. Andzelm, Chem. Phys. Lett. 523, 139 (2012). 24 P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett. 19, 155 (1992). 25 J. M. V. A. Koelman and P. J. Hoogerbrugge, Europhys. Lett. 21, 363 (1993). 26 S. Dutton, R. F. T. Stepto, and D. J. R. Taylor, Angew. Makromol. Chem. 240, 39 (1996). 27 N. Gilra, C. Cohen, and A. Z. Panagiotopoulos, J. Chem. Phys. 112, 6910 (2000). 28 R. S. Hoy, K. Foteinopoulou, and M. Kröger, Phys. Rev. E 80, 031803 (2009). 29 S. K. Venkataraman, L. Coyne, F. Chambon, M. Gottlieb, and H. H. Winter, Polymer 30, 2222 (1989). 30 R. Auhl, R. Everaers, G. S. Grest, K. Kremer, and S. Plimpton, J. Chem. Phys. 119, 12718 (2003). 31 Y. R. Sliozberg, R. A. Mrozek, J. D. Schieber, M. Kröger, J. L. Lenhart, and J. W. Andzelm, Polymer 54, 2555 (2013). 32 M. Kröger and S. Hess, Phys. Rev. Lett. 85, 1128 (2000). 10 M.

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