Prediction of Partition Coefficients of Plastic Additives between

Ind. Eng. Chem. Res. 2010, 49, 7263–7280

7263

Prediction of Partition Coefficients of Plastic Additives between Packaging Materials and Food Simulants Guillaume Gillet,†,‡,§ Olivier Vitrac,*,‡ and Ste´phane Desobry§ Laboratoire National de me´trologie et d’Essais, Centre Energie, Mate´riaux et Emballage, 29 aVenue Roger Hennequin, 78197 Trappes CEDEX, France, Institut National de la Recherche Agronomique, Joint Research Unit 1145 “Food Process Engineering” between INRA, Agroparistech and CNAM, 1 aVenue des Olympiades, 91300 Massy, France, and Nancy UniVersite´, LSGA-ENSAIA-INPL, 2 aVenue de la foreˆt de Haye, BP 172, 54505 VandoeuVre le`s Nancy, France

Partition coefficients, Ki,F/P, between liquids or food simulants, F, and amorphous regions of polymers, P, are important quantities to predict the sorption or desorption kinetics in various areas and in particular to predict the contamination of food by substances originating from food contact materials. This work extends an atomistic Flory-Huggins approach previously developed by us (Ind. Eng. Chem. Res 2009, 48, 5285-5301) to predict Ki,F/P values for large solutes such as antioxidants, light stabilizers, and surface agents. Two extensions were particularly considered. The first extension aims at determining by isobaric molecular dynamics (MD) simulation the contribution of translational entropy in liquids with increasing polarity (isopropanol, ethanol, methanol, ethyl acetate, water) for large and flexible solutes representative of plastics additives. It was found to be higher than the partial molar volume of such solutes, independent of the considered alcohol and satisfactory estimated by the volume accessible to a hydrogen probe. The validity of the coarse-graining approximation for large flexible solutes (octadecane and octacosane) was tested by computing the radial distribution function from MD simulations. A simple correction was proposed to account for the partial overlapping at the coarsegrained level between F molecules and large flexible segments of solutes. The second major improvement extends the whole methodology to water and water-ethanol mixtures. Intrinsic limitations of the Flory-Huggins approximation to handle hydrogen bond cooperativity were overcame by reweighting contact energies in water and by introducing tabulated nonideal properties of water-ethanol mixtures. All predictions agreed well with previously published partitioning data as well as those generated by this study. From experimental values and theoretical considerations, the possibility to predict the contamination of food emulsions with water-ethanol mixtures is finally discussed. 1. Introduction In Europe, all materials intended to be in contact with food must comply with the recently introduced framework regulation 2004/1935/EC1 specific for food contact materials, which enforces a safety assessment and risk management decision for all starting substances and possible degradation products coming from the material. For plastic materials, article 14 of the directive 2002/72/EC2 allows the use of mathematical modeling of desorption to demonstrate the compliance of these materials according to the specific migration limits (SML). Due to toxicological concern, among the 937 substances (including 340 monomers and 597 additives) positively listed in EU directives on plastics in contact with food,3 502 substances (including 230 monomers and 272 additives) are subjected to migration limits. As a result, the development of migration modeling is expected to contribute significantly to an increase of the safety of food products, while reducing significantly the cost of this safety. In other related areas, the recently introduced EU regulation on the safety of chemicals, so-called REACH directives,4 encourages the development of similar predictive approaches to assess the migration of chemicals into the environment or through the food chain. The principles of robust desorption modeling in presence of uncertainty on composition, contact condition and * Corresponding author. E-mail: [email protected]. † Laboratoire National de me´trologie et d’Essais, Centre Energie, Mate´riaux et Emballage. ‡ Institut National de la Recherche Agronomique, Joint Research Unit 1145 “Food Process Engineering”. § Nancy Universite´, LSGA-ENSAIA-INPL.

physicochemical properties has been detailed for monolayer and multilayer materials.5 In particular, probabilistic modeling6 has been applied to evaluate the contamination of 12 retailed food products by an ubiquitous antioxidant7 and to assess the exposure of French consumers to styrene originating from yogurt pots.8 It is argued that a similar methodology could be used by all involved stakeholders (packaging industry, food industry, regulation authorities, food safety agencies...) and updated according to the information available and the followed objective: compliance testing, tier assessment, sanitary survey.... Its dissemination and application to a wide range of substances, contact conditions (temperature, interactions with food) and materials (monolayer and multilayer materials, active materials) is currently limited by either the span of confidence intervals or safety margins mainly due to a large uncertainty on the partitioning of plastic additives with food or food simulants. Indeed, several predictive estimators of diffusion coefficients of additives including uncertainty estimates have been proposed for common polymers such as polyolefins9-13 but only rough approximations have been proposed to predict their chemical affinity for common food simulants. The current “practical guide” for the application of European directives14 advises that the concentrations at equilibrium are equal between the contact material and a fatty product whereas they are 103 times lower when the same material is in contact with an aqueous food. No recommendation is proposed to estimate partition coefficients between layered materials. The complications in measuring partitioning of plastic additives such as hindered phenolic antioxidants or hindered amine light stabilizers arise for their

10.1021/ie9010595  2010 American Chemical Society Published on Web 07/07/2010

7264

Ind. Eng. Chem. Res., Vol. 49, No. 16, 2010

low affinity for polar media (food simulant or polymer) and to their possible reactivity during migration or extraction.15 Alternative methods based on the prediction of activity coefficients are scarce and rarely applied to molecules as large as polymeric plastic additives. Group contribution methods and other semiempirical methods have been recently reviewed.16 The most extensive study was performed by Baner and Piringer17 and concerned volatiles and solutes with intermediate molecular weights. By mixing both the regular solution theory and a Flory-Huggins approach, a correlation model was proposed for partition coefficients between polyethylene materials and alcohols. Because the introduced correction factor was dependent on the size and the shape of the solute, the approach cannot however be extended reliably to plastic additives. Alternatives based on brute force calculations of excess chemical potentials in each phase by molecular simulations could be an alternative. For small solutes in liquid phases, several efficient techniques are available (see the review in18) involving either Gibbs Ensemble Monte Carlo19 or thermodynamics integration.20 In polymers, insertion21 or deletion22 techniques applied to bulky or stiff solutes fail to accommodate the shape of the polymer cavity to the solute shape due to a very slow relaxation of polymer segments at ambient temperature. An efficient free energy perturbation method of chemical potentials in polymers for bulky plastic additives has been devised only recently.23,24 As the method involves an explicit representation of entangled polymer segments and long-term molecular dynamics simulation, it remains highly computationally expensive and has been restricted to liquids and molten polymers at high temperature. To overcome previous limitations and with a goal of tailored calculations for a wide range of solutes, we recently devised an off-lattice generalized Flory-Huggins approach at atomistic scale25 to predict without any mathematical adjustment the partition coefficient of solutes, noted Ki,F/P, between the amorphous phase of a polyolefin, noted P, and a liquid simulant, noted F. As in the original Flory-Huggins approach, the interactions with P are accounted only with the neighbors in contact so that an explicit representation of the entangled polymer is not required. A similar approach has been applied to assess the activity coefficient in F. The total cost of the method was low since the distributions of interactions were calculated by using a full Monte Carlo sampling on molecule pairs and conformers and finally averaged over the number of neighbors. Since a full detailed forcefield at atomistic scale was applied, the method is expected to be applicable to any molecule (of solute i, of P or of F) even if they generate strong attractive interactions such as hydrogen bonding. The method was satisfactory applied to all molecules analyzed by Baner and Piringer,17 to series of flexible molecules (n-alkanes and n-alcohols) and to two prototypes of phenolic antioxidants, with a prediction error of a same magnitude order as the overall experimental error at least for solutes with molecular masses lower than 300 g · mol-1. Additionally, the simulations demonstrated that even if low polar molecules including plastic additives have by design a good solubility in P, they have also a significant chemical affinity for polar simulants such as small alcohols (methanol, ethanol). This effect, which was neglected in previous studies, was related to the higher number of possible rearrangements of F molecules when a solute much larger than F was introduced in the mixture. Within the Flory-Huggins approximation, this contribution is expressed as a change in translational entropy and it involves the number of F molecules reordered or exchanged when a single solute is inserted in the liquid. This number was denoted r-1 i,F . and implies that matter

can spread out randomly without any enthalpic constraint. For almost incompressible mixtures, we argued that a rough estimate of r-1 i,F would lie between the hardcore volume and the molar volume of i. For bulky additives with a melting point much higher than the common migration temperature (i.e., 40 °C), both bounds seem particularly unrealistic as they correspond respectively to a state where all the van-der-Waals surface of the solute is in contact with F and to an ordered state of i. Besides, our description applied to interacting liquids neglect nonideal behaviors due to the perturbation of the network of hydrogen bonds around i. The current work proposes to derive a more rigorous formulation for r-1 i,F for flexible and bulky solutes (beyond 250 g · mol-1) in polar liquids while keeping our molecular tailoring approach of excess chemical potentials and partition coefficients. The extensions were tested against new experimental partition coefficients of large hindered phenolic antioxidants and against isobaric molecular dynamics simulations in different liquids with increasing polarity and decreasing molecular size: isopropanol, ethanol, methanol and water. Since hydrolysis avoided direct measurements of Ki,F/P in water, different mixtures of ethanol in water were used instead. It is emphasized that the mixture 50:50 v/v of water and ethanol is also the new food simulant recommended by EU authorities to assess the migration of plastic additives into dairy products,26 after a significant amount of a photoinitiator (2-isopropylthioxanthone) for printing inks has been found in baby milk27 and of other packaged drinks.28 In presence of long-range Coulombic interactions, the initially proposed pair summation of intermolecular interactions was however expected to minimize significantly the contribution of hydrogen bonding: i) because only molecules in close contact could interact together and ii) because a single molecule could participate at maximum in two longlived hydrogen bonds through a head and tail arrangement pattern. A number of hydrogen-bonded neighbors ranged between 1 and 2 is acceptable for ethanol29 but is unrealistic for water, where a single molecule is involved in up to 6 hydrogen bonds.30 Our initial approach was extended to water by summing water-water pair interactions as there were involved in a 5 sites rigid model of water (i.e., tetrahedral model of water). A similar idea is used in the 5 sites transferable intermolecular potential function31 and it was shown that this geometrical description was be able to predict various properties of liquid water.32 The paper is organized as follows. Section 2 recapitulates the theoretical principles of calculation of Ki,F/P from a Monte Carlo (MC) sampling of pair contact interactions at atomistic scale. Section 3 details the experimental and isobaric molecular dynamics (MD) simulation conditions that were used to generate reference data and predictions. The computational cost associated to direct calculations by MD simulations of partial molar volumes of additives in F was significantly reduced by comparing their values with the volume enclosed to the surface accessible to an hydrogen atom. Section 4 summarizes the predictions of Ki,F/P in all tested conditions for large plastic additives. The improvements induced by new r-1 i,F estimates were analyzed against Ki,F/P values previously collected for homologous series of solutes17,33 and according F-solute radial distribution functions calculated for large flexible molecules. In particular, previous allegations of low partitioning of plastics additives in water and in water-ethanol mixtures were discussed according to a sensitivity analysis to nonidealities and according to experimental values generated by this study in water-ethanol mixtures and those extrapolated by Gandek34 in pure water. Finally, some suggestions are given on how current predictions


could be used to assess consumer exposure to ubiquitous hindered phenolic antioxidants and could be used to select an appropriate simulant of aqueous food products in the perspective of compliance testing of food packaging materials and plastics pipes used in water distribution systems.

K as prescribed in EU regulation. 〈〉 and 〈〉T are the average operator and the temperature ensemble average operator respectively. The latter was obtained by weighting the sampled distribution of energies pA+B(ε) with the Boltzmann factor, defined as exp(-ε/kBT):

2. Theory 2.1. Definition of Partition Coefficients. Since the diffusion of a solute i, such as a plastic additive, requires the cooperative motions of the polymer,12 the substance i is assumed to diffuse only in the amorphous phase of the polymer. The partition coefficient between F and the amorphous phase of the polymer, P, denoted Ki,F/P, is assessed as:25 Ki,F/P )

Ci,F | eq Ci,P | aeq

) (1 - c)

Cî,F | eq 1ˆ C | - biCi,P | t)0 ai i,P eq

(1)

where Ci,F|eq andCi,P|aeq are the concentrations at equilibrium in F and P respectively. {Cî,k|eq}k)F,P are the corresponding concentrations as experimentally assessed. The concentration in P is in particular associated with an extraction yield ai lower than 1. In addition, since large additives have a melting point much higher than the extrusion temperature of polyolefins, eq 1 takes also into account that a significant fraction bi of the initial additives content, Ci,P|t)0, could be crystallized and does not participate to the apparent equilibrium due to kinetic limitations. c is the crystallinity of the polymer (volume fraction of crystalline phase), which may vary according to the considered sample. 2.2. Partition Coefficient Estimates in Alcohols. According to Flory-Huggins theory,35-38 the activity coefficients relative to volume fractions in solute, denoted {γVi,k}k)F,P, can be calculated on a rigid lattice, whose mesh size is commensurable to the volume of the solute Vi. In this work, we used a similar approach without any lattice by sampling explicitly at atomistic scale the molecular pair interactions. Comparatively to previous approaches involving solubility coefficients, the superiority of this description comes from the sampling of the dispersion of interaction energies due to the different conformers and to accessibility constraints. This contribution on enthalpic terms, identified in many books as nonpositional entropy contribution, can therefore be explicitly calculated. In our initial formulation,25 the excess translational entropy in F was expressed within the cell theory of liquids: a bulky solute occupies r-1 i,F > 1cells of F, where the volume of a single cell was chosen equal to the volume of a single molecule of F in a pure liquid of F. By neglecting the sampling biases (see the full expression in eq 16 of ref 25), Ki,F/P was approximated as: where {χi,k}k)F,P is the

〈εA+B〉T )

∫ ∫

+∞ p (ε)e-ε/kBTεdε -∞ A+B +∞ p (ε)e-ε/kBTdε -∞ A+B

(3)

The sampling of pA+B(ε) was based on a large conformers set of A (seed molecule) and B (contact molecule) representative of their condensed state (up to 104 configurations) and based on all possible contacts of their van der Waals envelopes with spherical symmetric probability (up to 1011 configurations). As previously discussed,25 contacting was performed either by expansion or contraction to avoid an oversampling of cavities in large or flexible molecules. An infinite long polymer was simulated by replacing with an oligomer-like, where head and tail atoms were noncontact in order to mimic the hindrance associated to the connectivity of these atoms. Within the Flory’s “pair local energy approximation”, the optimal oligomer length represents the backbone bonds range of nonbond interactions within the polymer chain.39 In practice, χi,P depends weakly on the number of monomers40 and the optimal number of monomers was set to minimize the overestimation of noncontact “head and tail effects” while preserving the overall convergence of the sampling. As a rule of thumb, it was demonstrated that the overall bias was minimal when the test oligomer was commensurable to the solute.25 M The entropic term r-1 i,F was defined as the ratio Vi/VF , where M VF is the molar volume of pure F. The choice of Vi was highly critical as it affects significantly the quality of prediction via eq 2. By assuming that the mixtures i+F are compressible, Vi was estimated as the infinite dilution partial molar volume of i in F under isobaric and isothermal conditions, denoted VPi,F. As VPi,F values were not tabulated for plastics additives molecules in the literature, they were calculated from MD simulations. However, the assumption of a homogeneous compressibility close to i appeared questionable for large and highly flexible solutes, because the same cell (i.e., a volume commensurable to solute) can be simultaneously occupied by atoms of i and by F molecules due the persistence of a local order close to i. This deviation to mixing ideality is depicted in Figure 1 for a highly flexible solute dispersed in liquid methanol. The sketch represents two methanol molecules “enclosed” within the octacosane contour. Since the excess entropic contribution, r-1 i,F depends only on the rearrangement of the center of mass of molecules F and i regardless of the true shape of molecules,40 a general formulation of r-1 i,F was derived by correcting the number of “displaced” F molecules, expressed as the ratio VPi,F/VM F by the fractional number of F molecules overlapping i, denoted nFoverlap neighboring i: r-1 i,F )

Flory-Huggins parameter, kB is the Boltzmann constant, T is the absolute temperature. The enthalpies associated to all binary mixtures P+i, F+i, P+P and F+F were calculated as the product of the average number of neighbors, {zA+B}A,B)i,F,P, times the average enthalpy associated to a single pair {εA+B}A,B)i,F,P at the temperature T. The temperature was chosen comparable to the one used in reference data, commonly 313

7265

2

VPi,F VM F

overlap - nFneighboringi

(4)

For n-alkanes, nFoverlap neighboring i was derived from the ensembleaveraged pair radial distribution of F molecules around the solute i during isobaric MD simulations. Such results were used to validate a much less computationally expensive approximation, where VPi,F was replaced by the conformational ensemble average of the solute volume enclosed within the surface accessible to a hydrogen atom, denoted ViH. Indeed, it was

7266

Ind. Eng. Chem. Res., Vol. 49, No. 16, 2010 V -1 -1 ln γi,F ) (1 - φi) - ri,F φ - ri,F φ + [(χi,F1φF1 + 1+F2 1 F1 2 F2

( )

χi,F2φF2)(φF1 + φF2)] - χF1,F2

Vvdw i VFP2

φF1φF2 -

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

hihF2φF2

dχi,F2

- φiφF1φF2

∂χi,F1

- φiφF2 1

∂χi,F1

dhF2 ∂φF2 ∂φF2 Vdw ∂χi,F1 ∂χF1,F2 Vi φF1φF2 2 φiφF2 1 p ∂φF1 ∂φF2 VF2 ∂χF1,F2 ∂χi,F1,F2 Vvdw i φF2 1φF2 - φiφF2 1φF2 P ∂φF1 ∂φF1 V

-

F2

φiφF1φF2 2

∂χi,F1,F2 ∂φF2

- χi,F1,F2φF1φF2(1 - 2φi) (5)

where {φk}k)i,F1,F2 is the volume fraction in k in the mixture so that φi + φF1 + φF2 ) 1; hi ) φi/(φi + φF2) and hF2 ) φF2/(φi + φF2). At infinite dilution in solute i, φif0, the ternary Flory-Huggins interaction parameter χi,F1,F2 was reasonably supposed negligible. The entropic contribution was related to the size of the cavity to insert the solute i in -1 F1 and in F2, r-1 F1 and rF2 respectively. By noticing that φF1 ≈ (1 - φF2) when φi ≈ 0, the practical estimate of Ki,F,P was finally rearranged as a third-degree polynomial in φ2, while separating enthalpic and entropic (related to positions) contributions: Figure 1. Interpretation of packing of F molecules around a solute i at atomistic scale and at mesoscopic scale. As the translational entropy depends only on the position of the center of masses, it may be calculated at mesoscopic scale irrespectively to the true radial distribution is kept.

argued that close contacts between low polar solutes and F were performed mainly via hydrogen atoms of F rather than by any other heavy atoms. 2.3. Partition Coefficient Estimates in Water and in Water-Ethanol Mixtures. The estimate of γVi,F introduced in eq 2 was only valid in binary mixtures. The extension of Flory-Huggins theory to ternary mixtures, including a solute i and a simulant comprising two components F1 and F2, leads to a bilinear equation system.41-43 We used a formulation similar to43 by combining two coarse-graining processes that aims mainly at preserving the radial pair distributions around the solute i. Consistently with our sampling methodology of contact energies at atomistic scale, which forced contacts between van der Waals surfaces, enthalpic terms were renormalized at the scale of the van der Waals volume of each diffusant. The introduced approximation of solute-F interactions was particularly acceptable for almost spherical F molecules much smaller than the solute. By contrast, entropic terms were calculated as in binary mixtures at the P scale of ViH, chosen as a simple estimate of Vi,F 1+F2. In addition, the nonideal character of the water-ethanol mixtures was accounted by using partial molar volumes of F1 and F2, denoted VPF1 and VPF2 respectively, instead of their molar volumes. Since only the behavior of the mixture at infinite dilution of i was considered, ternary interaction terms were neglected. The activity coefficient associated to the volume fraction in i in F was finally expressed as:

Equation 6 is continuous in φF2, with values equal to r-1 i,F1 + χi,P - χi,F1 and r-1 + χ χ for binary mixtures i,F2 i,P i,F2 corresponding to φF2 ) 0 and φF2 ) 1 respectively. As high order terms depend only on enthalpic contributions of F1 and F2, they were strongly affected by the non ideal character of the water-enthanol mixture. The practical consequence was that ln Ki,(F1+F2)/P was a nonlinear function of the molar fraction in ethanol. To avoid an overestimation of these nonlinearities while keeping calculations tractable, the properties of water-ethanol were derived from literature data and models instead of being estimated from MD simulations. Flory-Huggins interaction parameter between water and ethanol, χwater,ethanol, was derived from the polynomial approximation of molar heat of mixing,44 m ∆Hwater+ethanol , via: χwater,ethanol ) xethanol(1 - φethanol)

m ∆Hwater+ethanol RT

(7)

where xethanol and φethanol are the molar fraction and the volume fraction in ethanol respectively. R is the universal gas constant. As depicted in Figure 2, χwater,ethanol was found to vary nonlinearly between -4.8 and -0.87 with an average value,


7267

χi,water ≈ 1 〉 + 〈zwater+i〉)〈εi+water〉T - 〈zi+i〉 an 〈εi+i〉T [(〈z 2kBT i+water 4〈zwater+water〉〈εwater+water〉T] (8) It is interesting to notice that an almost similar value can be derived with the conventional Flory-Huggins description on a 3D cubic lattice. On a 3 × 3 × 3 lattice, a seed molecule is in direct contact with 6 neighbors (lattice coordination number) and surrounded by 33 - 1 ) 26 molecules with possible longrange interactions. The discrepancy between a short-range sampling and a long-range sampling of interactions is about 6/26 ≈ 1/4.33. From these considerations, water in clusters appeared as a prototype of branched molecules, whose interactions on a lattice are discussed in ref 48. 3. Materials and Methods

Figure 2. Flory-Huggins interaction parameter between water and ethanol: (a) polynomial model derived from heat of mixing at 25 °C;35 (b) its derivative with Φethanol.

∫10χwater,ethanoldφethanol, about -2. By choosing F2 ) ethanol and a constant χwater,ethanol equal to its average, eq 6 becomes a second-degree polynomial in φethanol of which the second-degree term predicts a convex variation of ln Ki,(F1+F2)/P with φethanol. In particular, it is highlighted that the decrease of ln Ki,(F1+F2)/P with φethanol was expected to be higher in a mixture enriched in ethanol. In this work, additional nonlinearities were accounted by introducing for each water-ethanol mixture an estimate of ∂χwater,ethanol/∂φethanol (Figure 2b) and VPethanol values derived from experimental densities of water-ethanol mixtures.45 By choosing F1 ) water, the interaction between solute and water, χi,water could be inferred almost similarly as χi,ethanol. Nevertheless, hydrogen bond cooperativity could not be sampled, because our generalized Flory-Huggins formulation averaged interactions only between dimers of water and not on the likely structure of liquid water at room temperature involving tetramers.46 In other words, the current Monte Carlo method allowed only random contacts (i.e., with no preferred orientations) but not long-lived contacts as it should be assessed in the bulk. Indeed, the lifetime of hydrogen bonds in bulk water was estimated between 10-12 s and 20 × 10-12 s, while the lifetime of broken hydrogen bonds was about 10-13 s.47 To extend naturally our formulation without additional cost, pair interactions were assumed to be summed over a tetrahedral lattice and possible “dangling” hydrogen bonds were neglected. From the statistical point of view, the proposed summation mimicked the sampling of bulk water potential energy with a five-site water model such as the one involved in TIP5P forcefield.31 Formally, a factor 4 was therefore introduced for the term 〈zwater+water〉〈εwater+water〉T to account for the cooperativity of hydrogen bonding:

3.1. Polymer Formulation. Plastic additives were provided by CIBA (Basel, Switzerland) except for Erucamide, which was obtained from CRODA (Vimodrone, Italy). Their main characteristics are summarized in Table 1. Virgin high density polyethylene (HDPE) was obtained from ATOCHEM (Paris, France). Dried polymer flakes were formulated with the six tested additives during a first extrusion step prior to final processing at semi-industrial scale. Both extrusions were performed at 200 °C in a biscrew extruder (model BC-21 Clextral, France) and the material was finally calendered as 0.15 m wide ribbons. The final density and the crystallinity were of 940 kg m-3 and 70% respectively, with a melting point of 136 °C. 3.2. Experimental Determination of Partition Coefficients. Partition coefficients of six plastic additives were assessed by putting 10 cm3 of polymer samples cut in small 5 mm ×5 mm pieces in contact with 100 mL of simulant. Tested simulants included isopropanol (SIGMA-Aldrich, U.S.A.) and three mixtures of ethanol (MERCK, Germany) and water (Milli-Q water, MILLIPORE, U.S.A.) with volume fractions in ethanol of 95, 75, and 50%. The amount of nondesorbable substances, denoted bi in eq 1, was estimated from the intercept of desorption linear isotherm (expressing the concentration in P versus the concentration in F) obtained with different mass ratios of F:P (50, 35, 20, 10, 5). 3.3. Extraction of Additives from Polymer Materials. Initial concentrations, Ci,P|t)0, and residual concentrations after contact, Ci,P|eq, were assessed after pressurized solid-liquid extraction (model ASE 200, Dionex, USA) in a 75:25 v/v mixture of isopropanol-dichloromethane at 105 °C and 10.3 MPa. The protocol was similar to the one described by Coulier et al.49 and applied to 600 mg of samples mixed with sand in 11 mL extraction cells. Two static cycles of 15 min each were applied to achieve complete extraction. To prevent the degradation of additives during extraction, 100 µL L-1 of triethylphosphite (Sigma-Aldrich, U.S.A.) was added to solvents. It was verified that this methodology gave similar results than extraction by reflux during 40 h. 3.4. Concentration Measurements of Additives in Solvents and Simulants. The six additives were separated and detected either using high performance liquid chromatography (HPLC) associated to an UV diode array detector (DAD) and an evaporative light scattering detector (ELSD) in series or using gas chromatography (GC) associated to a flame ionization detector (FID). ESLD was only applied to Irganox PS802. GC was not applied to Irganox 3114 due to its insufficient volatility.

7268


Table 1. Main Characteristics of Experimentally Studied Plastics Additives additive

CAS number

Irganox 1076 Irganox PS802 Irgafos 168 Chimassorb 81 Tinuvin 326 Erucamide

technological function

Vivdw (Å3)

antioxidant antioxidant co-stabilizer UV-stabilizer UV-stabilizer surface agent

647 ( 12 850 ( 30 747 ( 3 344 ( 2 285 ( 3 424 ( 6

2082-79-3 00693-36-7 31570-04-4 01843-05-6 03896-11-5 00112-84-5

ViM (Å3)

melting range (°C) e

949 1235 1043 508 416 642

50-55 64-67 183-186 48-49 138-141 79

bi (%) 0.002 0 0.03 0.007 0.008 0

ViM was calculated at 313 K with densities (averaged or extrapolated at the good temperature) from the Chemical Abstracts database; melting ranges were obtained from technical datasheets of products; bi’s were estimated with results of Figure 3. Table 2. Linear Solvent Gradient of Our HPLC Method a

b

c

time (min)

solvent A (%)

solvent B (%)

solvent C (%)

0 2 30 32 37 42 53

40 40 0 0 0 40 40

60 60 100 0 0 60 60

0 0 0 100 100 0 0

a Solvent A: aqueous 10 mM NH4Ac (VWR, France) solution adjusted to pH 9.5 with 25% aqueous NH4OH (ACROS organics, Geel, Belgium) to which 500 µL/L n-hexylamine (SIGMA-ALDRICH, USA) was added. b Solvent B: acetonitrile (HPLC grade, CARLO ERBA, France) to which 500 µL/L n-hexylamine was added. c Solvent C: isopropanol to which 500 µL/L n-hexylamine was added.

The HPLC protocol was similar to the one described by GarridoLo´pez et al.50 The HPLC system consisted in a Waters 717 plus autosampler, a Waters 600 controller equipped with a thermostatted column compartment and an in-line degasser AF (Waters, U.S.A.). Separation was achieved on a Xterra C8 column (150 mm ×3.0 mm; 5 µm particles; Waters, U.S.A.) operated at 60 °C. The linear solvent gradient is described in Table 2. The flow rate was of 0.8 mL/min and the injection volume was of 30 µL. UV detection was performed between 200 and 300 nm with a PDA model 2996 (Waters, USA). Light scattering was assessed with a DDL 31 detector (Eurosep, France) using a nebulization temperature of 110 °C and an air gas flow of 1.5 L/min. The GC system comprised an AutoSystem XL (PerkinElmer, U.S.A.) used in on-column mode. The injector temperature was kept at 45 °C during 2 min followed by a linear increase up to 320 °C with a heating rate of 100 °C/min. Separation was carried out on a BPX5 capillary column (30 m × 0,25 mm × 0,25 µm; SGE Europe, France). After injection of 1 µL, the temperature was kept at 50 °C during 1 min, increased up to 330 at 15 °C/min and was at 330 °C for 5 min. Detector was permanently held to 340 °C. Because similar results were obtained with GC and HPLC methods, they are not distinguished in the remainder of this work. 3.5. Monte Carlo Sampling and MD Simulation Procedures. The main originality of our approach was to apply the same methodology and the same theoretical framework to sample the interactions between a solute and much larger polymer chains and to sample the interactions between a solute and much smaller F molecules. The tailoring and coverage of the whole methodology were reinforced by the use of an ab initio class II (i.e., including anharmonic and cross-coupling terms) forcefield (COMPASS forcefield, Accelrys, USA), which is optimized for condensed phases and which incorporates all the parameters required for the study of polymers,51 plastics additives and food simulant including in some extent water via a 3-site model.52 As only pair interactions between water molecules were explicitly considered, the geometric approximation of hydrogen bonding via a three-site model was argued to be appropriate.

As previously described,25 Monte Carlo sampling of pair interaction energies was performed by combining a proprietary pre- and postprocessing code, so-called Molecular Studio, developed for Matlab (Matworks)/Scilab and the classical MD simulation software Discover (Accelrys, U.S.A.). A large set of conformers (103 - 104) of each species was generated from long-term MD simulation at 313 K. Conformers of additives in polymer and conformers of polymer segments were generated in vacuum by neglecting intermolecular interactions. By contrast, conformers in F were extracted from isobaric MD simulations in periodic systems including 103 F molecules. Pair contact energies, {εA+B}A,B)i,F,P, relied on 107 samples for F molecules and up to 1011 samples for large plastics additives. Potential energies were calculated without applying any cutoff. The number of neighbors, {zA+B}A,B)i,F,P, was based on up to 104 trials of packing of contact molecules around a seed molecule. The contact surface was set to the van der Waals surface of each atom. Initial packing configurations in F including one to four solutes and 103 molecules were constructed at 313 K with the Amorphous Cell module (Accelrys, U.S.A.), using a compression box technique, as described by Rigby.52 The target density at this stage was setup to be 5-8% lower than that of the corresponding ideal mixture. Each cell was subsequently submitted to a long-term MD in the canonical NPT ensemble with an integration time step of 10-15 ns. Temperature and pressure were kept constant by coupling the simulated cell with an infinite external heat bath (Andersen thermostat) and with an isotropic pressure bath (Berensen barostat). MD simulations were equilibrated during 1 ns and repeated on 5 initial packing configurations. Van-der-Waals interactions were calculated at atomistic scale using a distance cutoff of 15.5 Å and an analytical spline tail correction of 5 Å (buffer width ) 2 Å). Coulombic interactions were calculated via a 3D Ewald summation technique with a lattice update width of 5 Å. The convergence of typical MD runs were tested by checking that simulated systems had a constant energy in the microcanonical (NVE) ensemble, that is without temperature and pressure baths coupling. Partial molar volumes, VPi,F, at T ) 313 K and a pressure p ) 105 Pa, were inferred from the volumes of equilibrated cells F+i, VF+i, containing ni ) 1-4 solute molecules dispersed among a constant number of F molecules, denoted nF. VPi,F, defined in eq 9, was calculated as the slope of the regression line of the mixture volume with the number of solutes i. VPi,F )

( ) ∂VF+i ∂ni

(9)

T,p,nF

Possible biases introduced by MD simulations were detected by comparing the intercept of the regression line with the theoretical volume of pure F. In addition, it was verified that doubling nF or changing the initial density did not affect the density at equilibrium.


7269

Figure 3. Desorption isotherms of typical additives in ethanol 95%: (a) linear scale; (b) log-log scale. Vertical and horizontal bars represent 95% confidence intervals due to experimental errors. The regression line with zero intercept and the normal regression line are plotted in continuous and dashed line respectively. Outliers are removed from the regression for Erucamide and Irgafos 168.

A cost-efficient estimate of VPi,F was tested by calculating the volume enclosed within the Connolly surface accessible to a hydrogen atom and denoted ViH. This surface was sampled for several conformers by rolling around the solute a spherical probe with a radius comparable to a hydrogen atom. The enclosed volume was finally calculated by Delaunay tessellation. 4. Results and Discussion 4.1. Typical Desorption Isotherms of Plastic Additives in Ethanol. Typical desorption isotherms of plastic additives including light stabilizers, an antioxidant and a lubricant (see Table 1) are plotted in Figure 3 for F ) ethanol 95%. Several equilibriums were carried out by varying the volume of F in contact. Contact times were longer than 3 months to ensure a thermodynamical equilibrium even for large additives and to mimic the long-term storage of food products. It is important to notice that all additives belonged to a same formulation of high density polyethylene and were consequently obtained in similar conditions. All isotherms appeared linear with an intercept close to 0. According to eq 1, the slope was proportional to the partition coefficient when bi was negligible. Estimates of bi were given by the intercept with Ci,F|eq and represented between 0 and 0.2% of the initial content in additive. These values were roughly correlated to the melting point of the considered additive (Table 1) and were used to derive accurate confidence intervals on partition coefficients including biases as defined in eq 1. Significant deviation to linearity occurred only for Erucamide. Its chemical affinity for F increased when the concentration in P at equilibrium increased. This effect was related to the poor compatibility of its unsaturated long chain carboxylic acid amide with polyethylene. The induced phase separation in polyethylene at high concentration is technologically used to modify the surface tension of polyolefins. The substance is commonly referenced as lubricant, antifogging or slip agent and is highly soluble in alcohols. In this case, only the first linear part of the isotherm was considered to determine partitioning. On the opposite, Irgafos 168, which is used as peroxide decomposer, had the lowest chemical affinity for ethanol. The significant discrepancy between concentration measurements of Irgafos 168

in particular in F was attributed to the high reactivity of its trivalent phosphorus. 4.2. Size of Cavities to Insert Plastic Additives in F. Figure 4a presents typical runs of isobaric 250 ps MD simulations at 313 K and 105 Pa to assess the volume at equilibrium of typical mixtures including 1-4 octacosane in 103 molecules of ethanol and methanol. The simulations started with initial configuration with a density lower than the one expected for ideal mixtures in targeted conditions of pressure and temperature. As a result, the systems underwent a contraction accompanied by a reduction in total potential energy caused by the progressive intensification of hydrogen-bonding attractive interactions. Volumes of mixtures at equilibrium were obtained by averaging asymptotic values during the last 50 ps. To avoid particular initial configurations, all simulations were repeated on five different initial cells. It is underlined that starting for more condensed states led to similar results but generated higher instabilities (e.g., temperature drifts) in the explicit integration scheme due to higher energy barriers to update the configurational space of the simulated system. The corresponding linear correlations between the volume of the mixture and the number of inserted octadecane are plotted in Figure 4b. The consistency of MD results appear as intercepts, which match the specific volume of the considered pure liquid simulants. Furthermore, although some fluctuations were noticed between five repetitions, a satisfactory linearity was achieved between one and four solute insertions. All calculated values are gathered in Table 3. It is worth noticing that VPi,F values obtained for a same solute and different simulants were closely related so that VPi,F would depend weakly on the considered simulant. This last assumption is analyzed in Figure 5 by comparing VPi,F against common simple descriptors of simulant and solute properties. From a geometrical point of view, when host molecules F were smaller, it could be thought that VPi,F of large additives should be smaller, with values possibly lower than their corresponding molar volumes, ViM. As shown in Figure 5a, such a behavior was never observed and all VPi,F values were significantly greater than ViM ones and were either independent of F or decreased when the size of F increased. This last effect

7270


Figure 4. Typical NPT MD runs at T ) 313 K and p ) 105 Pa in F ) methanol (empty symbols) and F ) ethanol (filled symbols) when O ni ) 1; ∆: ni ) 2; 0: ni ) 3; ]: ni ) 4 molecules of octacosane are added to 1000 F molecules: Total energy (a) and volume (b) are cumulative averages. Typical partial molar volumes are calculated as reference slopes (c).

was well identified with solutes including a BHT pattern (e.g., BHT itself or Irganox 1076) and was approximately correlated to the electrical dipole moment of F with an increase of VPi,F by 40% when isopropanol was replaced by water. These results showed that the volume of simulant molecules displaced by the insertion of a solute, r-1 i,F , was controlled by the volume of the inserted solute and by a nonideal effect related to the perturbation of the network of attractive interactions in F (see Figure 1 as a possible interpretation). In other words, the large external surface of the solute hindered locally the creation of attractive electrostatic forces, which were responsible for a higher molecular compaction and cohesion in pure F. Because all investigated solutes comprised only a few hydrogen bond donors, they were statistically in contact with the most abundant atoms in F, that is hydrogen atoms. It was therefore argued that VPi,F should be correlated with the volume enclosed within the solute external surface accessible to hydrogen atoms, ViH. A good agreement was observed for both small solutes and large additives (Figure 5b). In more polar simulants such as water, results obtained on BHT (Figure 5a) suggested that the perturbation region around large additives should be thicker than a single hydrogen atom and ViH values should underestimate r-1 i,F and consequently Ki,F/P. Due to the high computational cost of VPi,F calculations by MD simulations,

ensemble-averaged ViH values were privileged to assess partitioning via eq 22 or eq 6. This choice confirmed our first allegations,25 for which the large entropic contribution involved in the partitioning of plastic additives should be controlled by their large hardcore volume plus an interstitial volume that would be proportional to the surface exposed to F. In alcohols, using ViH led to an interstitial estimate that was on average 9% larger than the one derived using ViM as previously proposed as an upper bound of rF-1.25 For flexible solutes, the interstitial volume was thought to vary according to the configuration of solutes and according to considered interactions (e.g., without intermolecular interactions, in poor solvent). This possible effect was systematically tested for octadecane and octacosane. The corresponding gyration radii, ViH and VPi,F values are plotted in Figure 6. Without intermolecular interactions, large n-alkanes adopted mainly randomly coiled configurations, whereas more stretched configurations were privileged in ethanol and methanol. Surprisingly, although significant discrepancies appeared between sets of conformers, the ensemble-averaged ViH values were almost independent of the way conformers were sampled. In particular, such methodological results confirmed the principle of an inexpensive tailored sampling of volumes “in vacuum” (without intermolecular interactions). In details, the volume enclosed within the accessible surface to a hydrogen atom underestimated slightly the true volume of cavity because of additional excluded volumes due to the size of F molecules. There was, however, an exception for octacosane in methanol for which the cavity was smaller than ViH. This effect was related to a close packing between methanol molecules and large aliphatic chains close to their correlation length. 4.3. Solute-Simulant Radial Distribution Functions. Several plastics additives, among them plasticizers, surface agents, and antioxidants, include an aliphatic chain, which tends to “accommodate” the packing of small F molecules more efficiently than bulky rigid group (e.g., a BHT pattern). This effect was analyzed specifically as it counterbalanced the effect of solute in size in the translational entropic contribution. It was shown to be particularly significant for n-alkanes larger than tetradecane, where our initial formulation overestimated strongly experimental Ki,F/P values.25 Such deviations were inherent to our mean field aproximation, which assumed that two species could not occupy the same space below the length scale of a statistical segment. Since this length scale was chosen related to the solute volume to keep a similar chemical potential expression in P and in F, details smaller than the solute gyration radius were not captured. Fluctuations in F density around flexible solutes in methanol and ethanol are analyzed in Figure 7 by plotting the radial distributions of F-F and i-F pairs averaged over several cells and more than 100 independent time snapshots. Density correlations of small liquid alcohols, denoted gF,F(r), persisted over distances up to 6 molecular diameters (between 10 and 15 Å), that is on length scales close to the gyration radius of most of plastic additives. The first peak was related to a contact via a hydrogen bond (e.g., methanol or ethanol) or even a double hydrogen bond (e.g., ethanol), whereas following contacts involved hydrogen to another hydrogen or hydrogen to a carbon interaction. Around flexible solutes, the density in F molecules, gF,i(r), was by contrast always lower than in the bulk and increased regularly outward. For a molecule larger than its persistence length, such as octacosane, the solute center-of-mass was accessible to F molecules, but it was inacessible for a shorter solute such as octadecane. By assuming a mean field approximation of flexible solutes with a spherical symmetry,


7271

Table 3. Solute Partial Molar Volumes solvent F methanol

ethanol

isopropanol

water ethyl acetate

solute i

P Vi,F (Å3)

predicted F densitya (kg m-3)

experimental F density at 313K (kg m-3)b

BHT diphenylmethane D-limonene Irganox 1076 octacosane octadecane BHT diphenylmethane D-limonene Irganox 1076 octacosane octadecane BHT diphenylmethane D-limonene Irganox 1076 octadecane BHT BHT

408 ( 10 283 ( 3 280 ( 8 971 ( 1 776 ( 17 584 ( 11 384 ( 3 275 ( 2 282 ( 3 961 ( 8 909 ( 31 568 ( 17 390 ( 12 291 ( 12 277 ( 8 962 ( 1 564 ( 4 466 ( 23 411 ( 25

776 ( 3 773 ( 1 773 ( 2 770 ( 1 766 ( 3 774 ( 4 772 ( 1 771 ( 1 771 ( 1 771 ( 2 771 ( 1 771 ( 4 795 ( 2 798 ( 2 796 ( 2 796 ( 1 795 ( 1 943 ( 3 874 ( 4

772

0.9988 0.9999 0.9988 1.0000 0.9992 0.9993 0.9999 1.0000 0.9998 1.0000 0.9981 0.9994 0.9987 0.9974 0.9985 1.0000 0.9999 0.9957 0.9939

772

781

992 /

a Calculated from the intercept of the regression line of volumes against the number of solute molecules. data from the Chemical Abstracts Databases and from CRC Handbook.36

b

R2

Averaged or extrapolated at 313 K with

within the 4 Å around the center-of-mass and 7 Å after the center-of-mass. For octadecane, only the last position was available. There was however no difference in ethanol, where the first 4 Å and positions 7 Å and 11 Å were accessible. This behavior is illustrated on some simulation snapshots at atomistic scale in Figure 8 in the solute plane defined by its two main moment axes. The density in F was averaged along the endto-end distance of the solute in direction perpendicular to the inertia plane. The cavity to insert the solute appeared as an area with low density and following approximately its contours. The compaction around octacosane looked higher as the projected plane was thicker. Within the framework of our mean field approximation of F-i interactions, the average number of F molecules, nFoverlap neighboring i, which statistically overlapped the solute volume, was derived from the radial position of the inflection point of the following criterion: 4 C(r) ) πr3 - VM F 3

P Figure 5. Partial molar volumes of solutes i in different simulants, Vi,F , and their correlations with common properties of simulant and solute: (a) electrical dipole moment of F and (b) volume enclosed to the surface accessible to hydrogen atoms, VH i . Molar volumes of solutes appear as dotted horizontal lines in (a). R: i ) BHT; β: i ) diphenylmethane; γ: i ) D-limonene; δ: i ) Irganox 1076; ε: i ) octacosane; ζ: i ) octadecane.

the equivalent blob consists in the solute itself and in some F molecules in close contact. The excess equivalent potential mean force, PMFF neighboring i, to insert an F molecule close to the solute was estimated via the Boltzmann distribution law: PMFF neighboring i(r) ) -kBT log

gF,i(r) gF,F(r)

(10)

Corresponding normalized distributions of PMFF neighboring i are depicted in panels b and d of Figure 7. They showed that methanol molecules could be inserted without additional energy at two radial specific locations along the contour of octacosane:

∫g r

0

2 F,i(r)4πr dr

(11)

This criterion defined how the volume not occupied by the F molecules varied when the distance from the solute center-of-mass increased. Its first derivative is plotted in Figure 9 along with the cumulated number of F molecules around the solute. Beyond the inflection point, the solute influence on density ceased to be significant, and free voids around the solute and F molecules were almost independent. noverlap F neighboring i was defined as the number of molecules enclosed within the radius defined by the first inflection point. These radii are also reported in Figure 8. Up to 5.16 and 0.88 methanol molecules were found overlapping for octacosane and octadecane, respectively. These figures decreased down to 1.55 and 0.26 in ethanol due to the higher similarity between the packaging of ethanol molecules in the bulk and along the diffusant. Hence, ethanol molecules are mainly selforganized in dimers involving two upside-down hydrogen bonds.22 In isopropanol (results not shown) the corrections were significant only for the largest tested alkanes (beyond eicosane). In all cases, most of perturbed F molecules were mainly located close to elbows along the solute, where the possibilities of contact with solute were maximized. 4.4. Sampling of Pair Contact Energies Involving Bulky Plastic Additives. The sampling of intermolecular interactions between large hindered additives and themselves or with small

7272


Figure 6. Distribution of gyration radii (a, b) and volumes accessible to hydrogen atoms (c, d) for octadecane (a, c) and octacosane (b, d). Set A: conformers in vacuum; set B: conformers in methanol; set C: conformers in ethanol.

Figure 7. Pair distributions of F molecules: 2 with themselves, 9 with octadecane, b with octacosane. a, b: F ) methanol; c, d: F ) ethanol.

simulant molecules raised several difficulties that did not emerge as saliently for smaller or linear solutes.25 The main biases due to an overestimation of specific pair interactions are discussed in this subsection and are illustrated in Figure 10. Secondary polymeric antioxidants (or peroxide decomposer) are either phosphorous acid esters, which comprise Irgafos 168, or thioethers (also known as sulfides), which comprise Irganox PS802. Such molecules include a proton acceptor (a phosphorus

or sulfur atom), which is always located close to the center of mass because the symmetry of large substituents. Besides, the rigidity of either ether or ester bond was responsible for the creation of an open cavity, which tended to be oversampled by smaller contact molecules, as previously discussed.25 In the presence of a hydrogen bond donor simulant (alcohols, water), our standard pair contact algorithm, which proceeded by an expansion of both centers of mass until the van der Waals


7273

Figure 8. Density map of F in the main inertia plane of the solute i for: F ) methanol (a, b, c), F ) ethanol (d, e, f), i ) octadecane (c, f), i ) octacosane (a, b, d, e). Insets (a) and (d) detail the shape of the cavities when the solute is removed in (b) and in (e), respectively.

Figure 9. Calculation of nFoverlap neighboring i via the criterion C(r) defined in eq 11 for octadecane and octacosane in methanol (a, b) and ethanol (c, d).

surfaces of both molecules came in contact, tended to overestimate the contribution of the internal hydrogen bond with the simulant. This effect is depicted in Figure 10a for a conformer of Irganox PS802 in contact with methanol. The corresponding distribution of “thermalized” contact energies (see the ensemble average operator defined in eq 3) are plotted in Figure 10b. To

overcome such a difficulty, contraction methods were used instead of expansion ones and were combined with different directions of translation that were not collinear with centersof-mass of contact and seed molecules. In addition, replacing the sampling of many-molecules system by the sampling between a pair of molecules was shown to

7274


Figure 10. Biases associated to the sampling of contact energies involving bulky molecules: (a,b) Irganox PS802 and methanol, (c,d) Irganox 3114 and itself. The snapshots depicted in (a) and (c) have a minimal interaction energy. Sampled distributions of contact energies plotted in (b) and (d) are weighted with the corresponding Boltzmann factor at 313 K. Extraction and contraction methods are described in previous work.25

create artifacts when the molecular reference state of the system at the considered temperature was ordered. Such biases were particularly significant when the interactions between two large hindered phosphorous acid esters were sampled. Typical configurations of Irganox 3114 with very low potential energies are depicted in Figure 10c. Such configurations maximized the π-π interaction between their inner heterocyclic aromatic rings so that the left part distributions of pair contact energies at 313 K appeared highly noisy as plotted in Figure 10d. Such configurations were very unlikely as Irganox 3114 is crystalline at 313 K with a melting point ranging between 219 and 226 °C. The likely structure at 313 K of such hindered molecules would involve a stacking of convex regions by putting in contact its convex surface (bulging outward) with a concave one (bulging inward). Though they were associated to much lower contact energies, the few configurations, where two convex surfaces (e.g., Figure 10c) were in contact via their central aromatic rings (π stacking), were therefore discarded by introducing a posteriori a stacking weight, wA+B(ε) in eq 3, as follows: 〈εA+B〉T )

∫ ∫

+∞ p (ε)e-ε/kB · TwA+B(ε)εdε -∞ A+B +∞ p (ε)e-ε/kB · TwA+B(ε)dε -∞ A+B

(12)

Such additional refinements led to Flory interaction parameters in P, which were ranged between 0 and 3 (Table 4). Neglecting previous effects would yield much higher values, which would be unrealistic for plastics additives, which have by design a good solubility in P.

4.5. Predictions of Partitioning with Alcohols. Partition coefficients between alcohols and amorphous regions of polyethylene were calculated by means of eq 2. r-1 i,F was calculated from eq 4 using ViH values and the molar volume of F. The corrective term nFoverlap neighboring i was included for all molecules including an alkyl chain (e.g., n-alkanes, n-alcohols, Erucamide, Irganox 1076, Irganox PS802) on the basis of an extrapolation of the behavior observed from MD simulations octacosane and octadecane in the appropriate simulant. In practice, the correction was almost insignificant for segments including less than 16 carbons. Contact energies in F were derived conformers sampled from MD simulations at 313 K in F. They were finally “thermalized” by means of eq 3 at the corresponding experimental temperature, either 313 or 298 K. Corresponding χi,F and χi,P values are gathered in Table 4. Predictions were compared to reference data used in our previous work25 and to additional data generated by this study. The comparisons are plotted crudely in Figure 11 without applying any fitting procedure. The predictions were particularly satisfactory for plastic additives. For large n-alkanes, the predictions are poorly satisfactory in absence of correction but were highly confident when the overlapping effect was considered. Introduced new refinements improve drastically previous predictions of large molecules either flexible or not and made the predictions more robust for new solutes (e.g., surface agents, polymeric hindered antioxidants) and new food simulants (isopropanol). It is worth noticing that introduced corrections, highly nonlinear with chain length, were not derived from any chemical


7275

Table 4. Calculated Flory-Huggins Interaction Parameters solute decane undecane dodecane tridecane tetradecane pentadecane hexadecane heptadecane octadecane nonadecane eicosane docosane tetracosane octacosane decanol undecanol dodecanol tridecanol tetradecanol pentadecanol hexadecanol heptadecanol octadecanol nonadecanol eicosanol camphor diphenyl oxide diphenylmethane D-limonene D,L-menthol eugenol isoamyl acetate linalyl acetate phenylethyl alcohol BHT Chimassorb 81 Erucamide Irgafos 168 Irganox 1035 Irganox 1076 Irganox 245 Irganox 3114 Irganox PS802 stearic acid Tinuvin 326 a

χi,water

χi,methanol

χi,ethanol

/ / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / 30.45 ( 1.74c 33.84 ( 0.01a,c 30.13 ( 0.02a,c 44.17 ( 1.74a,c / 36.33 ( 0.12a,c / 45.36 ( 0.98c / / 30.97 ( 0.93a,c

5.89 ( 0.01 6.17 ( 0.005c 6.62 ( 0.001c 6.90 ( 0.005c 7.26 ( 0.01c 7.84 ( 0.02c 8.03 ( 0.05c 8.44 ( 0.05c 8.76 ( 0.05c 8.82 ( 0.05c 9.15 ( 0.02c 9.71 ( 0.02c 9.96 ( 0.6c 10.87 ( 0.7c 3.77 ( 0.05c 4.43 ( 0.05c 5.29 ( 0.01c 5.19 ( 0.04c 5.81 ( 0.03c 6.21 ( 0.02c 6.18 ( 0.06c 6.89 ( 0.06c 6.76 ( 0.03c 7.22 ( 0.06c 7.95 ( 0.02c -0.54 ( 0.005b -1.05 ( 0.005b 3.37 ( 0.02b 5.12 ( 0.03b 2.02 ( 0.03b 3.07 ( 0.02b 1.36 ( 0.02b 2.40 ( 0.03b 1.62 ( 0.02b 3.85 ( 0.003c 5.24 ( 0.15c 6.05 ( 0.01c 14.22 ( 0.1c 9.87 ( 0.61c 8.65 ( 0.13c 9.81 ( 1.09c 15.56 ( 0.14c 15.21 ( 0.14c 4.58 ( 0.03c 4.25 ( 0.09c c

χi,isopropanol

4.40 ( 0.01 4.65 ( 0.003c 5.06 ( 0.003c 5.31 ( 0.003c 5.64 ( 0.01c 6.16 ( 0.01c 6.35 ( 0.05c 6.72 ( 0.05c 7.01 ( 0.05c 7.09 ( 0.05c 7.38 ( 0.02c 7.91 ( 0.01c 8.21 ( 0.06c 9.08 ( 0.6c 2.49 ( 0.04c 3.03 ( 0.04c 3.71 ( 0.005c 3.77 ( 0.03c 4.35 ( 0.03c 4.80 ( 0.02c 4.92 ( 0.06c 5.62 ( 0.06c 5.71 ( 0.03c 6.04 ( 0.06c 6.68 ( 0.02c -0.81 ( 0.005b -0.78 ( 0.005b 2.73 ( 0.01b 3.80 ( 0.01b 1.55 ( 0.03b 2.59 ( 0.01b 0.71 ( 0.01b 1.77 ( 0.02b 1.30 ( 0.01b 3.12 ( 0.61c 4.40 ( 0.15a,c 5.38 ( 0.05a,c 12.77 ( 0.1a,c 9.49 ( 0.98c 7.14 ( 0.48a,c 8.48 ( 0.78c 14.54 ( 0.09c 14.15 ( 0.14c 3.78 ( 0.08c 4.50 ( 0.09a,c c

χi,P

3.28 ( 0.002 3.50 ( 0.001c 3.87 ( 0.001c 4.10 ( 0.01c 4.38 ( 0.02c 4.85 ( 0.01c 5.03 ( 0.04c 5.36 ( 0.04c 5.62 ( 0.04c 5.66 ( 0.04c 5.93 ( 0.01c 6.40 ( 0.01c 6.68 ( 0.05c 7.48 ( 0.5c 1.86 ( 0.04c 2.09 ( 0.04c 2.47 ( 0.001c 2.27 ( 0.03c 2.49 ( 0.03c 2.52 ( 0.01c 2.37 ( 0.05c 2.70 ( 0.05c 2.34 ( 0.02c 2.49 ( 0.05c 2.94 ( 0.01c -1.02 ( 0.005b -0.55 ( 0.005b 2.12 ( 0.01b 2.73 ( 0.01b 1.09 ( 0.03b 2.02 ( 0.04b 0.21 ( 0.01b 1.19 ( 0.02b 0.93 ( 0.01b 2.30 ( 0.1c 3.34 ( 0.15c 4.81 ( 0.10c 11.18 ( 0.1c 8.75 ( 0.1c 6.03 ( 0.1c 5.92 ( 0.1c 13.07 ( 0.06c 12.33 ( 0.1c 2.80 ( 0.02c 4.63 ( 0.1c c

0.3 ( 0.05c 0.30 ( 0.001c 0.31 ( 0.02c 0.30 ( 0.01c 0.23 ( 0.003c 0.24 ( 0.01c 0.37 ( 0.04c 0.28 ( 0.04c 0.44 ( 0.04c 0.55 ( 0.04c 0.68 ( 0.01c 0.11 ( 0.01c 0.39 ( 0.05c 0.35 ( 0.5c 0.79 ( 0.03c 0.77 ( 0.04c 0.69 ( 0.002c 0.68 ( 0.01c 0.40 ( 0.003c 0.46 ( 0.01c 0.54 ( 0.05c 0.47 ( 0.04c 0.70 ( 0.02c 0.72 ( 0.05c 0.76 ( 0.01c 0.67 ( 0.002b 0.42 ( 0.002b -0.05 ( 0.01b 0.04 ( 0.01b 0.56 ( 0.03b 0.43 ( 0.03b 0.49 ( 0.01b 0.20 ( 0.02b 0.62 ( 0.01b 0.51 ( 0.06c 1.12 ( 0.31c 1.75 ( 0.19c 0.12 ( 0.01c 1.38 ( 0.31c 0.15 ( 0.1c 2.44 ( 0.21c 2.57 ( 0.13c -0.31 ( 0.01c 1.56 ( 0.03c 1.49 ( 0.07c

Values used also in ternary mixtures. b Values at 298 K. c Values at 313 K.

specific interaction between small alcohols and aliphatic chains but from the coarse-graining process itself used to estimate the excess translational entropy. Indeed, after coarse-graining only the relative positions of center-of-mass (see Figure 1) were retained, and the excluded volume contribution tended to be overestimated as soon as interactions were described at a scale much larger than F molecules. For bulky solutes with almost spherical symmetry (e.g., BHT) or small solutes, this effect was negligible as the exposed surface to F was minimal. For asymmetric solutes and/or small liquids (e.g., methanol), the considered Flory-Huggins approximation lost by contrast some generalities, and the entropic term ceased to be invariant with the statistical blob size used to describe pair interactions. Mapping additional details from molecular simulations such as solute-F distribution functions solved this issue, while maintaining the complexity of the method at a reasonable level. As illustrated in Figure 8, flexible solutes in polar liquids could thus be captured at a coarse-grained level as spherical probes (i.e., blob) including a fraction of liquid molecules in it. The insufficient reliability of common experimental protocols, involving conventionally the long-time contact of small pieces of a plastic material with a large volume of F, made current predictions of partitioning, with computation time ranging

between few hours and one day, highly attractive and competitive for screening of possible contaminants of food products. As an illustration, our parallel implementations on a 3.2 GHz biquad core personal computer made it possible to assess accurately up to 20 activity coefficients per day (e.g., more than 10 partition coefficients). In particular, as low residual concentrations in P at equilibrium caused large uncertainties and biases in experimentally determined Ki,F/P values, the proposed theoretical method seemed particularly suitable to assess the partitioning of volatile substances encountered in polyolefins recycled for food contact or constituents of hot melt polyolefins used as adhesives. 4.6. Predictions of Partitioning with Ethanol-Water Mixtures. Prediction of partitioning in aqueous food simulant is both highly challenging and relevant for most food applications. Previous analyses neglected the presence of water in ethanol and ethanol 95% was therefore falsely considered as a pure substance. This subsection integrates explicitly the volume fraction of water in ethanol via eqs 6 and 8. The nonlinear behavior of ln Ki,F/P with the fraction of ethanol, φethanol, results from the virial expansion of the free energy for a ternary mixture. As expressed in eq 5, quadratic dependence of solute

7276


Figure 11. Comparisons between predicted (via eq 2) and experimental partition coefficients, Ki,F/P, between F ) alcohols and P ) amorphous regions of polyethylene: (a) F ) methanol, (b) F ) isopropanol, (c) F ) ethanol. Symbols are centered on published values. Corresponding confidence intervals are plotted as horizontal lines and are calculated with eq 1. They are therefore nonsymmetric, and they account for uncertainties of extraction yields and fraction of the initial i content, which do not participate to the apparent equilibrium. 1: decanea,b; 2: dodecanea,b; 3: tetradecanea,b; 4: hexadecanea,b; 5: octadecanea,b; 6: eicosaneb; 7: docosaneb; 8: tetracosaneb; 9: octacosaneb; R: dodecanolb; β: tetradecanolb; γ: hexadecanolb; δ: octadecanolb; a: camphora; b: diphenyl oxidea; c: diphenylmethanea; d: a a a a D-limonene ; e: DL-menthol ; f: eugenol ; g: isoamyl acetate ; h: linalyl acetatea; j: phenylethyl alcohola; A: BHTb; B: Chimassorb 81c,d; C: erucamidec,d; D: Irgafos 168c,d; E: Irganox 1035c,d; F: Irganox 1076b,d; G: Irganox 245c,d; H: Irganox 3114c,d; J: Irganox PS802c,d; K: stearic acidc,d; L: Tinuvin 326c,d. aExperimental values in methanol and ethanol from Baner et al. at 298 K.17 bExperimental values in methanol and ethanol from Vitrac et al. at 313 K.26 cExperimental values in ethanol 95% at 313 K, generated by this study. dExperimental values in isopropanol at 313 K, generated by this study.

chemical potential with φethanol was related to pair interactions between water and ethanol, χwater,ethanol, whereas a higher order dependence was caused either by ternary interactions, which were neglected, or by significant variations of χwater,ethanol with composition. Figure 12 plots the predicted partitioning of typical plastic additives in water-ethanol mixtures along with experimental values. The cubic trend was included or not to emphasize that the nonideal behavior of water-ethanol mixtures was needed to describesat least qualitativelyscollected experimental values and consequently that previously linear extrapolations used for regulatory purposes were supported neither by experimental evidence nor by theory. Sampled binary interaction parameters involved in eq 8 are tabulated in Table 4. Experi-

mental values in pure water were reported by Gandek et al.27 from desorption kinetics tough the experimental might be subjected to some partial hydrolysis. The quadratic and cubic models of partitioning predicted both a change in curvature at low fractions in ethanol. The solute chemical affinity for F tended either to decrease much more slowly or even to increase when more water was added to the mixture. As a result, the lowest chemical affinity for F could be obtained for a mixture, which might not be pure water. This condition would correspond to a fraction in ethanol between 0.1 and 0.3, where the behavior of the water-ethanol mixture deviates the most from the ideality. Indeed, the van-Laar m , is minimal for φethanol ≈ 0.3.35 This enthalpy, ∆Hwater+ethanol counterintuitive result was not possible to verify directly, because most of experiments were performed with a much higher fraction in ethanol. Only a global validation based on bounds values at φethanol ) 0 and at φethanol ) 0.95, along with their variations between 0.5 and 0.95, was possible. It is underlined that the latter range corresponds to values, where χwater,ethanol values were maximum (Figure 2a) and where the expected nonlinearities were therefore minimal. Assuming a constant χwater,ethanol value equal to its average in Figure 2a led to less accurate predictions but kept the ranking of solute according to their chemical affinity for F (Figures 12a and 12b). Adding nonidealities increased significantly the accuracy, while providing further insight on the effect of water fraction (c and d of Figure 12). A significant underestimation was only obtained for Erucamide, whose desorption isotherm exhibited also significant deviations to Henry’s law. For the other additives, it was shown that adding 5% v/v water in ethanol was accompanied by a sharp decrease in partitioning (Figure 12d). This trend would explain why predicted ln Ki,F/P values in pure ethanol would overestimate the experimental values in 95% ethanol mixtures (Figure 11c). The much less significant variation of ln Ki,F/P between φethanol ) 0.5 and φethanol ) 0.95 was well described. In pure water, the value of partitioning of Irganox 1076 was satisfactorily predicted, but the value of BHT was underestimated. 4.7. Compliance Testing of Materials in Contact with Dairy Products. From the sanitary point of view and because partitioning data in aqueous food simulants or in ethanol 50% are lacking in the literature, our theoretical framework and predictions could be used to check the reliability of some food simulant choices introduced in the EU regulation or to perform sanitary survey of plastic additives used in food packaging materials and plastics pipes. In this perspective, we outline simple calculations and magnitude orders associated to the migration of an ubiquitous antioxidant, such as Irganox 1076, into a typical dairy product packed in a HDPE bottle. From toxicological data, the maximum concentration authorized in food or food simulant, so-called specific migration limit, is set to 6 mg kg-1 (equivalently 6 mg L-1 for a mainly aqueous product) by the EU regulation 2002/ 72/EC. In agreement with EU regulation 2007/19/EC,19 the migration value can be estimated in ethanol 50% as an overestimate of the true migration in a real dairy product. According to Figure 12c, Ki,F/P is expected to be close to 10-3, which corresponds to an apparent partition coefficient -3 Kapp via eq 1 for a crystallinity of 70%. By i,F/P ) 3.3 × 10 considering a volume ratio between the packaging material and the food product of about LP/F ) 1/60 (e.g., bottle of 500 cm3, and 10 cm high and consisting in 0.1 mm thick) and a likely concentration in additive Ci,P|t)0 ) 103 mg L-1, the mass balance defined in eq 13, predicts a maximum


7277

Figure 12. Predicted and experimental partition coefficients, Ki,F/P, in F ) water-ethanol mixtures: (a,c) predictions without cubic terms (χwater,ethanol ) -2.04), (b,d) predictions with cubic terms. Insets of (c) and (d) are details of curves plotted in (a) and (b), respectively. aExperimental values generated by this study; bexperimental values generated by this study except values in pure water from ref 34; cexperimental value from ref 34; dvalues predicted with eq 5.

possible contamination of food at thermodynamical equilibrium of Ci,F|eq ) 2.8 mg L-1. Ci,F | eq )

Ci,P | t)0 1 1 + LP/F Kapp

(13)

i,F/P

The obtained value is about one-half of the specific migration limit. According to the expected diffusion coefficient of Irganox 1076 in HDPE,12 the equilibrium is reached in about one month at room temperature. The interested reader can find similar detailed derivations in refs 5 and 6. To protect the consumer, the so-determined value must overestimate the value that would be obtained in a real dairy product, with a typical fat content of fc ) 3% v/v (e.g., as in milk). In an oil in water emulsion, the apparent partition coefficient is a linear app app combination values in oil, Ki,fat/P and in water Ki,water/P : app app Kapp i,F/P ) (1 - fc)Ki,water/P + fcKi,fat/P

(14)

Assuming a value of Ki,fat/P close to 1 (similar chemical app affinity for fat and amorphous polyethylene) leads to Ki,fat/P ) 3.3. Choosing an underestimated value for pure water close to 2 × 10-6 (Figure 12c), one gets an underestimated partition app coefficient: Kapp i,F/P ≈ fcKi,fat/P, that is close to 0.10 and therefore much greater (30 times) than the one expected in 50% ethanol. Such considerations lead to an approximation of the contamination of 14 mg L-1, which exceeds the regulatory threshold. The discrepancy between estimates in ethanol 50% and in a true emulsion depends on how molecular interactions are calculated in a ternary mixture. In miscible water-ethanol mixtures, molecular contacts are averaged via eq 2, whereas in oil in water

emulsions, macroscopic occupation probabilities are averaged instead. Both points of views are not theoretically compatible and it remains unlikely that the recently proposed ethanol 50% may overestimate the partitioning with dairy products, except in cases of extremely low fat contents (e.g., 0.1% for the considered example). Similarly, the crude assessment of the contamination of food by phenolic antioxidants should be revised accordingly. 5. Conclusions Because most of previous methods of estimations of solute partition coefficients between packaging materials and liquids depended on the availability of experimental data, a general framework for their prediction have been developed in a series of two papers. As the relaxation of entangled polymer segments below their melting point is currently not tractable by molecular dynamics simulation as required in,19 it is underlined that the proposed approach is one of the rare predictive and tailored methods available for solutes with intermediate molecular weights (i.e., ranging between 100 and 1000 g mol-1). Activity coefficients in both phases are estimated through a generalized off-lattice Flory-Huggins formulation applied to semicrystalline materials, denoted P, and to interacting liquids simulating food products, denoted F, such as alcohols and water. The FloryHuggins approximation introduced three features, which participated to the attractiveness and satisfactory predictions of the method. First, the method samples independently contact energies between pair molecules and packing properties without requiring an explicit representation of the polymer or of the liquid simulant. The enthalpic contribution on excess chemical potential is accurately calculated for any molecule with technological interest (polymer, food simulant, additive or polymer

7278


residue) with an atomistic ab initio forcefield. Combined with a Monte Carlo sampling of both conformation and configurational spaces at atomistic scale, the nonpositional entropic contribution is naturally incorporated within the formulation. The enthalpic contribution is derived without additional cost for any close temperature (ca. ( 20 °C) by reweighting the sampled distribution with the corresponding Boltzmann factor. Second, the excess translational entropy, which was falsely neglected in previous studies applied to smaller solutes, is derived independently from geometrical considerations. Finally, the same multi-scale theory, combining atomistic scale calculations and coarse-grained approximations, can be applied indifferently to binary and ternary mixtures at solute infinite dilution. In well-characterized liquid mixtures (e.g., water-ethanol mixtures), excess solute chemical potential contributions in ternary mixtures can be straightforward extrapolated from the ones obtained in binary mixtures without repeating the sampling procedure. The application to bulky and flexible solutes showed however that refinements or extensions were required to apply the proposed off-lattice Flory-Huggins approach to migrants originating from food contact materials such as oligomers or plastic additives. The main difficulties arose from the sampling of interactions between molecules highly different in size and in shape. The configuration space was in addition significantly augmented by the presence of flexible molecules. Because the biases associated to the independent sampling of pair contact energies were analyzed in the first paper, this work focused on two specific extensions to increase mainly the accuracy and the robustness of predictions of partitioning with polar simulants. To validate introduced sophistications and assumptions, longterm isobaric molecular dynamics in F was used to provide conformers ensemble for flexible solutes, reference radial distributions and partial molar volumes. Independent validation was achieved by comparing with experimental values available in the literature and by values generated by this study. The first extension was devoted to improve the accuracy of the estimate of the size of cavity required to insert a large solute in F. In our previous work, the volume for the insertion of alkanes was guessed to be close to their hardcore volumes whereas it was envisioned close to or greater than their molar volumes for large hindered molecules. Isobaric MD simulations confirmed previous allegations by showing that different scenarios should be applied for flexible and stiff molecules. In addition, the volumes were found to be larger in water than in alcohols. To maintain the calculations tractable for a wide range of substances, a reliable estimate was proposed as the volume enclosed within the surface accessible to hydrogen atoms. Such a volume was shown to be almost independent of the considered conformational space (conformation in a gas phase or in a condensed phase) for most of tested solutes. In addition, pair distribution functions extracted from MD snapshots demonstrated that the local structure of liquid F molecules around large and flexible solutes were significantly modified so that an equivalent cell of liquid could be occupied by both solute atoms and F atoms. Our mean field approximation of the translational entropy was modified accordingly to account for the fractional number of F molecules overlapping the solute cavity. The correction was highly nonlinear: significant in methanol and/or for large alkyls chains involving more than 16 carbons. The second extension was less natural as it involved an intrinsic limitation of the Flory-Huggins approximation even when interactions were sampled at atomistic scale. Cooperative hydrogen bonding observed between more than two molecules

of water could not be sampled by pair interactions. Indeed, on the contrary to alcohols, the number of neighbors around one water molecule and related contact energies were correlated quantities. The ensemble average operator was modified accordingly to include the expected long lifetime tetrahedral coordination of liquid water at room temperature. This formal extension opened the way to the prediction of partitioning in water and in ethanol-water mixtures. In particular, it was found that a highly non linear variation of ln Ki,F/P should be expected with the fraction of ethanol. Such a behavior was related to the highly nonideal behavior of water-ethanol mixtures. To prevent additional artifacts, measured enthalpies of mixing of binary water-ethanol mixtures were introduced. A similar treatment could be envisioned for other mixtures used as food simulants, such as ethyl acetate - isopropanol proposed as a substitute of olive oil. As this work was initially motivated by the development of predictive methodologies to assess the contamination of food products by substances originating from plastic materials in contact, several technical recommendations were suggested to help the development of rapid screening methods in the perspective of sanitary surveys or regulatory issues. Future works will be focused on the validation of the current approach to other polymers (glassy polymers, copolymers, adhesives...) and to generate the first database of calculated activity coefficients and corresponding activity energies for plastic additives and residues. Acknowledgment We thank ACTIA and ANRT for their financial support. We also thank CTCPA for their technical help in assessing new experimental partition coefficients and Herve´ Perrier-Camby from the Institut National des Sciences Appliqueés in Lyon for processing the films. List of Symbols Roman symbols A,any molecule with a seed role B,any molecule with a contact role C(r),criterion defined in eq 11 Cî,k|eq,macroscopic concentration at equilibrium of solute i in k ) F,P as experimentally assessed (kg m-3) Ci,F|eq,concentration of solute i in liquid F at equilibrium (kg m-3) a Ci,P|eq ,concentration of solute i in the amorphous phase of polymer P at equilibrium (kg m-3) Ci,P|t)0,initial concentration of solute i in polymer P (i.e., before contact with L) (kg m-3) Ki,F/P,partition coefficient of solute i between F and the amorphous region of P app Ki,F/P ,apparent partition coefficient of solute i between F and P F,liquid (food simulant) Fk,component k ) 1,2 of mixture F. F1 and F2 were respectively chosen to be related to water and ethanol LP/F,volume ratio between P and F P,polymer PMFF neighboring i,excess equivalent potential mean force to insert a F molecule close to i R,universal gas constant (8.31 m2 kg s-2 K-1 mol-1) T,absolute temperature (K) Ubias|ε,weighting function of sampled pair contact energies ε VF+i,volume of equilibrated cells containing F and i molecules (m3) ViVdW,van der Waals volume of solute i (m3) VkM,molar volume of k ) i,F (m3)

Ind. Eng. Chem. Res., Vol. 49, No. 16, 2010 Vi,k,partial molar volume of i in k ) F,F1,F2 (m ) VkP,partial molar volume of k ) F1,F2 (m3) ViH,volume of solute i accessible to protons (m3) ai,extraction yield of solute i from polymer P (-). bi,relative amount of solute i, which is not “well-dispersed” in polymer P (-) c,polymer crystallinity (volume fraction of crystalline phase). fc,fat content expressed in volume fraction of F (-) gF,k(r),density in F molecules around k molecules (molecules · m3) hk,fraction of volume fraction of solute i and simulant F2 defined as hk ) φk/(φi + φF2) kB,Boltzmann’s constant (1.38 × 10-23 J K-1) i,solute index. nFoverlap neighboring i,fractional number of F molecules overlapping i nk,number of molecules of k ) i,F p,pressure (Pa). pA+B(ε),distribution of contact energies between A and B -1 ri,k ,number of k ) F,F1,F2 molecules to be rearranged when the solute i was displaced xethanol,molar fraction of ethanol zj+k,coordination number of the arrangement when j ) A,F,F1,F2,P,i is surrounded by molecules k ) B,F,F1,F2,P,I (J mol-1) P

3

Greek symbols V γi,k ,activity coefficient of solute i in k ) F, F1+F2,P m ∆Hwater+ethanol ,molar heat of mixing of water and ethanol εj+k,contact energy associated to the contact(s) between a molecule j ) A,F,F1,F2,P,i and a molecule k ) B,F,F1,F2,P,i (J mol-1) φk,volume fraction of k ) i,F1,F2 χk,j,Flory-Huggins parameter between k ) i,F1 and j ) F,F1,F2,F1+F2,P χi,F1,F2,ternary Flory-Huggins parameter between i, F1 and F2 Mathematical Operators 〈〉,averaging on the sampled configurational space 〈〉T,Boltzmann-weighted ensemble averaging on the sampled configurational space.

Literature Cited (1) European Community. Regulation (EC) No. 1935/2004 of the European Parliament and of the Council of 27 October 2004 on materials and articles intended to come into contact with food and repealing Directives 80/590/EEC and 89/109/EEC. Off. J. Eur. Union 2004, L338, 4-17. (2) European Community. Commission Directive 2002/72/EC of 6 August 2002 relating to plastic materials and articles intended to come into contact with foodstuffs. Off. J. Eur. Communities 2002, L220, 18-58. (3) European Commission. Substances listed in EU directives on plastics in contact with food. http://ec.europa.eu/food/food/chemicalsafety/foodcontact/ eu_substances_en.pdf; 2008. (4) European Community. Regulation (EC) No. 1907/2006 of the European Parliament and of the Council of 18 December 2006 concerning the Registration, Evaluation, Authorization and Restriction of Chemicals (REACH), establishing a European Chemicals Agency, amending Directive 1999/45/EC and repealing Council Regulation (EEC) No. 793/93 and Commission Regulation (EC) No. 1488/94 as well as Council Directive 76/769/EEC and Commission Directives 91/155/EEC, 93/67/EEC, 93/105/ EC, and 2000/21/EC. Off. J. Eur. Union 2006, L396, 1-849. (5) Vitrac, O.; Hayert, M. Design of Safe Packaging Materials Under Uncertainty. In New Trends Chemical Engineering Research; Berton, L. P., Ed.; Nova Science: New York, 2007; pp 251-292. (6) Vitrac, O.; Hayert, M. Risk Assessment of Migration from Packaging Materials into Foodstuffs. AIChE J. 2005, 51, 1080–1095. (7) Vitrac, O.; Challe, B.; Leblanc, J.-C; Feigenbaum, A. Contamination of Packaged Food by Substances Migrating from a Direct-Contact Plastic Layer: Assessment Using a Generic Quantitative Household Scale Methodology. Food Addit. Contam. 2007, 24, 75–94. (8) Vitrac, O.; Leblanc, J.-C. Exposure of Consumers to Plastic Packaging Materials: Assessment of the Contribution of Styrene from Yoghurt Pots. Food Addit. Contam. 2007, 24, 194–215. (9) Begley, T.; Castle, L.; Feigenbaum, A.; Franz, R.; Hirinchs, K.; Lickly, T.; Mercea, P.; Milana, M.; O’Brien, A.; Rebre, S.; Rijk, R.; Piringer,

7279

O. Evaluation of Migration Models that Might be Used in Support of Regulation for Food-Contact Plastics. Food Addit. Contam. 2005, 22, 73– 90. (10) Reynier, A.; Dole, P.; Feigenbaum, A.; Humbel, S. Diffusion Coefficients of Additives in Polymers. I. Correlation with Geometric Parameters. J. Appl. Polym. Sci. 2001, 82, 2422–2433. (11) Reynier, A.; Dole, P.; Feigenbaum, A.; Humbel, S. Diffusion Coefficients of Additives in Polyolefins. II. Effect of Swelling and Temperature on the D ) f(M) Correlation. J. Appl. Polym. Sci. 2001, 82, 2434–2443. (12) Vitrac, O.; Le´zervant, J.; Feigenbaum, A. Application of Decision Trees to the Robust Estimation of Diffusion Coefficients in Polyolefines. J. Appl. Polym. Sci. 2006, 101, 2167–2186. (13) Vitrac, O.; Hayert, M. Effect of the Distribution of Sorption Sites on Transport Diffusivities: a Contribution to the Transport of MediumWeight-Molecules in Polymeric Materials. Chem. Eng. Sci. 2007, 62, 2503– 2521, 9. (14) European Commission. Food contact materials. A practical guide for users of European directives. http://ec.europa.eu/food/food/chemicalsafety/ foodcontact/practical_guide_en.pdf (accessed September 2008). (15) Burman, L.; Albertsson, A.-C.; Ho¨glund, A. Solid-phase Microextraction for Qualitative and Quantitative Determination of Migrated Degradation Products of Antioxydants in an Organic Aqueous Solution. J. Chromatogr., A 2005, 1080, 107–116. (16) Piringer, O. G.; Baner, A. L. Partition Coefficients. In Plastic Packaging. Interactions with Food and Pharmaceuticals, 2nd ed.; WileyVCH: Weinheim, 2008; pp 89-121. (17) Baner, A. L.; Piringer, O. G. Prediction of Solute Partition Coefficients between Polyolefins and Alcohols using the Regular Solution Theory and Group Contribution Methods. Ind. Eng. Chem. Res. 1991, 30, 1506–1515. (18) Letcher, T. M. DeVelopments and Applications in Solubility; Royal Society of Chemistry: London, 2007. (19) Panagiotopoulos, A. Z. Direct Simulation of Fluid Phase Equilibria by Simulation in the Gibbs Ensemble: A Review. Mol. Simul. 1992, 9, 1– 23. (20) Van der Vegt, F. F. A.; Briels, W. J. Efficient Sampling of Solvent Free Energies in Polymers. J. Chem. Phys. 1998, 109, 7578–7582. (21) Widom, B. Some Topics in the Theory of Fluids. J. Chem. Phys. 1963, 39, 2808–2812. (22) Boulougouris, G. C.; Economou, I. G.; Theodorou, D. N. Calculation of the Chemical Potential of Chain Molecules Using the Staged Particle Deletion Scheme. J. Chem. Phys. 2001, 115, 8231–8237. ¨ zal, T.; van der Vegt, N. F. A. Fast-Growth (23) Hess, B.; Peter, C.; O Thermodynamic Integration: Calculating Excess Chemical Potentials of Additive Molecules in Polymer Structures Macromolecules. 2008, 41, 2283– 2289. ¨ zal, T.; Peter, C.; Hess, B.; van der Vegt, N. F. A. Modeling (24) O Solubilities of Additives in Polymer Microstructures: Single-Step Perturbation Method Based on a Soft-Cavity Reference State. Macromolecules 2008, 41, 5055–5061. (25) Gillet, G.; Vitrac, O.; Desobry, S. Prediction of Solute Partition Coefficients between Polyolefins and Alcohols using a Generalized FloryHuggins Approach. Ind. Eng. Chem. Res. 2009, 48, 5285–5301. (26) European Community. Commission directive 2007/19/EC amending Directive 2002/72/EC relating to plastic materials and articles intended to come into contact with food and Council Directive 85/572/EEC laying down the list of simulants to be used for testing migration of constituents of plastic materials and articles intended to come into contact with foodstuffs. Off. J. Eur. Union. 2007; Vol. L97, pp 50-69. (27) Bagnati, R.; Blanchi, G.; Marangin, E.; Zuccato, E.; Fanelli, R.; Davoli, E. Direct Analysis of Isopropylthioxanthone (ITX) in Milk by HighPerformance Liquid Chromatography/Tandem Mass Spectrometry. Rapid Commun. Mass Spectrom. 2007, 21, 1998–2002. (28) Sun, C.; Chan, S. H.; Lu, D.; Lee, H. M. W.; Bloodworth, B. C. Determination of isopropyl-9H-thioxanthen-9-one in Packaged Beverages by Solid-Phase Extraction Clean-Up and Liquid Chromatography with Tandem Mass Spectrometry Detection. J. Chromatogr., A. 2007, 1143, 162– 167. (29) Saiz, L.; Padro´, J. A.; Gua`rdia, E. Dynamics and Hydrogen Bonding in Liquid Ethanol. Mol. Phys. 1999, 97, 897–905. (30) Jorgensed, W. L. Transferable Intermolecular Potential Functions for Water, Alcohols, and Ethers. Application to Liquid Water. J. Am. Chem. Soc. 1981, 103, 335–340. (31) Mahoney, M. W.; Jorgensen, W. L. A Five-Site Model for Liquid Water and the Reproduction of the Density Anomaly by Rigid, Nonpolarizable Potential Functions. J. Chem. Phys. 2000, 112, 8910–8922. (32) Mahoney, M. W.; Jorgensen, W. L. Diffusion Constant of the TIP5P Model of Liquid Water. J. Chem. Phys. 2001, 114, 363–366.

