Modeling Polymer Glass Transition Properties from Empirical

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

pubs.acs.org/Macromolecules

Modeling Polymer Glass Transition Properties from Empirical Monomer Data with the SAFT‑γ Mie Force Field Amulya K. Pervaje,† Joseph C. Tilly,† David L. Inglefield, Jr.,§ Richard J. Spontak,†,‡ Saad A. Khan,† and Erik E. Santiso*,† Department of Chemical and Biomolecular Engineering and ‡Department of Materials Science and Engineering, North Carolina State University, Raleigh, North Carolina 27606, United States § Eastman Chemical Company, Kingsport, Tennessee 37660, United States Downloaded via KAOHSIUNG MEDICAL UNIV on November 20, 2018 at 11:42:15 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

S Supporting Information *

ABSTRACT: We apply a recently developed coarse-graining method to build models for polyester polyols, versatile polymers with applications in coatings, by combining models for the component monomers. This strategy employs the corresponding states correlation to the group-contribution SAFT-γ Mie equation of state [Mejia, A.; et al. Ind. Eng. Chem. Res. 2014, 53, 4131−4141] to obtain force-field parameters for the constituent monomer species. Results from simulations agree favorably with experimental values of mass density, glass transition temperature (Tg), and specific heat capacity change at Tg. Further simulations over a range of Mie parameters and polymer chemical compositions yield a correlation that relates the parameters directly to Tg. This correlation is validated by experimental data and can be used as a predictive tool within the tested parameter space to expedite the design of these coating materials.

1. INTRODUCTION Polymers synthesized from a variety of monomer species in specific combinations yield a wide range of thermomechanical properties. Modifying the chemical makeup and structure of polymers enables tunable properties.1−4 Polyester polyols constitute a technologically significant family of macromolecules that possess an assortment of chemical structures. Generally speaking, they contain multiple ester bonds along their backbone and are terminated by hydroxyl groups. They can be synthesized from diols and diacids, with a variety of different monomer identities and ratios possible, along with branching introduced via incorporation of triols.5 These resins are broadly employed in the development of coatings, adhesives, sealants, polyurethane foams, and elastomers.6−8 The identities of the monomers comprising the polyester polyol backbone greatly influence the properties of the polymer, as well as the properties of chemically modified versions of the polymer that might be needed for targeted applications.9−11 Because the glass transition behavior of polymers affects material processing and performance, being able to predict the glass transition temperature (Tg) for different possible monomer combinations and ratios would greatly expedite the design of polymers such as polyester polyols to satisfy application-specific requirements. The glass transition region and Tg are important in both amorphous and semicrystalline polymers since properties change dramatically as the polymer transitions from rigid and glassy below Tg to flexible and rubbery above Tg. The fundamental basis for the glass transition has been widely discussed in the literature.12−17 Experimental methods including differential scanning calorimetry (DSC) and dynamic © XXXX American Chemical Society

mechanical thermal analysis (DMTA) are commonly used to measure Tg.18 In DSC measurements, in particular, the onset and end of the glass transition region can be identified, along with the middle inflection value assigned to the Tg. Determination of Tg has practical significance especially with regard to polyester polyols, as this value can be compared across resins and related directly to material performance. The value of Tg relative to ambient and processing temperatures strongly affects the final material properties. Several empirical correlation approaches have been developed to predict Tg under different conditions. The Fox−Flory equation19,20 correlates Tg and molecular weight with a parameter related to free volume.21 The Fox equation, on the other hand, often permits accurate prediction of Tg for homogeneous polymer blends and random copolymers on the basis of the mass fractions and Tgs of the constituent species.21 Various other empirical expressions such as the Couchman equation, the Gordon−Taylor equation, and the Kwei equation build on the idea underlying the Fox equation with different functional forms and fit parameters.22,23 An equation proposed by Brostow et al.23 generalizes the Fox equation by extending the linear form of the Fox equation with a polynomial. This equation is effective for modeling binary blends and copolymers with errors generally less than 20 K for blends that deviate significantly from the Fox equation.23 A modified equation introduced by Lu and Weiss24 incorporates an interaction parameter especially relevant for miscible Received: August 10, 2018 Revised: October 23, 2018

A

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

Article

Macromolecules

