C G Pereira

R Lugo

J C Hemptinne

P Mougin

Guaiacol and its mixtures: new data and predictive models Part 1: Phase equilibrium

Keywords: phase equilibria, guaiacol, GC-PPC-SAFT, molecular simulation

In the present work, new experimental data of guaiacol mixture with methane were investigated. The results have been evaluated using several thermodynamic approaches.

Predictive calculations using the GC-PPC-SAFT (Group Contribution-Polar Perturbed Chain-Statistical Associating Fluid Theory) equation of state and Molecular Simulation using the AUA4 force field were performed. Data from literature for the binary systems of guaiacol with CO 2 , ethanol, octanol, acetone, butyl acetate and water were used to evaluate the thermodynamic models. The effect of the association scheme is discussed at length. Predictive phase equilibrium for systems containing small toxic compounds, such as hydrogen, carbon monoxide, hydrogen sulfide and ammonia were also performed. In GC-PPC-SAFT, two configurations of associative sites for guaiacol were considered. The predicted values showed to be consistent with new experimental data. The effect of conformational structure of guaiacol on phase equilibria was detected.

INTRODUCTION

The increase in the energy demand coupled with the gradual depletion of fossil fuels and the political and environmental requirements have motivated the search of new energy sources. Recently, the biomass conversion into hydrocarbon fuels has gained special attention due to its attractive characteristics in terms of availability and sustainability.

The use of bioresources in the production of energy and chemicals involves major environmental, economic and geostrategic implications. Indeed, creating new opportunities for bio-products in the field of energy, chemicals and materials should lead to a reduction of greenhouse gas emissions. The general objective behind this strategy is to ensure the protection of human beings and their environment while ensuring also the continuity of the economic activities with a reduced dependence on fossil energies. Unlike fossil fuels, bioresources are renewable, widely available and better distributed throughout the world, and their valorisation will certainly have a very positive impact at the regional level on the preservation of jobs and development of rural areas [START_REF] Pasty | Les débouchés non alimentaires des produits agricoles : un enjeu pour la France et l'Union Européenne[END_REF].

Chemical and petrochemical industries are now aiming at developing new green chemistry-based processes [START_REF] Dodds | Chemicals from Biomass[END_REF], [START_REF] Gallezot | Chapter 3: Process options for the catalytic conversion of renewables into bioproducts[END_REF], and lignocellulosic biomass (LCB) is called to play a major role as an alternative raw material. Indeed, the lignocellulosic material may be conditioned to obtain already "functionalised" molecules for specific applications, unlike hydrocarbons that require additional processing steps for the synthesis of specific chemicals.

The mixtures associated with the conversion of LCB include a large number of oxygenated molecules belonging to different families produced through the decomposition of the main constituents of LCB: cellulose, hemicellulose (polymers of sugars) and lignin (phenolic-derivatives polymers). The fast pyrolysis of LCB, for instance, produces solid, liquid (bio-oil) and gaseous fractions, in which the bio-oil is a source of several attractive compounds, including lighter components from different chemical families such as carboxylic acids, ethers or ketones and heavier compounds such as phenolic or furan derivatives. Recent studies have used guaiacol (or 2-methoxyphenol) as a model compound of bio-oil due to its molecular structure: it contains two oxygen-containing functional groups (phenolic -OH-and methoxy -OCH 3 -groups) that facilitate its repolymerization [4]- [6].

The processing and refining of these products is needed to improve their quality and to promote their use as a fuel. Several treatments can be considered: cracking, hydrocracking, decarbonylation, decarboxylation, hydrodeoxygenation, catalytic hydrotreatment, hydrogenation, others. However, in order to carry out studies on the technological and economic potential of LCB conversion process, it is necessary to develop predictive thermodynamic models with strong theoretical basis able to deal with multifunctional molecules where complex polar and associating interactions take place.

There have been several previous studies aiming at modelling systems with phenol derivatives. Lee et al [START_REF] Lee | Vapor-liquid equilibria for binary mixtures of carbon dioxide with 1,2-dimethoxybenzene, 2-methoxyphenol, or p-cresol at elevated pressures[END_REF] evaluated the phase equilibria of CO 2 + guaiacol and the data was correlated by different equations of state (EOS): the Soave-Redlich-Kwong, the Peng-Robinson, and the Patel-Teja EOS with the van der Waals mixing rules. They all require two binary interaction parameters for a good modeling and can therefore not be called predictive.

Hwang et al. [START_REF] Hwang | Isothermal vapor-liquid equilibria for mixtures composed of 1,2-dimethoxybenzene, 2-methoxyphenol, and diphenylmethane[END_REF] evaluated the phase equilibria of guaiacol + 1,2-dimethoxybenzene and guaiacol+ diphenylmethane with several versions of the UNIFAC. They were unable to predict the phase behavior of these systems, indicating that the interaction parameters of AC-OCH 3 binaries had to be re-evaluated. On the other hand, good data correlation were obtained using the Wilson, the NRTL, and the UNIQUAC models (not predictive). Espinosa et al. [9] applied the Group Contribution Associating Equation of State (GCA-EoS) to predict the phase equilibria of different phenol derivative systems; however, this work did not involve the mixtures evaluated here.

The GC-PPC-SAFT EoS has been used with some success on various mixtures of oxygenated compounds [START_REF] Nguyen-Huynh | Advances in Thermodynamics for Chemical Process and Product Design[END_REF]. The objective of this study is to evaluate the GC-PPC-SAFT model and molecular simulations using the AUA4 force field, to predict the phase equilibria of systems with guaiacol. The behavior of its mixtures with associative, non-associative, polar and apolar compounds was evaluated. Phase equilibria of mixtures of guaiacol with methane, carbon dioxide, ethanol, octanol, water, acetone, butyl acetate, n-hexadecane and other light components such as H 2, CO, H 2 S and NH 3 were evaluated in the present work. In order to evaluate the complexity of guaiacol, an analysis of its structure was also performed.

Experimental Section

Materials

For carrying out the phase equilibrium experiments, the following reagents were used: methane and carbon dioxide from Air Liquide (purity ≥ 99.998%), and guaiacol from SAFC (purity ≥ 98%), and guaiacol from Sigma Aldrich (purity ≥ 99.998%).

Apparatus and procedure

Experimental measurements of methane + guaiacol system at 443.15 K were performed in an equilibrium cell with a variable volume at constant pressure. The procedure and details of the equipment and its validation are provided in Pereira et al. [START_REF] Pereira | Phase Equilibria of systems containing aromatic oxygenated compounds with CH 4 , CO 2 , H 2 , H 2 S, CO and NH 3 : experimental data and predictions[END_REF].

Models used

In this work, the predictions of phase envelope were done using molecular-based GC-PPC-SAFT Equation of State as well as Monte Carlo Molecular Simulation with the AUA4 force field.

GC-PPC-SAFT

The GC-PC-SAFT EoS was originally proposed by Tamouza et al. [START_REF] Tamouza | Group contribution method with SAFT EOS applied to vapor liquid equilibria of various hydrocarbon series[END_REF] to model vapor-liquid equilibria of non-ideal mixtures. Later, Nguyen Huynh et al. [START_REF] Nguyenhuynh | Application of GC-SAFT EOS to polar systems using a segment approach[END_REF] extended the application of this model to polar-chain containing mixtures (the so-called PPC-SAFT EoS)