7280


(33) Vitrac, O.; Mougharbel, A.; Feigenbaum, A. Interfacial Mass Transport Properties Which Control the Migration of Packaging Constituents into Foodstuffs. J. Food Eng. 2007, 79, 1048–1064. (34) Gandek, T. P.; Alan Hatton, T.; Reaid, R. C. Batch Extraction with Reaction: Phenolic Antioxidant Migration from Polyolefins to Water. 2. Experimental Results and Discussion. Ind. Eng. Chem. Res. 1989, 28, 1036– 1045. (35) Flory, P. J. Thermodynamics of high polymer solutions. J. Chem. Phys. 1941, 9, 660–661. (36) Huggins, M. L. Theory of Solutions of High Polymers. J. Am. Chem. Soc. 1942, 64, 1712–1719. (37) Flory, P. J. Thermodynamics of High Polymer Solutions. J. Chem. Phys. 1942, 10, 51–61. (38) Huggins, M. L. Derivation of Molecular Relaxation Parameters of an Isomeric Relaxation. J. Chem. Phys. 1942, 46, 151–153. (39) Neyertz, S. Tutorial: Molecular Dynamics Simulations of Microstructure and Transport Phenomena in Glassy Polymers. Soft Mater. 2006, 4, 15–83. (40) Vitrac, O; Gillet, G. An Off-Lattice Flory-Huggins Approach of the Partitioning of Bulky Solutes between Polymers and Interacting Liquids. Int. J. Chem. Reactor Eng. 2010, A6, 1–32. (41) Flory, P. J. Principles of Polymer Chemistry; Cornell University: Ithaca, 1953. (42) Pouchly, J.; Zivny, A.; Solc, K. Thermodynamics Equilibrium in the System Macromolecular Coil-Binary Solvent. J. Polym. Sci., Part C: Polym. Lett. Polym. Symp. 1968, 23, 245–256. (43) Young, T.-H.; Chuang, W.-Y. Thermodynamic Analysis on the Cononsolvency of Poly(vinyl alcohol) in Water-DMSO Mixtures through the Ternary Interaction Parameter. J. Membr. Sci. 2002, 210, 349–359. (44) Boyne, J. A.; Williamson, A. G. Enthalpies of Mixing of Ethanol and Water at 25°C. J. Chem. Eng. Data 1967, 12, 318.

