Off-Lattice Flory–Huggins Approximations for the Tailored Calculation

Dec 22, 2016 - Citing Articles; Related Content. Citation ... For a more comprehensive list of citations to this article, users are encouraged to perf...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/IECR

Off-Lattice Flory−Huggins Approximations for the Tailored Calculation of Activity Coefficients of Organic Solutes in Random and Block Copolymers Phuong-Mai Nguyen,†,‡ Wafa Guiga,†,‡,§ Ahmed Dkhissi, and Olivier Vitrac*,†,‡ †

INRA, and ‡AgroParisTech, UMR 1145 Ingénierie Procédés Alimentaires, F-91300 Massy, France CNAM, UMR1145 Ingénierie Procédés Aliments, F-75003 Paris, France

§

S Supporting Information *

ABSTRACT: Predicting activity coefficients in polymers for arbitrary solute and polymers is of general interest for packaging, environmental, and membrane applications. We examined the extensions of the proposed off-lattice approach [Gillet et al. Ind. Eng. Chem.Res. 2009, 2010] for random and block copolymers. Based on a comparison with an explicit representation of random copolymers, the principles of a meanfield approximation using only the properties of homopolymers have been established. The approach was successfully validated against the Flory−Huggins coefficients collected by Fornasiero et al. [AIChE J. 2002] for 19 aromatic solutes in ethylene-vinyl acetate with acetylation rates varying from 33% to 100%. The role of substituents and isomerism is elucidated. Previously collected data and this study suggest that the real chemical affinity of substituted aromatic compounds for water could be strongly underestimated by previous oversimplified rules and that tailored methods, such as this one, would be preferable.

1. INTRODUCTION Thermodynamic data such as solubility, activity, or partition coefficients of organic substances are essential for assessing bioaccumulation and the contamination of food and of the environment, etc. In polymers, data have been tabulated mainly as solubility coefficients and similar cohesion parameters in several monographies.1−4 For nonvolatile or nonsolvent substances such as plastic additives the data are particularly scarce. Indirect methods involving residual concentration measurements in polymers or inverse gas chromatography are usually required. Predicting excess chemical potentials from molecular simulations and ab initio force fields offer attractive alternatives. Brute force calculations involve the partial difference of the free energy respectively to the number of solute particles.5,6 At infinite dilution, the statistics is poor and the difference needs to be repeated many times. Over the last years, several improved calculation strategies have been proposed to transpose the methods to large bulky solutes in polymers. As almost no void can accommodate large solutes, Widom’s insertion trials7 have to be replaced alternatively by a deletion method8−10 or by thermodynamic integration along a feasible path.11−13 Boulougouris14 recently extended common integration schemes by combining free energy perturbation over the number of molecules and thermodynamic integration over the thermodynamic state. Soft-core state reference ensemble15 offers another alternative for hindered solutes with almost spherical symmetry. All methods are, however, © 2016 American Chemical Society

insufficiently tailored for large sets of solutes and/or polymers. With this respect, off-lattice Flory−Huggins calculations at the atomistic scale offer more direct estimations.16−18 The derivations are based on a compressible version of the original mean-field Flory−Huggins lattice model,19,20 initially proposed by Bawendi and Freed.21 The main advantage is that they do not require any explicit representation of chain entanglements but only of polymer segments. The only cost is the optimal segment length, which maximizes the surface area of interactions with the considered solute, is not known a priori and needs to be determined by trial and error. As a substantial advantage, the athermal estimation of contact energies enables the reuse of sampled distributions to extrapolate the activity coefficient at slightly different temperatures (e.g., ±30 °C). However, the whole methodology suffers the same intrinsic limitations as the original Flory−Huggins theory.22,23 Thermodynamic effects are reduced to random pair contacts, which are assumed to be independent, not cooperative, and not separated by significant energetic barriers. Strong dissimilarities in size, shape, and flexibility as well as chiral centers along the chain reduce dramatically convergence rates (see ref 18 for details). Received: Revised: Accepted: Published: 774

September 22, 2016 December 20, 2016 December 22, 2016 December 22, 2016 DOI: 10.1021/acs.iecr.6b03683 Ind. Eng. Chem. Res. 2017, 56, 774−787

Article

Industrial & Engineering Chemistry Research

dispersed among two different polymer segments so-called A and B (i.e., ternary mixture i + A + B) with respective volume fractions {ϕk}k=A,B, FHT can be used to assess whether the solute prefers contacts with A or with B, or whether it can reduce the number of contacts between A and B. The standard theory assumes that all blobs are perfectly mixed and that the resulting mixture is incompressible. All blobs are commensurable and are constituted exclusively of one single species: i, A, or B. As demonstrated in eq (5) of the Supporting Information, the effective Flory−Huggins coefficient, χiPcopolymer associated with the volume fraction in copolymer ϕP = ϕA + ϕB is

The present work aims at extending the estimation approach of bulky solutes from apolar homopolymers16,17 to polar ones and to copolymers alternating apolar and polar segments with periodicity shorter than their persistence length. These effects cannot be captured by the mean-field lattice model without an additional assumption about the pair contact mixing rule, but they can be tested directly at the atomistic scale. In the next section, the main elements of the off-lattice generalized Flory− Huggins theory are presented for multicomponent mixtures as well as the set of reference data. Flory−Huggins coefficients inferred by Fornasiero et al.24 from experimental distribution coefficients between water and three ethylene vinyl acetates (EVA)25 were chosen as validation sets. Besides their analogies with monomers and additives used in printing inks, the 19 considered aromatic solutes are nonvolatile and capable of creating multiple weak hydrogen bonds with EVA. Section 3 details the simulation and the computation procedures for both homo and copolymers. The results obtained on homopolymers (polyethylene and polyvinyl acetate) and heteropolymers (EVA with arbitrary acetylation rates) are summarized in section 4. The comparison of results obtained on several isomers is particularly used to assess how multiple hydrogen bonds are seen at microscopic scale and in mean field theory. It is shown that, contrary to the assumption of Scatchard,26 the interactions between the nonelectrolyte and polymer segments are not independent of the composition. Section 5 summarizes our findings and recommendations to calculate activity coefficients and contains some comments important for the safety of materials intended to be in contact with food.

χicopolymer = ,P

hjk = ⟨εj + k⟩T ⟨zj + k⟩ + βjk +∞

where ⟨εj + k⟩T =

⎛ ⎛ r ⎞ r⎞ = ⎜1 − i ⎟ϕA + ⎜1 − i ⎟ϕB + χicopolymer ϕP 2 ,P RT rA ⎠ rB ⎠ ⎝ ⎝ ex

