Molecular dynamics simulation of O2 diffusion in polydimethylsiloxane

Molecular Simulation, 2013 http://dx.doi.org/10.1080/08927022.2013.830183

Molecular dynamics simulation of O2 diffusion in polydimethylsiloxane (PDMS) and end-linked PDMS networks Varun Ullala and Douglas E. Spearota,b*

Downloaded by [University of Arkansas Libraries - Fayetteville], [Douglas Spearot] at 12:55 05 September 2013

a

Department of Mechanical Engineering, University of Arkansas, Fayetteville, AR 72701, USA; bInstitute for Nanoscale Materials Science and Engineering, University of Arkansas, Fayetteville, AR 72701, USA (Received 19 April 2013; final version received 25 July 2013) Molecular dynamics simulations are used to compute diffusion coefficients for O2 molecules in polydimethylsiloxane (PDMS) and end-linked PDMS networks. The PDMS chains and penetrants are modelled using a hybrid interatomic potential which treats the Si and O atoms along the chain backbone explicitly while coarse-graining the methyl side groups and penetrants. In PDMS models with different molecular weights, diffusivity of the O2 penetrants is found to modestly decrease with an increase in chain length. To match typical experimental conditions, the end-linked PDMS networks are constructed with a PDMS to crosslinking (CL) molecule mass ratio of 5:1 or 10:1, demanding that the number of CL molecules exceeds the number of PDMS chains in each model. Despite end-linking, the presence of non-bonded CL molecules promotes increased O2 diffusivity in comparison with uncrosslinked PDMS. Temperature dependence is captured using the Williams –Landel – Ferry equation. Keywords: molecular dynamics; diffusion; polydimethylsiloxane; crosslinking

1.

Introduction

Polydimethylsiloxane (PDMS) is a silicon-based organic polymer belonging to the family of siloxanes with the chemical formula CH3[Si(CH3)2O]nSi(CH3)3. Its chemical inertness and non-toxic properties make it an ideal material for microelectromechanical systems (MEMS) applications, such as microfluidic pumps,[1] gas separation devices [2] and microcorrosion sensors.[3] In these applications, PDMS is cross-linked through chemical modification of the chain ends followed by mixing with a crosslinking (CL) agent capable of binding the PDMS chain ends together. For example, Sylgardw 184 silicone elastomer from Dow Corning includes vinyl-terminated PDMS chains and is mixed with a CL agent in a 10:1 (PDMS:CL molecule) mass ratio until cured.[4] CL transforms PDMS from a viscous gel to a solid elastomeric material through the formation of an end-linked network, with elastic properties that are dependent on the cure conditions. In gas separation [2] and sensing [3] applications, it is critical to understand the diffusion of small molecules through PDMS, including diffusion rates and mechanisms, as these nanoscale details fundamentally impact device performance. Over the past 20 years, the mechanisms by which atmospheric gases diffuse through polymers have been studied using molecular dynamics (MD) simulation. [5-13] One of the first works to study diffusion through polymers using MD simulation was carried out by Trohalaki et al. [5] They simulated the diffusion of 4 CO2 penetrants in a system of 25 polyethylene chains, each containing 20 united atom CH2 units. The computed

*Corresponding author. Email: [email protected] q 2013 Taylor & Francis

diffusion coefficients were larger than experiments; this was attributed to the presence of crystallites in the real polymer which were not considered in the simulation. Sok et al. [8] analysed the effect of penetrate size on the diffusion of atmospheric gases through PDMS. Their system consisted of five chains each with 30 PDMS monomer units and was modelled using the GROMOS interatomic potential. Calculated diffusion constants for CH4 and He penetrants were in good agreement with experiment. Charati and Stern [10] studied the influence of side groups on the diffusion of small molecular gases (He, O2, N2, CO2 and CH4) at 300 K through five different silicone polymers. Their simulations concurred with the experimentally observed trend of a decrease in the diffusion coefficient with an increase in size of the polymer side chains and the penetrant molecules. A similar MD study involving varying penetrant size was carried out by Jawalkar and Aminabhavi [11] using four different penetrants (N2, O2, CO2 and CH4). One of the significant advantages of molecular simulation is the ability to purposely manipulate the atomistic system to isolate the role of specific polymer mechanisms on diffusion. For example, Takeuchi and Okazaki [6] studied diffusion of O2 in two polyethylene models: one model included the full interatomic potential incorporating bond stretch, angle bending, dihedral torsion and non-bonded effects, while the second model was identical except that the dihedral potential was turned off. Elimination of the torsional potential increased the flexibility of the polyethylene chains and consequently the distribution and dynamics of the free volume in the


2

V. Ullal and D.E. Spearot