and applied it successfully to model the phase equilibrium behavior of different kind of chemical series, including alcohols, alkyl-benzenes, esters, ethers, ketones, aldehydes [START_REF] Nguyen-Huynh | Application of GC-SAFT EOS to polycyclic aromatic hydrocarbons[END_REF][15][START_REF] Nguyen-Huynh | Predicting VLE of heavy esters and their mixtures using GC-SAFT[END_REF] and small gas-containing systems [15,[START_REF] Nguyen-Huynh | Modeling Phase Equilibria of Asymmetric Mixtures Using a Group-Contribution SAFT (GC-SAFT) with a k(ij) Correlation Method Based on London's Theory. 2. Application to Binary Mixtures Containing Aromatic Hydrocarbons, n-Alkanes, CO 2 , N 2[END_REF][START_REF] Tran | Modeling VLE of H 2 + Hydrocarbon Mixtures Using a Group Contribution SAFT with a k(ij) Correlation Method Based on London's Theory[END_REF].

The

A mA A A A A       , (1)
where the first four terms relate to the non-polar interactions and follow the theory developed by Gross and Sadowski, the last incorporate the contribution of polar interactions that was obtained by extending the Gubbins and Twu [START_REF] Gubbins | Thermodynamics of Polyatomic Fluid Mixtures -1 Theory[END_REF] theory to chain molecules using the socalled "segment approach" of Jog & Chapman [START_REF] Jog | Application of Wertheim's thermodynamic Perturbation theory to dipolar hard sphere chains[END_REF], [START_REF] Jog | Application of Dipolar Chain Theory to the Phase Behavior of Polar Fluids and Mixtures[END_REF]. The details for the description of these terms can be found in the original papers.

Pure-Components parameters

In GC-PPC-SAFT, the segment parameters ( and ) and the chain parameter m of the molecule are calculated from group contribution parameters  k ,  k and R k using the following relations inspired by the Lorentz-Berthelot combining rules [START_REF] Tamouza | Group contribution method with SAFT EOS applied to vapor liquid equilibria of various hydrocarbon series[END_REF], [START_REF] Tamouza | Application to Binary Mixtures of a Group Contribution SAFT EOS (GC-SAFT)[END_REF]:

1 1 n groups groups k k k n n n k k        (2)      groups groups n k k n k k k n n 1 1   (3) n groups k k k1 m n R    (4)
where n k is the number of groups k in the molecule that is made of n groups different groups.

In the case of polar compounds (water, aromatic compounds, alkanols, ...) additional parameters shall be considered: for dipolar compounds (water, alkanols), the dipole moment  and dipole fraction   m x p .



, and for quadrupolar compounds (aromatic hydrocarbons), the quadrupole moment Q and its corresponding quadrupolar fraction   m x Q p . . These parameters are only relevant to polar and quadrupolar groups.

Finally, association interactions are considered using a specific term that requires additional information and parameters: an association scheme and a self-association energy

( AB 
) and volume ( AB  ). Self-association parameters are required for water and alkanols, and more generally, for all molecules containing both electron-donor and electron-acceptor sites.

All the pure compounds parameters for methane, carbon dioxide, ethanol, octanol, water and small molecules (H 2 , H 2 S, NH 3 , CO) were already determined in previous works, by the simultaneous regression of vapor pressure and pure saturated liquid density. These parameters are presented in Table 1.

For octanol and guaiacol, the segment parameters were calculated using group contributions using the equations ( 2) -( 4). The group parameters defined in Table 2 were used.

(Insert Tables 1 and2)

Mixture parameters

For mixtures containing compounds of different structures and polarities, combining rules must be defined. To improve the reproduction, two interaction parameters can be considered:

dispersive binary interactions parameters and cross-association parameters.

-Dispersive Binary Interactions Parameters: PC-SAFT EoS can be applied to mixtures using the van der Waals one-fluid model and the modified Lorentz-Berthelot mixing rules, through the potential parameters  ij and  ij between segments of molecules i and j.

  jj ii ij ij ε ε k = ε  1 (5) 2 ii jj ij σ +σ σ=    (6)
In order to keep the predictive characteristic of the GC-PPC-SAFT model, a predictive k ij was here considered (Eq 7, [15]) for all systems with guaiacol.

j i j i ij J J J J k    2 1 ( 7
)
where J i is the pseudo-ionization energy of the molecule i.

-Association Interactions Parameters: PC-SAFT EoS can be applied to mixtures using

The combination rule for the cross-association parameters is as follows:

2 B A AB        (8) A B B A         (9)
They were used without any correction in this work. Note the use of the indices  and  rather than i and j. This is to clearly label that the parameters refer to the groups  and  rather than to compounds. The small molecules of Table 1 can be considered single group molecules.

Molecular Simulation

To simulate CH 4 , CO 2 , ethanol, H 2 , CO, H 2 S and NH 3 , the molecular models used were those of Möller et al. [START_REF] Möller | Prediction of thermodynamic properties of fluid mixtures by molecular dynamics simulations: methane-ethane[END_REF], Harris and Young [START_REF] Harris | Carbon Dioxide's Liquid-Vapor Coexistence Curve And Critical Properties as Predicted by a Simple Molecular Model[END_REF], Darkrim et al. [START_REF] Darkrim | Monte Carlo simulations of nitrogen and hydrogen physisorption at high pressures and room temperature. Comparison with experiments[END_REF], Perez-Pellitero et al. [37], Kristof and Liszi [START_REF] Kristof | Effective Intermolecular Potential for Fluid Hydrogen Sulfide[END_REF] and Eckl et al. [START_REF] Eckl | An optimised molecular model for ammonia[END_REF], respectively. Details of these force fields are given in Table 3.

(Insert Table 3)

To simulate guaiacol, the AUA4 force field was used. To take advantage of the transferability feature of this force field, Lennard-Jones parameters and electrostatic charges distribution were directly taken from the parameterization of alkanes [START_REF] Ungerer | Optimization of the anisotropic united atoms intermolecular potential for n-alkanes[END_REF], aromatics [START_REF] Ahunbay | Optimized Intermolecular Potential for Aromatic Hydrocarbons Based on Anisotropic United Atoms. III. Polyaromatic and Naphthenoaromatic Hydrocarbons[END_REF],

alcohols [START_REF] Ferrando | Transferable Force Field for Alcohols and Polyalcohols[END_REF] and ethers [START_REF] Ferrando | A Transferable Force Field to Predict Phase Equilibria and Surface Tension of Ethers and Glycol Ethers[END_REF]. For intramolecular energy calculation, united-atoms separated by more than 3 bonds interact with a Lennard-Jones potential with the same parameters as was used for intermolecular energy calculation. Intramolecular electrostatic energy was computed using the dipole-dipole methodology introduced by Ferrando et al. [START_REF] Ferrando | Transferable Force Field for Alcohols and Polyalcohols[END_REF]. Two dipoles were introduced in the guaiacol molecule: the first one was defined by the 3 electrostatic charges of the ether function, and the second one by the 3 electrostatic charges of the alcohol function. A charge belonging to a given dipole interact with all the charges of the second dipole through the Coulomb's law. Bonded and non-bonded parameters of this force field are summarized in Tables 4 and5, respectively.

(Insert Tables 4 and5)

The choice of a mixing rule for cross-interactions between a center of force belonging to a small gas molecule and a center of force belonging to guaiacol was not obvious. From a previous study [START_REF] Pereira | Phase Equilibria of systems containing aromatic oxygenated compounds with CH 4 , CO 2 , H 2 , H 2 S, CO and NH 3 : experimental data and predictions[END_REF], it was shown that the cross-interaction between a center of force of a phenolic compound and a center of force of the various gases studied was better represented using the classical Lorentz-Berthelot rules, except for CH 4 which was better represented with the Kong rules [START_REF] Kong | Combining rules for intermolecular potential parameters. II. Rules for the Lennard-Jones (12-6) potential and the Morse potential[END_REF]. Thus, in this work either of them was adopted according to the gas molecule considered .

Pure guaiacol was simulated in the NVT Gibbs Ensemble [START_REF] Panagiotopoulos | Direct determination of phase coexistence properties of fluids by Monte Carlo simulation in a new ensemble[END_REF]. In this ensemble, the two phases in equilibrium were introduced in two separate boxes without an explicit interface.

When employed to study pure component equilibrium, this ensemble is applied at constant total volume of the two phases in order to respect the phase rule, the temperature and the total number of molecules being also imposed. A total number of 300 molecules was used. The average properties were computed during a production run lasting for 50 million Monte Carlo steps, one step corresponding to a single Monte Carlo move. Before each production run, a preliminary run lasting for 30 to 70 million Monte Carlo steps was carried out to achieve equilibrium. In the case of Lennard-Jones interactions, a spherical cut-off equal to half of the simulation box was used while the classical tail correction was employed [START_REF] Allen | Computer Simulation of Liquids[END_REF]. For long-range electrostatic energy, the Ewald summation technique was used, with a number of reciprocal vectors k equal to 7 in all three space directions and a Gaussian width α red equal to 2 in reduced units. The different Monte Carlo moves and their corresponding attempt probabilities used during the simulations were molecular translation (20%), molecular rotation (20%), regrowth with configurational bias [START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications[END_REF] (20%), transfer of a molecule between boxes with pre-insertion bias [START_REF] Mackie | Vapour-Liquid Phase Equilibria Predictions of Methane-Alkane Mixtures by Monte Carlo Simulation[END_REF] (39.5%), and volume change (0.5%). The amplitude of translations, rigid rotations and volume changes was adjusted during the simulation to achieve an acceptance ratio of 40% for these moves. The simulations were carried out using the GIBBS software [START_REF] Ungerer | Applications of Molecular Simulation in the Oil and Gas Industry[END_REF] jointly developed by IFP Energies nouvelles and the Laboratoire de Chimie Physique (CNRS-Université Paris-Sud). Note that vapor pressure was calculated using the Virial equation in the vapor phase, and the molar vaporization enthalpy h vap with the following relationship:

vap liq vap liq vap vap a vap liq vap liq U U V V h N P N N N N             (8)
where N a is the Avogadro number, P vap the pressure calculated in the vapor phase, and <U i >, <N i > and <V i > the average potential energy, total number of molecules and volume of phase i, respectively.

The simulations of binary mixtures gas + guaiacol were carried out in the NPT Gibbs Ensemble [START_REF] Panagiotopoulos | Phase equilibria by simulation in the Gibbs ensemble[END_REF]. A total number of 700 molecules was used. The equilibrium run lasted for 50 to 200 million Monte Carlo steps, according to the temperature simulated (the lowest temperatures requiring a longer equilibrium step). The production run lasted for 50 to 100 million Monte Carlo steps. Other simulation parameters were identical to those previously described for pure compound simulations.

Results and parameterization

Experimental Results

For the system methane + guaiacol there are no experimental data from literature. Thus, measurements were carried out within this work (see Table 6), and the results are compared with models in section 5.3 below.

(Insert Table 6)

Molecular Description

Several researchers debate on the strong inter and intramolecular HBs of phenolic compounds such as with guaiacol (2-methoxyphenol) [START_REF] Brune | Correlation between Adsorption and Small Molecule Hydrogen Bonding[END_REF]- [START_REF] De Heer | Poly Methoxy Phenols in Solution: O-H Bond Dissociation Enthalpies, Structures, and Hydrogen Bonding[END_REF]. This is because HBs affect the total enthalpy of specific interactions between solute and solvent, and then their physical behavior.

It is known also that pure guaiacol can actually exist as a mixture of two [START_REF] Varfolomeev | Pairwise Substitution Effects, Inter-and Intramolecular Hydrogen Bonds in Methoxyphenols and Dimethoxybenzenes. Thermochemistry, Calorimetry, and First-Principles Calculations[END_REF] or even four [START_REF] Dorofeeva | Molecular structures of 2-methoxyphenol and 1,2-dimethoxybenzene as studied by gasphase electron diffraction and quantum chemical calculations[END_REF] conformational isomers (see Figure 1). The formation of intramolecular HBs has a pronounced effect on the molecular structure and properties of molecules. However, in infinite dilution solutions, practically all molecules form intramolecular HB [START_REF] Dorofeeva | Molecular structures of 2-methoxyphenol and 1,2-dimethoxybenzene as studied by gasphase electron diffraction and quantum chemical calculations[END_REF].

(Insert Figure 1)

According to Heer et al. [START_REF] De Heer | Poly Methoxy Phenols in Solution: O-H Bond Dissociation Enthalpies, Structures, and Hydrogen Bonding[END_REF] phenols are known to form intermolecular hydrogen bonds that are able to have hydrogen-bond-accepting properties. Nevertheless, when they present an ortho methoxy group, they may create intramolecular hydrogen bonds as well, that are not linear. Two possibilities of conformers were evaluated considering the lowest energy conformations: one with an intramolecular hydrogen-bond and another with no hydrogenbond. In both cases, the hydroxyl and the methoxy group are in the plane of the phenyl ring.

The enthalpy difference between the two conformations, 4.4 kcal/mol, was defined as the intramolecular hydrogen-bond enthalpy. Several experiments were performed with the different isomers and the authors demonstrated the presence of strong intramolecular HB in 2-methoxyphenol. They also evaluated the strength of the HBs in the different conformations. The difference in energy of the conformers represents the energy of intramolecular HBs. This shows the complexity observed when studying systems with 2-methoxyphenol (guaiacol).

In order to better understand this effect, the liquid phase structure of pure guaiacol was investigated using molecular simulation. A Monte Carlo simulation of pure guaiacol was carried out in the monophasic NPT ensemble at 323 K, 0.1 MPa, and the equilibrated configurations were used to analyze hydrogen bond formation in the system. Figure 2 shows the intermolecular radial distribution function (rdf) between two hydroxyl oxygen atoms (noted O(H)), and between an ether oxygen atom (noted O) and a hydroxyl oxygen atom.

(Insert Figure 2)

The O(H)-O(H) rdf clearly exhibits a peak located at 3 Å which reveals the existence of an intermolecular hydrogen bond network in the system. It is interesting to see that the O(H)-O rdf does not exhibit such a peak. In other words, the hydroxyl hydrogen atom forms hydrogen bond only with the hydroxyl oxygen atom and not with the ether oxygen atom. This can be explained by the fact that the electrostatic partial charge of the ether oxygen is lower than the charge of the hydroxyl oxygen. A commonly used geometric criterion to define the presence of a hydrogen bond is the oxygen-oxygen distance that must range between the beginning of the rdf (2.5 Å) and the first local minimum (3.5 Å), and simultaneously the O(H)-H--O(H) angle that must range between 140 and 180° [START_REF] Jorgensen | Optimized intermolecular potential functions for liquid alcohols[END_REF], [START_REF] Chen | Monte Carlo Calculations for Alcohols and Their Mixtures with Alkanes. Transferable Potentials for Phase Equilibria. 5. United-Atom Description of Primary, Secondary, and Tertiary Alcohols[END_REF], [START_REF] Ferrando | Transferable Force Field for Alcohols and Polyalcohols[END_REF]. With these geometric criteria, we found a monomer fraction of 0.79 (fraction of molecules not involved in a hydrogen bond) and an average number of hydrogen bond per molecule of 0.23. This is considerably lower than the values obtained for phenol with the same force field (average number of hydrogen bond of 1.95 and a monomer fraction close to 0 at room temperature) [START_REF] Ferrando | Transferable Force Field for Alcohols and Polyalcohols[END_REF]. This phenomenon could be explained by analyzing intramolecular radial distribution functions. Figure 3 shows the intermolecular distance histogram between the ether oxygen atom O and the hydroxyl hydrogen H atom.

(Insert Figure 3)

It can be seen that all the O-H distances range from 2.1 to 2.4 Å. That means that all guaiacol molecules in the bulk are in anti-syn conformation as defined in Figure 1. Such a proximity of the ether oxygen atom and the hydrogen atom can result in steric hindrance, which could explain the difficulty to form intermolecular hydrogen bond compared to phenol for example. Based on this analysis, different conformational schemes of association can be used, and two of them were considered on the GC-PPC-SAFT predictions to be evaluated in this study.

Guaiacol pure-Components parameters

Auger [START_REF] Auger | Equilibres de phases des fluides oxygénés issus du traitement thermochimique de la biomasse lignocellulosique[END_REF] has investigated this compound with the GC-PPC-SAFT model, but it appears that the behavior depends strongly on the association scheme that is selected. The association term in the model cannot describe explicitly the intra-molecular bonds. These can only be taken into account by the fact that they des-activate possible bond-formation sites. This is the approach considered here. Figure 2 shows the groups used: 1(OH); 4(CH_aro); 2(C_aro); 1(-O-); 1(CH 3 ); 1(aromatic ring). The parameters proposed by this author are given in Table 2. The association scheme is that which is shown as scheme 2 in Figure 4. The global deviation using this parameterization is 18.9% on vapor pressure and 5.37% on molar volumes. Using a predictive approach, we can also try to compute the vapor pressure of guaiacol using the published parameters by Nguyen-Huynh et al [START_REF] Nguyen-Huynh | Advances in Thermodynamics for Chemical Process and Product Design[END_REF], also provided in Table 2. In this case, the average deviation on vapor pressure is 23% (2.5% on volume). In this work, we didn't consider this to be accurate enough and tried to improve by refitting some parameters.

As will be seen below, results are very sensitive to the association scheme. For this reason, two association schemes were considered in this work: i) Scheme 1 with 3 associating sites, one positive and one negative on the group (OH) and one negative on the (O) of the ether function (notice that this scheme implicitly takes into account the steric hindrance effect of intramolecular association, i.e. less associating sites in the molecule); and ii) Scheme 2 with 6 associating sites, one positive and two negative on the group (OH), two negative on the group (O) of the ether function and one negative on the aromatic ring. The two schemes are represented in Figure 4.