⎛ r⎞ = ⎜1 − i ⎟ϕP + χic, Popolymer ϕP 2 rP ⎠ ⎝

where ri is the partial molar of solute i; {rj}k=A,B and {ϕk}k=A,B are respectively the partial volumes and the volume fractions associated with segments A and B. The effective partial volume

(

ϕA

+

rA χcopolymer i,P

ϕB rB

−1

)

(3) ε

( ) ) exp(− ) dε

j+k ∫−∞ pr(εj + k) exp − k T ε dε +∞

∫−∞ pr(εj + k

B εj + k k BT

is the temperature

ensemble averaged contact energy inferred from the distribution of pairwise contact energies, pr(ε), with kB being the Boltzmann constant, between j and k; ⟨zj+k⟩ is the average number of test molecules k around the seed molecule j, and βjk is the expected covariance between ⟨εj+k⟩T and ⟨zj+k⟩. It is worth noticing that ⟨εj+k⟩T = ⟨εk+j⟩T, whereas ⟨zj+k⟩ ≠ ⟨zk+j⟩ for k ≠ j. The binary Flory−Huggins coefficient χjk is defined accordingly as the dimensionless excess mixing enthalpy:

(1)

of the copolymer AB is given by rP = (ϕA + ϕB)

(2)

Equation 2 leads to several paradoxes already mentioned by Krause:31 (i) the effective solute−polymer coefficient does not depend on segmental composition (i.e., different values of ϕA and ϕB would lead to similar results); (ii) it cannot depend on more than two segment types; (iii) quadratic corrections (see eq (4) in the Supporting Information) do not have physical meaning. 2.3. Microscopic Averaging of Pair Interactions for Polynary Mixtures. Equation 2 considers that only energies per volume (assuming uniform contact types) have to be averaged. By following the same idea of Simha32 to reach a theory of chain copolymerization reactions (see discussion in note 18 of ref 32), we will assume alternatively that contact energies have to be averaged at the microscopic scale. In this description, each solute i is assumed to be at infinite dilution and surrounded by many types of molecules rather than by a preferential species. By following the pairwise decomposition of contact energies at the molecular scale proposed in refs 16−18, the enthalpy of interaction between two species j and k (being either a solute or a polymer segment), hjk, is defined as

2. THEORY AND INTERPRETATION OF REFERENCE DATA 2.1. Effective Flory−Huggins Coefficients for Copolymers. As shown in the Supporting Information, effective Flory−Huggins coefficients for copolymers, denoted χcopolymer , i,P enable to consider solute-copolymer mixtures as binary systems. By replacing any binary copolymer AB by an effective homopolymer, denoted P. The expression of the solute excess chemical potential becomes accordingly (see eq 2 in the Supporting Information): Δμi

1 (χ + χi B − χAB ) 2 iA

.

χjk (T ) =

can be The following subsections will review how theoretically and experimentally assessed. It is shown, in particular, that mean field lattice approximations and local descriptions of pair contacts lead to dissimilar conclusions for random copolymers. 2.2. Mean-Field Lattice Approximation in Perfectly Mixed Ternary Systems. The classical Flory−Huggins theory (FHT)27−30 is the first mean field approximation offering a reasonable estimation of the trade-off between attractive forces generated by the contact of united atom groups and the entropy of mixing. In the special case where a solute i is

hjk (T ) + hkj(T ) − hjj(T ) − hkk (T ) 2kBT

(4)

Equation 4 assumes random contacts between the molecule/ segment j and the molecule/segment k. The next section discusses the generalization of eq 4 when the solute i is surrounded by polymer segments of ν types. The number of effective pair contacts between i and k, {⟨zi+k⟩}k=1,...,ν, must be weighted accordingly by the probabilities to find simultaneously a segment of kth type in the neighborhood and in contact with solute i. The latter are denoted {pij}j=i,1,...,ν or {pji}j=i,1,...,ν according to whether the solute i or the molecule/ segment j is chosen as seed. At thermodynamic equilibrium, 775

DOI: 10.1021/acs.iecr.6b03683 Ind. Eng. Chem. Res. 2017, 56, 774−787

Article

Industrial & Engineering Chemistry Research

At high acetylation rates, the amount of absorbed water was larger than the amount of absorbed solutes so that a ternary mixture between P, i, and water was required. The authors consider that the interactions between the solute i, the copolymer PVAc, and water, denoted with subscript w, obeyed a microscopic model similar to eq 5. The effective Flory− Huggins coefficient was expressed accordingly as r r χieff = χicopolymer ϕP + i χi ,w ϕw − i χw,copolymer ϕw ϕP ,swollen P ,P rw rw P

they are both equal. The effective excess enthalpy accounts for all possibilities of pair contacts between i and k = 1,...,ν segments but also on possible pair contacts between segments, when they are in contact with i: χicopolymer = ,P



pii pik hik + pii pki hki − pii pik hii − pii pki hkk (pii + pii )kBT

k = 1...ν

=

∑ k = 1...ν

pik χi , k −

∑ j = 1,...,ν k = 1,...,ν j nopt i,P , χi,P estimates converged to a similar value, but uncertainty was much higher due to the difficulty to sample “good” contacts between i and P as well as between P and P. The effect of solute size on nopt i,P values is reported in Figure 3. An apparent correlation emerged

reference χi,P values could not be directly measured and required a thermodynamical interpretation of apparent partition coefficients between PVAc and water. Because of the absorption of water by PVAc, three pair interactions needed to be considered: i+PVAc, i+water, and PVAc+water in the microscopic model (6). Finally, as the experimental temperature was below the melting point of 11 of tested solutes, fugacity corrections between solid and molten states were applied to enable comparison with our amorphous reference states. The corrections are plotted as vertical downward arrows. Without applying any fitting procedure, calculated and experimentally derived χi,PVAc values fall globally along the first bisector. Fugacity corrections improved dramatically the predictions for nitrophenol and nitroaniline groups. Predictions reproduce the observed two classes of mixing properties: solutes associated with exothermic mixing with PVAc (negative net heat of sorption: χi,P < 0) and those with endothermic mixing (χi,P > 0). In addition, calculations at the atomistic scale confirmed that isomers did not share the same behavior. This ability motivates the use of atomistic calculations and generic molecular force fields. Alternative methods to represent nonideal mixtures such as group contribution methods and among them, self-associating fluid theory (SAFT), require important sophistications (e.g., explicit identification of interaction sites) and proper fitting against experimental or simulation results. The applicability of SAFT methods to

Figure 3. Correlation between PE chain lengths nPE and van der Waals volumes of solutes. van der Waals volumes are obtained by rolling a ball of 0.1 Å on the atom surface of 500 conformers. The internal volume to the contact surface with the ball was calculated subsequently by tessalation with a step of 0.1 Å. Symbols indicate solute melting temperature range: ▲, Ti,m ≤ 25 °C, ■, 25 °C ≤ Ti,m ≤ 100 °C and ●, Ti,m ≥ 100 °C. 780

DOI: 10.1021/acs.iecr.6b03683 Ind. Eng. Chem. Res. 2017, 56, 774−787

Article

Industrial & Engineering Chemistry Research

between water and PVAc. With low melting points (193 and 230 K, respectively) and relatively high boiling points (518 and 388 K), both solutes fit well with the Flory description at the atomistic scale. The good accuracy of their cohesive energy determinations confirmed this feature (see Figure 1). Though all 19 solutes have a higher chemical affinity for PVAc than for water, benzophenone along with naphthols has the highest (partition coefficient of 1782). Fornasiero et al.25 reported, however, a χi,PVAc value of 0.75, which makes benzophenone not fully soluble in dry PVAc. Benzophenone presented several difficulties not only for the interpretation of experimental data also for molecular modeling in the Flory approximation. It was the largest studied solute with a melting point of 320 K close to the set temperature of 298 K. It required also the largest number of VAc monomers to enable the sampling of pair contact energies (nopt i,PVAc = 5). According to Table 5 in ref 25, it is the only solute for which the experimental chemical affinity for the polymer was not monotonously decreasing with acetylation rate. Theory (see eq 5, section 4.4, and Figure 7) predicts conversely a minimum with acetylation rate but not a maximum as reported in experimental partitioning values. 4.4. χi,P Values for EVAc. In this study, χi,P values for EVAc were calculated according to two formulations: from a meanfield model (macroscopic) reusing χi,P values previously calculated for homo PE and PVAc and from a very computationally expensive microscopic (atomistic) model. The interests of the explicit representation of the random copolymer are multiple. On the one hand, the proposed microscopic model can validate or invalidate eq 5 without using experimental data and it can give further insights on the nonlinear behavior of energy of mixing based on atomistic details (random sequence and tacticity of VAc monomers). On the other hand, only the microscopic model can estimate the spreading of the mixing energy based on the variations of the arrangement of EVAc segments. The macroscopic approach estimates, indeed, only the first moment (average) of the energy of mixing. Macroscopic Model. At infinite solute dilution, pair contacts were weighted in eq 5 by the volume fraction of vinyl acetate segments (by assuming the mixture incompressible), defined as

copolymers is discussed in ref 58 and recent advances to describe hydrogen-bonding in interacting liquids are presented in ref 59. Prediction errors tended to be lower for endothermic mixtures and larger for exothermic ones. The deviations observed in exothermic mixtures had two possible origins: uncertainty in reference χi,PVAc coefficients (highly dependent on the robustness of analyses made by Fornasiero et al.25) and possible sampling biases in pairs, where electrostatic interactions (dipole−dipole and weak hydrogen bonds) were dominating. Inferring reference χi,P values from the partitioning between a ternary mixture (i+PVAc+water) and a binary one (PVAc+water) is particularly challenging as water is dominating most of the interactions with polar solutes (see eq 6). Equilibrated PVAc absorbs indeed 7% w/w of water. From the computational point of view followed in this work, as electrostatic interactions were randomly sampled without preferential arrangement or consideration of their lifetimes, some sampling biases or errors could be expected. The Boltzmann averaging procedure guarantees the absence of significant biases of pair interaction energies in amorphous phases. The convergence was, nevertheless, slow and required intensive sampling. It has been shown that increasing the number of samples, conformers, or stereoisomers did not reduce the deviation to experimental data. Nonreducible deviations were consequently thought to be intrinsic to the formulation itself or to its approximation. As an example, the Flory approximation can generate significant biases: all states are assumed accessible and not separated by significant energy barriers. This assumption is valid for polymers at rubber state but it requires corrections at the glassy state (see the discussion in ref 60). Dry PVAc as considered for χi,PVAc values reported in this work is glassy at room temperature, whereas PVAc used in the experiments of Fornasiero et al.25 was significantly plasticized by water. Similarly, PVAc was reported to be experimentally oversaturated in water with an excess of 3% w/w due to “additional sorbed water [that] fills microscopic cavities in the polymer matrix”. On the contrary, Flory approximation associated with eq 6 assumes that water was perfectly mixed with polymer segments (no water cluster). Finally, surface absorption was estimated about 0.5% of the total sorption (estimations for acetophenone, 2-nitroaniline, 2-nitrophenol according to Fornasiero et al.24) but it was neglected in our calculations. The deviations are discussed with more details for the two solutes leading apparently to the largest prediction errors: nicotine and benzophenone. Nicotine is fully soluble in water and might be also distributed between PVAc and possible water clusters. Due to the experimental procedure (initial dissolution in water), the original partition coefficients estimated by Fornasiero et al.24 for nicotine were subjected to an error up to 50% and should be considered with caution as mentioned by the authors themselves. As pyridine and benzoic acid, nicotine was, indeed, partly dissociated in the conditions of the experiment (1%, 100%, and 6%, respectively). These effects were not considered in our simulations, where all atoms were at their fundamental state. For all these mentioned reasons, the χi,PVAc value of −4.27 for nicotine looks particularly low. Our calculations predicted also a negative energy of mixing but with a higher value of −0.87, which looks more realistic. In particular, nicotine and pyridine were accordingly estimated with very similar χi,PVAc values, in good agreement with their similar partitioning

ϕVA (wVA ) =

wVA /ρPVAc wVA /ρPVAc + (1 − wVA )/ρPE

(11)

where wVA is the weight fraction in vinyl acetate (VAc) monomer; ρPVAc and ρPE are the densities of amorphous PVAc (at rubbery state) and amorphous PE. The following values were retained for this study: 1189 kg·m−3 and 850 kg·m−3, respectively. The Flory−Huggins coefficient between PVAc and PE segments, χPVAc,PE, was determined similarly as for solutes using the values of nP, which were optimal for each considered solute. This approach led to an average χPVAc,PE value of 7 ± 1. The data of the macroscopic model were compared crudely with the values of χi,EVAc estimated by Fornasiero et al.25 for acetylation rates of 45% w/w (denoted EVAc45) and 33% w/w (denoted EVAc33). The comparisons are plotted in Figures 5 and 6, respectively. Since our calculations were not independent of calculations plotted in Figure 4 with full acetylation, the same biases were expected to be present and possibly amplified. From the experimental point of view, reference values χi,EVAc were easier to extract by Fornasiero et al.25 In particular, equilibria were faster to achieve, and they 781

DOI: 10.1021/acs.iecr.6b03683 Ind. Eng. Chem. Res. 2017, 56, 774−787

Article

Industrial & Engineering Chemistry Research

substituted ones. Without hydroxyl moiety, the predictions of nitrobenzene were almost without significant error for all acetylation rates. Polar moieties seemed to cause the largest source of deviation but the picture was more complex because all considered solutes included a polar functional group: amino, carbonyl, carboxyl, hydroxyl, nitro. Independent studies on the sorption of phenolic compounds61 and of monosubstituted benzoic acids62 invoked a possible complexation between hydroxyl groups and acetate groups. The main argument was the nonlinear dependence of water−PVAc partition coefficients with the vinyl acetate content. This effect was so specific that it was observed with phenolic compounds but not with 1-chloro4-nitrobenzene.61 Nonlinear behaviors cannot be captured with the mean-field eq 5, which assumes that all pair-contacts are equiprobable. Recent examination of vapor pressures of various substituted benzoic acids63 suggested that steric effects would be dominant in ortho positions, whereas resonance effects would be prevailing for meta and, to a lesser extent, para positions. Microscopic Model. The microscopic model was used to establish the nature of interactions and their impact on χi,EVAc values when the random segments of vinyl acetate and of ethylene are commensurable to solute size. For example, the approach makes it possible to compare repeated patterns with similar acetylation rate but with different lengths for polar/ apolar segments, such as AAEA, EAAA (E = ethylene segment, A = vinyl acetate segment). Without homogenization, there is no a priori approximation on the nature and location of preferred contacts. They are all sampled equivalently and weighted according to the value of their respective contact energy via the Boltzmann factor. In agreement with results depicted in Figures 3 and 4, the following series were studied: nitrophenols, nitrotoluenes, acetophenone, benzonitrile, benzyl alcohol, nitrobenzene, and benzophenone (11 solutes). These rather hydrophobic substances presented all possible variations on substituents and on their positions. The simulations were carried-out on pseudo-oligomers (head and tail atoms not in contact) for all possible combinations of E and A. The χi,EVAc values for the 11 solutes are depicted in Figure 7 as a function of the vinyl acetate content and compared with their homogenized estimates from eq 5. Microscopic calculations were highly consistent with their mean-field approximations. The values were only slightly above the homogenized estimates. The high value of χPVAc,PE was responsible for the negative curvature of the mean-field model. This property was reproduced at microscopic scale with minima of χi,EVAc occurring near a common acetylation rate of 0.7−0.8. Only 3-nitrophenol, 4-nitrophenol, and benzophenone did not present any local minimum. Their mixtures were all exothermic in pure PVAc. Sampling of electrostatic and polar interactions with vinyl acetate groups required intensive sampling over large sets of stereoisomers and conformers. Pseudo-oligomers enabled only discrete acetylation rates, and oligomers lengths were needed to be varied around the optimal nopt value to reproduce a P pseudocontinuous acetylation rate scale. For the largest solute, benzophenone, the results were dramatically influenced by solute size and illustrated the general difficulty. The optimal nopt P value should be adjusted to assess the interactions with copolymers presenting monomers of different sizes. This effect is exemplified in Figure 8 by plotting the ratio of the van der Waals surfaces between benzophenone and the tested pseudooligomer. Keeping this ratio constant imposed to increment np,

Figure 5. Comparison of χi,EVAc45 predictions from eq 6 with values estimated from experiments of Fornasiero et al.25 The continuous line represents the first bisector. Symbols indicate solute melting temperature range: ▲, Ti,m ≤ 25 °C; ■, 25 °C ≤ Ti,m ≤ 100 °C; and ●, Ti,m ≥ 100 °C.

Figure 6. Comparison of χi,EVAc33 predictions from eq 6 with values estimated from experiments of Fornasiero et al.25 Symbols indicate solute melting temperature range: ▲, Ti,m ≤ 25 °C; ■, 25 °C ≤ Ti,m ≤ 100 °C; and ●, Ti,m ≥ 100 °C.

were less hindered by the mutual sorption of water and solute. Water sorption was found accordingly about 0.45% in EVAc45 and almost negligible in EVAc33. The same issues as reported in Figure 4 for PVAc were also apparent for EVAc45 (Figure 5) and in a lesser extent for EVAc33 (Figure 6). The macroscopic model (eq 5) showed that dominating electrostatic interactions in VAc rich-segments were mitigated by van der Waals ones in ethylene groups. Most of χi,EVAc were thus positive or close to 0 in EVAc33. Correcting fugacity effects due to solid reference states significantly improved predictions, with values that coincided almost with experimental ones. The final prediction errors were lower in EVAc33 than in PVAc. The largest error was for 3-nitrophenol. More globally, the prediction errors were higher for isomers of nitrophenols and nitrotoluene substituted in meta and para positions, whereas they were fully acceptable for ortho782

DOI: 10.1021/acs.iecr.6b03683 Ind. Eng. Chem. Res. 2017, 56, 774−787

Article

Industrial & Engineering Chemistry Research

Figure 7. Comparisons of χi,PVAc values calculated from mean-field approximation (6) (continuous line) with the microscopic approach representing explicitly the copolymer at prescribed acetylation rates ωVA. The effects of segment length nP are also indicated.

theory supporting eq 5 predicts this trend when the differences in χi,P values are low between the two homopolymers. Experimental partition coefficients between EVAc and water (see Table 5 of Fornasiero et al.25) exhibited the highest chemical affinity for the polymer one and the lowest relative dependence with acetylation rate. The partition coefficient with EVAc33 and PVAc was estimated to only 1.3. By expressing activity coefficients as exp (1 + χi,P) (see ref 16 for justification), microscopic results give partition coefficients close to 3. Configurations of minimum energies (the likeliest ones) brought further insight of electrostatic interactions between aromatic solutes and vinyl acetate groups. For a same acetylation rate, the number of hydrogen-bonds (see classification in ref 64) varied according to the distance and the accessibility of hydrogen acceptors. For solutes of similar size, the number of hydrogen bonds significantly affected pair contact energies, in return. These effects are illustrated in Figure 9 on 3-nitrophenol, benzonitrile, and nitrobenzene. The role of isomerism was particularly identifiable on nitrophenols. Meta and para positions of hydrogen donor groups reduced intramolecular interactions while increasing the probability of intermolecular H-bonds. Pair contact energies with random copolymers were maximized when the vinyl acetate sequence was matching the distance between hydrogen donor groups. As up to simultaneous five weak-to-moderate H-bonds were identified, the complexation proposed in ref 62 made sense. The either linear or quadratic dependence with acetylation rate, first reported in ref 61, could be interpreted in the light of the number of H-bonds created between aromatic solutes and EVAc. The relationship is linear when one single H-bond (e.g., 2-nitrophenol) is created on average, and quadratric when this number is close to 2 (e.g., 3- and 4-nitrophenols). 4.5. General Comment for the Estimation of Partition Coefficients between Polyolefins and Water. Fornasiero et al.24 reported partition coefficients between EVAc and water for the 19 studied solutes but did not report any value for the most hydrophobic polymer: polyethylene (i.e., with zero

Figure 8. van der Waals surface ratios between segments of PVAc and benzophenone. The results are plotted as a function of segment length nP and acetylation rate ωVA.

when acetylation rate was increased from pure PE to pure PVAc. Consistent with results presented in Figure 2, a ratio of 2 was chosen. This value maximizes the contact surface between an essentially planar solute and a linear polymer segment shorter than its persistence length (i.e., only one-half of the solute surface can at maximum be covered by a rigid linear segment, which cannot be fully folded). In practice, np values lower than optimal ones led to poor convergence and to an artificial spread of χi,EVAc values. For benzophenone, nopt P values laid between 2 and 6 for acetylation rates varying from 0 and 60%. It increased dramatically beyond with the fractal dimension of polymer segments. By following these rules the estimated Flory−Huggins coefficient of benzophenone presented pseudoplateau with acetylation rate. The mean-field 783

DOI: 10.1021/acs.iecr.6b03683 Ind. Eng. Chem. Res. 2017, 56, 774−787

Article

Industrial & Engineering Chemistry Research

Figure 10. Partition coefficients K̂ i,PE/water at 298 K estimated from eq 12 (filled symbols) vs their theoretical solubility limit (see text). Experimental values are plotted as empty diamond symbols. Lower bounds corresponding to EU rules for fatty (Ki,PE/water = 1) and aqueous (Ki,PE/water = 1000) contacts are plotted as horizontal lines. The grayed area corresponded to the domain min(103, K̂ i,PE/water ≈ − ln 1/γvi,w − χi,PE − 1) with 0 ≤ χi,PE ≤ 7.05. Symbols indicate solute melting temperature range: ▲, Ti,m ≤ 25 °C; ■, 25 °C ≤ Ti,m ≤ 100 °C; and ●, Ti,m ≥ 100 °C.

Figure 9. Configurations of minimum energies of three typical solutes with two EVAc sequences of length n = 5 and presenting the same acetylation rate ωVA ≈ 0.82. Hydrogen bonds are plotted as dashed lines.

103. Only benzophenone presented a value greater than 100, which could be close to the suggested lower bound. All others fell far below with values even lower than unity for nitroanilines and nitrophenols substituted in meta and para positions, as well as for very polar nicotine and pyridine. K̂ i,PE/water values distributed along a regression line with a general equation ln v K̂ likely i,PE/water ≈ − ln 1/γi,w − 3.86 and corresponding to a likely χi,PE value of 2.86, denoted χlikely i,PE . Since χi,PE values were ranged between 1.22 to 7.05 and since no significant negative χi,PE values are expected in an apolar polymer, a robust upper bound of Ki,PE/water is expected by choosing χi,PE = 0 in eq 12: K̂ i,PE/water < min (103, γvi,w/exp(1)). This upper bound would give the minimum concentration in water after a prolonged contact between amorphous PE and water. For a semicrystalline polyethylene (with c usually ranged between 25% and 90%), the upper bound becomes min(103, (1 − c)γvi,w/exp(1)). The lower bound is more difficult to establish as there is a balance between solubility in water and the enthalpy of mixing with PE segments. One possible choice is to consider that solutes are ideal in water (γvi,w = 1), that is K̂ i,PE/water ≈ exp(−3.86). Among the tested solutes, only pyridine behaves as an ideal solute. This approximation fails to encompass 3-nitrophenol, which remains outside bonds because of its very high χi,PE value of 7.05. From our set of values, the optimal lower bound was: K̂ i,PE/water > γvi,w/ exp(8.05) ((1 − c)γvi,w/exp(8.05) if crystallinity effects are included). The likely interval within suggested lower and upper bounds is depicted as gray surface in Figure 10. To assess the reliability of our estimates and safety margins, the experimental values of partition coefficients of three solutes (acetophenone, benzonitrile, and benzoic acid) are also depicted. To ensure that benzoic acid was not dissociated during the experiment, the pH was decreased down to 1.9 (significantly below its pKa of 4.1) by adding hydrochloric acid. Experimental values confirmed very well the magnitude orders of the predictions. Average prediction errors were about 30%. Despite their low solubility

acetylation rate). Polyethylene−water partitioning is of general interest for many packaging applications and environmental issues. This highly formulated polymer is a possible source of many contaminants, including phenolic antioxidants, ultraviolet absorbers such as benzophenones, and resin precursors such as acetophenone used in adhesives or printing inks. Very few partition coefficients Ki,P/F have been reported between plastics and water. As quoted by Landra and Satyro65 “the amount of hydrocarbons in the water phase is usually severely underestimated”. For the assessment of the safety of food contact materials, the recent practical guide for food contact materials published by the European Commission66 suggests to use a lower limit value of 1000 for all lipophilic substances distributed between any thermoplastic and aqueous food or food simulant. This lower bound has been proposed as a worst-case value to overestimate the relative affinity of plastic additives and monomers for water. Among the 19 substances, studied by Fornasiero et al.,24 distributed between EVAc33 (the least polar polymer tested by authors) and water, only one single substance falls into this category: benzophenone with a partition coefficient of 1356 ± 30. All other substances have partition coefficients of much lower values, ranging from 1 to 599. Comparatively, polyethylene can offer reliable estimates for the most apolar material. Partition coefficients between amorphous PE and water, K̂ i,PE/water, were estimated at infinite dilution as ln K̂ i ,PE/water ≈ ln γiv, w − χi ,PE − 1

(12)

where χi,PE values were calculated by this work (see Table 3) and solute activity coefficients related to volume fractions, γvi,w, were taken from Table 3 of Fornasiero et al.25 The results are plotted in Figure 10 versus solubility limit, idealized as 1/γvi,w. This approximation of the solubility limit assumes a reversible linear isotherm in water. All 19 partition coefficients K̂ i,PE/water were severely overestimated by the recommended bound of 784

DOI: 10.1021/acs.iecr.6b03683 Ind. Eng. Chem. Res. 2017, 56, 774−787

Article

Industrial & Engineering Chemistry Research

liquid at the considered temperature. As a result, the main deviations were observed for polar solutes distributed between water and polar polymers. Despite these limitations, both experimental determinations and calculations demonstrated that the community studying packaging materials tends to underestimate the real chemical affinity of substituted aromatic solutes for water. Distribution of such solutes between water and the polymer is not as favorable to the polymer, as it is usually thought, in particular, to estimate migration issues. The chemical affinity for the polymer depends not only on the number and type of polar groups but also on the type of isomers. In the particular case of EVAc, it has been shown that meta and para isomers were more prone to find a pattern match that maximizes the number of hydrogen bonds than ortho isomers. In the specific context of food contact materials, the present study suggests that the migration rates of substances such as benzephenones, identified as endocrine disrupting chemicals,68 could be underestimated and should encourage a revision of existing evaluation rules beyond the affirmation “not soluble in water”. This large solute was, in particular, shown to be the less sensitive to the chain polarity and to structure of the polymer. In complement with the proposed approach for homo- and heteropolymers, solvent implicit methods using Molecular Density Functional Theory69 could emerge as powerful and tailored approach to estimate solubilities in water in a few minutes, while relying on atomistic force-fields or quantum calculations.

limits in water, the results confirmed that these three solutes had almost a similar chemical affinity for PE and water. This effect was associated with relatively high χi,PE values around 3 (to be compared with χlikely i,PE ≈ 2.86). By noticing that for food risk assessment, the concentration in the food at thermodynamical equilibrium, denoted Ceq i,F, needs to be overestimated and that the food-to-packaging volume ratio, denoted LF/P, is much larger than unity, additional recommendations can be drawn. Hence, for monolayer materials in contact with food, Ceq i,F reads (see ref 67 for a general discussion): Cieq ,F =

1 Cit,=P 0 K i ,P/F + L F/P

(13)

where Cti,P= 0 is the initial concentration in the polymeric material in contact. ̂ Using K̂ likely i,PE/water in eq 13 instead of Ki,PE/water leads instead to a relative error on Ceq i,F equal to likely K̂ i ,PE/water − K̂ i ,PE/water K̂ i ,PE/water + L F/P

This error remains acceptable while it is positive (i.e., overestimation of Ceq i,F). The risk of underestimation was maximal with 2-nitrotoluene with a concentration underestimated by ca. −17% for LF/P = 100 (volume ratio corresponding to yogurt pot). It could be therefore suggested v that ln K̂ likely i,PE/water ≈ − ln 1/γi,w − 3.86 appears safe for all conditions such that LF/P ≫ 10, while minimizing the risk of a false positive (i.e., too large overestimation of Ceq i,F).



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.6b03683. Mean-field approximation of FHT in perfectly mixed polynary systems; Table S1, studies solutes and key properties; Figures listed S1 and S2 as mentioned in the text (PDF)

5. CONCLUSIONS AND RECOMMENDATIONS Flory−Huggins approaches based on atomistic descriptions of pair contact energies can be used to generate reliable estimates of excess-chemical potentials and Flory−Huggins coefficients in random and block copolymers. In this work, we compared a mean-field approach, requiring only the Flory−Huggins coefficients of homopolymers with a more extensive approach representing explicitly random copolymers. As both approaches led to very similar results even for solutes as large as benzophenone, we confirm the validity of eq 5 for homogenizing pair contact-interactions at the microscopic scale and the possibility to compute rapidly the excess quantities for copolymers, including two or more patterns. The only limitation is that the length of the largest uniform sequence must be smaller than its persistence length, so that the concept of microscopic mixing instead of a macroscopic one holds. When the lengths of uniformly repeated patterns are not random, contact probabilities, pi,k, must be chosen accordingly. For random copolymers, they coincide with their volume fractions. The direct validation against experimental data of the proposed approach was satisfactory. Ethylene-vinyl acetate (EVAc) and aromatic solutes offered a proper challenge test as they combined a polymer and aromatic solutes with different shapes and polar moieties. For nonvolatile substances, reference Flory−Huggins coefficients were inferred from liquid−polymer partition coefficients followed by an identification of pair interactions: solute−liquid, solute−polymer interactions. The Flory approximation was again considered at this stage, and it required a particular attention when the polymer and the liquid were interacting together or when the pure solute was not



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel.: (+33) 169-935063. Fax: (+33) 169-935-024. ORCID

Olivier Vitrac: 0000-0001-7787-5962 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors would like to acknowledge the French National Research Agency for funding this research as part of the collaborative project SAFEFOODPACK DESIGN (reference ANR-10ALIA-009).



REFERENCES

(1) Barton, A. F. M. CRC Handbook of Solubility Parameters and Other Cohesion Parameters, 2nd ed.; CRC Press: Boca Raton, 1991; p 768. (2) Bicerano, J. Prediction of Polymer Properties, 3rd ed.; CRC Press: New York, 2002; p 784. (3) Hansen, C. M. Solubility Parameters: A User’s Handbook, 2nd ed.; CRC Press: Boca Raton, 2012; p 544. 785

DOI: 10.1021/acs.iecr.6b03683 Ind. Eng. Chem. Res. 2017, 56, 774−787

Article

Industrial & Engineering Chemistry Research (4) van Krevelen, D. W.; te Nijenhuis, K. Properties of Polymers: Their Correlation with Chemical Structure; their Numerical Estimation and Prediction from Additive Group Contributions; Elsevier Science: Oxford, 2009; p 1030. (5) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1991; p 408. (6) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed.; Elsevier Science: Bodmin, 2001; p 638. (7) Widom, B. Some Topics in the Theory of Fluids. J. Chem. Phys. 1963, 39, 2808. (8) Boulougouris, G. C. Calculation of the Chemical Potential beyond the First-Order Free-Energy Perturbation: From Deletion to Reinsertion†. J. Chem. Eng. Data 2010, 55, 4140. (9) Boulougouris, G. C. On the Estimation of the Free Energy, From a Single Equilibrium Statistical Ensemble, via Particle Reinsertion. J. Phys. Chem. B 2012, 116, 997. (10) De Angelis, M. G.; Boulougouris, G. C.; Theodorou, D. N. Prediction of Infinite Dilution Benzene Solubility in Linear Polyethylene Melts via the Direct Particle Deletion Method. J. Phys. Chem. B 2010, 114, 6233. (11) Hess, B.; Peter, C.; Ozal, T.; van der Vegt, N. F. A. Fast-Growth Thermodynamic Integration: Calculating Excess Chemical Potentials of Additive Molecules in Polymer Microstructures. Macromolecules 2008, 41, 2283. (12) Hess, B.; van der Vegt, N. F. A. Predictive Modeling of Phenol Chemical Potentials in Molten Bisphenol A−Polycarbonate over a Broad Temperature Range. Macromolecules 2008, 41, 7281. (13) Köddermann, T.; Reith, D.; Arnold, A. Why the Partition Coefficient of Ionic Liquids Is Concentration-Dependent. J. Phys. Chem. B 2013, 117, 10711. (14) Boulougouris, G. C. Multidimensional direct free energy perturbation. J. Chem. Phys. 2013, 138, 114111. (15) Ö zal, T. A.; Peter, C.; Hess, B.; van der Vegt, N. F. A. Modeling Solubilities of Additives in Polymer Microstructures: Single-Step Perturbation Method Based on a Soft-Cavity Reference State. Macromolecules 2008, 41, 5055. (16) Gillet, G.; Vitrac, O.; Desobry, S. Prediction of Solute Partition Coefficients between Polyolefins and Alcohols Using a Generalized Flory-Huggins Approach. Ind. Eng. Chem. Res. 2009, 48, 5285. (17) Gillet, G.; Vitrac, O.; Desobry, S. Prediction of Partition Coefficients of Plastic Additives between Packaging Materials and Food Simulants. Ind. Eng. Chem. Res. 2010, 49, 7263. (18) Vitrac, O.; Gillet, G. An Off-Lattice Flory-Huggins Approach of the Partitioning of Bulky Solutes between Polymers and Interacting Liquids. Int. J. Chem. React. Eng. 2010, 8, No. 6580.2094, DOI: 10.2202/1542-6580.2094. (19) Flory, P. J. Thermodynamics of High Polymer Solutions. J. Chem. Phys. 1941, 9, 660. (20) Huggins, M. L. Solutions of Long Chain Compounds. J. Chem. Phys. 1941, 9, 440. (21) Bawendi, M. G.; Freed, K. F. Systematic corrections to FloryHuggins theory: Polymer-solvent-void systems and binary blend-void systems. J. Chem. Phys. 1988, 88, 2741. (22) Sariban, A.; Binder, K. Critical properties of the Flory−Huggins lattice model of polymer mixtures. J. Chem. Phys. 1987, 86, 5859. (23) Schweizer, K.; Curro, J. Microscopic theory of the structure, thermodynamics, and apparent χ parameter of polymer blends. Phys. Rev. Lett. 1988, 60, 809. (24) Fornasiero, F.; Olaya, M. M.; Esprester, B.; Nguyen, V.; Prausnitz, J. M. Distribution coefficients and diffusivities in three polymers for nineteen aqueous nonvolatile solutes. J. Appl. Polym. Sci. 2002, 85, 2041. (25) Fornasiero, F.; Olaya, M. M.; Wagner, I.; Brüderle, F.; Prausnitz, J. M. Solubilities of nonvolatile solutes in polymers from molecular thermodynamics. AIChE J. 2002, 48, 1284. (26) Scatchard, G. Equilibria in Non-electrolyte Solutions in Relation to the Vapor Pressures and Densities of the Components. Chem. Rev. 1931, 8, 321.

(27) Flory, P. J. Thermodynamics of High Polymer Solutions. J. Chem. Phys. 1942, 10, 51. (28) Huggins, M. L. Some Properties of Solutions of Long-chain Compounds. J. Phys. Chem. 1942, 46, 151. (29) Flory, P. J. Statistical Thermodynamics of Liquid Mixtures. J. Am. Chem. Soc. 1965, 87, 1833. (30) Flory, P. J. Fifteenth Spiers Memorial Lecture. Thermodynamics of polymer solutions. Discuss. Faraday Soc. 1970, 49, 7. (31) Krause, S. Note on a paper by kilb and bueche on solution and fractionation properties of graft polymers. J. Polym. Sci. 1959, 35, 558. (32) Simha, R.; Branson, H. Theory of Chain Copolymerization Reactions. J. Chem. Phys. 1944, 12, 253. (33) Dudowicz, J.; Freed, K. F.; Douglas, J. F. Theory of competitive solvation of polymers by two solvents and entropy-enthalpy compensation in the solvation free energy upon dilution with the second solvent. J. Chem. Phys. 2015, 142, 214906. (34) Stockmayer, W. H.; Moore, L. D.; Fixman, M.; Epstein, B. N. Copolymers in dilute solution. I. Preliminary results for styrene-methyl methacrylate. J. Polym. Sci. 1955, 16, 517. (35) Shiomi, T.; Suzuki, M.; Tohyama, M.; Imai, K. Dependence of miscibility on copolymer composition for blends of poly(vinyl chloride-co-vinyl acetate) and poly(n-butyl methacrylate-co-isobutyl methacrylate). Macromolecules 1989, 22, 3578. (36) Krause, S.; Smith, A. L.; Duden, M. G. Polymer Mixtures Including Random Copolymers. Zeroth Approximation. J. Chem. Phys. 1965, 43, 2144. (37) Kilb, R. W.; Bueche, A. M. Solution and fractionation properties of graft polymers. J. Polym. Sci. 1958, 28, 285. (38) Ten Brinke, G.; Karasz, F. E.; MacKnight, W. J. Phase behavior in copolymer blends: poly(2,6-dimethyl-1,4-phenylene oxide) and halogen-substituted styrene copolymers. Macromolecules 1983, 16, 1827. (39) Kambour, R. P.; Bendler, J. T.; Bopp, R. C. Phase behavior of polystyrene, poly(2,6-dimethyl-1,4-phenylene oxide), and their brominated derivatives. Macromolecules 1983, 16, 753. (40) Takeshita, H.; Shiomi, T.; Suzuki, T.; Sato, T.; Miya, M.; Takenaka, K.; Wacharawichanant, S.; Damrongsakkul, S.; Rimdusit, S.; Thongyai, S.; Taepaisitphongse, V. Miscibility in blends of linear and branched poly(ethylene oxide) with methacrylate derivative random copolymers and estimation of segmental χ parameters. Polymer 2005, 46, 11463. (41) Lindvig, T.; Hestkjær, L. L.; Hansen, A. F.; Michelsen, M. L.; Kontogeorgis, G. M. Phase equilibria for complex polymer solutions. Fluid Phase Equilib. 2002, 194−197, 663. (42) Lindvig, T.; Economou, I. G.; Danner, R. P.; Michelsen, M. L.; Kontogeorgis, G. M. Modeling of multicomponent vapor−liquid equilibria for polymer−solvent systems. Fluid Phase Equilib. 2004, 220, 11. (43) Díaz, I.; Díez, E.; Camacho, J.; León, S.; Ovejero, G. Comparison between three predictive methods for the calculation of polymer solubility parameters. Fluid Phase Equilib. 2013, 337, 6. (44) Oh, J. S.; Bae, Y. C. Liquid-liquid equilibria for binary polymer solutions from modified double-lattice model. Polymer 1998, 39, 1149. (45) Oh, S. Y.; Bae, Y. C. Phase Equilibrium Calculations of Ternary Liquid Mixtures with Binary Interaction Parameters and Molecular Size Parameters Determined from Molecular Dynamics. J. Phys. Chem. B 2010, 114, 8948. (46) Hildebrand, J. H. Solubility. XII. Regular Solutions 1. J. Am. Chem. Soc. 1929, 51, 66. (47) Hildebrand, J. H.; Scott, R. L. Regular solutions; Prentice-Hall: Englewood Cliffs, N.J., 1962; p 180. (48) Sun, H. COMPASS: An ab Initio Force-Field Optimized for Condensed-Phase ApplicationsOverview with Details on Alkane and Benzene Compounds. J. Phys. Chem. B 1998, 102, 7338. (49) van der Vegt, N. F. A.; Briels, W. J.; Wessling, M.; Strathmann, H. Free energy calculations of small molecules in dense amorphous polymers. Effect of the initial guess configuration in molecular dynamics studies. J. Chem. Phys. 1996, 105, 8849. 786

DOI: 10.1021/acs.iecr.6b03683 Ind. Eng. Chem. Res. 2017, 56, 774−787

Article

Industrial & Engineering Chemistry Research (50) Prausnitz, J. M.; Lichtenthaler, R. N.; de Azevedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria, 3rd ed.; Prentice Hall PTR: NJ, 1999. (51) Lyman, W. J. Estimation of physical properties. In Environmental exposure from chemicals; Neely, W. B., Blau, G. E., Eds.; CRC Press: Boca Raton, 1985; Vol. 1, pp 13−48. (52) Domalski, E. S.; Evans, W. H.; Hearing, E. D. Heat Capacities and Entropies of Organic Compounds in the Condensed Phase. J. Phys. Chem. Ref. Data 1984, 13, 1. (53) Domalski, E. S.; Hearing, E. D. Heat Capacities and Entropies of Organic Compounds in the Condensed Phase Volume II. J. Phys. Chem. Ref. Data 1990, 19, 881. (54) Domalski, E. S.; Hearing, E. D. Heat Capacities and Entropies of Organic Compounds in the Condensed Phase. Volume III. J. Phys. Chem. Ref. Data 1996, 25, 1. (55) RSC ChemSpiderThe free chemical database (Royal Society of Chemistry). www.chemspider.com (access date: 10/10/2013). (56) White, C. M. Prediction of the boiling point, heat of vaporization, and vapor pressure at various temperatures for polycyclic aromatic hydrocarbons. J. Chem. Eng. Data 1986, 31, 198. (57) Hilal, S. H.; Karickhoff, S. W.; Carreira, L. A. Prediction of the Vapor Pressure Boiling Point, Heat of Vaporization and Diffusion Coefficient of Organic Compounds. QSAR Comb. Sci. 2003, 22, 565. (58) Economou, I. G. Statistical Associating Fluid Theory: A Successful Model for the Calculation of Thermodynamic and Phase Equilibrium Properties of Complex Fluid Mixtures. Ind. Eng. Chem. Res. 2002, 41, 953. (59) Fouad, W. A.; Wang, L.; Haghmoradi, A.; Asthagiri, D.; Chapman, W. G. Understanding the Thermodynamics of Hydrogen Bonding in Alcohol-Containing Mixtures: Cross-Association. J. Phys. Chem. B 2016, 120, 3388. (60) Kadam, A.; Karbowiak, T.; Voilley, A.; Bellat, J.-P.; Vitrac, O.; Debeaufort, F. Sorption of n-hexane in amorphous polystyrene. J. Polym. Sci., Part B: Polym. Phys. 2014, 52, 1252. (61) Maurin, M. B.; Dittert, L. W.; Hussain, A. A. Thermogravimetric analysis of ethylene-vinyl acetate copolymers with Fourier transform infrared analysis of the pyrolysis products. Thermochim. Acta 1991, 186, 97. (62) Maurin, M. B.; Dittert, L. W.; Hussain, A. A. Mechanism of diffusion of monosubstituted benzoic acids through ethylene-vinyl acetate copolymers. J. Pharm. Sci. 1992, 81, 79. (63) Booth, A. M.; Bannan, T.; McGillen, M. R.; Barley, M. H.; Topping, D. O.; McFiggans, G.; Percival, C. J. The role of ortho, meta, para isomerism in measured solid state and derived sub-cooled liquid vapour pressures of substituted benzoic acids. RSC Adv. 2012, 2, 4430. (64) Grabowski, S. J. Hydrogen BondingNew Insights; Springer: Dordrecht, 2006; Vol. 3, p 524. (65) Landra, C.; Satyro, M. A. Mutual Solubility of Water and Hydrocarbons. J. Chem. Eng. Data 2016, 61, 525. (66) Hoekstra, E. J.; Brandsch, R.; Dequatre, C.; Mercea, P.; Milana, M. R.; Störmer, A.; Trier, X.; Vitrac, O.; Schäfer, A.; Simoneau, C. Practical guidelines on the application of migration modelling for the estimation of specific migration. In support of Regulation (EU) No 10/2011 on plastic food contact materials. EUR 27529 EN; Joint Research Centre: 2015; DOI: 10.2788/04517. (67) Vitrac, O.; Hayert, M. Risk assessment of migration from packaging materials into foodstuffs. AIChE J. 2005, 51, 1080. (68) Muncke, J. Endocrine disrupting chemicals and other substances of concern in food contact materials: An updated review of exposure, effect and risk assessment. J. Steroid Biochem. Mol. Biol. 2011, 127, 118. (69) Sergiievskyi, V. P.; Jeanmairet, G.; Levesque, M.; Borgis, D. Fast Computation of Solvation Free Energies with Molecular Density Functional Theory: Thermodynamic-Ensemble Partial Molar Volume Corrections. J. Phys. Chem. Lett. 2014, 5, 1935.

787

DOI: 10.1021/acs.iecr.6b03683 Ind. Eng. Chem. Res. 2017, 56, 774−787