system. Upon elimination of the rotational barrier, O2 diffused five times faster.[6] Further investigation of the influence of the torsional barrier on diffusion was carried out by Boyd and Krishna Pant [7] who showed that raising the torsional barrier completely halted the diffusion of O2 in polyethylene. These simulations also imply that significant error in the computed diffusion coefficient can be realised if artificial system constructions are used; significant care is taken in this study to accurately build both uncrosslinked and crosslinked PDMS models. The above studies provide a very clear understanding of the mechanisms associated with small molecule diffusion through polymers and the role of polymer attributes (side group size, torsional flexibility, etc.) on diffusion. However, these previous diffusion studies consider only uncrosslinked homogeneous polymers. In typical MEMS applications, PDMS is mixed with a CL molecule to create an end-linked network. Furthermore, the mixing is carried out with a specified PDMS:CL molecule mass ratio [4] to ensure that a sufficient number of PDMS chain ends are reacted, meaning that many unreacted CL molecules will be present within the cured PDMS elastomer. Therefore, the objective of this study is to study O2 diffusion (as a representative atmospheric penetrant) through end-linked PDMS with experimentally relevant chain lengths and mixing ratios. Both CL [14 – 22] and diffusion have been studied independently using MD simulations; the unique aspect of this study is that the combined influence of CL and mixing ratio on small molecule diffusion is explored. Following Heine et al. [17], this study uses silanolterminated PDMS chains and the silane-terminated CL molecule tetrakis(dimethylsiloxy)silane, as shown in Figure 1. The end-linking reaction occurs by forming a bond between a terminal SiH group on the CL molecule to the nearest OH end of the PDMS chain, both of which are modelled using the united-atom approximation. Note that crosslinked PDMS is often prepared experimentally using vinyl-terminated PDMS chains.[4] The chemical details of the end-linking procedure are not expected to appreciably

influence the diffusion of O2 molecules through the PDMS network because the CL molecules are of comparable size in both cases.

2.

Atomistic simulation of PDMS systems

2.1 Interatomic potential MD simulation is a computational technique by which the time evolution of an interacting set of atoms is determined via the numerical solution of a prescribed set of equations of motion.[23] When modelling complex molecular systems, a simplifying approach is used which involves the grouping of multiple atoms into a single point mass, thereby reducing the computational cost of the simulation. This modelling technique is known as the united-atom model or coarse-grained approach.[24] In this study, the PDMS chains are modelled using a hybrid interatomic potential developed by Frischknecht and Curro [25] which treats the Si and O atoms along the chain backbone explicitly while coarse-graining the methyl side groups. Bond and angle potentials are specified using harmonic forms while dihedral torsion interactions are modelled using a periodic function, U bond ðbÞ ¼ kb ðb 2 b0 Þ2 ;

ð1Þ

U angle ðuÞ ¼ ka ðu 2 u0 Þ2 ;

ð2Þ

U torsion ðfÞ ¼ kd ½1 þ cosðnfÞ;

ð3Þ

where kb , ka and kd are the stiffness constants for bond, angle and dihedral interactions, respectively; b, u and f are the instantaneous values of bond lengths, bond angles and dihedral angles, respectively; b0 and u0 are the equilibrium bond length and bond angles, respectively, and n controls the period of the dihedral torsion function. Non-bonded interactions are modelled using the following equations: sab 12 sab 6 U LJ ðrÞ ¼ 41 2 ab vdw r r ð4Þ when a or b ¼ CH3 ; U IIvdw ðrÞ

sab 9 sab 6 ¼ 1ab 2 23 r r

ð5Þ

when a and b ¼ Si or O; U coulombic ðrÞ ¼ Figure 1. (Colour online) (a) Silanol-terminated PDMS chain and (b) silane-terminated CL molecule. Silicon atoms are in orange, oxygen atoms are in red and methyl united atoms are in white. The ends of the PDMS chain are modelled using an OH united atom, while the ends of the CL molecule are modelled using a SiH united atom.

qa qb ; 4 p 10 r

ð6Þ

where 1ab and sab are the Lennard– Jones (LJ) parameters, r is the separation distance between a pair of atoms, qa and qb are partial charges of the atoms and 10 is the permittivity of free space. Non-bonded interactions occur between atoms in different chains and between atoms further apart


Molecular Simulation than third nearest neighbours within the same chain. Nonbonded interactions within the siloxane backbone use a softer LJ 9-6 form (Equation (5)) to better capture the flexibility of the PDMS chain while all other non-bonded pair interactions use a standard LJ 12-6 form (Equation (4)).[25] The O2 penetrants in this study are also modelled using a coarsegrained approach with united-atom parameters extracted from Sok et al. [8]. Cross-interactions between atomic species are computed through standard mixing rules.[23]