(Insert Figure 4)

In order to improve the predictions, the association parameters as well as chain length contribution of the OH and the O groups were regressed to fit pure compound data (vapor pressure and liquid volume). The predicted values were compared with experimental data from literature (DIPPR validated [START_REF] Rowley | Design Institute for Physical Properties Data[END_REF]), in terms of vapor pressure, liquid volume and enthalpy of vaporization for the pure compound, and also vapor-liquid equilibrium (VLE) for the system CO 2 + guaiacol [START_REF] Lee | Vapor-liquid equilibria for binary mixtures of carbon dioxide with 1,2-dimethoxybenzene, 2-methoxyphenol, or p-cresol at elevated pressures[END_REF]. Table 7 presents these parameters adjusted considering the vapor pressure and density of guaiacol.

(Insert Table 7)

All other group parameters originate from Table 2, including those of the aromatic ring. Not unexpectedly, the larger the number of sites, the lower the values of the association parameters.

With these parameters, the vapor pressure, liquid volume and vaporization enthalpy of guaiacol were calculated accurately using GC-PPC-SAFT. Figures 5, 6 and7 shows the deviations between the predicted values (GC-PPC-SAFT and molecular simulation) and experimental data [START_REF] Rowley | Design Institute for Physical Properties Data[END_REF]. Table 8 presents the values obtained by Molecular Simulation, and Table 9 presents the deviations for each model. It is observed that for the pure compound, the two approaches (MS and GC-PPC-SAFT) yield rather low deviations. The schemes 1 and 2

provided almost identical results: in fact, almost no difference is observed between the pure component results of the two association schemes. As will be seen later, the differences appear for the mixtures with hydrogen-bonding compounds.

(Insert Tables 8 and9)

(Insert Figures 5, 6 and 7)

Small deviations were also observed in molecular simulation predictions: < 5% for vapour pressure and vaporization enthalpy and only 1.3% in terms of liquid volume deviations.

Discussion: mixtures computations

The difficulty to work with guaiacol is related to its molecular structure and strong and complex interactions (inter-and intramolecular hydrogen bonds -HBs), which can occur in antagonistic manner: aiding or opposing each other [START_REF] Brune | Correlation between Adsorption and Small Molecule Hydrogen Bonding[END_REF]- [START_REF] Varfolomeev | Pairwise Substitution Effects, Inter-and Intramolecular Hydrogen Bonds in Methoxyphenols and Dimethoxybenzenes. Thermochemistry, Calorimetry, and First-Principles Calculations[END_REF] (more details are discussed in section 4.2). We will investigate systems that could help us understand the phenomena at hand. Therefore, we shall focus on systems with water (extremely HB), with ethanol and noctanol (moderately HB), acetone and, butyl acetate (only cross HB), n-hexadecane (no HB).

We shall also have a look at some gases as methane and CO 2 . At the end, we shall also look at some other systems for which no data exist.

For the predictions of phase equilibria of guaiacol systems with GC-PPC-SAFT, the predictive binary interaction parameters (k ij ) were used (eq. 7). For that, the pseudo-ionization energy of each component (J i ) was considered and applied in Eq. ( 3). The pseudo-ionization energy of guaiacol (7.76) was assumed to be equal to the pseudo-ionization energy of benzene (i.e. not computed from group contributions). For the other components, the pseudoionization energies from literature [15], [START_REF] Nguyen-Huynh | Modeling Phase Equilibria of Asymmetric Mixtures Using a Group-Contribution SAFT (GC-SAFT) with a k(ij) Correlation Method Based on London's Theory. 2. Application to Binary Mixtures Containing Aromatic Hydrocarbons, n-Alkanes, CO 2 , N 2[END_REF], [START_REF] Tran | Modeling VLE of H 2 + Hydrocarbon Mixtures Using a Group Contribution SAFT with a k(ij) Correlation Method Based on London's Theory[END_REF], [START_REF] Nguyen-Huynh | Modeling Liquid-Liquid and Liquid-Vapor Equilibria of Binary Systems Containing Water with an Alkane, ans Aromatic Hydrocarbon, an Alcohol or a Gas (Methane, Ethane, CO 2 or H 2 S), Using Group-Contribution Polar Perturbed-Chain Statistical Associating Fluid Theory[END_REF], [START_REF] Nguyen-Huynh | Modélisation thermodynamique de mélanges symétriques et asymétriques de composés polaires oxygénés et/ou aromatiques par GC-SAFT[END_REF] were used. For ammonia, it was found in a previous work that J NH3 = 0 [START_REF] Pereira | Phase Equilibria of systems containing aromatic oxygenated compounds with CH 4 , CO 2 , H 2 , H 2 S, CO and NH 3 : experimental data and predictions[END_REF], [START_REF] Pasty | Les débouchés non alimentaires des produits agricoles : un enjeu pour la France et l'Union Européenne[END_REF], which means that k ij for the system NH 3 + guaiacol is zero.

Systems with water and hexadecane

In a first instance, we can evaluate the behavior of the model on binary systems with water and with n-hexadecane. Interestingly, both systems guaiacol + water (Figures 8 and9)

and guaiacol + n-hexadecane (Figure 10) show liquid-liquid phase splits.

(Insert Figures 8 and 9)

Figures 8 and 9 both show the very nice behavior of scheme 1 for computing the water + guaiacol phase diagram. The average deviations on water concentration in the organic phase is 7%, while that of guaiacol on the aqueous phase is 29%. On the other hand, scheme 2 is unable to describe the liquid-liquid phase split, and the bubble pressure is clearly underestimated.

The reason for the strong difference in behavior between the two association schemes is most probably related to the large number of possible cross-associations that are offered by scheme 2. The more cross-association, the better the solubility behavior.

On the opposite extreme, it is possible to observe a liquid-liquid equilibrium between a fully non-associating alkane, n-hexadecane, and guaiacol. As observed in Figure 10, here both association schemes provide identical, quite reasonable results. Please note that these calculations are fully predictive, and that no SAFT model is capable of describing correctly liquid-liquid critical phenomena. We therefore focus on the quality of the phase composition away from the critical end point.

(Insert Figure 10)

Systems with moderately hydrogen bonding systems

Here we investigate the case of guaiacol with alcohols: ethanol and n-octanol, and with acetates or esters: acetone and butyl acetate.

The results were different for systems containing alcohol. Figures 11 and12 show the phase envelopes for the systems ethanol + guaiacol and octanol + guaiacol with predictions using GC-PPC-SAFT at different temperatures. For these systems, the predictions using Scheme 1 provided systematically better results.

(Insert Figure 11)

For the system ethanol + guaiacol (Figure 11), a deviation of 2.5% (disregarding data at 1.9 × 10 -3 MPa) was observed for the predicted values obtained when the Scheme 1 was used in GC-PPC-SAFT. The Scheme 2 tends to overestimate the solubility of ethanol in guaiacol (~ 7%). This overestimation was also observed for the system octanol + guaiacol.

