In Silico Screening of Green Plasticizers for Poly(vinyl chloride

Article Cite This: Macromolecules XXXX, XXX, XXX−XXX

pubs.acs.org/Macromolecules

In Silico Screening of Green Plasticizers for Poly(vinyl chloride) Hüsamettin D. Ö zeren,† Marcel Balçık,† M. Göktuǧ Ahunbay,*,† and J. Richard Elliott‡ †

Department of Chemical Engineering, Istanbul Technical University, Istanbul 34469, Turkey Department of Chemical and Biomolecular Engineering, University of Akron, Akron, Ohio 44325, United States

‡

Macromolecules Downloaded from pubs.acs.org by UNIV OF TEXAS AT DALLAS on 03/13/19. For personal use only.

S Supporting Information *

ABSTRACT: Phthalate derivative plasticizers used in poly(vinyl chloride) (PVC) processing have been a subject of concern because of their possible toxicity. Hence, there is a growing interest toward new, nontoxic, “green” plasticizers. In this work, the performances of biobased plasticizers including esters of succinic, levulinic, oleic, and adipic acids were compared in reference to the conventional plasticizer bis(2-ethylhexyl) phthalate (aka dioctyl phthalate or DOP). For this purpose, molecular dynamics (MD) simulations were used to determine polymer/plasticizer interactions and to predict thermomechanical properties of polymer mixtures. The variation of glass temperatures (Tg) of the systems was investigated, and the stabilities of the polymer/ plasticizer mixtures were compared through the Flory−Huggins solubility parameter. The mechanical properties were investigated through nonequilibrium MD simulations. Young’s modulus and yield strength values were predicted through stress−strain curves. The results suggest that succinic acid derivatives have the potential to replace phthalate derivatives due to their good solubility in PVC and their effectiveness in reducing the Tg.

■

INTRODUCTION

Plasticizers are used in formulations to produce relatively inexpensive materials. Thus, the plasticizers to replace phtalate derivatives must be inexpensive. The challenge to design “green” plasticizers for industrial applications starts from synthesis because these compounds are relatively expensive and have a more complex synthesis process than phthalate derivatives. Besides the plasticization properties and cost constraints, there are other concerns about performance of plasticizers in industrial applications. They should be substantially soluble in PVC since some applications may require up to 35 wt % plasticizer. On the other hand, the plasticizers should not migrate through the polymer over short time periods or the surface properties, and the long-term stability of the polymer may suffer. These constraints often establish a trade-off in the molecular size of the plasticizer. Among the phthalate-free alternatives, biobased “green” plasticizers are of particular interest because they are nontoxic to human health as well as renewable and biodegradable.7 To

Poly(vinyl chloride) (PVC) is one of the most widely used industrial plastics. Its mechanical properties and processability are tuned for targeted applications by the use of additives called plasticizers. Plasticizers facilitate processability and increase flexibility and toughness of the final product by lowering glass transition temperature and modifying the mechanical properties of the polymer. They are generally organic esters, among which phthalate derivatives are typical in the PVC industry. PVC accounts for 80−90% of world plasticizer consumption.1 The use of phthalates allows efficient and low-cost processes to manufacture PVC products with customized mechanical properties for a wide variety of industrial applications. However, significant concerns on the use of phthalates have arisen in recent years regarding environmental and health issues. Studies show that the hormone disruptor effect of phthalate derivatives is critical in early childhood development,2−5 motivating the PVC industry to shift toward the use of phthalate-free alternatives, especially in childhood products such as toys and care products. Relatedly, research into phthalate-free plasticizers has also increased for medical applications.6 © XXXX American Chemical Society

Received: October 8, 2018 Revised: March 1, 2019

A

DOI: 10.1021/acs.macromol.8b02154 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

where electrostatic charges were calculated using density functional theory (DFT) to obtain an accurate electrostatic distribution to be used with the AMBER force field. In addition to these studies, PVC plasticization was investigated with various plasticizers through molecular simulations. Palacios et al.23 studied PVC with several plasticizer pairs and investigated the plasticization effect on PVC. He used specifically glycerol and isosorbide for esters−PVC interaction with the MOPAC-PM6 force field and observed the thermodynamic interactions such as enthalpy ratio, free energy, and total energy ratio. The compatibility and performance of the plasticizers were predicted through the enthalpy and Gibbs free energy changes. In this study, various biobased molecules have been investigated to assess their potential to replace phthalate derivatives for PVC by employing atomistic scale simulations. The candidate molecules shown in Figure 1 were selected

date, some phthalate-free plasticizers have already been developed and used commercially for PVC in industrial applications. Various sebacate derivatives and acetyl tributyl citrate are commercially used. They provide superior low temperature flexibility properties which make them excellent choices for inks and coatings in food packaging applications. Furthermore, some epoxidized vegetable oils, which are produced from soybean oil or fatty acid ester, are commercially available. Some noncommercial green plasticizer candidates for PVC are also reported in the literature. Firlotte et al.8 have modeled 1,5-pentanediol dibenzoate, and they found that it is a good candidate compared to commercial plasticizers. Erythropel et al.9,10 studied the influence of of plasticizer geometry and alkyl groups on PVC. Their works were based on designing green plasticizers such as succinic, fumarate, and maleate acid derivatives. They worked on several aspects of plasticizer including mechanical tests, biodegradation, and Tg values. Glucose derivatives and other compounds were also tested as green plasticizers in the literature.11−13 In a recent study, macrodiol that was synthesized by the alcoholysis of highmolecular-weight poly(propylene carbonate) was reported to be a potential candidate to replace DOP as plasticizer.14 To overcome these challenges and to replace current commercial plasticizers with green plasticizers, identification of the important functional groups and their polymeric interactions is important. Molecular simulation is one of the efficient techniques to screen and evaluate the phthalate alternatives. Their compatibility with PVC can be estimated in terms of their solubility values in solution. Their effects on glass transition temperature and mechanical properties can be determined by molecular dynamics. In the past three decades, molecular simulation methods have been applied for amorphous atomic, molecular, ionic, and polymeric materials.15 Molecular dynamics (MD) and Monte Carlo (MC) are proven methods to investigate structure−property relationships of polymeric materials. Nonpolar and polar polymers have been studied, and their properties such as glass transition temperature (Tg) have been obtained with MD simulations.16,17 Han et al. applied MD simulations successfully to nonpolar polymers to estimate glass transition temperature by reproducing the pressure−volume−temperature (PVT) behavior of the polymers. Boyd et al. also studied PVT behavior and partial charge distribution with polyethylene, which has a polar contribution.18 These polar groups lead to stronger electrostatic interactions between the polymer chains, inducing a partial charge distribution. PVC is also a mildly polar polymer. PVC has also been being studied in the literature to develop an accurate force field. Smith et al.19 have integrated the OPLS-AA force field to amorphous atactic PVC polymer in melt simulation. Bond and angle parameters have been characterized separately for varying tacticities. The cell and chain size are important for polar polymers to model the polymer accurately. Abu-Sharkh20 has worked with Smith’s model and also modeled the CH2 and CHCl groups as rigid units in the polymer. The rigid unit simulations resulted in nearly same Tg as Smith’s model. PVC has also been modeled with force fields such as COMPASS21 and AMBER.22 Fermeglia et al. studied the compatibility of PVC blends through the integration of Flory−Huggins theory into polymer simulation. Shaikh et al. have studied hydration effects and antifouling properties of the membrane that has had PVC side in a copolymer membrane. They modeled neat PVC polymer

Figure 1. Selected biobased plasticizers.

among the esters of the top value added building blocks from biomass that were identified by the US Department of Energy (DOE).24 It is anticipated that these building blocks can be derived in sufficient quantities to supply a commodity industry like plasticizers. The approach of this study depends on MD simulations to determine polymer/plasticizer interactions and observe plasticization effects via the lowering of Tg. The solubility of the biobased plasticizers was evaluated with Flory−Huggins theory for compatibility of the plasticizer with amorphous PVC, and Tg’s of the polymer−plasticizer mixtures and neat polymers were obtained via pressure−volume− temperature (PVT) behavior of the polymer. The mechanical properties of the systems also were investigated for the systems through uniaxial MD simulations. The stress−strain curves B


Article

Macromolecules

Construction of Simulation Cells. The Amorphous Cell module of Accelrys Materials Studio 5.5 (San Diego, CA) was used in creating all simulation cells in this study. Amorphous PVC chains were constructed with equal fractions of the meso and racemic diads, distributed randomly along the chain to reproduce the structure of commercial PVC.33,34 Three PVC chains, each composed of 150 monomers, were then packed into the simulation cell at the experimental density (1.381 g/ cm3). The selected plasticizers were then added in the polymer with weight fractions of 10%, 20%, and 30%. The well-established 21-step slow decompression method was applied to each simulation cell to equilibrate the structures with selected force fields.35 The highest pressure (Pmax) in the procedure was chosen to be 200 bar. Details of the procedure are given in the Supporting Information. Systems were equilibrated at 300 K to compare their densities with experimental results. The LAMMPS molecular dynamics simulation package36 was used for the MD simulations. The Nosé−Hoover thermostat37,38 and barostat38,39 were used when temperature and pressure were required to be kept constant, respectively. The cutoff distances for both LJ and Coulombic interactions were set to 12 Å. The contribution of LJ and Coulombic interactions beyond the cutoff distance was calculated via tail correction and Ewald summation, respectively.40 Determination of Glass Transition Temperature. Experimentally, the glass temperature, Tg, is one of the most important descriptors in assessing the effectiveness of a plasticizer effect the on the chain mobility. In this study, Tg values were estimated via isotropic−isobaric cooling, where the polymer−plasticizer systems were cooled from 450 to 150 K with intervals of 20 K at 1 atm. MD simulations in the isotropic NPT were performed for 1.5 ns, amounting to a total of 24 ns. The last 200 ps of the simulations was used to obtain average specific volumes. The piecewise regression method, also known as broken-stick regression, was applied to determine Tg values: Two lines were fitted to our resulting temperature vs specific volume curve for the glassy and rubbery phase for each sample.16 The breakpoint of these two lines was identified as the Tg. Miscibility of the Polymer−Plasticizer Systems: Flory−Huggins (FH) Theory. FH theory provides a feasible theoretical framework for our work. The theory describes polymer−polymer or polymer−molecule interactions through the dimensionless FH interaction parameter (χFH). Solubility could also be estimated using the Hildebrand approach,41 but our analysis accounts for system specific interactions by evaluating the mixing energy directly through simulation results. System specific interactions are important when the polar mixtures are considered. The FH theory has been previously integrated successfully into MD simulations and applied to evaluate the solubility parameters of pharmaceuticals.42,43 Moreover, Caddeo et al. successfully applied the FH theory to a polymer−solvent mixture.44 The Flory−Huggins interaction parameter is dimensionless and temperature-dependent, and in the original theory it is concentration-independent, as shown in eq 3:

were demonstrated with the stress behavior of the systems at an applied constant strain rate.

■

METHODS The all-atom approach was adopted in modeling of the PVC chains to take into account polymer tacticity. Two different force fields, Generalized Amber Force Field (GAFF)25 and OPLS-AA Force Field,26 were considered in prediction of properties of amorphous PVC. A geometric combining rule was used for both force fields. The electrostatic charges that were used with the GAFF force field were adapted from a previous study by Shaikh,22 wherein charges were calculated using density functional theory (DFT). The charge parameters are listed in the Supporting Information. The electrostatic charges for the OPLS-AA force field are directly taken from the study of Hirvi et al.27 Coulombic interactions were included via the Ewald summation method. For atoms closer than three covalent bonds, nonbonded interactions were neglected in both force fields. For atoms separated by three covalent bonds, however, nonbonded interactions are taken into account, yet with a scaling factor, as they still have bonded interactions through the dihedral term. This scaling factor was 0.833 for the LJ interactions 0.5 for the electrostatic interactions in the GAFF force field, whereas it was 0.5 for both the LJ and electrostatic interactions in the OPLS-AA force field. Intramolecular bond and angle interactions of PVC were calculated using harmonic stretching and harmonic bending, respectively. The torsion potential energy was calculated with the Fourier equation, given in eq 1, when the GAFF force field is used. Etor =

∑

K i[1 + cos(niϕ − di)]

i = 1, m

(1)

where Ki is the torsion constant, ϕ is the torsion angle, ni is a constant integer dependent on i, and di is the phase (shift) angle. However, the OPLS-AA force field uses the harmonic dihedral equation in eq 2. 3

Etor =

∑ Sn(cos ϕ)n n=0

(2)

where Sn is the torsional energy constant. In OPLS-AA, contributions of torsions involving H and Cl atoms as end atom were neglected. The molecular models of the plasticizers were constructed based on the Transferable Potentials for Phase Equilibria (TraPPE) force field.28−32 The TraPPE force field provides accurate predictions of thermophysical properties of various esters, especially for liquid density. The 1−4 nonbonded Lennard-Jones (LJ) interactions were excluded in the TRAPPE force field while Coulombic interactions were scaled by a factor of 0.5. The parameters for the plasticizer structures are listed in the Supporting Information. The solubility of plasticizers in PVC is directly determined by polymer−plasticizer interactions; meanwhile, polymer− plasticizer interactions are dependent on their molecular structures. Hence, accurate modeling of both intermolecular and intramolecular interactions is crucial. In intermolecular interactions, the geometric combining rule was applied as a reasonable approximation. For intramolecular interactions within the polymer, combining rules dictated by the abovementioned force fields were applied.

χFH =

1 ΔHm 1 kT Vm ϕ1ϕ2

(3)

where ϕ1 and ϕ2 are the volume fractions of the mixing molecules, ΔHm is the enthalpy change upon mixing per unit C


Article

Macromolecules volume, and Vm is the molar volume the PVC monomer (42 cm3/mol). For a given plasticizer, χFH can be inferred from the ΔHm vs ϕ diagram. It relates also to the cohesive energy density (CED) that is defined as the energy of vaporization per unit volume: CED =

Ecoh V

simulations represented ideal systems, but experimental systems could have many microscopic or even macroscopic defects. In other earlier works, molecular dynamics simulations were also applied to predict elastic properties of polymer− carbon nanotube composites.48,49 The mechanical properties such as Young’s modulus and yield strength are derived from a stress−strain curve (SSC) which indicates the spatial deformations of a material when a stress is applied. Hossain et al.50 examined polyethylene in various regimes of temperature, stress rate, and chain length. Their results show the same deformation trends but different yield strength values. In our simulations, Young’s modulus was extracted by linear regression of the first 4% of the SSC data. Hossain et al.50 showed that the value of Young’s modulus can increase with the increase in the strain rate. Hence, the value of Young’s modulus was obtained via uniaxial MD deformation. It should be noted that MD simulations are expected to yield higher values due to higher strain rates than those used in experiment since the simulations are performed on the nanoscale. Furthermore, PVC becomes brittle under a high magnitude of stress, which makes it even more difficult to predict the experimental yield points of these materials, where the transition from elastic to nonelastic occurs. Therefore, the offset method can be applied to determine the yield point of these types of materials. Because the yield strength of the simulated systems is difficult to predict, a 0.2% offset is applied, wherein the linear line of Young’s modulus is shifted parallelly to the x-axis by 0.2% to identify the yield point.51 The deformation tests were performed via uniaxial NPT-MD simulations using LAMMPS, in which plasticizer incorporated PVC systems were deformed under tensile deformation with constant strain rate. The work of Hossain et al.50 clearly showed that tensile deformation needs higher cell size to reduce the fluctuations in SSC data. Therefore, the cells were constructed to contain 12 PVC chains having 150 monomers each, with plasticizer concentrations of 10, 20, and 30 wt %. The resulting cubic simulation cells had a cell length between 60 and 65 Å. Constant stress was applied to each principal axis of the simulation cell in deformation simulations, while the cell dimension in the other two directions were allowed to change at 1 atm and 300 K. This pressure was set to enable the two other directions to independently respond to the applied stress. Moreover, the mechanical test of neat PVC was repeated at 450 K to investigate the effect of temperature on the mechanical properties.

(4)

The term ΔHm in eq 3 can be expressed as the difference between the energy of the mixed and of the unmixed states per unit volume as in eq 5. E − x1E1 − x 2E2 ΔHm = total Vm Vm

(5)

where E1 and E2 are the individual mixings energy of each species. The energies of E1 and E2 refer to the pure component molar energies, where we neglect the contribution of excess volume to the mixture enthalpy. E E xV E xV ΔHm = total − 1 1 1 − 2 2 2 Vm Vm V1 Vm V2 Vm = CEDm − CED1ϕ1 − CED2 ϕ2

(6)

By combining eqs 3 and 6, and evaluating ΔE for a range of plasticizer loadings, one can estimate χFH through the equation χ12 =

Vref ΔE 1 kT Vm ϕ1ϕ2

(7)

In general, a positive value of the Flory−Huggins interaction parameter is considered indicative of the immiscibility of highmolecular-weight polymer blends, but actually the critical value of χFH obeys eq 8: χFH,c =

1 ijj 1 jj + 2 jk r1

1 r2

yz zz zz {

(8)

where r1 and r2 are the chain lengths of the polymer blend. However, the r1 value can be taken as 1 for the plasticizer in a polymer−plasticizer solution. In this case, χFH,c can be estimated around 0.5 ± 0.05. In this study, the FH theory was applied to estimate the miscibility of plasticizers in PVC at 300 and 450 K over a range of compositions. The PVC−DOP system was considered as a reference to assess the performance of bioplasticizers, where the benzene group of the phthalate plasticizer was constructed using the nine-site TraPPE model with three point charges inside the ring.31 Prediction of Polymer Mechanical Properties. Previous studies in the literature showed that MD simulations can be used to predict the mechanical properties of polymers and can yield results in agreement with the experimental data. Brown and Clarke45 performed detailed studies of uniaxial deformation at various strain rates of an amorphous polyethylene-like polymer system. The simulated elastic deformation captures the experimental deformation behavior. Li et al.46 have performed MC simulation to determine the mechanical behavior of polyethylene-like polymer. However, early works of elongation simulation overestimated Young’s modulus, shear modulus, bulk modulus, and yield strength. Fan et al.47 simulated polycarbonate, and the calculated values of Young’s and shear moduli were about twice as high as the experimental data. They explained this disagreement by reasoning that the

■

RESULTS AND DISCUSSION MD simulations were performed for PVC bulk polymer models to determine the validity of the selected force fields by reproducing the change in the specific volume of the polymer as a function of temperature (υ vs T) and by estimating its glass transition temperature (Tg). The comparison of the generated υ−T curves and Tg values with experimental data and previous computational works in the literature is shown in Figure 2 and summarized in Table 1. It can be seen in Figure 2 that the slope of the υ−T curve predicted with the AMBER force field changes near the experimental Tg, and predicted specific volume values are consistent with the experimental data. On the other hand, the OPLS-AA force field overpredicts the specific volumes compared to experimental data, even though obtained T g value is relatively consistent. In comparison, the overprediction reported by Abu-Sharkh20 and Hirvi et al.27 for the OPLS-AA force field is lower. Note D


Article

Macromolecules

Table 2. Change in the Glass Transition Temperature (Tg) of PVC with Plasticizer Content Tg (K) plasticizer/loading bis(2-ethylhexyl) phthalate (DOP) exptl53 this work dioctyl adipate (DOA) methyl oleate (MO) dimethyl succinate (DMS) diethyl succinate (DES) dibutyl succinate (DBS) bis(2-ethylhexyl) succinate (DEHS) methyl levulinate (ML) ethyl levulinate (EL) butyl levulinate (BL) 2-ethylhexyl levulinate (EHL)

Figure 2. Variation of the PVC specific volume with temperature. Comparison of the GAFF and OPLS force fields to the literature data: AMBER (gray ■), OPLS (red ◀), OPLS (blue ●),27 OPLS (yellow ▲),20 AMBER (black ◆),22 experimental (green ▶).52

52

experimental this work−GAFF this work−OPLS simulation−OPLS27 simulation−OPLS20 simulation−AMBER22

Tg (K)

density (g/cm3)

% error

343−353 345 352

1.381 1.385 1.268 1.302 1.316 1.362

0.289 −8.911 −6.068 −4.939 −1.395

20 wt %

30 wt %

319 320 318 297 319 303 301 297 330 327 318 317

271 281 303 292 300 297 288 279 316 303 298 297

243 263 276 283 283 278 263 256 300 296 280 277

for the biobased plasticizer candidates are also shown in Table 2. The assessment of the results in Table 2 suggests that as a replacement plasticizer DEHS may be as effective as DOP, since it is the only biobased plasticizer candidate that can reduce the Tg to a value lower than DOP does, for all three loadings. It is also of interest to estimate the effect of the number of carbon atoms in the substituents of ester groups on Tg. Figure 3 shows the computed Tg values of succinic acid derivatives as

Table 1. Comparison of Densities at 300 K and Tg Results with Previous Works system

10 wt %

also that Shaikh et al.22 had slight inconsistency in specific volume even though their model was very similar to ours, which may be due to the geometric combining rule in the GAFF model, instead of the Lorentz−Berthelot rule. The force field validation of the selected biobased plasticizer was conducted via density calculation at 300 K. Each plasticizer, which was constructed individually in Material Studio, was packed in a simulation cell at the experimental density. Sizes of the resulting simulation cells were 30 and 35 Å. The systems were then equilibrated via MD simulation in the NVT ensemble first at 600 K and then at 300 K for 1.5 ns at each temperature. Finally, the densities of the cells were calculated via NPT-MD simulations at 1 atm and at 300 K for 2 ns, of which the last 200 ps was used to estimate the densities. The results are compared to the experimental density data in Table 1. The glass temperature decrease as a function of plasticizer loading was first evaluated for the reference GAFF-PVC/DOP system. The Tg values for the polymer containing 10, 20, and 30 wt % DOP are presented in Table 2. Although the error in the estimated specific volumes increases with increasing DOP loading in solution, the decreasing trend in the Tg is well reproduced. This increase in error is probably due to MD timescales on the order of nanoseconds that correspond to a very fast cooling rate. This prevents individual atoms to pack in the most efficient way, and consequently, leads to Tg values that are higher than the experimental ones measured over much longer time scales.20 The changes in Tg of the polymer

Figure 3. Effect of carbon length in the substituents of ester groups on Tg of PVC for succinic derivatives for three different loadings (10, 20, and 30 wt %).

a function of the number of carbon atoms in the substituent group. Succinic acid derivatives and levulinic acid derivatives have similar backbone structures. However, the difference is that levulinic acid has only one ester group and one methyl group connected to the backbone. This makes the levulinic acid derivatives less effective to interfere with PVC chains. For a deeper understanding of the effect of carbon number in the backbone on Tg, DOA and DEHS were compared. DOA has four and DEHS has two carbon atoms in their backbones, respectively, with the same ester group. The results show that increasing the carbon atoms in the backbone increases the E


Article

Macromolecules flexibility of the structure since adipate-based polyesters have lower Tg,54,55 but yet it failed to lower the Tg as expected. This might be due to the linear substituent group instead of a branched one in the structure. Consequently, the backbone with lower carbon number and longer carbon chain in the substituents of the ester group in the biobased plasticizer is more effective as plasticizer for PVC polymer. Flory−Huggins Solubility Parameter. The Flory− Huggins interaction parameter, χ12, is estimated by fitting it to eq 7 using the polymer−plasticizer interaction energy at three different plasticizer loadings at 300 and 450 K. The fitting implies that the χ12 value for each plasticizer is independent of the plasticizer loading. However, the plots of χ12 value with respect to volume fraction of PVC in the Supporting Information show that the χ12 parameter varies with PVC volume fraction at 300 K for all plasticizers except DOP and succinates. Nevertheless, for the simulated systems, all χ12 values are obtained negative; this suggests that the molecules exhibit mutual attraction to each other,56−58 and the increasing value shows that the polymer−plasticizer mixture is more thermodynamically miscible.59 As an example, Figure 4

Figure 5. Flory−Huggins diagram of succinic acid with different chain lengths at 300 K.

Figure 4. Flory−Huggins diagram of DOP at 300 and 450 K.

shows the FH diagram for the DOP/PVC system at 300 and 450 K. While the individual FH diagrams for the rest of the plasticizers are provided in the Supporting Information, it is instructive to compare the behavior of the succinic acid derivatives with different chain lengths at 300 K, as shown in Figure 5. It can be seen that in Figure 4 (and in those provided in the Supporting Information) the free energy decreases with increasing temperature, which hereby confirms that PVC− plasticizer systems exhibit lower critical solution temperature (LCST) type mixture behavior.59 Furthermore, the comparison in Figure 5 shows intermolecular energy decreases when the carbon length increases. Hence, χ12 decreases with increasing chain length of the end group of plasticizer as also given in Figure 6, where the same colors represent the same type of biobased molecule with a different substituent of ester. It is clear that χ12 values of DEHS, DOA, and EHL are close to each other while that of MO is slightly lower compared to DOP. Furthermore, the values of χ12 point out those plasticizers that have high molecular weight are less soluble in polymers. On the other hand, plasticizers having low

Figure 6. χ12 parameters of PVC−plasticizer systems at (a) 300 K and (b) 450 K.

molecular weight can easily migrate from the final product, resulting in a trade-off between solubility and migration tendency. Hence, to design an efficient plasticizer, it will be necessary to understand the relationship between solubility and mobility. As DOP is soluble, as well as an efficient plasticizer, the calculated χ12 values suggest that DEHS, DOA, and EHS have also similar solubility in PVC, and they would have the potential to be effective plasticizers. Stress−Strain Curves. Three different strain rates of 108, 9 10 , and 1010 s−1 were used to observe the strain rate effect on neat PVC. We selected the SSC data of 108 s−1 up to 2% due to pronounced severe nonlinearity in the stress−strain curve. For the other strain rates, we followed the same procedure mentioned above. These rates were applied in uniaxial MD simulations, and stress−strain responses are plotted in Figure F


Article

Macromolecules

This behavior in strain softening at higher rates has also been reported.63−65 It should be noted that a rigid PVC sample would be fractured when it is loaded experimentally, but our simulation is not designed for this particular behavior. The predicted elastic moduli of the neat PVC sample at the investigated strain rates are also close to experimental values.66 Therefore, it can be concluded that uniaxial MD simulations can predict mechanical properties of neat PVC in reasonable agreement with experiments. The influence of temperature on the mechanical behavior of PVC was also investigated by conducting same MD simulations of uniaxial tension at various temperatures. Figure 9 shows that SSC data taken from MD deformation tests using

7. The simulations showed that yield strength strongly and linearly depends on logarithmic strain rate while the elastic

Figure 7. Stress−strain response of the neat PVC system at 300 K under three different strain rates (108, 109, and 1010 s−1).

modulus is less sensitive, consistent with the report by Yang and Qu.60 Figure 8 shows that the computed Young’s modulus and yield strength vary linearly from 3.3 to 3.7 with logarithmically

Figure 9. Thermal effect on the stress−strain behavior of PVC at 150, 300, and 450 K

a strain rate of 108 s−1 at 150, 300, and 450 K. The change in the temperature shows a very clear effect on the PVC stress− strain response. A decrease in the elastic modulus is seen with increasing temperatures, as expected. The sub-Tg SSC data follow the typical stress−strain curve behavior. The yield strength and softening regime are more distinguishable at lower temperatures. When the temperature is above Tg, these points are not easily distinguished, which means that the material is flexible and easily processable. The influence of the DEHS plasticizer, which is suggested to be the most effective plasticizer candidate according to computed Tg values, was also investigated by conducting the same MD simulations of uniaxial tension for the three different DEHS loadings of 10, 20, and 30 wt % at a strain rate of 109 s−1 and 300 K. As can be seen in Figure 10, increasing the DEHS loading decreased elastic modulus significantly, which indicates that the material is plasticized and gains more flexibility. These results follow the same trend as the computed Tg values, yet the computed moduli are close to glassy regions. Experimental studies showed that the moduli were linearly dependent on logarithmic strain rate, even in slow strain rates.61 Although computed moduli decrease significantly, they are higher than the experimental values because of the unavoidable extremely higher tensile strain rates in the simulations. This can be seen in the comparison of the simulated and experimental strength data. However, Kendall et al.63 reported that even plasticized PVC exhibits glassy polymer

Figure 8. Elastic modulus (black ■) and yield strength (red ●) derived under different loading rates for the neat PVC system at the temperature of 300 K.

increasing strain rates from 108 to 1010 s−1. This linear− logarithmic relationship for strength has been observed experimentally.61 Even though the magnitudes of used strain rates are higher by far, the strength values are just slightly higher than the experimental ones. The SSC data of 108 and 109 s−1 follow the same trend compared to experimental mechanical test.62,63 Furthermore, as shown in Figure 7, the strain rate of 1010 s−1 has strain softening while other rates have strain hardening. The increased strain softening of neat PVC is expected to be due to the higher rate; it can be expounded by adiabatic heat stored in the system, where the high rate does not allow any heat release during deformation. G


Article

Macromolecules

■

Force field functions and parameters for PVC and biobased plasticizers; 21-step slow decompression equilibration method; specific volume vs temperature plots of plasticizer−polymer systems; Flory−Huggins diagrams of plasticizer−polymer systems (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail [email protected]. ORCID

M. Göktuǧ Ahunbay: 0000-0003-3595-6419 J. Richard Elliott: 0000-0002-1159-4786 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the Scientific and Technological Research Council of Turkey (TUBITAK) through grant # 114M178. Helpful discussions with Roger Avakian of PolyOne Corp are also gratefully acknowledged.

Figure 10. Effect of DEHS content in PVC (wt %) on elastic modulus at a strain rate of 109 s−1 and 300 K.

behavior experimentally at higher strain rates. This also explains the higher computed moduli of PVC/DEHS systems due to the shifting of α-transition to the higher temperature.

■

■

CONCLUSION The plasticization effect of biobased plasticizer candidates and their solubility in PVC were investigated at the atomistic level through MD simulations. PVC was modeled with GAFF and OPLS-AA force fields, and the validation of these force fields was conducted through the comparison of the specific volume and Tg with experimental results. Flory−Huggins interaction parameters (χ12) of PVC/ plasticizer systems were calculated by interpolating the FH diagrams at 300 and 450 K. The solubility comparison of biobased plasticizers was demonstrated by calculated χ12 with respect to simulated DOP. As expected, increasing the number of carbon atoms in the substituent groups of esters decreased intermolecular energy significantly and FH interaction parameter (χ12). Furthermore, increasing temperature also decreased the magnitude (made less negative) of the χ12 parameter. A more negative value of χ12 indicates that plasticizer is more soluble in the polymer. According to this, the selected plasticizers were clearly soluble in PVC except methyl oleate, which is shown to interact weakly with PVC. This could lead to phase separation in the PVC/MO mixture. The results are very promising and suggest that biobased plasticizers have the potential to replace phthalate derivatives. Especially, succinic acid derivatives show good compatibility, and they are more effective on Tg compared to levulinic derivatives. One of the succinic acid derivatives, DEHS, has a higher potential to be used in industrial applications as a replacement for DOP. In summary, it has been demonstrated in this study that purely atomistic simulations are potent tools to design polymer additives such as plasticizers. As a future work, this work can be extended to understand the relationship between solubility and mobility/migration tendency of the plasticizer candidates.

■

REFERENCES

(1) Knight, D. J. EU Regulation of Chemicals: REACH. Rapra Rev. Rep. 181 2006, 16 (1). (2) Kovacic, P. How Dangerous Are Phthalate Plasticizers? Integrated Approach to Toxicity Based on Metabolism, Electron Transfer, Reactive Oxygen Species and Cell Signaling. Med. Hypotheses 2010, 74 (4), 626−628. (3) Latini, G. Monitoring Phthalate Exposure in Humans. Clin. Chim. Acta 2005, 361 (1−2), 20−29. (4) Krauskopf, L. G. How about Alternatives to Phthalate Plasticizers? J. Vinyl Addit. Technol. 2003, 9 (4), 159−171. (5) Chiellini, F.; Ferri, M.; Morelli, A.; Dipaola, L.; Latini, G. Perspectives on Alternatives to Phthalate Plasticized Poly(Vinyl Chloride) in Medical Devices Applications. Prog. Polym. Sci. 2013, 38 (7), 1067−1088. (6) Gimeno, P.; Thomas, S.; Bousquet, C.; Maggio, A.-F.; Civade, C.; Brenier, C.; Bonnet, P.-A. Identification and Quantification of 14 Phthalates and 5 Non-Phthalate Plasticizers in PVC Medical Devices by GC−MS. J. Chromatogr. B: Anal. Technol. Biomed. Life Sci. 2014, 949, 99−108. (7) Kermanshahi, A.; Cooper, D. G.; Maric, M.; Nicell, J. A. Towards the Development of Green Plasticizers. 13th Annu. Green Chem. Eng. Conf., 2009. (8) Firlotte, N.; Cooper, D. G.; Marić, M.; Nicell, J. A. Characterization of 1,5-Pentanediol Dibenzoate as a Potential “Green” Plasticizer for Poly(Vinyl Chloride). J. Vinyl Addit. Technol. 2009, 15 (2), 99−107. (9) Erythropel, H. C.; Dodd, P.; Leask, R. L.; Maric, M.; Cooper, D. G. Designing Green Plasticizers: Influence of Alkyl Chain Length on Biodegradation and Plasticization Properties of Succinate Based Plasticizers. Chemosphere 2013, 91 (3), 358−365. (10) Erythropel, H. C.; Shipley, S.; Börmann, A.; Nicell, J. A.; Maric, M.; Leask, R. L. Designing Green Plasticizers: Influence of Molecule Geometry and Alkyl Chain Length on the Plasticizing Effectiveness of Diester Plasticizers in PVC Blends. Polymer 2016, 89, 18−27. (11) Jia, P.; Zhang, M.; Hu, L.; Zhou, Y. Green Plasticizers Derived from Soybean Oil for Poly(Vinyl Chloride) as a Renewable Resource Material. Korean J. Chem. Eng. 2016, 33 (3), 1080−1087. (12) Shi, G.; Cooper, D. G.; Maric, M. Poly(ε-Caprolactone)-Based “green” Plasticizers for Poly(Vinyl Choride). Polym. Degrad. Stab. 2011, 96 (9), 1639−1647. (13) Kastner, J.; Cooper, D. G.; Marić, M.; Dodd, P.; Yargeau, V. Aqueous Leaching of Di-2-Ethylhexyl Phthalate and “Green” Plasticizers from Poly(Vinyl Chloride). Sci. Total Environ. 2012, 432, 357−364.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.macromol.8b02154. H


Article

Macromolecules (14) Li, W.; Qin, J.; Wang, S.; Han, D.; Xiao, M.; Meng, Y. Macrodiols Derived from CO2-Based Polycarbonate as an Environmentally Friendly and Sustainable PVC Plasticizer: Effect of Hydrogen-Bond Formation. ACS Sustainable Chem. Eng. 2018, 6 (7), 8476−8484. (15) Abrams, C.; Site, L. D.; Kremer, K. Multiscale Computer Simulations for Polymeric Materials in Bulk and Near Surfaces, 2002. (16) Han, J.; Gee, R. H.; Boyd, R. H. Glass Transition Temperatures of Polymers from Molecular Dynamics Simulations. Macromolecules 1994, 27, 7781−7784. (17) Boyer, R. F. Effect of Plasticizers on Some Physical Properties of Polymers. Tappi 1951, 34, 357−362. (18) Boyd, R. H.; Gee, R. H.; Han, J.; Jin, Y. Conformational Dynamics in Bulk Polyethylene: A Molecular Dynamics Simulation Study. J. Chem. Phys. 1994, 101 (1), 788−797. (19) Smith, G. D.; Jaffe, R. L.; Yoon, D. Y. Conformations and Order in Atactic Poly(Vinyl Chloride) Melts from MolecularDynamics Simulations. Macromolecules 1993, 26 (2), 298−304. (20) Abu-Sharkh, B. Glass Transition Temperature of Poly (Vinylchloride) from Molecular Dynamics Simulation: Explicit Atom Model versus Rigid CH 2 and CHCl Groups Model. Comput. Theor. Polym. Sci. 2001, 11, 29−34. (21) Fermeglia, M.; Coslanich, A.; Ferrone, M.; Humar, T.; Paneni, M. S.; Pricl, S.; Paglione, E. PVC/LDPE Incompatible Blends for Industrial Applications: A Computational Multiscale Approach. AIChE Annu. Meet. Conf. Proc. 2004. (22) Shaikh, A. R.; Rajabzadeh, S.; Matsuo, R.; Takaba, H.; Matsuyama, H. Hydration Effects and Antifouling Properties of Poly(Vinyl Chloride-Co-PEGMA) Membranes Studied Using Molecular Dynamics Simulations. Appl. Surf. Sci. 2016, 369, 241−250. (23) Suárez Palacios, O. Y.; Narváez Rincón, P. C.; Camargo Pardo, M.; Corriou, J.-P. Methodology To Predict PVC Plasticization Using Molecular Simulation by Pairs. Ind. Eng. Chem. Res. 2013, 52 (43), 15094−15103. (24) Werpy, T.; Petersen, G. Top Value Added Chemicals from. Biomass. U. S. Dept. Energy 2004, 1, 76. (25) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25 (9), 1157−1174. (26) Jorgensen, W. L.; Maxwell, D. S.; Tiradorives, J. Development and Testing of the OPLS All Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (27) Hirvi, J. T.; Pakkanen, T. A. Molecular Dynamics Simulations of Water Droplets on Polymer Surfaces. J. Chem. Phys. 2006, 125 (14), 144712. (28) Maerzke, K. A.; Schultz, N. E.; Ross, R. B.; Siepmann, J. I. TraPPE-UA Force Field for Acrylates and Monte Carlo Simulations for Their Mixtures with Alkanes and Alcohols. J. Phys. Chem. B 2009, 113 (18), 6415−6425. (29) Kamath, G.; Robinson, J.; Potoff, J. J. Application of TraPPEUA Force Field for Determination of Vapor-Liquid Equilibria of Carboxylate Esters. Fluid Phase Equilib. 2006, 240 (1), 46−55. (30) Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. United-Atom Description of n -Alkanes. J. Phys. Chem. B 1998, 102 (97), 2569−2577. (31) Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 6. United-Atom Description for Ethers, Glycols, Ketones, and Aldehydes. J. Phys. Chem. B 2004, 108 (45), 17596−17605. (32) Wick, C. D.; Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 4. United-Atom Description of Linear and Branched Alkenes and Alkylbenzenes. J. Phys. Chem. B 2000, 104 (33), 8008−8016. (33) Smith, G. D.; Ludovice, P. J.; Jaffe, R. L.; Yoon, D. Y. Conformations of 2,4-Dichloropentane and 2,4,6-Trichloroheptane and a Force Field for Poly(Vinyl Chloride) Based upon Ab Initio Electronic Structure Calculations. J. Phys. Chem. 1995, 99 (1), 164− 172.

(34) Smith, G. D.; Jaffe, R. L.; Yoon, D. Y. Conformations and Order in Atactic Poly(Vinyl Chloride) Melts from MolecularDynamics Simulations. Macromolecules 1993, 26 (2), 298−304. (35) Larsen, G. S.; Lin, P.; Hart, K. E.; Colina, C. M. Molecular Simulations of Pim-1-like Polymers of Intrinsic Microporosity. Macromolecules 2011, 44 (17), 6944−6951. (36) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1. (37) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31 (3), 1695− 1697. (38) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81 (1), 511−519. (39) Hoover, W. G. Constant-Pressure Equations of Motion. Phys. Rev. A: At., Mol., Opt. Phys. 1986, 34 (3), 2499−2500. (40) Ewald, P. P. Die Berechnung Optischer Und Elektrostatischer Gitterpotentiale. Ann. Phys. 1921, 369 (3), 253−287. (41) Hildebrand, J. H. Is There a “Hydrophobic Effect”? Proc. Natl. Acad. Sci. U. S. A. 1979, 76 (1), 194. (42) Huynh, L.; Grant, J.; Leroux, J.-C.; Delmas, P.; Allen, C. Predicting the Solubility of the Anti-Cancer Agent Docetaxel in Small Molecule Excipients Using Computational Methods. Pharm. Res. 2008, 25 (1), 147−157. (43) Kasimova, A. O.; Pavan, G. M.; Danani, A.; Mondon, K.; Cristiani, A.; Scapozza, L.; Gurny, R.; Möller, M. Validation of a Novel Molecular Dynamics Simulation Approach for Lipophilic Drug Incorporation into Polymer Micelles. J. Phys. Chem. B 2012, 116 (14), 4338−4345. (44) Caddeo, C.; Mattoni, A. Atomistic Investigation of the Solubility of 3-Alkylthiophene Polymers in Tetrahydrofuran Solvent. Macromolecules 2013, 46 (19), 8003−8008. (45) Brown, D.; Clarke, J. H. R. Molecular Dynamics Simulation of an Amorphous Polymer under Tension. 1. Phenomenology. Macromolecules 1991, 24 (8), 2075−2082. (46) Li, J.; Mulder, T.; Vorselaars, B.; Lyulin, A. V.; Michels, M. A. J. Monte Carlo Simulation of Uniaxial Tension of an Amorphous Polyethylene-like Polymer Glass. Macromolecules 2006, 39 (22), 7774−7782. (47) Fan, C. F.; Cagin, T.; Chen, Z. M.; Smith, K. A. Molecular Modeling of Polycarbonate. 1. Force Field, Static Structure, and Mechanical Properties. Macromolecules 1994, 27, 2383−2391. (48) Griebel, M.; Hamaekers, J. Molecular Dynamics Simulations of the Elastic Moduli of Polymer-Carbon Nanotube Composites. Comput. Methods Appl. Mech. Eng. 2004, 193 (17−20), 1773−1788. (49) Han, Y.; Elliott, J. Molecular Dynamics Simulations of the Elastic Properties of Polymer/Carbon Nanotube Composites. Comput. Mater. Sci. 2007, 39 (2), 315−323. (50) Hossain, D.; Tschopp, M. A.; Ward, D. K.; Bouvard, J. L.; Wang, P.; Horstemeyer, M. F. Molecular Dynamics Simulations of Deformation Mechanisms of Amorphous Polyethylene. Polymer 2010, 51 (25), 6071−6083. (51) Rottler, J.; Robbins, M. O. Yield Conditions for Deformation of Amorphous Polymer Glasses. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2001, 64 (5), 8. (52) Walsh, D.; Zoller, P. Standard Pressure Volume Temperature Data for Polymers; CRC Press: 1995. (53) Heydemann, P.; Guicking, H. D. Specific Volume of Polymers as a Function of Temperature and Pressure. Kolloid Z. Z. Polym. 1963, 193 (1), 16. (54) Park, J. W.; Im, S. S. Phase Behavior and Morphology in Blends of Poly(L-Lactic Acid) and Poly(Butylene Succinate). J. Appl. Polym. Sci. 2002, 86 (3), 647−655. (55) Penning, J. P.; St. John Manley, R. Miscible Blends of Two Crystalline Polymers. 1. Phase Behavior and Miscibility in Blends of Poly(Vinylidene Fluoride) and Poly(1,4-Butylene Adipate). Macromolecules 1996, 29 (1), 77−83. (56) Silva, M. A.; De Paoli, M. A.; Felisberti, M. I. Flory-Huggins Interaction Parameter of Poly(Ethylene Oxide)/PolyI


Article

Macromolecules (Epichlorohydrin) and Poly(Ethylene Oxide)/Poly(EpichlorohydrinCo-Ethylene Oxide) Blends. Polymer 1998, 39 (12), 2551−2556. (57) Silva, M.; Mano, V.; Pacheco, R. R. J.; Freitas, R. F. S. Miscibility Behavior of Polyacrylamides Poly (Ethylene Glycol) Blends: Flory Huggins Interaction Parameter Determined by Thermal Analysis. J. Mod. Phys. 2013, 4 (7), 45−51. (58) Men’shikov, E. A.; Bol’shakova, A. V.; Yaminskii, I. V. Determination of the Flory-Huggins Parameter for a Pair of Polymer Units from AFM Data for Thin Films of Block Copolymers. Prot. Met. Phys. Chem. Surf. 2009, 45 (3), 295−299. (59) Frezzotti, D.; Ravanetti, G. P. Evaluation of the Flory-Huggins Interaction Parameter for Poly(Styrene-Co-Acrylo-Nitrile) and Poly(Methyl-Methacrylate) Blend from Enthalpy of Mixing Measurements. J. Therm. Anal. 1994, 41 (6), 1237−1243. (60) Yang, S.; Qu, J. Computing Thermomechanical Properties of Crosslinked Epoxy by Molecular Dynamic Simulations. Polymer 2012, 53 (21), 4806−4817. (61) Retting, W. Das Mechanische Verhalten von Thermoplasten Bei Mittlerer Und Hoher Deformation Bis Zum Bruch. Rheol. Acta 1969, 8 (3), 259−267. (62) Rostam, S.; Ali, A. K.; AbdalMuhammad, F. H. Experimental Investigation of Mechanical Properties of PVC Polymer under Different Heating and Cooling Conditions. J. Eng. 2016, 2016, 1−5. (63) Kendall, M. J.; Siviour, C. R. Rate Dependence of Poly(Vinyl Chloride), the Effects of Plasticizer and Time-Temperature Superposition. Proc. R. Soc. London, Ser. A 2014, 470 (2167), 20140012. (64) Arruda, E. M.; Boyce, M. C.; Jayachandran, R. Effects of Strain Rate, Temperature and Thermomechanical Coupling on the Finite Strain Deformation of Glassy Polymers. Mech. Mater. 1995, 19 (2−3), 193−212. (65) Kendall, M. J.; Froud, R. F.; Siviour, C. R. Novel Temperature Measurement Method & Thermodynamic Investigations of Amorphous Polymers during High Rate Deformation. Polymer 2014, 55 (10), 2514−2522. (66) Povolo, F.; Schwartz, C.; Hermida, É . B. Stress Relaxation of PVC below the Yield Point. J. Polym. Sci., Part B: Polym. Phys. 1996, 34 (7), 1257−1267.

J


In Silico Screening of Green Plasticizers for Poly(vinyl chloride

Recommend Documents