2.2 PDMS model construction The systems modelled in this study contain 50 PDMS chains with molecular weights of approximately 5000, 10,000, 20,000 or 40,000 g/mol (the entanglement molecular weight of PDMS is 24,500 g/mol [26]). In uncrosslinked PDMS models, the chains are terminated with three methyl groups while in the case of end-linked PDMS models, the chains are terminated by one OH united atom and two methyl groups.[17] The number of CL molecules in each system is determined from the mass ratio of the components, as shown in Table 1. PDMS chains and CL molecules are created using a configurational bias Monte Carlo (CBMC) approach (using the code Towhee [27]) in an constant number of atoms, volume and temperature (NVT) ensemble at a temperature of 300 K in a large simulation cell corresponding to a density of less than 10% of the experimental density of PDMS.[26] Three independent models of each PDMS system are created with different initial random chain conformations. The CBMC simulations are run for 250 steps to randomise each system. Although typical Monte Carlo simulations employ a much larger number of steps, the purpose of using CBMC in this study is simply to generate an accurate initial structure for each PDMS chain and CL molecule with respect to the equilibrium bond, angle and dihedral configurations specified in the interatomic potential. After the CBMC chain construction is complete, each atomic configuration is imported into the MD code LAMMPS [28] and is subjected to a thermodynamic equilibration in the NVT ensemble at 300 K for 1 ns with non-bonded interactions (Equations (4) – (6)) turned off. This step allows the PDMS chains and CL molecules (if present) to mix thoroughly in the underdense simulation Table 1. Number of CL molecules required to meet a given PDMS:CL mass ratio.

Pure PDMS 10:1 mass ratio 5:1 mass ratio

5000 g/mol

10,000 g/mol

20,000 g/mol

40,000 g/mol

0 75 150

0 151 303

0 303 606

0 607 1215

3

cell. Next, an equilibration in the NVT ensemble at 300 K for 10 ps is carried out where a ‘soft’ purely repulsive potential is switched on for the non-bonded interactions to push overlapping atoms and chains away from each other. [17] Finally, the appropriate 12-6, 9-6 and Columbic nonbonded interactions (Equations (4) –(6)) are turned on and each model is equilibrated in the constant number of atoms, pressure and temperature (NPT) ensemble at 300 K and 1 atm for 1 ns. During this step, the simulation cell naturally contracts, attaining the equilibrium density for each PDMS system.

3.

Formation of crosslink bonds in PDMS

3.1 CL algorithm After the PDMS models are equilibrated, CL is carried out to create the end-linked network. Several approaches have been proposed in the literature to incorporate crosslink bond formation into atomistic simulations.[14 –22] A common conclusion of these studies is that the allotment of a relaxation time between crosslink bond formations is imperative to avoid the creation of large internal strains which promote unnatural bond, angle and dihedral distributions in the final crosslinked configuration. Most previous studies focused on crosslink formation in epoxy. For example, Varshney et al. [21]. conducted an MD study on epoxy networks with an emphasis on structural properties and volume shrinkage as a function of degree of curing. Comparisons between the total energy of the resulting crosslinked topology and computational time were made between four different CL approaches. The generation of a relaxed topology is ensured through a multistep relaxation process during the CL procedure. Li and Stratchan [22] conducted MD studies on epoxy to understand the role of the CL procedure on the internal strain of the resulting topology and the role of chemical kinetics on the properties of epoxy. Both studies [21,22] accounted for the change in charge distribution in the local vicinity of crosslinks. PDMS is different from epoxy in that the crosslink bonds form at the ends of the chemically modified PDMS chains. Heine et al. [17] conducted MD simulations of endlinked network formation in PDMS using the united-atom framework. Their study used a dynamic CL approach based on a cut-off distance criterion. The relaxation of the crosslinked topology was achieved via the use of a modified bonding potential which was harmonic at distances below a specified cut-off and linear at larger distances, U xlink bond ðrÞ ¼ 81 2 < 2 kbond ðr 2 r 0 Þ ; : 1 kbond ðr cut 2

r # r cut

2 r0 Þ þ 2

1 2 kbond ðr

2 r cut Þ;

r . r cut ;

ð7Þ


4


Figure 3. (Colour online) Distribution of the OZSiZO angle after CL showing the influence of bonding potential and thermodynamic relaxation between crosslink formation steps.

final crosslinked system is equilibrated for an additional period of 1 ns in the NPT ensemble. This CL algorithm is summarised in Figure 2.

3.2 Figure 2. Schematic of the CL algorithm implemented in this study (modified from Varshney et al. [21]).