(45) Lide, D. R., Ed. CRC Handbook of Chemistry and Physics; CRC: Boca Raton, 2005; pp15-41. (46) Cruzan, J. D.; Braly, L. B.; Liu, K; Brown, M. G.; Loeser, J. G.; Saykally, R. J. Quantifying Hydrogen Bond Cooperativity in Water: VRT Spectroscopy of the Water Tetramer. Science 1996, 271, 59–62. (47) Keutsch, F. N.; Saykally, R. J. Water clusters: Untangling the Mysteries of the Liquid, one Molecule at a Time. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 10533–10540. (48) Yang, J.; Peng, C.; Liu, H.; Hu, Y.; Jiang, J. A Generic Molecular Thermodynamic Model for Linear and Branched Polymer Solutions in a Lattice. Fluid Phase Equilib. 2006, 244, 188–192. (49) Coulier, L.; Kaal, E. R.; Tienstra, M.; Hankemeier, Th. Identification and Quantification of (Polymeric) Hindered-Amine Light Stabilizers in Polymers using Pyrolysis-Gas Chromatography-Mass Spectrometry and Liquid Chromatography-Ultraviolet Absorbance Detection-Evaporative Light Scattering Detection. J. Chromatogr., A 2005, 1062, 227–238. ´ .; Tena, M. T. Experimental Design Approach (50) Garrido-Lo´pez, A for the Optimisation of Pressurised Fluid Extraction of Additives from Polyethylene Films. J. Chromatogr., A 2005, 1099, 75–83. (51) Fried, J. R. Computational Parameters. In Physical Properties of Polymers Handbook, 2nd ed.; Mark, J. E., Ed.; Springer: Berlin, 2006; pp 59-82. (52) Rigby, D. Fluid Density Predictions using the COMPASS Force Field. Fluid Phase Equilib. 2004, 217, 77–87.

ReceiVed for reView July 6, 2009 ReVised manuscript receiVed May 3, 2010 Accepted June 4, 2010 IE9010595

Prediction of Partition Coefficients of Plastic Additives between

Recommend Documents