This phenomenon, which corresponds to an underestimation of the bubble pressure, can easily be explained by stronger cross-association compared to that in scheme 1. A GC description of ethanol proposed by [START_REF] Rowley | Design Institute for Physical Properties Data[END_REF] was also evaluated, but it did not provide any further improvement with respect to the molecular description used to obtain Figure 11.

(Insert Figure 12)

The acetone + guaiacol system (Figure 13) is peculiar in the sense that it shows a slightly negative deviation from ideality: the bubble pressure data are slightly smaller than the straight line connecting the pure component vapor pressures. Please note that guaiacol is solid at 290.15K, but an extrapolation of the DIPPR vapor pressure correlation provides the low value of 6.3Pa; also, a 6.5% overestimation of the acetone vapor pressure by GC-PPC-SAFT is observed, which happens to be the largest deviation of the full temperature range. This deviation is quite acceptable for a predictive model. Considering that cross-association generally yields such negative deviations from ideality, it is quite unexpected that the guaiacol association scheme with the most hydrogen-bonding sites yields a larger vapor pressure. Yet, we should recognize that the behavior of such systems is driven by a balance between selfassociation (which will yield positive deviation from ideality) and cross-association (that yields negative deviation). In addition, we can also observe that both guaiacol association schemes have only one proton-donor, and that acetone has one proton acceptor site. Hence, increasing the number of proton acceptor sites on the guaiacol results in an increase of selfassociation without increasing the cross-association. This clearly explains why the scheme 2 leads to a more positive (or less negative) deviation from ideality compared to scheme 1.

We observe again that scheme 1 is closer to the true behavior. If we want to improve the representation of such a system, the cross-association parameters must probably be adjusted.

(Insert Figure 13)

In the guaiacol + butyl-acetate system that is shown in Figure 14, it is impossible to visualize whether a negative or a positive deviation from ideality exists. Yet, the same conclusion as for the acetone + guaiacol mixture can be drawn.

(Insert Figure 14)

Systems with gases

The mixtures with methane or CO 2 are shown in the Figures 15 and16. For methane, we can observe the same phenomenon as with hexadecane: both schemes provide identical results.

(Insert Figure 15)

An example for the phase envelopes for the CO 2 + guaiacol system using GC-PPC-SAFT and molecular simulations are presented in Figures 10. The same phenomenon is observed as with acetone, that is that scheme 2, with a large number of sites, results in larger bubble pressures (underestimation of solubilities). This can be explained in the same way as for acetone: the more guaiacol can self-associate (scheme 2), the larger its tendency to feature positive deviations from ideality. This is why the bubble pressures are larger.

Espinosa et al [9] also evaluated the CO 2 + guaiacol system at the same conditions of temperature. Their results were very similar to our results for predictions with GC-PPC-SAFT using Scheme 1. The difference between the GCA EoS used by Espinosa et al [9] and the one used in this work is in the additional terms for the formation of chain and multipolar interactions.

All the investigated systems with their standard deviations are presented in Table 10.

(Insert Figure 16)

Final discussion on the predictive quality of the GC-PPC-SAFT model with guaiacol

According to Table 10, the predictions made with GC-PPC-SAFT showed good agreement with experimental data, and deviations were lower than those obtained with molecular simulation predictions (Figure 16). Even though binary interactions parameters k ij were used in the GC-PPC-SAFT model, these parameters were completely predicted and the values obtained were small, keeping the predictive nature of the model.

In all cases, it is observed that scheme 1 provides better predictions.

(Insert Table 10)

The analysis that is shown here illustrates following interesting features:

-Both for pure guaiacol as for mixtures with non-associating species, the association scheme has no effect on the final result. This was exemplified by mixtures with methane and with hexadecane. In all these cases, the model was capable to reproduce, both qualitatively as quantitatively the vapor-liquid and liquid-liquid phase envelopes.

-For mixtures with compounds that can only cross-associate, i.e. when it contains only a proton acceptor, the effect of having a large number of proton acceptor sites on the guaiacol molecule (scheme 2) only results in increasing self-association without allowing for an increased cross-association (both schemes feature only one proton donor site). As a consequence, the mixture shows a stronger positive deviation from ideality for scheme 2 than for scheme 1.

-For mixtures where both compounds contain proton-donors and proton-acceptors (alcohols, and water), the fact of increasing the number of proton acceptors on the guaiacol molecule results in an increase of cross-association. As a consequence, the mixture behavior shows a decrease in its positive deviation from ideality. This can be observed on the bubble curves of mixtures with both for n-octanol and ethanol, and has an even stronger effect on mixtures with water where the scheme 2 is unable to catch the liquid-liquid equilibrium.

In all comparisons with experimental data, scheme 1 showed better predictive results (Table 10). As a conclusion, we recommend using the smaller number of association sites to represent the molecule of guaiacol, i.e, scheme 1.

Cross-over phenomenon

Another interesting point is shown in the calculations: the cross-over behavior of the CO 2 + guaiacol system indicated by literature is also predicted by GC-PPC-SAFT predictions (Figure 17). It is observed that at low pressure, increasing temperature results in a larger solubility of guaiacol in the supercritical CO 2 (due to the increase in guaiacol vapor pressure), while at high pressure, the opposite is true (because higher temperature means a lower CO 2 density which in turn yields a lower solubility). This behavior is common for systems in supercritical conditions. Similar phenomena are observed in systems as 1,2 dimethoxybenzene + CO 2 and p-cresol + CO 2 [START_REF] Lee | Vapor-liquid equilibria for binary mixtures of carbon dioxide with 1,2-dimethoxybenzene, 2-methoxyphenol, or p-cresol at elevated pressures[END_REF], other systems [START_REF] Raeissi | Application of double retrograde vaporization as an optimizing factor in supercritical fluid separations[END_REF][START_REF] Lee | Solubilities of chlormezanone, metaxalone and methocarbamol in supercritical carbon dioxide[END_REF][START_REF] Leitão | Anacardium occidentale L. leaves extraction via SFE: Global yields, extraction. kinetics, mathematical modeling and economic evaluation[END_REF][START_REF] Pereira | Extraction of indole alkaloids from Tabernaemontana catharinensis using supercritical CO 2 +ethanol: an evaluation of the process variables and the raw material origin[END_REF][START_REF] Wolbach | Using molecular orbital calculations to describe the phase behavior of cross-associating mixtures[END_REF].

(Insert Figure 17)

Other Systems

During the processing and treatment of bio-oil, other important compounds are produced, and in the absence of data for systems with toxic compounds, we performed predictive calculations for the mixtures of guaiacol + toxic gases (H 2 S, NH 3 , H 2 , CO). Figures 18-21 present the predictions for these systems using GC-PPC-SAFT and molecular simulation. It is possible to observe that the model is consistent in the reproduction of phase envelopes. For H 2 S, H 2 and CO the predictions using GC-PPC-SAFT, lead to higher solubilities of these gases in guaiacol than predicted by molecular simulation. The opposite was observed for NH 3 . Figure 22 provides the same type of information in terms of Henry constants (here, it was not possible to compare with molecular simulation).

(Insert Figures 18 to 22)

It is observed that for H 2 S and NH 3 , the Henry constant is much smaller than for CO and H 2 , this means that the solubility is much larger. In addition, the slope of these two sets of gases is opposite. This was already observed in our previous paper [71] where solubility was investigated in phenol-and furan type solvents: a gas that features association sites will dissolve better, and its solubility decreases with temperature (increasing Henry constant); on the opposite, the solubility of very light, non hydrogen-bonding gases is driven by free volume available in the solvent [72], which increases with temperature: their solubility increases with temperature, or their Henry constant decreases.

Conclusions

In this work, new experimental data of the system guaiacol-methane were presented and compared with predictive models (GC-PPC-SAFT and Molecular Simulation).

Despite the lack of experimental data, an examination of the different systems makes it possible to apprehend the predictive capacity of the model.

The use of a SAFT-type model implies the definition of an association scheme. The molecular structure of guaiacol was discussed in that perspective: it appears that intramolecular association may de-activate some proton-donor or acceptor sites. In order to mimic this, two possible schemes have been proposed and analyzed. They indicated no difference for the pure compound properties or for mixtures with non-associating species. In the presence of species featuring proton donors and proton acceptors, the larger the number of sites, the more cross-association is observed. Yet, in the presence of species that only contain proton acceptors, the scheme with the larger number of sites shows more self-association and therefore a larger positive deviation (or less negative) from ideality.

As a general conclusion, for the guaiacol, the lower number of association sites (scheme 1) provided better agreement with experimental results, including liquid-liquid equilibrium with water.

Using the model that was thus developed, a number of observations could be made related to the cross-over phenomenon in supercritical solvents and related to the solubility of gases in guaiacol.

Acknowledgements

The authors wish to acknowledge the financial support of the Tuck foundation chair for thermodynamics for biofuels, and to CNPQ (246524/2012-0) for financial support..

Acknowledgment

The authors are grateful to the Tuck foundation and to CNPQ for financial support. Figure 13: Bubble pressure of the guaiacol + acetone binary system at 290.15K using GC-PPC-SAFT with the two guaiacol schemes. Experimental data from [START_REF] Weissenberger | On the Knowledge of Aromatic Molecular Compounds. XVI. Bivalent Phenols and Their Ethers[END_REF]. Ideal mixture assumption between the extrapolated vapor pressures of the two pure components. 6).

Figure 16: Phase diagram of the system CO 2 + guaiacol at 393.15 K using GC-PPC-SAFT with the two guaiacol schemes and Molecular Simulation. Experimental data from [START_REF] Lee | Vapor-liquid equilibria for binary mixtures of carbon dioxide with 1,2-dimethoxybenzene, 2-methoxyphenol, or p-cresol at elevated pressures[END_REF].

NOMENCLATURE

  k ij = binary interaction parameter ε = segment-segment interaction energy  = diameter of a segment m = number of segments μ= dipole moment x p μ m = number of dipoles on the molecule Q= quadrupole moment x p Q m = number of quadrupoles on the molecule n sites = number of association sites ε AB /k = association energy κ AB = association volume R 1-4 = chain fraction T C = critical temperature P C = critical pressure ω= acentric factor

Figure 1 :

 1 Figure 1: Conformational isomers of guaiacol [56].

Figure 2 :

 2 Figure 2: Intermolecular radial distribution function O(H)-O(H) (solid line) and O(H)-O (dashed line) in liquid guaiacol at 323.15 K, 0.1 MPa.

Figure 3 :

 3 Figure 3: Intermolecular histogram O-H in liquid guaiacol at 323 K, 0.1 MPa.

Figure 4 :

 4 Figure 4: Association schemes of for guaiacol.

Figure 5 :

 5 Figure 5: Comparison of vapor pressure data of guaiacol from literature [60] and calculated using Molecular Simulation (□) and GC-PPC-SAFT using Scheme 1 (∆) and Scheme 2 (○).