where the spring constant kbond is set to 35.01 kcal/mol ˚ .[17] Using this approach, structural and r cut is 1.71 A inhomogeneities resulting from stretched bonds, which could adversely affect the natural bond, angle and dihedral distributions in the system, were mitigated. In this study, a CL approach similar to the one used by Heine et al. [17] is adopted. A bond between one arm of a CL molecule and one end of a silanol-terminated PDMS chain (Figure 1) is formed with unit probability if the appropriate atoms are within a chosen reaction distance. ˚ and Equation (7) is This distance is specified to be 6.5 A used to describe bonded interactions.[17] The chosen reaction distance is larger than the sum of the radii of the SiH and OH united atoms promoting the CL reaction to occur. This procedure is carried out in the NPT ensemble at a pressure of 1 atm and a temperature of 300 K. Every 10 ps during the simulation, all potential reactive pairs are identified and crosslinks are generated between eligible pairs. The topological information is updated and the system is relaxed for 10 ps alleviating unnatural changes in the spatial coordinates of atoms in the local region of the newly formed bonds. This is repeated until the system achieves an 80% crosslink density (since there are 50 PDMS chains in the models, this corresponds to the formation of 80 bonds at the chain ends) upon which the

Algorithm validation

To illustrate the importance of relaxation time, a sample system consisting of 50 short silanol-terminated PDMS chains (n ¼ 10) and 25 CL molecules is created using Towhee. The end-linking process continues until a crosslink density of 50% is achieved. Four different simulations are carried out: (1) CL with the harmonic bond potential (Equation (1)) without relaxation allowed between crosslink formations, (2) CL with the modified bond potential (Equation (7)) without relaxation, (3) CL with the harmonic bond potential allowing relaxation between crosslink formations and (4) CL with the modified bond potential with relaxation. To quantify the impact of the CL procedure, the bond, angle and dihedral distributions are computed and plotted as a histogram after the 50% crosslink density is achieved, as shown in Figure 3 for the OZSiZO angle which has an equilibrium value of 107.828.[25] Figure 3 shows the distribution for the OZSiZO angle in the case in which no relaxation is allowed between successive bond formations and the distribution in the case in which relaxation is allowed for 10 ps between crosslink bond formations. An optimal bin size minimising the error in the shape of the distribution is determined separately for each set of data based on the number of data points and the distribution of OZSiZO angles. Upon comparison of the two histograms, the presence of outliers is evident for the case in which CL without relaxation is carried out, especially in the case of the unmodified harmonic bond potential (Equation (1)). When relaxation of the structure is allowed, angles generated using Equation (7) display a


Molecular Simulation tighter distribution around the equilibrium OZSiZO angle. This indicates that using the modified bond potential with relaxation promotes the generation of a less distorted PDMS end-linked network. Furthermore, the influence of relaxation appears to be more important than the bond potential modification made in Equation (7). Histogram analysis conducted on bond and dihedral distributions for the same system yields similar results. As a second validation step, the glass transition temperature, Tg, for uncrosslinked and crosslinked PDMS models is computed using MD simulation. The glass transition temperature is the temperature at which a transition occurs from a soft ‘rubbery’ state to a solid ‘glassy’ state during cooling. Several authors have shown that MD simulations can be used to extract general trends in Tg, although exact comparison with experiment is not possible due to the timescale associated with MD simulation.[29] Two observations which are consistent with experiments [30] are made from these calculations: (1) a very modest increase in Tg is observed with PDMS chain molecular weight and (2) glass transition temperature is independent of CL. The lack of dependence of Tg on CL was observed experimentally in PDMS by Andrady et al. [30] who attributed this effect to the flexibility of the PDMS chains. Even though CL is known to introduce constraints on the mobility of the chains, in PDMS it fails to reduce the flexibility of the chain backbones, primarily due to the fact that crosslink bonds are formed only at the chain ends.

4.

Diffusion results

4.1 Diffusion coefficient calculation Simulations of diffusion in this study are carried out at temperatures ranging from 225 to 400 K in 25 K increments in an NPT ensemble with the system pressure maintained at 1 atm. All of these temperatures are above the glass transition temperature for PDMS.[30] A total of 100 united-atom O2 penetrants are used as the diffusing species; penetrants are introduced into the system at random (overlaps with PDMS chains are avoided) and the entire system is equilibrated for 500 ps. Finally, an NPT production run is carried out for 1 ns to study diffusion mechanisms and to extract diffusion coefficients for O2 through uncrosslinked PDMS and endlinked PDMS networks. A common method to compute diffusion coefficients using MD simulation is to track the mean-squared displacement (MSD) of the diffusing species.[23] Then, in the long time limit, the diffusion coefficient, D, is proportional to the slope of the MSD with respect to time, D¼

1d ½rðtÞ 2 rð0Þ2 ; 6 dt

ð8Þ

where the angle brackets indicate that the MSD is averaged over the number of penetrants in the simulation and rð0Þ

5

and rðtÞ are the positions of the penetrants at initial and later times during the simulation, respectively. As shown in Figure 4, the initial transient portion of the data is removed (first 30 ps) and linear fits are carried out over the last 70 ps of each simulation. Furthermore, the time origin of the MSD calculation is reset every 100 ps during the production run. Resetting of the time origin is analogous to carrying out 10 separate simulations to obtain the MSD of the diffusing species and hence allows better statistics in the computed average. In this study, the role of temperature on diffusion is captured using the Williams – Landel – Ferry (WLF) equation,[31] log ðDÞ ¼ y0 2