Several studies have examined modeling Tg through the use of MD simulations.40−51 Because Tg constitutes an example of a second-order phase transition,52 a plot of a first-derivative property (such as specific volume or energy from experimental data) versus temperature exhibits an abrupt change in slope at Tg. Therefore, simulations follow this same methodology by routinely identifying Tg at abrupt changes in the slope of specific volume, density, or energy when plotted against temperature.40−51 This straightforward method involves fitting two lines, one in the lower temperature range and one in the higher temperature range, to locate Tg at the intersection of the two lines. Recently, a protocol for quantifying Tg and its uncertainty from MD simulations has been developed by Patrone et al.53 This method involves fitting a hyperbolic function to density data over the entire temperature range, thereby avoiding the arbitrary delineation of the lower and higher temperature regions, which ultimately affects Tg identification. For this reason, we employ this approach in our analysis of the simulations performed in this work. Previously, Rossi et al.54 have used the MARTINI force field to develop a coarse-grained model of a thermoset coating based on polyester−polyol chemistry. The MARTINI force field is a method designed to build coarse-grained models by mapping several heavy atoms to one coarse-grained bead.55−57 Each coarse-grained bead needs to be specified as an interaction site type with interaction levels previously parametrized within the context of the MARTINI force field.55 Parametrization in MARTINI facilitates reproduction of thermodynamic data such as oil/water partitioning coefficients.55 The polyester polyol modeled by Rossi et al.57 uses different MARTINI bead types to represent parts of the polyester polyol backbone, and bond lengths and angles are matched to atomistic simulations. In their work, the predicted Tg of a thermoset coating formed by reacting a cross-linker with a polyester polyol composed of neopentyl glycol, phthalic acid, and adipic acid lies within 30 K of the experimental Tg value.57 Here, we focus on polyester polyols possessing different chemistries by using a “top-down” coarse-graining approach58 based on the corresponding states correlation26 to the groupcontribution SAFT-γ Mie EoS.27−29 This thermodynamicsbased strategy permits direct calculation of polymer density, Tg and heat capacity. Section 3 provides details regarding experimental methods, while section 4 describes the model and methods used in the computational modeling of four polyester−polyol resins, as well as in exploratory simulations within the Mie model parameter space. Section 5 discusses the results generated from this work and presents a correlation between Mie model parameters and Tg values obtained from the exploratory simulations. Section 6 concludes this work by summarizing the implications of this study. A comparison of the model’s structural properties with experimental data and complementary results from united atom MD simulations are provided in the Supporting Information (section S3).

systems with strong specific interactions. A limitation of these empirical equations is that the Tg of each species comprising a copolymer or blend must be known a priori. Taking a different approach, Wiff and Altieri25 have devised a semiempirical method to predict Tg based on the chemical structure of repeat units in linear polymers, random copolymers, and cured reactive oligomers.25 Using a database of 190 polymers, they have successfully tabulated group-contribution parameters for a wide range of chemical structures to reliably predict Tg. However, their approach requires an experimental polymer database to parametrize the chemical structure-based contributions. Moreover, their study does not include parameters for chemical structures in the polyester polyols considered in this study. Here, we relate the glass transition directly to the monomer chemistry and composition of the polymer using coarse-grained molecular simulations. As a result, we obtain a correlation between the Tg, polymer composition, and thermodynamic parameters of the monomers.