Figure 6 :

 6 Figure6: Comparison of liquid volume data of guaiacol from literature[START_REF] Rowley | Design Institute for Physical Properties Data[END_REF] and results using Molecular Simulation (□) and GC-PPC-SAFT using Scheme 1 (∆) and Scheme 2 (○).

Figure 7 :

 7 Figure 7: Comparison of vaporization enthalpy data of guaiacol from literature [60] and results using Molecular Simulation (□) and GC-PPC-SAFT using Scheme 1 (∆) and Scheme 2 (○).

Figure 8 :

 8 Figure8: Phase diagram of the system water + guaiacol at 313.6, 334.4, 353.9 and 372 K using GC-PPC-SAFT with the two guaiacol schemes. Experimental data from[START_REF] Hahn | Determination of the VLE in the systems water -guaiacol and butyl acetateguaiacol[END_REF].

Figure 9 :

 9 Figure9: Guaiacol + water liquid-liquid phase equilibrium at 0.1 MPa, using GC-PPC-SAFT with the Scheme 1 in the predictions (solid line). The use of Scheme 2 didn't reveal any liquid-liquid phase split. Experimental data (*) from[START_REF] Balina | Kalininskij Gosudarstvennyj Universite, Investigation of the hydrazine -guaicol -water system[END_REF] .

Figure 10 :

 10 Figure 10: Phase diagram n-hexadecane + guaiacol system at 0.1 MPa, using GC-PPC-SAFT with the two guaiacol schemes. Experimental data from [63].

Figure 11 :

 11 Figure 11: Phase diagram of the system ethanol + guaiacol at 290.15 K using GC-PPC-SAFT with the two guaiacol schemes. Experimental data from [64].

Figure 12 :

 12 Figure 12: Phase diagram of the system octanol + guaiacol at 433.15, 448.15 and 463.15 K using GC-PPC-SAFT with the two guaiacol schemes. Experimental data from [8].

Figure 14 :

 14 Figure14: Phase behavior description of the butyl-acetate + guaiacol binary system at 101.3 kPa using GC-PPC-SAFT with the two guaiacol schemes. Experimental data from[START_REF] Li | Isobaric (vapour + liquid) equilibria of binary systems containing butyl acetate for the separation of methoxy aromatic compounds (anisole and guaiacol) from biomass fast pyrolysis oil[END_REF].

Figure 15 :

 15 Figure 15: Phase diagram of the system methane + guaiacol 443.15 K using GC-PPC-SAFT with the two guaiacol schemes and Molecular Simulation. Experimental data (this work, Table6).

Figure 17 :

 17 Figure17: Cross-over behavior for the CO 2 + guaiacol system using predicted values from GC-PPC-SAFT, using Scheme 1.

Figure 18 :

 18 Figure 18: Phase diagram of the system H 2 S + guaiacol at 373.15 K and 573.15 K using GC-PPC-SAFT (Scheme 1) and Molecular Simulation.

Figure 19 :

 19 Figure 19: Phase diagram of the system NH 3 + guaiacol at 373.15 K and 573.15 K using GC-PPC-SAFT (Scheme 1) and Molecular Simulation.

Figure 20 :

 20 Figure 20: Phase diagram of the system H 2 + guaiacol at 573.15 K using GC-PPC-SAFT (Scheme 1) and Molecular Simulation.

Figure 21 :

 21 Figure 21: Phase diagram of the system CO + guaiacol at 573.15 K using GC-PPC-SAFT (Scheme 1) and Molecular Simulation.

Figure 22 :

 22 Figure 22: Henry constant of gases in guaiacol.

,

  where x i is the molar liquid fraction of the solvent (1), except for (c), b Disregarding the value at 1.9 × 10 -3 MPa, c relative deviation expressed in bubble data at isobaric conditions (101.3 kPa).

Figure 1 .

 1 Figure 1.

Figure

  Figure(s)

  Figure 5.

  Figure 6.

  Figure 7 .

  Figure 9.

Figure 11 .

 11 Figure 10.

  Figure 13.

  Figure 14.

  Figure 15.

  -SAFT-Scheme 1 GC-PPC-SAFT-Scheme 2 GCA EoS-Spinosa et al 2005

Figure 16 .

 16 Figure 16.

  Figure 18.

  Figure 19.

  Figure 20.

  Figure 21.

Figure 22 .

 22 Figure 22.

Table 3 :

 3 Molecular model for small gases[START_REF] Möller | Prediction of thermodynamic properties of fluid mixtures by molecular dynamics simulations: methane-ethane[END_REF] -[START_REF] Eckl | An optimised molecular model for ammonia[END_REF].

	Molecule	Center of force		Position		Lennard-Jones Parameters	Charge (e)
			X	Y	Z	ε/k (K)	σ (Å)
	CH 4	CH 4	0	0	0	149.920	3.7327	0
		C	0	0	0	28.129	2.7570	0.6512
	CO 2	O1	1.149	0	0	80.507	3.0330	-0.3256
		O2	-1.149	0	0	80.507	3.0330	-0.3256
		H 2	0	0	0	36.700	2.9600
		q-	0	0	0	-	-	-0.9328
	H 2
		q+(1)	0.3705	0	0	-	-	0.4664
		q+(2)	-0.3705	0	0	-	-	0.4664
		C	0	0	0	47.676	3.2490	0.05376
	CO
		O	1.1283	0	0	30.820	3.4052	-0.05376
		S	0	0	0	250.000	3.7300	0.4000
		H1	0.9639	0.9308	0	-	-	0.2500
	H 2 S
		H2	-0.9639	-0.9308	0	-	-	0.2500
		q-	0	0.1862	0	-	-	-0.9000
		N	0	0	0.0757	182.900	3.3760	-0.9993
		H1	0.9347	0	-0.3164	-	-	0.3331
	NH 3
		H2	-0.4673	0.8095	-0.3164	-	-	0.3331
		H3	-0.4673	-0.8095	-0.3164	-	-	0.3331

Table 4 :

 4 Non-bonded parameters of the AUA4 force field for guaiacol.

	Center of force	 /k (K)	 (Å)	δ (Å)	Charge (e)
	CH 3	120.15	3.607	0.216	+0.223
	CH arom	89.4	3.246	0.407	0
	C arom (-OH)	37.7	3.246	0	+0.265
	C arom (-O)	37.7	3.246	0	+0.223
	O (-CH 3 )	59.69	2.991	0	-0.446
	O (OH)	125.01	3.081	0.010	-0.700
	H (OH)	0	0	0	+0.435

Table 5 :

 5 Bonded parameters of the AUA4 force field for guaiacol.

				Parameter
	Groups	Bond	Angle	Torsion
		lenght r 0	θ 0	k bending	a i (K)
		(Å)	(deg)	(K)
	CH 3 -O	1.360	-	-	-
	CH x arom =CH x arom (x=0 or 1)	1.400	-	-	-
	C arom -O	1.430	-	-	-
	O-H	0.945	-	-	-
	CH 3 -O-CH arom	-	112	69000	-
	CH x arom =CH x arom =CH x arom (x=0 or 1)	-	120	rigid	-
	CH x arom =C arom -O (x=0 ou 1)	-	120	rigid	-
	C arom -O-H	-	108.5 61000	-
	CH x arom -CH x arom -CH x arom -CH x-arom (x=0 or 1)	-	-	-	rigid
	CH x arom -CH x arom -CH x arom -O (x=0 ou 1)	-	-	-	rigid
					a 0 = 845.65; a 1 =0;
	CH x arom -CH x arom -O-H	-	-	-
					a 2 = -845.65; a 3 = 0
	O-CH x arom -CH x arom -CH x arom	-	-	-	rigid
	O-C arom -C arom -O	-	-	-	rigid
					a 0 =1631.27; a 1 = 2651.77;
					a 2 = 868.86; a 3 = -2184.14;
	CH arom -CH arom -O-CH 3	-	-	-	a 4 = -3028.37; a 5 = -3065.12;
					a 6 = 1926.79; a 7 = 1723.03;
					a 8 = -497.96

Table 6 :

 6 Experimental data of systems methane + guaiacol. Uncertainties (U) on the measurements are as follows: U(T) = 0.05 K, U(P) = 0.05MPa, u(x 1 ) = 0.002

	Systems	T (K)	P (MPa)	x 1
		443.5	4.967	0.03421
		443.5	9.988	0.06633
	Methane (1) + guaiacol (2)
		443.5	14.874	0.10439
		443.5	24.929	0.20514
		443.5	34.916	0.26742

Table 7 :

 7 Associative Parameters for guaiacol.

				Parameters
		Group
			ε AB /k (K)	κ AB	R
	Scheme 1	OH	812.17	0.2429	0.4040
		O	2298.84	0.0051	0.4516
	Scheme 2	OH	534.50	0.1466	0.4102
		O	2303.49	0.0050	0.4511

Table 8 :

 8 Results of vapor pressures, density and heat of vaporization of guaiacol obtained from Molecular Simulation.

	Temperature (K)	Vapor Pressure (MPa)	Density (kg/m 3 )	Heat of vaporization (kJ/mol)
	450	0.051	967.3	48.6
	500	0.173	922.9	45.4
	550	0.460	870.0	41.2
	600	0.986	802.3	36.3

Table 9 :

 9 Relative deviations from experimental data for correlated values of vapor pressure, liquid volume and vaporization enthalpy of guaiacol using GC-PPC-SAFT and Molecular Simulation (in %).

			GC-PPC-SAFT	Molecular Simulation
		Temperature range (K)	Scheme 1	Scheme 2
	Vapor Pressure	361-478	2.15	2.17	4.8
	Liquid Volume	313-478	0.36	0.38	1.3
	Vaporization Enthalpy	305-627	4.5	4.5	4.8

Table 10 :

 10 Relative deviations of the solubility for systems with guaiacol.

			Relative Deviation (%) a
	System
		Temperature (K)	GC-PPC-SAFT	Molecular	Ref
			Scheme 1	Scheme 2	Simulation
	CO 2 + guaiacol	323.15	11.9	10.4	-	[7]
		353.15	12.5	9.4	-
		393.15	10.3	4.7	56.4
	Methane + guaiacol	443.15	20.2	19.5	-	This work
	Ethanol + guaiacol	290.15	6.2 (or 2.5 b )	9.1 (or 7.1 b )	-	[64]
	Octanol + guaiacol	433.15	1.7	3.2	-	[8]
		488.15	1.3	3.7	-
		463.15	0.6	2.4	-
	Acetone + guaiacol	290.15	50 c	71 c	-	[64]
	Butyl acetate + guaiacol (401.2 -473.5) d	26	36	-	[65]
	Water + guaiacol	313.6	7.8	>100	-	[61]

(*): the association scheme shows the number of electron-donor (-) and electron acceptor (+) sites attributed to each group. '0' means no sites. (*) See original paper for details.