c1 ðT 2 T g Þ ; c2 þ ðT 2 T g Þ

ð9Þ

where c1 and c2 are material constants. The WLF equation provides a more accurate description of diffusion in polymers than the Arrhenius equation as temperatures approach glass transition.[32] In situations in which free volume availability is the critical factor which determines molecular mobility, the WLF model is appropriate to describe temperature dependence of the process.[33 –35] The material constants, c1 and c2 , are sometimes considered universal for a wide range of polymers.[31] However, Angell [35] argues that although these parameters are fairly robust, with c1 being more universal than c2 , it is best to fit the WLF equation to a given set of data. Parameter y0 provides a shift necessary to fit the WLF function to a specific phenomenon, such as viscosity in Urakawa et al. [36] Specifically, Urakawa et al. studied the viscosity of low-molecular weight polystyrene; the WLF equation provided an excellent model for the self-

Figure 4. (Colour online) Average MSD of 100 O2 penetrants during diffusion through one configuration of an uncrosslinked PDMS model with 50 chains with molecular weight of 5000 g/ mol. The transient portion of the MSD data is removed and linear fits are carried out at each temperature.

6


diffusivity of the molecular chains as the temperature was cooled towards Tg.


4.2 Role of chain length on O2 diffusion in uncrosslinked PDMS Figure 5 shows the temperature dependence of the O2 diffusion coefficient through uncrosslinked PDMS models with different molecular weights; Figure 5(a) shows the raw data (Equation (8)) while Figure 5(b) shows the log of the diffusion coefficient to illustrate the WLF model fits. Solid lines represent the best fit of the WLF equation to the data while perfect Arrhenius behaviour is represented in Figure 5(b) with a dashed line. Error bars in Figure 5(a)

illustrate ^ 1 standard deviation in D computed at each temperature from the three independent realisations. Figure 5 clearly shows that D increases with an increase in temperature and that D modestly decreases with an increase in chain length. The variation in D for models with different molecular weights is more pronounced at higher temperatures, as shown in Figure 5(a). The reduction in O2 diffusivity as chain length increases can be attributed to the presence of physical entanglements in the uncrosslinked PDMS models. Entanglements restrict the local movement of the PDMS chains, restricting the necessary evolution of the free volume in the system necessary for small molecule diffusion to occur. A previous MD study of diffusion of atmospheric penetrants through uncrosslinked PDMS systems consist-

Figure 5. (Colour online) Temperature dependence of the diffusion coefficient for O2 diffusion through PDMS showing the influence of molecular weight. The WLF fit captures the role of temperature. The Arrhenius relationship (valid well above Tg) is illustrated with a dashed line.


Molecular Simulation ing of 50 chains with n ¼ 10 repeating units (much shorter than that used in this study) reported a diffusion coefficient for O2 of 5.38 £ 1025 cm2/s at 300 K.[12] In this study, for the uncrosslinked PDMS model with a molecular weight of 5000 g/mol (n ¼ 66), the diffusion coefficient computed at 300 K is 2.06 £ 1025 cm2/s. This discrepancy is due to the fact that the molecular weights considered in this study are much larger than those used in the previous work. Other MD studies of O2 diffusion in PDMS [10,11] have reported diffusion coefficients between 1.8 £ 1025 and 3.7 £ 1025 cm2/s at 300 K, in agreement with the current MD simulations. Figure 5 also shows that the WLF equation provides an excellent fit to the diffusion data. The variation of the averaged WLF parameters as a function of molecular weight for uncrosslinked PDMS models is shown in Figure 6. No definitive trend is observed for parameters c1 and c2 , while a very modest decrease is observed in y0 over the range of molecular weights considered. For the four different molecular weight models, c1 values range from 5.65 to 6.07, which are relatively close to the value of 6.1 reported by Angell for PDMS.[35] The parameter c2 varies from 72.51 to 81.17 K and y0 parameter ranges from 2 8.36 to 2 8.83, with the smallest negative value corresponding to the 5000 g/mol uncrosslinked PDMS model. The authors are not aware of the reported c2 and y0 values for O2 diffusion through PDMS in the literature.

4.3 Role of end-linking and mixing ratio on O2 diffusion Figure 7 shows the influence of end-linking and PDMS:CL molecule mass ratio on O2 diffusion for each molecular weight considered in this study. MD simulations show that

7