2. THEORETICAL BACKGROUND In this work, we model Tg in molecular simulations with interaction parameters obtained from monomer data. We obtain interaction parameters using the corresponding states correlation26 to the group-contribution SAFT-γ Mie equation of state (EoS),27−29 which is based on the statistical associating fluid theory framework derived from thermodynamic perturbation theory and used in several different forms.30,31 The SAFT-γ Mie EoS describes the thermodynamic properties of particles that are bonded to form chains and that interact with each other via the Mie potential (ϕMie) given by ÅÄÅ ÑÉ ÅÅij σ yz λr ij σ yz λa ÑÑÑ Mie Å ϕ (r ) = CεÅÅjj zz − jj zz ÑÑÑ ÅÅk r { k r { ÑÑÑÖ (1) ÅÇ where r represents the distance between particles, σ and ε are parameters quantifying particle size and interparticle interaction energies, respectively, and the exponents λr and λa are related to the nature of the repulsive and attractive interactions, respectively. The particles are beads that are bonded together to represent an entire molecule. The constant C can be expressed in terms of the exponents as C=

λr ijj λr yzz j z λr − λa jjk λa zz{

λa /(λr − λa)

(2)

Because the SAFT-γ Mie EoS inherently depends on molecular interaction potential parameters, it can be used to determine the potential in eq 1 by fitting experimental thermodynamic properties such as mass density and vapor pressure. Molecules ranging from carbon dioxide,32 water,33 and hydrocarbons34,35 to perfluorohydrocarbons,36 polystyrene,37 and surfactants38 have all been successfully modeled utilizing the SAFT-γ Mie formalism, yielding good descriptions of thermodynamic, interfacial and structural properties.35 The corresponding states correlation26 to the SAFT-γ Mie EoS provides a shortcut to estimate the Mie potential parameters from minimal experimental data. This correlation relates the Mie parameters to the critical temperature (Tc), acentric factor, and liquid density at a reduced temperature of T/Tc = 0.7 (where T denotes absolute temperature) and allows for expedient estimation of interaction potentials for a broad range of molecules.39

3. EXPERIMENTAL SECTION 3.1. Resin Synthesis. We synthesized four representative polyester−polyol resins and used them as the basis for the experimental validation of our coarse-grained models. These polymers were synthesized at Eastman Chemical Co. and used as provided for this study. The monomers incorporated in these resins included neopentyl glycol (N), trimethylolpropane (T), isophthalic acid (I), and 1,4-cyclohexanedicarboxylic acid (C). These monomers, B

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

Article

Macromolecules commonly used in the production of polyester polyols, were selected because of their commercial relevance and the valuable thermomechanical properties they impart to their corresponding polymers.2,59 A chemical depiction of these monomers, as well as the overall polyester polyol polycondensation (step-growth polymerization) synthesis scheme, is presented in Figure 1. The four polyester−polyol resins

polymers by measuring the mass and volume displacement of each resin in water using a Denver Instrument APX-100 mass balance and a 25 mL graduated cylinder.

4. MODELS AND METHODS 4.1. Monomer Data and Mie Parameters. We propose a coarse-grained model specifically tailored for polyester polyols on the basis of their constituent monomers (although the approach could be extended to include other polymers). The data required for the corresponding states correlation26 to the SAFT-γ Mie EoS are available from online databases60−62 for all the diol and diacid monomers considered here. For the triol monomer, we split the trimethylolpropane molecule into one 1-butanol fragment and two methanol fragments, thereby constructing a trifunctional surrogate, as depicted in Figure 2. In this trifunctional model, the masses of the beads are appropriately adjusted to account for the loss of hydrogen in connecting the trimethylolpropane fragments.

Figure 1. Polyester−polyol resin synthesis from hydroxyl or carboxylic acid monomers. Single-letter abbreviations used in this study are displayed below each monomer. are designated by a resin number (PP1 to PP4), followed by the single-letter abbreviations included in Figure 1 of the specific monomers incorporated: PP1NI, PP2NTI, PP3NIC, and PP4NTIC. Table 1 displays the monomer composition and molecular weight distribution parametersnumber-average molecular weight (Mn), weight-average molecular weight (Mw), and z-average molecular weight (Mz)of each experimental polyester−polyol resin. Molecular weight characterization was performed by size exclusion chromatography conducted in tetrahydrofuran (THF) with polystyrene (PS) standards. 3.2. Calorimetry Measurements. We conducted differential scanning calorimetry (DSC) on a TA Instruments Q2000 DSC unit with a liquid nitrogen cooling system to determine the Tg of each resin and obtain the specific heat capacity changes with temperature. The DSC was calibrated with heat flow and cell constant calibration, enthalpy/temperature calibration with indium, and heat capacity calibration with sapphire. We used four different temperature ramp rates in the DSC experiments (5, 10, 15, and 20 K/min) with two DSC trials per ramp rate for each resin and an additional DSC run at 10 K/min for each resin to measure the specific heat capacity over an extended temperature range. The temperature interval examined in the eight DSC experiments ranged from 253 to 393 K, whereas the temperature range in the heat capacity measurements at 10 K/min extended from 173 to 503 K. A heat/cool/heat procedure was followed in each experiment, and the value of each Tg was discerned from the inflection point of the heat flow during the second heating ramp. In addition, the onset and end temperatures of the glass transition, as well as the specific heat capacities, were likewise ascertained from the second heating cycle with the TA Instruments Universal Analysis tool. 3.3. Density Measurements. Because the resins are more dense than water at ambient temperature, we obtained the densities of the

Figure 2. Construction of the coarse-grained trifunctional surrogate for trimethylolpropane: (a) atomistic representation, (b) 1-butanol (yellow) and methanol (purple) fragments, and (c) coarse-grained model.

We calculated the Mie parameters for each molecule using the corresponding states correlation26 in conjunction with the critical temperature Tc, the acentric factor ω, and the liquid density ρ at 0.7Tc as input parameters. The m parameter, which defines the number of beads chosen to represent one molecule, was set to the lowest value allowed by the corresponding states correlation26 given the acentric factor. In using the m parameter to represent a molecule composed of m beads bonded together, there is no specific mapping of each bead to a specific chemical fragment in the molecule. In SAFT, molecules are represented as groupings of different beads without specific reference to connectivity, and therefore the topology of the resulting coarse-grained model does not necessarily correspond to that of the original molecule. This prevents mapping the coarse-grained model back to the original atomistic model without additional information. Empirical data and Mie model parameters for the constituent molecules of the polyester polyols are listed in Table 2. The Mie model parameters were calculated from the empirical data

Table 1. Experimental Compositions and Average Molecular Weights of the Polyester−Polyol Resins Investigated in This Work alcohol content (wt %)

acid content (wt %)

specimen designation

N

T

I

C

Mn (g/mol)

PP1NI PP2NTI PP3NIC PP4NTIC

100 93.3 100 93.35

0 6.7 0 6.65

100 100 50 50

0 0 50 50

1300 1300 1300 1300

C

± ± ± ±

70 70 100 110

Mw (g/mol) 1900 2200 1900 2200

± ± ± ±

30 30 50 50

Mz (g/mol) 2600 3300 2700 3500

± ± ± ±

80 140 150 210

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

Article

Macromolecules Table 2. Thermodynamic Parameters for Each Monomer and the Resulting Mie Parametersa Tc (K) neopentyl glycol isophthalic acid 1,4-cyclohexanedicarboxylic acid 1-butanol methanol

60

677 100760 88961 56360 51360

ω

ρ at 0.7Tc (mol/cm3)

m

ϵ (kcal/mol)

0.00800760 0.00738360 0.00501262 0.09694961 0.02267960

3 5 5 3 3

1.18 1.55 1.34 0.86 0.76

60

0.726 1.06160 1.03661 0.58860 0.56660

(594 (780 (674 (433 (382

K) K) K) K) K)

σ (Å)

λr

λa

3.86 3.50 3.73 3.56 2.67

48.0 49.6 45.7 29.3 27.1

6 6 6 6 6

The parameters ε, σ, λr, and λa specify the Mie interaction potential between beads calculated from ref 26. Superscripts in the table data indicate the source references.

a

using the corresponding states correlation equations developed by Mejiá et al.26 We constructed polymer models by connecting monomer models with harmonic bonds. To account for the loss of water from the reaction between monomers responsible for creating ester bonds, we adjusted the masses of the two beads bonded from different monomers to have 18 g/mol less mass at each ester bond. This allowed the polyester polyol coarse-grained model to match the polymer mass while retaining the properties of the monomers in the force field. 4.2. Molecular Weight Distributions. We considered two molecular weight distributions for each polymer system: a “simple” distribution that is either monodisperse or bidisperse and a “complex” polydisperse distribution optimized to match the first three moments of the experimental molecular weight distributions of the resins. To construct the complex distribution, we employed a genetic algorithm63 with initial guesses from the Schulz−Zimm function.64 The simple distribution was monodisperse for the linear resins (PP1NI and PP3NIC) and bidisperse for the branched resins (PP2NTI and PP4NTIC), whereas the complex model incorporated a polydisperse distribution with up to 100 distinct molecular weights. Histograms showing the number of molecules of each molecular weight in both models are depicted in Figure 3 for the four resins. The black lines represent the simple distributions, whereas the red lines identify the complex distributions. In the simple distributions, the predicted value of M n resided between the M n and M w values of the corresponding experimental resin. The complex distributions matched the experimental values of Mn, Mw, and Mz within experimental uncertainty. Additional details regarding the development of the different distributions, as well as comparisons with experiment, are available in the Supporting Information (section S1). 4.3. Molecular Dynamics Simulations. We modeled the polyester polyol chains as sequences of coarse-grained Mie monomers joined by harmonic bond potentials with a force constant of 300 kcal/(mol Å2) and an equilibrium bond length defined as the σ interaction parameter between neighboring beads. To be consistent with the SAFT theory (which assumes tangential spheres), the bond constant was chosen to be as high as possible while still maintaining energy conservation at the time step size used in the simulations. No angle- or dihedral-bonded potentials were included for consistency with the original SAFT-γ Mie formalism. Nonbonded interactions were considered as pairwise Mie potentials with the parameters obtained from the monomer properties provided in Table 2. These interactions were used for all bead pairs except those directly connected by a bond. Unlike-pair interaction parameters were obtained from the standard mixing rules specified in the corresponding states correlation, viz.26

σAB =

εAB =

1 (σAA + σBB) 2

(3)

εAA εBB

λAB = 3 +

(4)

(λAA − 3)(λBB − 3)

(5)

The λ in eq 5 represents the repulsive exponent λr, while the attractive exponent (λa) was set to 6 in all the simulations. The cutoff radius for the nonbonded potential was set at rc = 20 Å, with a switching function65 S(r) between rl = 18 Å and rc to avoid discontinuities at the cutoff distance. The interaction potential in the range rl < r < rc was weighted so that ϕtotal(r ) = ϕMie(r )S(r )

(6)

where ϕMie(r) is the Mie potential from eq 1 and ÄÅ É ÄÅ É ÄÅ É ÅÅ r − r ÑÑÑ3 ÅÅ r − r ÑÑÑ4 ÅÅ r − r ÑÑÑ5 l Ñ l Ñ l Ñ Å Å Å Å Ñ Å Ñ Å ÑÑ S(r ) = 1 − 10ÅÅ Ñ + 15ÅÅ Ñ − 6ÅÅ ÅÅÇ rc − rl ÑÑÑÖ ÅÅÇ rc − rl ÑÑÑÖ ÅÅÇ rc − rl ÑÑÑÖ (7)

We formulated initial configurations for the simulations by first creating coarse-grained structures for the different polymer molecules. Then, by using the packmol tool,66 we populated the simulation box by specifying the number of each unique molecule to be added to generate the (simple or complex) molecular weight distribution of a particular resin. The data files for the polymer structures used in the simulations are included in the Supporting Information. Molecular dynamics simulations were performed with the large-scale atomic/molecular massively parallel simulator (LAMMPS).67 The time step in all the simulations was 4 fs, which was sufficiently small to ensure energy conservation during trial simulations in the microcanonical (NVE) ensemble at the temperatures studied. Periodic boundary conditions were maintained in all directions. In LAMMPS, all the simulations employed the isobaric−isothermal (NPT) ensemble with the canonical (NVT) ensemble included during equilibration, which incorporated the Nosé−Hoover style thermostat with a time constant of 0.4 ps and the Nosé− Hoover style barostat with a time constant of 4 ps. For each polyester−polyol resin, we conducted five sets of simulations: three replicates of conventional MD simulations with the simple molecular weight distribution, parallel tempering with the same simple distribution, and conventional MD simulations with the complex distribution. Beyond the simulations of the four polyester−polyol resins, 30 additional sets of simulations were completed to explore the Mie model parameter space to discern if the Mie model parameters and polymer composition could be correlated to macroscopic properties, as described in a later section (section 4.3.2). 4.3.1. Polyester−Polyol Resins Simulations. The simulations of the four experimental polyester polyols were D

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

Article

Macromolecules

the Polymatic algorithm,68 which relies on a compressionrelaxation scheme at different temperatures to ensure adequate equilibration. The resulting configuration constituted the starting point for further equilibration and production simulations, as described in the following sections. Conventional Molecular Dynamics. In the conventional MD simulations, after the initial equilibration described above, we cooled the resulting configuration to 200 K at a simulation time rate of 10 K/200 ps and stored intermediate configurations every 25 K. After cooling, we conducted 13 MD simulations in the NPT ensemble at 1 atm and temperatures ranging from 200 to 500 K in steps of 25 K using the cooled states as starting points. Each of these simulations lasted for 10 ns, and properties were averaged over the last 4 ns of simulation time. No significant differences among the average properties calculated over each half of the final 4 ns period were discernible. Parallel Tempering. In parallel tempering, multiple system replicas were simulated simultaneously at different temperatures, with random attempts to swap the temperatures of adjacent simulations.69 The attempts, intended to enhance sampling of conformational space, were accepted on the basis of a Metropolis criterion.70 In our parallel-tempering simulations, the initial configurations were obtained from equilibration at 500 K. We then conducted independent parallel simulations for 2 ns at 64 different temperatures between 200 and 500 K. After this period, parallel tempering continued for an additional 8 ns, over which time switching was attempted every 4 ps (1000 timesteps). The frequency of configuration switching and the temperature spacing were selected to ensure good sampling. The production period was restricted to the last 4 ns of the simulations. More information about the parallel-tempering simulations is available in the Supporting Information (section S2). 4.3.2. Exploring the Mie Parameter Space. Following the study of the four representative polyester polyol resins, we performed simulations spanning the Mie potential parameter space in the vicinity of their original parameters. These simulations did not correspond to any experimental resins. Instead, they were intended to elucidate the effect of each model parameter and identify a correlation to predict the Tg of other similar resins. We created 30 simulations following a design of experiment (DoE) protocol71 using the SAS JMP Pro 12 software package.72 The simulations proceeded with the same methodology as the conventional MD simulations but were limited to a smaller temperature range extending from 200 to 400 K. Table 3 lists the parameters varied in these exploratory simulations as well as the values of each parameter. The values in Table 3 represented reasonable extremes for the parameters based on those corresponding to the four polyester−polyol resins considered in the previous sections. We only considered the diol and diacid molecules in setting the extremes for the interaction parameters, and not the parameters from the branching agent. For each parameter, we considered a small value (colored blue in Table 3) and a large value (colored red). For some of the parameters, depending on the composition of the polymers, an intermediate value (colored yellow) was also included. The specific parameter combinations employed for the 30 exploratory simulations are listed in Table 4. The values of each parameter in Table 4 are indicated by the color-coding introduced in Table 3. A box left blank indicates that the parameter was not relevant.

Figure 3. Predicted simple (black) and complex (red) molecular weight distributions of model resins.

performed over a temperature range spanning from 200 to 500 K. We conducted conventional MD simulations as well as parallel-tempering simulations. For both methods, the model systems were equilibrated to 500 K and 1 atm using energy minimization, followed by the 21-step equilibration protocol of E

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

Article

Macromolecules

simulation was set as close to 5000 as possible. This ensured that different simulations were of similar size despite differences in chain lengths. 4.3.3. Simulated Tg Values. We obtained Tg values from the simulations by fitting a hyperbolic function to the density− temperature curve generated by each set of simulations. The fitting function was given by53

Table 3. List of Parameters Varied in the DoE Simulations and the Values Used

ρ(T ) = ρ0 − a(T − T0) − b /0(T , T0 , c)

(8)

where ρ0, T0, a, b, and c are adjustable parameters, and /0(T , T0 , c) =

1 (T − T0) + 2

(T − T0)2 + ec 4

(9)

where e is the exponential function. In accordance with the recommendation by Patrone et al.,53 the simulated Tg was set to T0, since no significant difference exists between the hyperbola center and the intersection of the tangents of the hyperbola. The glass transition regions in the simulations also followed the convention of Patrone et al.53 This region spanned from T0 − ΔT to T0 + ΔT, with ΔT provided by ΔT =

e c/2(27̂ − 1) 7̂ (1 − 7̂ )

(10)

with 7̂ = 0.9.53 We established that the system size for the temperature range studied did not affect the value of Tg and that the sharpness of the transition made pooling analysis53 unnecessary. 4.4. Specific Heat Capacities. We calculated the specific heat capacity (C p ) in the simulations from enthalpy fluctuations in the NPT ensemble, as indicated by eq 11: Cp =

Our simulations considered the scenarios analogous to the original four resins, with monomer 1 alternating with monomers 2 and 3. The Mie parameters of the difunctional monomers 1, 2, and 3 varied across different simulations on the basis of our DoE. The composition of the polymers in the DoE was selected so that monomer 2 was present at an equal ratio with monomer 3 if monomer 3 was present. The length of a molecule was dictated by the number of monomer 1 units in the chain, with monomer 2 or 3 separating monomer 1 occurrences. Following the procedure for the “simple” molecular weight distribution described earlier, the simulations without branching possessed only one molecular type (linear), while those with branching exhibited two molecular types (one linear and one branched). The branching element was the same trimethylolpropane introduced in our initial polyester polyol models, and only one was permitted in each branched chain, which consisted of three branches of approximately equal length. All the molecules were end-capped by monomer 1 at each terminus. The molecular weights of the chains were established by the values of the molecular weights of each bead and the chemical makeup of the chains. We used the same scheme described in section 4.1 to connect beads of different monomers and account for the mass loss associated with forming ester bonds. In the cases with branching, the ratio of linear molecules to branched molecules was defined by the parameter specifying the number fraction of linear chains. In the DoE, Mn ranged from 201 to 3476 g/mol, and the total number of polymer molecules was scaled so that the overall number of beads in the coarse-grained representation of each

(⟨H2⟩−⟨H ⟩2 )NPT kBT 2

(11)

where H represents the specific enthalpy and kB is the Boltzmann constant.

5. RESULTS AND DISCUSSION 5.1. Polyester−Polyol Resins. 5.1.1. Mass Density. The results provided in Figure 4 display the mass density values from one replicate of the simple distribution simulated by conventional MD, the simple distribution simulated with parallel tempering, and the complex distribution simulated by conventional MD. The two additional simulation replicates of the simple MD for each resin overlap with the first and are omitted from Figure 4 for clarity. Although results from the simple distribution densities do not depend on whether conventional MD or parallel tempering MD was used, the densities predicted from the complex distributions are generally higher than those from the simple distributions. This difference is most clearly seen in the case of PP1NI, wherein the density from the complex distribution is significantly higher than that from the simple distributions. These higher densities are attributed to improved packing of polydisperse particles relative to monodisperse particles, previously established for hard spheres.73 In our polymer models, bonded neighbors pack closer than nonbonded neighbors based on the equilibrium bond distance (σ) used and the Mie potential form. Longer linear molecules will have relatively more bonded neighbors than shorter molecules, leading to closer packing with the inclusion of longer F

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

Article

Macromolecules Table 4. Parameter Combinations Used in the DoE Studya

a

The columns indicate the parameters that were varied; the colored boxes represent the values of the parameters used for each simulation according to the color scheme introduced in Table 3.

densities without any additional need to adjust bond lengths. This result may be attributed to the similarity between monomer and polymer densities at the same temperatures. 5.1.2. Glass Transition Behavior. The glass transition temperatures obtained from the DSC experiments and the five simulations/resin are provided in Figure 5, and their values are listed in Table 6. The average values are included in Figure 5 for direct comparison. In calculating the average Tg values, data from the different cooling rates in the experiments are considered together, as they are fairly similar. In addition to the average Tg and the average width, the minimum onsets and maximum end temperatures of the glass transition region across all the simulations and experiments are identified in Figure 5. As in the case of the mass density, the simulation and experimental results for Tg are in reasonably good agreement, with the transition regions overlapping in all cases. The

molecules. Perhaps the effect of polydispersity is exaggerated in the linear PP1NI resin complex model because the distribution, presented in Figure 3a, contains many more linear molecules with molecular weights between 3000 and 5500 g/mol than the other resins. This skewness of the molecular weight distribution might affect molecular packing. Included in Figure 4 are the experimental density values measured at ambient temperature (296 K). Close examination confirms favorable agreement between the experimental and simulation densities, with the exception of PP1NI complex MD, as discussed above. In Table 5, we list the values of the experimental densities and selected simulation densities. The standard errors in the simulations are on the order of 10−5 g/ mL and are omitted from the table. The agreement between simulated and experimental densities suggests that, for polyester−polyols, using the monomer density to set the sigma (σ) parameter is sufficient to reproduce polymer G

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

Article

Macromolecules

Figure 4. Temperature dependence of mass density obtained from simulations and experimental measurements of each resin. Values of Tg for selected simulations are also included. Figure 5. Tg values and the glass transition regions from DSC and simulations of four polyester−polyol resins. Stars identify experimental Tg values, whereas orange lines correspond to the experimental glass transition regions. Triangles indicate Tg values from simulations and the corresponding glass transition regions are identified by blue bars. Summary values from experiments and simulations (expressed as averages) are also included with horizontal lines indicating Tg region extremes.

Table 5. Experimental and Simulation Densities near Ambient Temperature and 1 atm for the Four Polyester− Polyol Resins Studied resin designation

simple model density (g/mL)

complex model density (g/mL)

PP1NI PP2NTI PP3NIC PP4NTIC

1.18 1.17 1.07 1.06

1.27 1.18 1.08 1.08

exptl density (g/mL) 1.16 1.18 1.10 1.11

± ± ± ±

0.02 0.02 0.03 0.03

ration of the branching monomer, trimethylolpropane, has little, if any, discernible effect on Tg at the concentration used. The simulated Tg values are consistent across the five simulation sets per resin. The spread in the Tg values within one resin is