(**) the association scheme shows the number of electron-donor (-) and electron acceptor (+) sites attributed to each group. '0' means no sites.


Chih-Hung Chen

Tristan Cambonie

Veronique Lazarus

Matteo Nicoli

Antonio Pons

Alain Karma

Crack Front Segmentation and Facet Coarsening in Mixed-Mode Fracture

Keywords: numbers: 62.20.Mk, 46.50.+a, 46.15.x

HAL is

Crack propagation is a main mode of materials failure. Understanding and controlling this complex phenomenon continues to pose both fundamental and practical challenges. While quasi-static planar crack growth with a tensile stress normal to the fracture plane (mode I) is well-understood, geometrically much more intricate crack patterns can form in varied conditions [1]. A few examples include thermal or drying stresses that can cause cracks to oscillate and branch [2,3], or re-organize into complex three-dimensional patterns [4][5][6], nonlinear elastic effects that can induce crack front instabilities even in mode I [7], or the superposition of mode I and a shear stress parallel to the crack front (mode III). This mixedmode I+III fracture is observed in a wide range of engineering and geological materials to produce arrays of daughter cracks, which are shaped as tilted facets and form by a geometrically complex crack front segmentation process [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23].

Recent theoretical progress has been made to characterize the crack-front instability leading to segmentation [24,25] and to describe the propagation of daughter-crack arrays [26]. However, theory and experiments have not produced a consistent picture. Griffith's energetic criterion [27] predicts that planar crack growth is possible when the elastic energy release rate

G = 1 2µ (1 -ν)K 2 I + K 2 III , (1)
exceeds a critical material-dependent threshold G c , where K I and K III are the mode I and mode III stress intensity factors (SIF), respectively, which characterize stress divergences near the crack front, µ is the shear modulus and ν is Poisson's ratio. Phase-field simulations of brittle mixed-mode I+III fracture have revealed that planar growth is linearly unstable against helical defor-mations of the crack front, which couple in-plane and out-of-plane perturbations and develop nonlinearly into facets [24]. A subsequent linear stability analysis in the framework of linear elastic fracture mechanics (LEFM) [25] has predicted that this helical instability should occur when K III /K I exceeds a threshold

K III K I c = (1 -ν)(2 -3ν) 3(2 -ν) -4 √ 2 (1 -2ν) , (2)
which only depends on Poisson's ratio. However, crack front segmentation is experimentally observed for K III /K I values much smaller than this threshold [8,23], or even vanishingly small [22]. This apparent disagreement between linear stability analysis and experiment raises the question of whether LEFM and phase-field modeling are adequate theories to describe crack propagation in mixed-mode I+III brittle fracture. Also poorly understood is "facet coarsening", the progressive increase of facet width and spacing with propagation length from the parent crack. Phase-field modeling [24] and experiments [28] suggest that coarsening may be a self-similar process, but its precise mechanism and dependence on mode mixity are not well understood.

In this letter, we investigate both facet propagation and coarsening by mixed-mode I+III fracture experiments that allow us to visualize in-situ complex crack morphologies during quasi-static propagation, thereby providing much more detailed geometrical information on crack front evolution than conventional post-mortem fractography. Moreover, we carry out phase-field simulations of those experiments that allow us to relate experimental observations to LEFM theory. The results help resolve the puzzling discrepancy between linear stability analysis and experiments with regards to facet formation and shed new light on the coarsening process.

Experiments are carried out using plexiglas beams and a traditional three or four point bending setup [29]. To introduce some amount of mode III, the initial planar notch in the sample is tilted at an angle from the mode I central plane of symmetry [19,30]. A special procedure is used to initiate a sharp crack with a straight front [29]. The corresponding values of the SIF for each angle and hence K III /K I have been obtained by finite element calculations, which show that K III /K I varies between approximately 0.1 and 0.5 when the notch angle varies be-tween 15 • and 45 • , where zero angle corresponds to pure mode I loading. This range was selected because it contains the linear instability threshold (K III /K I ) c ≈ 0.39 predicted by Eq. ( 2) for Poisson's ratio of plexiglass ν ≈ 0.38. Finite element calculations also show [19] that K III /K I is reasonably constant away from sample edges, thereby allowing us to investigate crack propagation at constant K III /K I along a wide section of the parent crack inside the sample. Several beams were broken by fatigue in the bending set-up [29]. The advantage of this cyclic type of loading is that the crack advance (i) is quasi-static, while leaving the crack path unchanged in comparison to the one obtained under monotonical increasing loading [START_REF] Linning | Mixed-Mode Fatigue and Fracture (ESIS 14)[END_REF] and (ii) controlled by the number of cycles so that complex crack morphologies can be observed in-situ at different stages of crack growth. Observations were made using a Leica binocular or a Keyence numerical microscope by transparency.

Examples of experimental images are shown in Fig. 1(a)-(g) for K III /K I values of 0.3 and 0.5 corresponding to initial notch angles of 30 • and 45 • , respectively. Those images reveal several important features. Firstly, facets have a finger-shape with curved tips and flat sides that is consistent with the shape predicted by phasefield simulations (Fig. 1(i) and Movie 1 of [29]). Secondly, facets form for values of K III /K I both below and above the linear stability threshold (K III /K I ) c ≈ 0.39. Within optical resolution, only energetically favored type A facets are observed to emerge from the parent crack with a well-defined tilt angle θ from the original fracture plane. Thirdly, facets coarsen by elimination of other facets leading to an increase of both facet width and facet spacing along the array with increasing propagation length. Coarsening is clearly visible from top views in Fig. 1(b) and in the sequence Fig. 1(c)-(e), which moreover shows that surviving facets maintain the same angle while overgrowing others. Additional views are given in [29].