O2 diffusion in the end-linked PDMS models is consistently higher than that in the uncrosslinked PDMS models at all temperatures and molecular weights considered. Furthermore, Figure 7 shows that the endlinked PDMS models with a mass ratio of 5:1 show higher O2 diffusivity than the 10:1 models at all temperatures and molecular weights considered. Recall from Table 1 that for each end-linked PDMS model, more CL molecules are included than are necessary to exactly achieve the 80% crosslink density criterion; the crosslinked model with a mass ratio of 5:1 has more CL molecules than the 10:1 model. Thus, it is hypothesised that the observation of enhanced diffusivity in end-linked PDMS systems is attributed to the presence of free or non-bonded CL molecules in these models, which due to their small size (relative to the PDMS chains) are significantly more locally mobile, facilitating the necessary evolution in the system free volume for O2 diffusion to occur. It is also interesting to note that deviation from Arrhenius behaviour at low temperatures (225 and 250 K) is reduced in the end-linked PDMS models as compared with the uncrosslinked PDMS samples. This observation can also be attributed to the high local mobility of the CL molecules as compared with that of the PDMS chains. To provide quantitative insight into the details of the end-linked PDMS network, Figure 8 presents a bar chart with the number of bonds that are created for each CL molecule (from 0 to 4) during the CL procedure (Figure 2) for mass ratios of 5:1 and 10:1 and for the four different molecular weights examined in this study. In Figure 8, the zero on the x-axis represents the number of free CL molecules present after 80% crosslink density is achieved. The largest PDMS end-linked model (molecular weight of 40,000 g/mol with a mass ratio of 5:1) has 1131 free CL

Figure 6. (Colour online) Relationship between chain length (molecular weight) and averaged WLF parameters for uncrosslinked PDMS models.


8


Figure 7. (Colour online) Temperature dependence of the diffusion coefficient for O2 diffusion through crosslinked PDMS models with 5:1 and 10:1 mass ratios and for molecular weights of (a) 5000 g/mol, (b) 10,000 g/mol, (c) 20,000 g/mol and (d) 40,000 g/mol. The WLF equation captures the role of temperature.

molecules, compared with 22 free CL molecules in the smallest model (molecular weight of 5000 g/mol with a mass ratio of 10:1). The large number of free CL molecules in PDMS models with the highest molecular weights results in enhanced O2 diffusivity. Only 4.65% of the CL molecules form two or more bonds in the largest model, while the percentage rises to 59.25% in the smallest case. Thus, a dense network is generated in the smaller models with significantly fewer non-bonded CL molecules. Similar to physical entanglements, crosslinks reduce the local mobility of the PDMS chains and accordingly retard the evolution of the free volume in the polymer network required for O2 diffusion.

4.4 Effect of chain length on diffusion in end-linked PDMS Figure 9 shows the influence of chain length on O2 diffusivity for constant crosslink density and mass ratio. Recall, in the case of uncrosslinked PDMS (Figure 5(b)), the chain length had a measurable influence on the diffusion of O2; however, in the end-linked PDMS models, the

delineation between PDMS samples with different molecular weights is negligible in the 5:1 mass ratio samples and very minimal in the 10:1 mass ratio samples. The high mobility of the CL molecule, relative to the PDMS chains, combined with their excessive presence (particularly in the 5:1 mass ratio models) masks the influence of molecular weight on O2 diffusivity. Deviation from perfect Arrhenius behaviour is observed in the crosslinked systems; however, it is not as severe as observed in the uncrosslinked PDMS samples. Crosslinked systems with a mass ratio of 10:1 show a greater deviation from Arrhenius behaviour than the 5:1 model, as the temperature is cooled towards the glass transition temperature. Figure 10 shows a plot of the variation of the average WLF parameters for crosslinked PDMS models with different molecular weights. For O2 diffusion through crosslinked PDMS systems with 5:1 mass ratio, a fit to the WLF equation provides c1 values ranging from 5.06 to 5.39, c2 values ranging from 76.72 to 95.66 K and the shift parameter y0 varying from 2 8.19 to 2 7.71 for the four different molecular weights studied. For the 10:1 mass ratio models, c1 values range from 5.24 to 6.01, c2 values


Molecular Simulation

Figure 8. (Colour online) Crosslink distribution for models with (a) 5:1 and (b) 10:1 mass ratios display the number of bonds added to each CL molecule to illustrate the underlying network structure.

range from 62.14 to 95.72 K and the y0 parameter varies from 2 8.90 to 2 7.91. Similar to the conclusions drawn from Figure 6, no definitive trends in the variation of c1 and c2 WLF parameters are apparent as a function of molecular weight of the PDMS chains. This insensitivity is attributed to the excessive number of free CL molecules in these systems, which mitigate any dependence of chain length. Furthermore, for the CL models, the shift y0 is also independent of chain length. However, clear trends emerge as a function of mass ratio. Parameter c1 is consistently lower in the 5:1 mass ratio models than in the 10:1 mass ratio models. In addition, when compared with Figure 6, c1 values for end-linked PDMS models are all lower than those for uncrosslinked PDMS models, leading to the deviation from Arrhenius behaviour observed in Figure 7. The y0 parameter is consistently higher in the 5:1 mass ratio models than in the 10:1 mass ratio models. Again, when compared with Figure 6, the majority of y0 values for end-linked PDMS models are higher (less negative) than those for uncrosslinked PDMS models.