Simulations were carried out using a phase-field model of brittle fracture that, like gradient damage models [6,[START_REF] Bourdin | [END_REF], regularizes stress-field divergences on a process zone scale ∼ ξ around the crack front. All energy dissipation takes place on a characteristic timescale τ [33]. As shown by an asymptotic analysis of the phase-field model in the limit where ξ is much smaller than all other dimensions [34], fracture in this model is governed by standard crack propagation laws assumed in the LEFM theoretical framework, namely Griffith's criterion and vanishing mode II SIF [35]. Since we are primarily interested in modeling crack evolution in a region away from the experimental sample boundaries where K III /K I is approximately uniform [19,28], we carried out simulations in a rectangular slab geometry of length D x , width D y and height D z , defined in Fig. 2(b), with the origin defined at the center of the slab. We impose fixed displacements at y = ±D y /2, u y (x, ±D y /2, z) = ±∆ y (mode I) and u z (x, ±D y /2, z) = ±∆ z (mode III), periodic boundary conditions in z that allow us to model a periodic array of daughter cracks infinite in z [24]. We use a "treadmill" that adds a strained (y, z) layer at x = D x /2 and removes a layer at x = -D x /2 when the crack has advanced by one lattice spacing. This allows us to simulate crack propagation lengths much longer than D x (a D x ), thereby modeling propagation in a slab infinitely long in x [29]. We also choose D x ≥ 2.5D y to eliminate the influence of the two end-boundaries of the slab (x = ±D x /2) on the central region of the slab (|x| D x ) where the average crack front position is maintained by the treadmill. Standard expressions of linear elasticity are used to relate ∆ y and ∆ z to the SIF [29] and therefore to K III /K I and G/G c where G c ≈ 2γ (twice the surface energy) is known in the phase-field model [33,34]. All simulations are performed with ν = 0.38 of plexiglass. We simulated both quasi-static propagation, where the elastic field is relaxed at each time step of crack advance, and dynamic propagation by solving the full elastodynamic equations. Both sets of simulations yielded similar results for the range G/G c ≤ 1.5 where the ratio of the crack propagation speed to the shear wave speed v/c ≤ 0.3 is small enough to neglect inertial effects [29].

We first carried out simulations to check quantitatively the theoretical prediction of Eq. ( 2). For this purpose, we slightly perturbed the planar parent crack with a small amplitude helical perturbation of the form δx front + iδy front = A 0 e -ikz , where δx front and δy front indicate the x and y components of deviations of the front from the reference planar crack, respectively, and k = 2π/D z fits one wavelength D z = Λ of the perturbation in the periodic domain in z. The stability of planar crack propagation is then determined by tracking the amplitude of the perturbation that grows or decays exponentially in time [29] if propagation is unstable, as illustrated in Fig. 2(a)-(d), or stable, respectively. Simulations were carried out by increasing K III /K I in small steps to determine the threshold (K III /K I ) c , and repeating this procedure for increasing values of D y /Λ to quantify finite size effects. Fig. 2(e) shows that (K III /K I ) c increases monotonously with D y /Λ and approaches a value reasonably close to the prediction (K III /K I ) c ≈ 0.39 of Eq. ( 2) in the large system size (D y /Λ 1) limit. Consistent with the result of Fig. 2(e), an examination of strain fields shows that finite size effects becomes negligible when D y /Λ ≥ 2 [29]. We conclude that LEFM theory (Eq. ( 2)) and phase-field modeling predict similar linear instability thresholds in the large system size limit, and that facets are experimentally observed well below this threshold.

Next, in order to explore the nonlinear character of the bifurcation from planar to segmented crack front, we measured experimentally the facet tilt angle θ extracted from three-dimensional maps of post-mortem fracture surfaces obtained using a profilometer as detailed in [28]. The angle θ is plotted versus K III /K I in Fig. 3(a). Furthermore, we investigated computationally the propagation of periodic arrays of A facets, where coarsening is suppressed by choosing D z = Λ due to the periodic boundary conditions along z. In this geometry, we tracked the steady-state branch of propagating solutions by decreasing K III /K I starting from values above the linear instability threshold to values below this threshold, as low as 0.07 to span the entire experimental range of mode mixity. For each K III /K I value, we allowed the facet to relax to a new stationary shape and tilt angle, as illustrated in Fig. 3(b) for a simulation where K III /K I was decreased from 0.5 to 0.07. The computed tilt angles are compared to experimental results in Fig. 3(a) with the corresponding facet shapes shown in Fig. 3(c). Both the facet shapes, which gently curve at their extremities in the yz plane due to elastic interactions between neighboring facets, and the tilt angles are in good quantitative agreement with experimental observations within measurement errors. Fig. 3(a) also shows that computed tilt angles are weakly dependent on system size (D y /Λ) and fall below the prediction of a simple theory, which assumes that facets are shear-free [16,24]. Those results demonstrate that propagating segmented front solutions exist over the entire range of K III /K I investigated experimentally, including values less than (K III /K I ) c . We conclude that the bifurcation from planar to segmented front is strongly subcritical, with bistability of planar and segmented crack growth for K III /K I < (K III /K I ) c as illustrated schematically in Fig. 3(d).

To characterize coarsening in phase-field simulations, we investigated the stability of periodic array of facets by repeating the above series of simulations with two facets (D z = 2Λ). This geometry is motivated by the striking similarity between the coarsening behavior of facets in the present experiments (Fig. 1(a)-(g)) and coarsening of curved fronts in other interfacial pattern forming systems, in particular viscous fingering [36] and dendritic crystal growth [37,38]. In those systems, it is wellestablished that coarsening of finger arrays is associated with a spatial period-doubling linear instability of the array leading to elimination of one of every two fingers in the array by exponential amplification of small perturbations. Results of simulations illustrated in Fig. 4(a) show that arrays of facets exhibit a similar spatial period doubling instability driven by elastic interactions between facets. This instability yields an increase (decrease) of the SIF and hence the energy release rate at the tips of leading (lagging) facets. The amplification rate of instability is obtained by computing the difference of x-tip position ∆x tip (t) between leading and lagging facets, which grows exponentially in time starting from an infinitesimal perturbation, ∆x tip (t) ≈ ∆x tip (0)e ωv0t/Λ , where v 0 and Λ are the initial facet growth velocity and spacing, respectively. The slopes of semi-log plots of ∆x tip (t)/Λ versus v 0 t/Λ in Fig. 4(b) yield values of ω that increase markedly with K III /K I , showing that a larger mode III component leads to a faster elimination rate of facets. Coarsening, clearly visible in Fig. 1(b) and other experimental views [29], was quantified experimentally by analyzing post-mortem fracture surfaces [28]. The results show that the relation between the mean facet spacing Λ and the crack propagation length a is approximately linear, with a mean slope β ≡ dΛ/da increasing with K III /K I (inset of Fig. 4(b)). To relate the coarsening rates in phase-field simulations and experiments, we derive a simple evolution equation for the average array spacing Λ based on dynamical mean-field picture as previously done for dendritic arrays [37]. The coarsening rate β ≡ dΛ/da ≈ ∆Λ/∆a where ∆Λ is the change of array spacing due to elimination of one of every two facets along the array or ∆Λ ≈ Λ, while ∆a is the distance that the facets propagated during the elimination process. Since elimination occurs via exponential amplification of small perturbations, facets will propagate an average distance ∆a ∼ Λ/ω during this process, yielding the prediction β ∼ ω, or β = Cω where C is a constant prefactor of order unity. The comparison in the inset of Fig. 4(b) shows that this simple theory is able to predict reasonably well the increase of the coarsening rate with K III /K I up to the value of the constant prefactor C = 0.198 determined from a global best fit to the experimental data for all K III /K I values.

The present results reconcile the prediction of linear stability analysis (Eq. ( 2)) with experimental observations by showing that the bifurcation from planar to segmented crack growth is strongly subcritical; facet arrays exist as fundamental crack propagating solutions of LEFM for a range of K III /K I values extending below the instability threshold. They further show that coarsening is driven by a spatial period doubling instability of facet arrays with a growth rate that depends on mode mixity. The reasonably good quantitative agreement between simulated and observed morphologies suggests that LEFM is an adequate theory to describe complex geometrical features of both brittle and fatigue cracks in mixed mode I+III fracture. While the present results show that the subcritical propagation of segmented cracks is theoretically possible, they do not identify the mechanism and scale of subcritical facet formation. As suggested by a recent LEFM analysis, materials imperfections may contribute to this process [39]. However, this scenario, and even more fundamentally the ability of LEFM to model subcritical facet formation, remain to be explored both computationally and experimentally.

The To initiate a true crack with a sharp tip and a smooth straight front (Fig. 1, right), we pushed a wedge (razor blade) into the slit, quasi-statically, under controlled slowly increasing force, using a tensile machine. The total length d of the crack (slit+sharp crack) is d ∼ W/3 in both cases. The residual stresses introduced by the slit manufacturing, were relaxed by heating the samples at 90-95 • C during 10 hours. Their annihilation is checked using polarizers. To introduce some amount of mode III, the initial planar notch in the sample is tilted at an angle Γ 0 from the mode I central plane of symmetry, Γ 0 varying between 15 • and 45 • , where zero angle corresponds to pure mode I loading. Larger and smaller samples were loaded in 3 and 4 point bending setups, respectively. The loading frequency was f = 5 Hz and K min /K max = 0.1. The amplitude of the stress intensity factor ∆K ∼ 0.5 MPa•m 1/2 has been chosen well below the brittle fracture threshold K c ∼ 1 MPa•m 1/2 , to avoid brittle fracture, while being large enough to ensure propagation [1].

Sets of experiments were performed on series of similar beams by stopping the propagation at different stages of the propagation. Figure 2 gives some in-situ views of the corresponding crack morphologies, acquired with a numerical Keyence microscope. The coalescence rate β cannot be quantified on the in-situ samples due to optical distortions induced by the plexiglas. For this purpose, other samples were broken completely (Fig. 3 (a)-(c)) and three dimensional profilometer maps of their fracture surfaces were done. Quantification of θ (Fig. 3 text) has been done performing some statistical post-treatment over several facets of these profiles. The procedure is explained in detail in [2].

The more clear zones in Fig. 3(a)-(c) correspond to fragments left behind by partial breaking of the zones located between adjacent facets. Indeed, beside propagating in the x-direction, the facets propagate also in the lateral z-direction. As they interact, they curve to form the well-known [3,4] "en-passant S-shape" pattern (see Fig. 3(d)). Since the symmetry of this interaction is sensible to any perturbation, generally only one out of two crack tips connects to the adjacent facet, leaving on one of the fracture surface, a partially broken fragment. These fragments are visible in white on Fig. 3(a)-(c) and are lacking on the complementary fracture surface.

SIMULATIONS

We conducted large-scale simulations of mixed-mode I+III fracture using the phase-field model originally proposed by Karma, Kessler and Levine (KKL) [5], which has been used to solve fracture propagation problems over the last decade [6][7][8][9][10]. The model introduces a scalar order parameter φ to distinguish between intact and broken states of the material. The total energy of the system is given by the functional [5,11]

E = ρ 2 ∂ u ∂t 2 + κ 2 (∇φ) 2 + g (φ) (e strain -e c ) dV, (1)
where ρ is the mass density; g (φ) = 4φ 3 -3φ 4 is conventionally chosen with the properties that g (0) = 0, g (1) = 1 and g (1) = g (0) = 0 [5,11]; u = (u x , u y , u z ) represents the displacement field; e strain = λ (ε ii ) 2 /2 + µ (ε ij ) 2 is the strain energy density and ε ij = (∂ i u j + ∂ j u i ) /2 is the strain tensor of a linear elastic material with i = 1, 2 and 3 corresponding to x, y and z, respectively, λ and µ are the Lamé coefficients. The function g decreases monotonously from the intact state (φ = 1) to the fully broken state (φ = 0) and accounts for elastic softening at large strain. When the strain energy density e strain exceeds the threshold e c , the broken state becomes energetically favored. Energy dissipation occurs in the process zone of size ξ = κ/ (2e c ) around the crack tip, where φ varies smoothly between 0 and 1. In the phase-field model, the Griffith's criterion G c is given by G c = 2 √ 2κe c 1 0 φ 1g (φ) [5,11]; the iso-surfaces of φ = 0.5 are conventionally defined as the fracture surfaces. The equations of motion for the phase-field φ and displacement field u are derived variationally from the energy functional [5,11]:

∂φ ∂t = -χ δE δφ , (2)
ρ ∂ 2 u i ∂t 2 = - δE δu i . (3)
We studied crack-front instabilities of mixed mode I+III fracture propagation along the x-axis in a rectangular slab of length D x , width D y and thickness D z with center at the origin. The mixed I+III loading was imposed by fixed displacements at y = ±D y /2, u y (x, ±D y /2, z) = ±∆ y (mode I) and u z (x, ±D y /2, z) = ±∆ z (mode III), and periodic boundary conditions in the z-direction. To allow full relaxation of the crack front to a stationary state, simulations were carried out on a "treadmill" that adds a strained (y, z) layer at x = D x /2 and removes a layer at x = -D x /2 when the crack has advanced by one lattice spacing. This allows us to simulate crack propagation lengths much longer than D x (a D x ), thereby modeling propagation in a slab infinitely long in x. The crack dynamics is controlled by two key parameters K III /K I and G/G c , where K I = 2µ∆ y 2 (1 + λ/µ) /D y and K III = 2µ∆ z / D y are the stress intensity factors of a semi-infinite planar crack and G = (1ν) K 2 I + K 2 III / (2µ) is the corresponding energy release rate, where ν is Poisson's ratio. Simulations were carried out with the range G/G c ≤ 1.5 where the ratio of the crack propagation speed to the shear wave speed v/c ≤ 0.3, which is small enough to neglect inertial effects. The slab length D x was chosen greater than 2.5 of the slab width D y to eliminate the influence of the two end-boundaries of the slab (x = ±D x /2) on the central region of the slab (|x| D x ) where the average crack front position is maintained by the treadmill. As a result, the planar crack propagation speed becomes independent on the system length, as shown in Fig. 4.

We rewrite the phase-field model in a dimensionless form by measuring length in units of the fracture process zone scale ξ and time in units of the characteristic dissipation time scale τ = 1/ (µχ). The remaining simulation parameters are chosen as follows, e c /µ = 1/2, ν = 0.38 and cτ /ξ = 2.76, where the shear wave speed c = µ/ρ. Large scale simulations of the order of 10 7 -10 8 grid points were performed using graphics processing units (GPUs) with the CUDA parallel programming language.

The energy functional of the phase-field model was discretized in space on a cubic grid of uniform mesh size ∆ = 0.3ξ (Fig. 5). To write the discretized energy in compact form, we define the superscripts , m, and n to refer to the gridpoint at x = ∆, y = m∆, and z = n∆, the subscripts {i, j} = {1, 2, 3} to refer to the corresponding {x, y, z} components of the displacement fields, and use the standard Kronecker delta defined by δ ij = 1 if i = j and δ ij = 0 if i = j. We discretize ∂u i /∂x j , (∂u i /∂x j )

2 and (∂φ/∂x j ) 2 on the gridpoint ( , m, n) by the following approximations:

∂u i ∂x j ( , m, n) ≈ 1 2∆ u +δ1j ,m+δ2j ,n+δ3j i -u -δ1j ,m-δ2j ,n-δ3j i , ∂u i ∂x j 2 ( , m, n) ≈ 1 2∆ 2 u +δ1j ,m+δ2j ,n+δ3j i -u ,m,n i 2 + u ,m,n i -u -δ1j ,m-δ2j ,n-δ3j i 2 , and ∂φ ∂x j 2 ( , m, n) ≈ 1 2∆ 2 φ +δ1j ,m+δ2j ,n+δ3j -φ ,m,n 2 + φ ,m,n -φ -δ1j ,m-δ2j ,n-δ3j 2 .
Accordingly, the strain energy density e strain on gridpoint ( , m, n) becomes e ,m,n strain =

1 4∆ 2 (λ + 2µ) H ,m,n 1 + H ,m,n 2 + H ,m,n 3 + λ I ,m,n 1,2 + I ,m,n 2,3 + I ,m,n 3,1 +µ J ,m,n 1,2 + J ,m,n 2,3 + J ,m,n 3,1
.

where we have defined the functions

H ,m,n i ≡ u +δ1i,m+δ2i,n+δ3i i -u ,m,n i 2 + u ,m,n i -u -δ1i,m-δ2i,n-δ3i i 2 , I ,m,n i,j ≡ u +δ1i,m+δ2i,n+δ3i i -u -δ1i,m-δ2i,n-δ3i i u +δ1j ,m+δ2j ,n+δ3j j -u -δ1j ,m-δ2j ,n-δ3j j , and
J ,m,n i,j ≡ u +δ1j ,m+δ2j ,n+δ3j i -u ,m,n i 2 + u ,m,n i -u -δ1j ,m-δ2j ,n-δ3j i 2 + u +δ1i,m+δ2i,n+δ3i j -u ,m,n j 2 + u ,m,n j -u -δ1i,m-δ2i,n-δ3i j 2 + u +δ1j ,m+δ2j ,n+δ3j i -u -δ1j ,m-δ2j ,n-δ3j i u +δ1i,m+δ2i,n+δ3i j -u -δ1i,m-δ2i,n-δ3i j .
The discretized form of the energy functional Eq. ( 1) on the cubic grid becomes

E ∆ 3 ≈ ,m,n 3 i=1 ρ 2 ∂ t u ,m,n i 2 + 3 i=1 κ 4∆ 2 F ,m,n i + g φ ,m,n e ,m,n strain -e c , (4)
where, in addition to the functions defined above,

F ,m,n i ≡ φ +δ1i,m+δ2i,n+δ3i -φ ,m,n 2 + φ ,m,n -φ -δ1i,m-δ2i,n-δ3i 2 .
From Eq. ( 2), Eq. (3) and Eq. ( 4), we obtain the equations of motion

∂φ ,m,n ∂t = -χ ∂ E/∆ 3 ∂φ ,m,n , (5)
which are integrated in time using a Verlet scheme with a timestep size ∆t/τ = 0.001.

To investigate the onset of instability of mixed mode I+III crack propagation, we carried out simulations starting from a parent planar crack in yz plane for x ≤ 0. As in [10], the planar crack front was initially perturbed by a helical perturbation, δx front + iδy front = A 0 e -ikz , where δx front and δy front indicate the x and y components of deviations of the front from the reference planar crack, respectively. The amplitude A (t) = A 0 e σt/τ was then found to grow (σ > 0) or decay (σ < 0) exponentially in time (Fig. 6).

The theoretical instability threshold predicted by Eq. (2) of the main text is only valid in the large system size limit [12]. To characterize finite size effects, simulations with D y /Λ between 0.5 and 4 were carried out. The results show that finite size effects reduce significantly (K III /K I ) c when D y is comparable to Λ (Fig. 2(e) of the main text). For such system sizes, the strain fields along the y = ±D y /2 boundaries are spatially varying along z and hence strongly influenced by the facets. However for D y /Λ ≥ 2, the strain fields become nearly independent of z as expected in the large system size limit D y Λ (Fig. 7). In addition, the results of simulations with Λ/ξ between 30 and 60 show that the crack-front instability is independent of Λ in the limit Λ ξ.

1 FIG. 1 :

 11 FIG. 1: (Color online). In-situ microscope images (a)-(g) of fatigue cracks in plexiglass at different stages of crack advance in mixed mode I+III loading depicted schematically in (h) and corresponding example of crack-front segmentation in phase-field simulation (i). KIII /KI ≈ 0.3 in (a)-(e) and ≈ 0.5 in (f)-(g); (a), (b) and (f) are experimental views from a direction approximately perpendicular to the plane of the parent crack with facets propagating downwards, while views (c), (d), (e) and (g) are views with the crack propagation direction out of the page. Views (c), (d) and (e) correspond to different stages of crack advance increasing from (c) to (e). Broken (pristine) regions of the samples appear in black (light blue) or darker (lighter) grey depending on the viewing direction. The bar scale is 1 mm in all images. The red dashed lines in (a) highlight the curved fronts of two facets as guide to the eye; curved tips are clearly visible in (f). (i) Snapshots of phase-field fracture surfaces (φ = 1/2 surfaces) at different stages of crack advance increasing from top to bottom, showing that energetically favored A facets [18] propagate ahead of B facets eventually outgrowing them completely. Simulation parameters are G/Gc = 1.5, KIII /KI = 0.5, and box dimensions Dx = 307ξ, Dy = 100ξ and Dz = 200ξ.

1 FIG. 2 :

 12 FIG. 2: (Color online). Snapshots phase-field simulations illustrating the destabilization of planar crack growth for KIII /KI = 0.4. The crack propagation length a increases from (a) to (d) and both the crack front (blue lines) and its inplane and out-of-plane projections (red lines) are shown. (e) Plot of linear instability threshold (KIII /KI ) c versus Dy/Λ. Planar growth is unstable (stable) above (below) the filled circles, where error bars reflect the uncertainty in stability threshold resulting from the fact that KIII /KI was increased in finite steps in the simulations. KIII /KI values corresponding to the top (bottom) of each error bar were simulated and found to yield unstable (stable) propagation. In all simulations, G = 1.5Gc, Dx = 230ξ and Dz = Λ = 60ξ.

1 FIG. 3 :

 13 FIG. 3: (Color online). (a) Comparison of facet tilt angles obtained from experiments and simulations, where red and blue arrows indicate the instability thresholds of planar crack propagation for Dy/Λ = 1 and Dy/Λ = 2, respectively (see Fig. 2(e)), and theoretically predicted assuming shearfree facets (dashed line) [16, 24]. (b) Snapshots of a phasefield simulation for Dy/Λ = 1 demonstrating the subcritical nature of the bifurcation from planar to segmented crack propagation. A propagating segmented front solution for KIII /KI = 0.5 was used as initial condition (θ = 31 • ). The facet continuously rotated towards a lower angle in response to the decrease in KIII /KI and then reached its steady state (θ = 11.2 • ) after propagating a distance a = 5.5Λ (see Movie 2 of [29]). (c) Out-of-plane and in-plane (inset) crack-front projections. In all simulations, Dx = 154ξ, Dy = Dz = 60ξ, Λ = 60ξ and G = 1.5Gc. (d) Schematic diagram of subcritical bifurcation recapitulating the experimental and simulations results with solid (dashed) lines representing stable (unstable) solutions.

1 FIG. 4 :

 14 FIG. 4: (Color online). (a) Illustration of spatial period doubling instability in a phase-field simulation for KIII /KI = 0.5; out-of-plane and in-plane projections of crack fronts at different times are plotted in the top panel and the bottom panel, respectively (see Movie 3 of [29]). (b) Semi-log plot of difference of tip positions along the propagation x-axis between leading and lagging facets versus scaled time for different KIII /KI . Inset: coarsening rate β versus KIII /KI obtained from experiments and phase-field simulations. In all simulations, Dx = 307ξ, Dy = 60ξ, Dz = 120ξ, Λ = 60ξ and G = 1.5Gc.

1 1FIG. 1 :

 11 FIG. 1: Left hand side: 3 Point Bending setup. Right hand side: Lateral view of the initial slit made by micro-milling and of the additional sharp crack introduced by pushing under a controlled manner, a razor blade in the initial slit.

  Front and bottom views are given in the first and second row respectively. Each column corresponds to the same sample. The beam size is L = 100 mm, W = 10 mm, b = 10 mm; columns (a)-(c) correspond to Γ 0 = 30 • (K III /K I ≈ 0.3) and column (d) to Γ 0 = 20 • (K III /K I ≈ 0.18). The development of rotated facets that coalesce during propagation are clearly visible. One can observe i) comparing columns (a) to (c), that the rotation angle of those facets is approximately constant during propagation, for a given value of K III /K I and ii) comparing columns (a)-(c) with (d), that the rotation angle of these facets decreases with decreasing K III /K I .

FIG. 3 :

 3 FIG. 2: In-situ numerical microscope images of partially broken samples. Each column corresponds to one sample, the first row being a front view and the second row being a bottom view as indicated by cartesians axes corresponding to Fig. 1(h) of the main text. KIII /KI ≈ 0.3 in (a)-(c) and ≈ 0.18 in (d).In the front views, the initial slit and the facets appear in black. In the bottom views, the initial larger slit is dark and the facets appear in white. The bar scale is 1 mm in all images.

5 FIG. 4 :

 54 FIG. 4: Plot of pure mode I crack propagation speed v (scaled by the shear wave speed c) versus G/Gc. Phase-field simulations show that the crack speed becomes independent of the simulation box length when Dx/Dy ≥ 2.5 (inset). Simulation parameters are Dy = 60ξ and Dz = 30ξ.

1 FIG. 5 :K 35 K 25 FIG. 6 :

 1535256 FIG.5: Spatial discretization of phase-field model on a cubic grid of uniform mesh size ∆. Gridpoint labeled ( , m, n) represents a unit cubic of material of mass ρ∆ 3 at location ( ∆, m∆, n∆).

1 FIG. 7 :

 17 FIG. 7: Finite size effects: (a) (Color online). Images of strain field εyz of a slice along x behind the facet tips (the corresponding slice position is illustrated in the top panel) from simulations with Dy/Λ = 1, 2 and 4. (b) Plots of strain field εyz of the sliced regions at y = Dy/2 versus z (scaled by Λ). The strain fields on boundaries are strongly influenced by facets when Dy/Λ = 1 (red line) nearly constant along z for both Dy/Λ = 2 (blue dashed line) and Dy/Λ = 4 (green double dashed line), where boundary effects become negligible. In all simulations, KIII /KI = 0.35, Dx/Dy = 3.75 and Dz = Λ = 30ξ.

ACKNOWLEDGMENTS

The research at Northeastern University was supported by Grant No. DE-FG02-07ER46400 from the U.S. Department of Energy, Office of Basic Energy Sciences and a seed grant from the Massachusetts Green High Performance Computing Center. The research at University Paris Sud benefited of financial support from ANR GeoSMEC (2012-BS06-0016-03). We thank posthumously the late F. Buchholz, who has broken the largest samples, we have further exploited here. We also thank L. Auffray, D. Bonamy, V. Doquet, J.-C. Eytard, R. Pidoux, A. Tanguy for their help in the experiments.