9

Figure 9. (Colour online) Temperature dependence of the diffusion coefficient for O2 diffusion through crosslinked PDMS models with (a) 5:1 and (b) 10:1 mass ratios for varying molecular weights.

Figure 10. (Colour online) Relationship between averaged WLF parameters for various chain lengths at a weight ratio of 5:1 and 10:1 for crosslinked PDMS systems.

5. Conclusions In this study, MD simulations are used to study the steadystate diffusion of O2 penetrants in uncrosslinked and endlinked PDMS networks. This study uses silanol-terminated


10


PDMS chains and tetrakis(dimethylsiloxy)silane as the CL molecule. The number of PDMS chains is fixed at 50 and the number of CL molecules is varied to satisfy 5:1 or 10:1 by mass PDMS:CL molecule mixing ratios, chosen specifically to match experimental conditions.[4] By tracking the trajectories of the O2 penetrants in the system and subsequently computing their MSD, diffusion coefficients are obtained. Interestingly, MD simulations predict that O2 diffusion coefficients in end-linked PDMS models are higher than those in uncrosslinked PDMS models. The presence of a large number of non-bonded CL molecules in the system, as a result of the specified mixing ratio, leads to the enhanced O2 diffusivity. Quantitative analysis is carried out to illustrate the fraction of CL molecules that do not participate in the end-linking process; the fraction of free CL molecules varies significantly as a function of PDMS chain length and mixing ratio. The CL molecule has a molecular weight which is at least an order of magnitude smaller than the PDMS chains. Thus, this molecule has increased mobility within the end-linked PDMS network, compared with the PDMS chains, facilitating the evolution in the free volume necessary for small molecule diffusion through PDMS. This is emphasised in Figure 7, where for a given molecular weight the diffusivity can be arranged in a descending order: crosslinked system with weight ratio 5:1 (most CL molecules), crosslinked system with weight ratio 10:1 and uncrosslinked PDMS. The WLF equation provides a good fit to the MD simulation data. At temperatures well above glass transition, Tg, the diffusion of the O2 penetrants through PDMS and end-linked PDMS networks can be approximated as Arrhenius; however, as the simulation temperature is cooled towards Tg, diffusion of the O2 penetrants deviates from Arrhenius behaviour. The role of chain length and mixing ratio on the WLF parameters is presented in Figures 6 and 10. No clear trend is observed between chain length and the c1 or c2 WLF parameters, regardless of the presence or absence of CL. A modest trend is observed in the shift parameter, y0 , reflecting the role of chain length on diffusion in the uncrosslinked PDMS models. This trend is masked by the mobility of the CL molecules in the end-linked PDMS networks, as shown in Figure 9. On the other hand, clear trends are observed between mixing ratio and the WLF parameters. These observations can ultimately be used to predict diffusion rates of O2 molecules in PDMS samples as a function of mixing ratio. Acknowledgements This work was supported by the National Science Foundation under grant CMMI#0800718. MD simulations were carried out on cyberinfrastructure resources at the University of Arkansas which are supported in part by the National Science Foundation

under Grants ARI#0963249, MRI#0959124 (Razor), EPS #0918970 (CI TRAIN) and a grant from the Arkansas Science and Technology Authority.

References [1] Pooran R, Tung S, Kim JW. Application of a PDMS microsieve for the patterning of flagellar motors in a microfluidic system. MEMS. 2005;7:37–41. [2] Nour M, Berean K, Griffin MJ, Matthews GI, Bhaskaran M, Sriram S, Kalantar-Zadeh K. Nanocomposite carbon-PDMS membranes for gas separation. Sens Actuators B. 2012;161:982– 988. [3] Pan F, Spence H, Spearot D, Huang A. Nano-particle polymer composite MEMS corrosion. Proceedings of the 6th IEEE International Conference on Nano/Micro Engineered and Molecular Systems. 2011:1144–1148. [4] Dow Corning Corporation. Materials Safety Data Sheet for Sylgardw 184 Silicone Elastomer Base and Curing Agent ; 2010. [5] Trohalaki S, Rigby D, Klockowski A, Mark JE, Roe RJ. Estimation of diffusion coefficients for CO2 in polyethylene by molecular dynamics simulation. Polym Prepr. 1989;30:23– 24. [6] Takeuchi H, Okazaki K. Molecular dynamics simulation of diffusion of simple gas molecules in a short chain polymer. J Chem Phys. 1990;92:5643–5652. [7] Boyd RH, Krishna Pant PV. Molecular packing and diffusion in polyisobutylene. Macromolecules. 1991;24:6325–6331. [8] Sok RM, Berendsen HJC, van Gensteren WF. Molecular dynamics simulation of the transport of small molecules across a polymer membrane. J Chem Phys. 1992;96:4699–4704. [9] Tamai T, Tanaka H, Nakanishi K. Molecular simulation of permeation of small penetrants through membranes. 1. Diffusion coefficients. Macromolecules. 1994;27:4498 –4508. [10] Charati SG, Stern SA. Diffusion of gasses in silicone polymers: molecular dynamics simulations. Macromolecules. 1998;31: 5529–5535. [11] Jawalkar SS, Aminabhavi TM. Molecular dynamics simulations to compute diffusion coefficients of gases into polydimethylsiloxane and poly{(1,5-naphthalene)-co-[1,4-durene-2,20 -bis(3,4-dicarboxyl phenyl)hexaflouropropane diimide]}. Polym Int. 2007;56:928– 934. [12] Sudibjo A, Spearot DE. Molecular dynamics simulation of diffusion of small atmospheric penetratesin polydimethylsiloxane. Mol Simulat. 2011;37:115 –122. [13] Spearot DE, Sudibjo A, Ullal VN, Huang A. Molecular dynamics simulations of diffusion of O2 and N2 penetrants in polydimethylsiloxane-based nanocomposites. J Eng Mater Technol. 2012;134 (021013):1–8. [14] Shy LY, Leung YK, Eichinger BE. Computer simulation of gelation. Macromolecules. 1985;18:983–986. [15] Schulz M, Frisch HL. Microphase and macrophase separation in irreversible reacting chemical binary mixtures. J Chem Phys. 1994; 101:5013 –5016. [16] Doherty DC, Holmes BN, Leung P, Ross RB. Polymerization molecular dynamics simulations. I. Cross-linked atomistic models for poly(methacrylate) networks. Comput Theor Polym Sci. 1998;8:169–178. [17] Heine DR, Grest GS, Lorenz CD, Tsige M, Stevens MJ. Atomistic simulations of end-linked poly(dimethylsiloxane) networks: structure and relaxation. Macromolecules. 2004;37:3857–3864. [18] Tsige M, Lorenz CD, Stevens MJ. Role of network connectivity on the mechanical properties of highly cross-linked polymers. Macromolecules. 2004;37:8466–8472. [19] Tsige M, Stevens MJ. Effect of cross-linker functionality on the adhesion of highly cross-linked polymer networks: a molecular dynamics study of epoxies. Macromolecules. 2004;37:630– 637. [20] Komarov P, Tsung CY, Ming CS, Khalatur PG, Reineker P. Highly cross-linked epoxy resins: an atomistic molecular dynamics simulation combined with a mapping/reverse mapping procedure. Macromolecules. 2007;40:8104–8113. [21] Varshney V, Patnaik SS, Roy AK, Farmer BL. A molecular dynamics study of epoxy-based networks: cross-linking procedure

Molecular Simulation

[22] [23] [24] [25]


[26] [27] [28] [29]

and prediction of molecular and materials properties. Macromolecules. 2008;41:6837 –6842. Li C, Stratchan A. Molecular simulations of crosslinking process of thermosetting polymers. Polymer. 2010;51:6058–6070. Allen MP, Tildesley DJ. Computer simulation of liquids. Oxford: Clarendon; 2000. Torrens IM. Interatomic potentials. New York, NY: Academic Press; 1972. Frischknecht AL, Curro JG. Improved united atom force field for poly(dimethylsiloxane). Macromolecules. 2003;36:2122–2129. Sperling LH. Introduction to physical polymer science. Hoboken, NJ: Wiley; 2006. Martin MG, MCCCS Towhee. http://towhee.sourceforge.net/algor ithm/cbmc.html LAMMPS. http://lammps.sandia.gov Han J, Gee RH, Boyd RH. Glass transition temperatures of polymers from molecular dynamics simulations. Macromolecules. 1994;27: 7781–7784.

11

[30] Andrady AL, Llorente MA, Mark JE. Some dynamic mechanical properties of unimodal and bimodal networks of poly(dimethylsiloxane). Polym Bull. 1991;26:357–362. [31] Williams M, Landel R, Ferry J. The temperature dependence of relaxation mechanisms in amorphous polymers and other glassforming liquids. J Am Chem Soc. 1955;77:3701–3707. [32] Gabbott P. Principles and applications of thermal analysis. Ames: Blackwell Publishing; 2007. [33] Ferry JD. Viscoelastic properties of polymers. New York, NY: Wiley; 1980. [34] Lomellini P. Williams–Landel–Ferry versus Arrhenius behavior: polystyrene melt viscoelasticity revised. Polymer. 1992;33: 4983–4989. [35] Angell CA. Why C1¼16-17 in the WLF equation is physical and the fragility of polymers. Polymer. 1997;38:6261–6266. [36] Urakawa O, Swallen SF, Ediger MD, Von Meerwall ED. Selfdiffusion and viscosity of low molecular weight polystyrene over a wide temperature range. Macromolecules. 2004;37:1558–1564.