Olivier Bonnefoy
email: olivier.bonnefoy@emse.fr

Frédéric Gruy

Jean-Michel Herri

[ O Bonnefoy

J.-M Herri

Van der Waals interactions in systems involving gas hydrates

Keywords: Methane gas hydrate, Hamaker constant, Agglomeration, Dielectric response function, Van der Waals interaction potential

Van der Waals interactions in systems involving gas

hydrates.

Introduction

Context Gas hydrate clathrates are crystalline compounds composed of host gas molecules trapped in a lattice of water molecules. In industrial contexts, the host gas is commonly methane, ethane or carbon dioxide but clathrates may form with other gases like light alcanes (up to C 4 ), hydrogen sulfide, dioxygen, dinitrogen and some rare gases (argon, neon, krypton, . . . ). Literature is abundant on gas hydrates and the novice can safely refer to Sloan's book [START_REF] Sloan | Clathrate Hydrates of Natural Gases[END_REF] for a detailed presentation of these crystals. Among the industrial contexts where they appear, we shall cite : hydrate plugs obstructing oil-or gas-pipelines, energy resources in form of methane hydrates trapped in permafrost or submarine sediments, natural gas transport in form of a slurry or gas separation by fractionated crystallization. The intelligent exploitation of hydrates in these contexts require a considerable amount of scientific data and the research efforts of many nations (United States, Russia, France, Germany, Great-Britain, Japan and China, to cite the most active ones) greatly contribute to this goal. From this perspective, it seems to us that the knowledge of the Hamaker constant of systems involving gas hydrates is a key point, in particular for those who wish to predict the agglomeration behaviour of these systems. Since we do not expect the Hamaker constant to depend too much from the host gas nature, we restrict our study to methane hydrates.

Interactions between molecules Before going into detail, let us remember some basic concepts about the interactions between molecules. When electrostatic interactions are absent or negligible, molecular interactions are of two natures. Firstly, we find the short range, repulsive interactions due to the electronic clouds that cannot overlap. Different models describe this interaction. Among them, let us quote the hard sphere potential, exponential decay potential and power law potential. Secondly, we have the long range, attractive interactions. Three of them are known and they all have a potential that varies as the inverse sixth power of the intermolecular distance.

• The orientation interaction, or Keesom interaction, tends to correlate the relative orientation of two polar molecules. The Keesom potential [START_REF] Keesom | [END_REF] depends on the molecule dipole moments µ i :

u (r) = -1 (4πǫ 0 ) 2 . µ 2 1 µ 2 2
3kT.r 6 where ǫ 0 = 8.854 × 10 -12 C 2 .J -1 .m -1 is the vacuum permittivity and k = 1.38054 × 10 -23 J.K -1 is the Boltzmann constant.

• The induction interaction, or Debye interaction, is observed between an apolar molecule and a polar molecule. The latter creates an electric field that induces a dipole moment on the apolar but polarisable molecule. The Debye potential [3] is

u (r) = -1 (4πǫ 0 ) 2 .
α 1 µ 2 2 + α 2 µ 2 1 r 6

• The dispersion interaction, or London interaction is observed between two apolar molecules and was first explained by Eisenschitz & London [4]. The time fluctuations of the electronic cloud density, produce a transient dipole moment, which average is zero. At every moment, however, this dipole moment generates an electric field that acts on the second molecule and induces a dipole moment. For two molecules of polarisability α i and first ionisation frequency ν i , the London potential [5] is

u (r) = -1 (4πǫ 0 ) 2 .
3α 1 α 2 .hν 1 ν 2 2(ν 1 + ν 2 ).r 6 where h = 6.626 × 10 -34 J.s is the Planck constant.

Retardation effect We should note the existence of a retardation effect of the dispersion interaction. When two apolar molecules are separated by a relatively long distance, the time required for the electric field to travel (t = 2 * L/c) can become comparable to the fluctuation period of the dipole itself 1 . In this case, the field that comes back on the first molecule finds that the instantaneous dipole direction changed and produces a smaller attraction. Casimir & Polder [6] showed that the retardation effect is negligible for distances below ≈ 10 nm and that, for distances above ≈ 100 nm, this effect creates a dependence in -1/r 7 instead of -1/r 6 .

Van der Waals interaction potential The Van der Waals interaction(s) is the generic name for the three attractive interactions. Its potential is the sum of the corresponding potentials

u (r) = -C/r 6 ( 1
)
where C is the Van der Waals constant, positive and equal to the sum of the three contributions

C = 1 (4πǫ 0 ) 2 . µ 2 1 µ 2 2 3kT + α 1 µ 2 2 + α 2 µ 2 1 + 3α 1 α 2 .hν 1 ν 2 2(ν 1 + ν 2 ) (2)
The table 1 gives these constants for different pairs of molecules. Data are from Israelachvili [START_REF] Israelachvili | Intermolecular and Surface Forces[END_REF].

Hamaker approach : theory

This Hamaker approach to calculate the interactions between macroscopic bodies is also known as the microscopic approach since it focuses at the molecular structure scale.

Hypothesis

In 1937, Hamaker [START_REF] Hamaker | [END_REF] proposed a method to calculate the interaction force between two macroscopic bodies (1 and 2). He made two assumptions

• The retardation effect of the dispersion interaction is neglected, whatever the distances are,

• The interaction potential between two molecules keeps the same form even if other molecules surround them. This is often called additivity assumption.

Hamaker used a potential in 1/r 6 to describe the Van der Waals interaction forces between two molecules and the expression (2) to calculate the Van der Waals constant C. Then he summed up the pair-potentials between molecules in body (1) and molecules in body (2).

Interaction potential

The Hamaker method leads to expressions in form of a product of two functions A and f . The function f depends only on the shape of the bodies and on the separation distance d ≥ 0. The

Tadmor [START_REF] Tadmor | [END_REF] proposes a expression for f in the case of two spherical shells, one spherical shell and a sphere as well as one spherical shell and a half space. These expressions may be of interest for example in the case of water droplets undergoing crystallisation into hydrate.

For a sphere in a cylindrical pore, the work of Bhattacharjee & Sharma [11] reviews the existing calculation methods. Zeman & Wales [START_REF] Zeman | Synthetic Membranes : HF and UF Uses[END_REF] as well as Papadopoulos & Kuo [13] proposed numerical methods to solve this question. However, they require prohibitive computing time. Sharma uses the two Hamaker assumptions to simplify these expressions. He studies the limiting case of a particle near the revolution axis and near the pore wall. This study is of major interest to evaluate the interaction potential between a particle and the wall of a porous medium. In the particular case of an infinite cylindrical pore containing a spherical particle near the wall, we have the following approximation

f (d,geometry) = λ 3 3 [(1 -η(1 -λ)) 2 -λ 2 ] . (1 -η(1 -λ)) 2 .(1 + η(1 -λ)) -λ 2
where λ = Rpore R particle and η = r Rpore-R particle with R pore , R particle and r representing respectively the pore radius, the particle radius and the distance between the particle center and the revolution axis of the pore.

Hamaker constant according to Hamaker

When bodies (1) and (2) are separated by vacuum, the Hamaker constant

A H is A H 12 = π 2 .C 12 .ρ 1 .ρ 2 (4)
where C 12 is the Van der Waals constant for the pair of molecules (1,2) and ρ i the molecular density2 constituting the body (i). When bodies (1) and ( 2) are separated by a medium (3), we use again the additivity assumption and the Hamaker constant A H becomes

A H 132 = A H 12 + A H 33 -A H 13 -A H 23 (5)
The Hamaker constant depends only on the intrinsic properties of the materials. It usually ranges between 10 -21 and 10 -19 J. The attractive or repulsive force is simply the derivative of the interaction potential with respect to the separation distance between the two bodies.

Hamaker approach : calculation

We consider the system hydrate/water/hydrate and calculate its Hamaker constant. The use of formulae ( 4) and ( 5) requires the knowledge of Van der Waals constants for water/vacuum/water but also water/vacuum/hydrate and hydrate/vacuum/hydrate as well as the corresponding molecular densities. For the sake of simplicity, we introduce the notion of hydrate "molecule" which structure is known (structure sI or sII) and composition equally (46 water molecules and 8 gas molecules for sI and, respectively, 136 and 24 for sII).

We begin to study two hydrate molecules (1) and (2) in the vaccum, separated by a distance L. Let us consider a molecule A of water or gas in the molecule [START_REF] Sloan | Clathrate Hydrates of Natural Gases[END_REF]. All other molecules I of (1) and all other molecules J of (2) each exert on A a force deriving from the Van der Waals interaction potential. The total force that acts on molecule A is then

F A = I∈(1)-A F AI + J∈(2)
F AJ So, the hydrate molecule (1) is submitted to an attractive force from (2)

F = A∈(1) F A
We shall note here that the half of the terms vanishes since the sum of the forces exerted by molecules of (1) on molecules on (1) is zero (internal forces). Consequently,

F = A∈(1) J∈(2)
F AJ Given the interaction potential between two molecules A and J, separated by a distance r, u AJ(r) = -CAJ r 6 , we can calculate the force exerted by J on A :

F AJ(r) = -∂u AJ(r) ∂r . AJ AJ = -6.C AJ AJ 7 .

AJ AJ

The calculations have been performed numerically using data in table 3 and make it possible to compute the intensity of the attractive force between two hydrate molecules as a function of their separation distance. From that, we can deduce the interaction potential and, by formal analogy with the equation (1), find the Van der Waals constant C hh . In fact, this "constant" varies with distance and tends to a finite limit Chh = 3.500×10 -74 J.m 6 for long3 distances (above 30 nm). The figure 1 represents the Van der Waals constant C hh as a function of the center-to-center distance of two hydrate molecules 4 . When d=1.5 nm, we can see that it is 20 times bigger than its asymptotic value. We can also calculate the Van der Waals constant between an hydrate molecule and a water molecule. For long distances, this constant tends to Chw = 6.930 × 10 -76 J.m 6 . For short distances, the figure 2 represents its variation. When the distance is above 1 nm, the value of C hw does not differ from its asymptotic value by more than 38%.

We can now calculate the Hamaker constant for hydrates in vacuum. With the asymptotic values, we have

A H hh = π 2 .C hh .ρ 2 h
where ρ h is the hydrate molecular density. Since the methane hydrate lattice parameter is 1.203 nm, the density is 912.6 kg/m 3 and the molecular density ρ h is equal to 5.744 × 10 26 m -3 . Consequently, the Hamaker constant is 1.140 × 10 -19 J or5 approximatively 30.2 kT. When vacuum is replaced by liquid water, the Hamaker constant shall be calculated by

A H hwh = A H hh + A H ww -2A H hw The constant A H hw is equal to A H hw = π 2 .C hw .ρ h .ρ w that
A H hwh = 4.589 × 10 -21 J ≈ 1.217 kT
In the hypothetical case of an hydrate without host gas (methane gas in our case), we would have Chw = 6.394 × 10 -76 and Chh = 2.9412 × 10 -74 J.m 6 and then, the constant would be A H hwh = 3.928 × 10 -20 J. Noting θ the occupancy rate of the hydrate cages, the Hamaker constant can be written

A = [3.928 -3.469θ] × 10 -20 J
4 Lifshitz approach : theory

Notations

The Hamaker method suffers from two drawbacks : the additivity assumption is strong and the retardation effect is not taken into account. That is the reason why Lifshitz developed an other theory, called macroscopic theory. This theory was then extended in 1961 by Dzyaloshinskii, Lifshitz and Pitaevskii [15] (known as "DLP" theory). This theory considers the medium as continuous and characterized by interconnected parameters

• the refractive index n (w) and the absorption coefficient k (w) (both real).

The latter measures the fraction of the wave energy that will be dissipated in heat; it is zero for transparent materials

• the complex index of refraction N (w) = n (w) + i.k (w) • the complex dielectric constant ǫ (ω) such that ǫ (ω) = N 2 (w)
• the real part ǫ ′ (ω) and imaginary part ǫ ′′ (ω) of the complex dielectric constant

ǫ (ω) = ǫ ′ (ω) + i.ǫ ′′ (ω)
The complex dielectric constant ǫ (ω) of a continuous medium is the linear response, in terms of electric displacement, to an electromagnetic field oscillating at the pulsation ω 6 . It contains all information on translational, rotational and vibrational relaxations and on the electronic relaxations of the molecules of the considered medium. The imaginary part ǫ ′′ (ω) is always positive and represents the absorption of the medium, that is to say the energy dissipation of a propagating wave. The graphical representation of ǫ ′′ (ω) is called the absorption spectrum. The real part ǫ ′ (ω) quantifies transmission properties of the medium. Its graphical representation is called transmission spectrum. When ω = 0, the complex dielectric constant is equal to the static dielectric constant ǫ 0 .

In the theory of linear response [START_REF] Landau | Statistical Physics[END_REF], the definition of the complex dielectric constant ǫ can be extended to complex variables. For ω in the complex ensemble, ǫ(-ω) = ǫ(ω) where the bar symbolizes the conjugated complex. For the purely imaginary pulsation iω, this relationship implies ǫ (iω) = ǫ (iω) and consequently that the imaginary part of

ǫ (iω) = ǫ ′ (iω) + i.ǫ ′′ (iω) is zero. Hence, ǫ (iω) = ǫ ′ (iω)
. This shows the reality of ǫ (iω) that is called dielectric response function. Its knowledge is indispensable to the application of the Lifshitz theory.

Hamaker constant according to Lifshitz

The Lifshitz theory is complete and exact in that it takes into account the multi-body interactions and the retardation effect. However, its complexity prevents any application to complex geometries. In what follows, we will study the particular case of two parallel semi-infinite plates separated by a very small distance d.

We have previously seen in equation ( 3) that the Hamaker method leads to an approximated expression of the interaction potential7

U H (d) = -A H .f (d,geometry)
The Lifshitz theory gives the exact expression of this potential. We call Hamaker constant in the sense of Lifshitz the quantity A L such that the exact potential U is

U (d) = -A L .f (d,geometry)
where f is the above defined function. Rigorously, A L depends on the distance d. This is thus not a constant and, for that reason, it is sometimes called Hamaker coefficient. However, when distance d is very small, this dependence disappears and A L tends to a value A L0 that we call non retarded Hamaker constant in the sense of Lifshitz. The coefficient A L (d) can be calculated with the following equations as well as its limit value A L0 . Using the field quantum theory and the Feynman diagrams technique, Dzyaloshinskii, Lifshitz and Pitaevski proposed the following formula

A L (d) = -3kT 2 ∞ n=0 ´ ∞ rn x. ln D (x) .dx
where the prime on the sum means that the term n = 0 has to be multiplied by

D (x) = 1 -∆ 13 .∆ 23 .e -x . 1 -∆ 13 .∆ 23 .e -x ∆ j3 = x.ǫ j(iωn) -x j .ǫ 3(iωn)
x.ǫ j(iωn) + x j .ǫ 3(iωn)

∆ j3 = x -x j x + x j x j = x 2 -1 - ǫ j(iωn) ǫ 3(iωn) .r 2 n r n = 2d. 2πkT.n h.c . √ ǫ 3(iωn)
ω n = 2πkT h .n When the distance d is very small, all the r n terms can be approximated by zero. Consequently, x 1 = x 2 = x and ∆ 13 = ∆ 23 = 0. The Hamaker constant becomes

A L0 = -3kT 2 ∞ n=0 ´ ∞ 0 x. ln D (x) .dx
where the function D is now

D (x) = 1 -∆ 13 .∆ 23 .e -x ∆ j3 = ǫ j(iωn) -ǫ 3(iωn) ǫ j(iωn) + ǫ 3(iωn)
As noticed by Hough & White [START_REF] Hough | [END_REF], the function D is comprised between 0 and 1 strictly. So, it is possible to develop this term in a series of powers and integrate it term by term to obtain

A L0 = 3kT 2 ∞ n=0 ´∞ s=1 (∆ 13 .∆ 23 ) s s 3
In case of truncation of the sum on s just after the first term, the error is below 15% since ∆ 13 .∆ 23 < 1. In general, this error is even lower than 5%. For this reason, the following expression is often written

A L0 ≈ 3kT 2 ∞ n=0 ´∆13 .∆ 23
However, the use of the computer is anyway required and it is not more difficult to entirely calculate the sum on s.

Because of the sampling (summation on n), the knowledge of the permittivities is not necessary on the whole spectrum. Notably, their evolution within the frequency range comprised between 0 and ν 1 = 3.578 × 10 13 s -1 at 0 • C8 does not play any role on the interaction potential. More precisely, the main contributions come from the zero-frequency term and the ultra-violet range (between 280 and 400 nm).

5 Lifshitz approach : calculation

Dielectric data for basic materials

To calculate the Hamaker constant according to the Lifshitz theory, we need the expression of the dielectric response function and the static dielectric constant for the following materials : methane gas molecules, liquid water molecules and water molecules in hydrate structure. For the latter, we will use the dielectric response function of water molecules in ice structure, as will be justified in section 5.2.1. Additionnally, we give the dielectric data for quartz and dodecane.

Liquid water

The dielectric response function ǫ 3(iω) of liquid water, which is the intermediate medium in our geometric configuration, can be calculated by different methods. We use here a common way, based on a parametric representation of the dielectric behaviour of materials, using a model of damped oscillators proposed by Ninham & Parsegian [18]

ǫ (ω) = 1 + d 1 -i.hω.τ + j f j ω 2 j -i.g j hω -(hω) 2
The second term of the sum, called Debye relaxation term, models the molecule behaviour in micro-waves. The relaxation time is noted τ and the oscillator strength d. In the third term of the sum , f j , ω j and g j represent respectively the oscillator strength, its resonance frequency and its bandwidth. This equation9 is valid for any number, real or complex, and in particular for the purely imaginary pulsation iω

ǫ (iω) = 1 + d 1 + hω.τ + j f j ω 2 j + g j hω + (hω) 2
Parameters for water, consigned in the table 4, are given by Parsegian [START_REF] Parsegian | Van der Waals Interactions[END_REF] and modified by Roth & Lenhoff [START_REF] Roth | [END_REF].

Concerning the static dielectric constant of liquid water, Klein & Swift [21] proposed the correlation ǫ 0 = 88.045-0.4147T +6.295×10 -4 T 2 +1.075×10 -5 T 3

where the temperature T is in Celsius. In what follows, we will take the value ǫ 30 = 88.045. The dielectric response function for liquid water is displayed in the figure 5.

Free gas

For the molecule of free gas (f), we use the classical model of a simple harmonic oscillator with an absorption frequency ν e = ωe 2π , which means that the gas has only one absorption line. We have then

ǫ f (ω) = 1 + n 2 -1 1 -ω ωe 2
where n is the refractive index in the visible range, approximately equal to ǫ (ωvis ). Hence, the dielectric response function is

ǫ f (iω) = 1 + n 2 -1 1 + ω ωe 2
According to Israelachvili [START_REF] Israelachvili | Intermolecular and Surface Forces[END_REF], the absorption frequency is about ν e = 3 × 10 15 Hz that is a pulsation ω e = 1.885 × 10 16 rad.s -1 . For methane, the index of refraction in the visible range is n = 1.000444 [START_REF]CRC Handbook of Chemistry and Physics[END_REF]. The static dielectric constant of free methane is ǫ f 0 = 1.000944 at 0 • C and 760 mm Hg [START_REF]CRC Handbook of Chemistry and Physics[END_REF].

Ice : parametric method

Elbaum & Schick [START_REF] Elbaum | [END_REF] have used experimental data of Kislovskii [24], Daniels [25] and Seki [26] to propose two parametric representations of the ice dielectric constant. The parameters given in the table 5 correspond to the formula

ǫ (ω) = 1 + j f j ω 2 j -ihωg j -(hω) 2
In the following, by Elbaum-Seki's representation, we will mean the parametric representation proposed by Elbaum & Schick using Kislovskii's data [24] for IR-range and Seki's ones for UV-range [26]. At 0 • C, ice is slightly birefringent; its static dielectric constant is equal to 106.4 along the c-axis and to 91.6 perpendicularly to the latter [START_REF] Israelachvili | Intermolecular and Surface Forces[END_REF]. The average is 1 3

(ǫ || + 2ǫ ⊥ )=96.5.
In what follows, we take the value ǫ i0 = 96.5.

Ice : direct calculation

The parametric representation of the complex dielectric constant is based on an oscillators model. This model is sometimes abusively simple and can lead to inaccurate representation of the material properties. Since the final result is very sensitive to the precision of these data, we will use another way. So, we hope, we will introduce as few errors as possible between the experimental data (ice complex refraction index compiled by Warren [27] over a broad wavelength spectrum) and the final output we want, that is to say the complex dielectric constant for imaginary frequencies.

Kramers-Kronig relationship Generally speaking, the imaginary and real parts of a response function are linked by the Kramers-Kronig transformations. This transformation has a mathematical origin and hence is not an approximation but a strict equality. It states

ǫ ′ (ω) = 1 + 2 π . ∞ 0 x.ǫ ′′ (x) x 2 -ω 2 .dx ǫ ′′ (ω) = -2ω π . ∞ 0 ǫ ′ (x) -1 x 2 -ω 2 .dx
On the other side, Landau & Lifshitz [START_REF] Landau | Electrodynamics of continuous media[END_REF] have established an other relation that looks similar

ǫ (iω) = 1 + 2 π . ∞ 0 x.ǫ ′′ (x)
x 2 + ω 2 .dx

We shall note here that the integrals are not ordinary ones. It is the principal part of Cauchy integral, which means that the discontinuity in ω is circumvented in the integration. The symbol ∞ 0 dx used for the sake of simplicity should, rigorously, be replaced by lim δ→0 ω-δ 0 dx + +∞ ω+δ dx . Practically, δ must be above 10 2 rad/s to avoid convergence problems and below 10 8 rad/s to ensure an accurate result. In what follows, we take δ = 10 5 rad/s. Adaptation of the Kramers-Kronig formula Practically, the experimental data of the ice refractive index of Warren [27] are not sufficient to calculate ǫ ′′ (x) on the whole spectrum but only in the interval Ω 1 = [ω 1 ; ω 4 ]. We may write

ǫ ′ (ω) = 1 + 2 π . ω1 0 x.ǫ ′′ (x)
x 2 -ω 2 .dx + ω4 ω1

x.ǫ ′′ (x)

x 2 -ω 2 .dx + ∞ ω4 x.ǫ ′′ (x)
x 2 -ω 2 .dx

Let us imagine that we want to know ǫ ′ (ω) for ω in the interval Ω 2 = [ω 2 ; ω 3 ] ⊂ Ω 1 such that ω2 ω1 2 is very high compared to 1 and ω3 ω4 2 is very small compared to 1. In this case, the first and third integrals may be approximated and we have

ǫ ′ (ω) ≈ 1 + 2 π . -K (ω1) ω 2 + I (ω) + Q (ω4)
where K and Q are defined by

K (t) = t 0 x.ǫ ′′ (x) .dx and Q (t) = ∞ t ǫ ′′ (x)
x .dx. Since ω 1 and ω 4 are known, K (ω1) and Q (ω4) are two constants . The integral I is equal to

I (ω) = ω4 ω1 x.ǫ ′′ (x) x 2 -ω 2 .

dx

The figure 3

represents the residual part R (ω) = π 2 ǫ ′ (ω) -1 -I (ω) calculated directly from Warren data for 200 values of ω ∈ Ω 2 = [ω 2 ; ω 3 ] = [ω 1 .D; ω 4 /D].
The calculations have been performed with D=50. No obvious dependence to ω can be observed. We deduce that

K (ω 1 ) ω 2 is negligible compared to Q (ω4)
. In what follows, we take K (ω1) = 0 and Q (ω4) = 0.07867.

We can now evaluate the dielectric constant on the imaginary axis ǫ (iω) . For that, we write

ǫ (iω) = 1 + 2 π . ω1 0 x.ǫ ′′ (x)
x 2 + ω 2 .dx + ω4 ω1

x.ǫ ′′ (x)

x 2 + ω 2 .dx + ∞ ω4 x.ǫ ′′ (x)
x 2 + ω 2 .dx and, as it has been done above, give the approximation

ǫ (iω) ≈ 1 + 2 π . K (ω1) ω 2 + ω4 ω1
x.ǫ ′′ (x)

x 2 + ω 2 .dx + Q (ω4) ≈ 1 + 2 π . J (ω) + Q (ω4)
where the integral J is defined by the relation

J (ω) = ω4 ω1
x.ǫ ′′ (x)

x 2 + ω 2 .dx
As illustrated on the figure 4, the calculation on the frequency domain of interest gives a result very close to that obtained with the Elbaum-Seki parametric representation [START_REF] Elbaum | [END_REF].

The dielectric response function of ice is compared to liquid water one in the figure 5.

Other materials

For other materials, like quartz and dodecane, Hough & White [START_REF] Hough | [END_REF] proposed the following representation

ǫ (iω) = 1 + C IR 1 + ω ωIR 2 + C UV 1 + ω ωUV 2
where the parameters are given in the table 6. The pulsations ω IR and ω UV are expressed in rad/s.

Dielectric data for hydrate

Dielectric response function

To calculate the hydrate dielectric response function, we will do the following three assumptions 1. The hydrate ( 1) is considered as the superposition of a water structure (e) and of gas (g). For any frequency, its dielectric permittivity ǫ 1 is given by the mixing rule of Clausius-Mosotti that expresses the additivity of the total polarisabilities

ǫ 1 -1 ǫ 1 + 2 = ǫ e -1 ǫ e + 2 + ǫ g -1 ǫ g + 2 2.
A water molecule (e) in the hydrate crystalline lattice has the same polarisability 10 as a water molecule in the ice crystalline lattice (i). The Clausius-Mosotti formula states that

α 4πǫ 0 = ǫ -1 ǫ + 2 . 3V 4π
Noting ǫ i the dielectric constant of the ice (i), we have

ǫ e -1 ǫ e + 2 = ǫ i -1 ǫ i + 2 . ρ e ρ i 3.
A gas molecule (g) in a hydrate cavity has the same polarisability as a free gas molecule (f) at atmospheric pressure and 0 • C. Thanks to Clausius-Mosotti formula, we have

ǫ g -1 ǫ g + 2 = ǫ f -1 ǫ f + 2 . ρ g ρ f
The sum of the two contributions allows us to calculate the hydrate dielectric response function. We can give an interpolation of this result. For ω ∈ 2.248 × 10 14 ; 4.255 × 10 16 rad.s -1 , it is correctly approximated by

ǫ 1(iω) = 8 n=0 a n .Log n 10 (hω)
The coefficients are in the table 7 The figure 5 gives a graphical summary of the dielectric response functions of liquid water, ice and hydrate.

Static dielectric constant

For water ice, the density is ρ i = 917 kg.m -3 . The water density in hydrate being ρ e = 46×18×10 -3 NA.(12.03×10 -10 ) 3 =790.6 kg.m -3 , we have ǫ e0 =16.29 for the static dielectric constant of water in hydrate. In the same logic, for gas, with ρ g = 8×16×10 -3 NA.(12.03×10 -10 ) 3 =122.1 kg.m -3 and ρ f = P M RT =0.7139 kg.m -3 , we have ǫ g0 =1.171 for the static dielectric constant of gas (g) in hydrate. Consequently, with the help of the mixing rule given above, we can calculate the hydrate static dielectric constant ǫ 10 = 17.46.

Hamaker constant

With all these data, we performed the numerical computations and found the Hamaker constant to be A L0 = 4.8647 × 10 -20 J when the intermediate medium is vacuum and A L0 = 8.2520 × 10 -21 J when it is liquid water. Expressed in multiples of kT, we obtain respectively 12.90 and 2.188 kT. Finally, the non retarded Hamaker constant in the sense of Lifshitz is

A L0 hwh = 8.2520 × 10 -21 J ≈ 2.188 kT
6 Hybrid approach : theory

General philosophy

As already mentioned, the full application of the DLP theory has two advantages (the retardation effect and the multi-body interactions are taken into account) but an analytical expression is available only for simple geometries. Inversely, the application of the Hamaker method is applicable for complex geometries but needs the additivity hypothesis and neglects the retardation effect. To circumvent these problems, we will make an hybridation 11 . The interaction potential will be approximated by

Ū(d) = A L0 .f (d,geometry) .Γ (d,geometry)
We must give some comments on this hybridation.

• The function f (d,geometry) is calculated with the Hamaker approach. The table 9 shows its value for a given geometrical configuration. For other geometries, even complex ones, its value can be easily calculated (Cf.

Tadmor [START_REF] Tadmor | [END_REF]).

• The constant A L0 is the non retarded Hamaker constant in the sense of Lifshitz, computed in the case of two parallel, semi-infinite plates. This way, the multi-body interactions and the presence of the intermediate (liquid water in our case) are taken into account.

• The function Γ (d,geometry) expresses the retardation effect. When the distance tends to zero, this function Γ tends to 1.

Retardation effect

To take into account the retardation effect in this hybrid approach, it is usual to modify the London equation [START_REF] Sloan | Clathrate Hydrates of Natural Gases[END_REF]. So, the interaction potential between two molecules takes the form u (r) = -C r 6 .g (r) . The function g tends to 1 when the intermolecular distance r tends to zero. The literature gives many expressions for the corrective function g (r) . For two infinite half-spaces, we can obtain the exact expression using the Lifshitz theory and a formal analogy. It is the same for two half spheres (Cf. Langbein [START_REF] Langbein | Theory of Van der Waals attraction[END_REF]). These expressions are however complicated and various approximations have been proposed : Casimir & Polder [6] used quantum electrodynamics to find the exact expression for two identical neutral atoms where only one excited state contributes to the London energy. It is the most accurate expression. Overbeek [START_REF] Th | Colloid Science[END_REF] simplified this expression and gave piece-wise functions on two domains. Schenkel & Kitchener [START_REF] Schenkel | [END_REF] proposed a unique function for the range λ ≥ 8 nm. Anandarajah & Chen [32] proposed a unique function valid between 0 and 1000 nm. These authors used then these corrective functions g to calculate the Γ function for different geometries.

Two parallel semi-infinite plates

For two parallel semi-infinite plates separated by a distance d, Overbeek [33] used his function g to find the retarded interaction potential d'

Γ o (d) = 1.01 -1.12ad + 0.696a 2 d 2 -0.229a 3 d 3 + 0.031a 4 d 4 Γ o (d) = 0.490 ad - 0.102 a 2 d 2
respectively for d ≤ 3 2a and for d > 3 2a . The constant a is equal to a = π λ with λ = 100 nm. Anandarajah & Chen [32] used their own corrective function g and came to the expression :

Γ ac (d) = 1 - 4d c - 18d 2 c 2 - 12d 3 c 3 + 12d 2 c 4 .(d + c) 2 . ln 1 + c d c = bλ
2π with λ = 100 nm and b = 3.1.

A sphere and a semi-infinite plate

A sphere of radius r, separated by a distance d of a semi-infinite plate, is submitted to an interaction potential given in the equations ( 7), ( 8) and ( 9) of Chen & Anandarajah [34], who use their corrective function g.

Two spheres

Two spheres of radius r 1 and r 2 , which centers are at a distance h are submitted to an interaction potential given in the equations ( 18), ( 19) and (20) of Chen & Anandarajah [34]. In their article, this expression is compared to the exact expression that Clayfield et al. [35] proposed for two spheres using the corrective function of Casimir & Polder [6]. Langbein [START_REF] Langbein | Theory of Van der Waals attraction[END_REF] applied the Lifshitz theory to this geometric configuration and obtained the exact solution. However it requires heavy calculations and numerical convergence problems make it difficult to use. The Lifshitz theory can also be applied indirectly using the Derjaguin approximation [36] or the one of Papadopoulos-Cheh [37]. The article of Thennadil [38] proposes a complete review of their calculation methods and their respective advantages.

Hybrid approach : calculation

As an example for this hybrid approach, we consider two spherical hydrate particles. According Chen & Anandarajah [34], two spheres of radius r 1 and r 2 placed at a center-to-center distance h are submitted to the interaction potential

U (h) = A L0 . [U 0 (h) + U 1 (h)] U 0 (h) = -1 6 . 2r 1 r 2 d 1 d 3 + 2r 1 r 2 d 2 d 4 + ln d 1 d 3 d 2 d 4 U 1 (h) = -1 30c 4 h
. 4cr 1 r 2 (4r 2 1 + 4r 2 2 + 2h 2 + 7hc + 9c 2 ) + c 4 (c + 5h) ln

d 2 d 4 d 1 d 3 + 4 i=1 b i ln e i d i
with following notations : The figures 6 and 7 show the absolute value of the interaction potential as well as the attractive force between two equal spheres, when the radius varies.

d 1 = h -r 1 -r 2 d 2 = h + r 1 -r 2 d 3 = h + r 1 + r 2 d 4 = h -r 1 + r 2 e i = d i + c b i = (-1) i+1 e 4 i [4e i -5(h + c)] -

Discussion

Concerning hydrate, we find a value of ǫ 10 = 17.46 for its static dielectric constant. This value is 70% smaller than the value of 58 found experimentally by Makogon [START_REF] Makogon | Hydrates of hydrocarbons[END_REF]. However, these values are both credible since they both are between the static dielectric constant of liquid water (88.045) and of gas (1.171), dodecane (1.991) or quartz (4.29). With Makogon's value, the non retarded Hamaker constant is A L0 = 7.03 × 10 -21 J instead of 8.25 × 10 -21 J, that is a difference of only 15%.

We shall cite here the work of Camargo & Palermo [START_REF] Camargo | Proc. 4th Int. Conf. on Gas Hydrates[END_REF] who use Israelachvili's formula [START_REF] Israelachvili | Intermolecular and Surface Forces[END_REF] to get a first estimate of the Hamaker constant. They found a value of 5.2 × 10 -21 J, which is used in Yang's work et al. [START_REF] Yang | Micromechanical measurements of hydrate particle attractive forces, 15th Symp. on Thermophysical Properties[END_REF]. Qualitatively, our results valid their estimate a posteriori.

Concerning the ice dielectric response function, which is the base to calculate hydrate's one, we have seen that our result, obtained by application of the Kramers-Kronig transformations, is slightly different from the result obtained using the parametric representation that Elbaum & Schick [START_REF] Elbaum | [END_REF] built from Seki's data [26]. If we use this dielectric response function, the Hamaker constant would be A L0 = 8.36 × 10 -21 J instead of 8.25 × 10 -21 J that is a difference of the order of 1%. Here again, the robustness of our result seems reassessed.

If we use Makogon's value [START_REF] Makogon | Hydrates of hydrocarbons[END_REF] for the ice static dielectric constant and the parametric representation of Elbaum & Schick [START_REF] Elbaum | [END_REF] for ice, then the Hamaker constant would be A L0 = 7.13 × 10 -21 J instead of 8.25 × 10 -21 J that is a difference of the order of 14%.

We can extend our work to systems containing a third material M . The table 8 shows the Hamaker constants A L0 for Hydrate-Water-M or Hydrate-Vacuum-M systems. The dielectric response functions ǫ (iω) are from Hough & White [START_REF] Hough | [END_REF] for dodecane, quartz, sapphire and teflon. For the three metals (gold, silver and copper), data are from Parsegian & Weiss [START_REF] Parsegian | [END_REF].

Finally, we can predict the most probable configurations. For instance, the Hamaker constant of the system Quartz-Hydrate-Water (6.509 × 10 -21 J) is higher than that of the system Quartz-Water-Hydrate (2.890 × 10 -21 J). It means that quartz attracts more water than hydrate and thus that the stable situation is an hydrate particle separated from quartz by a liquid film. Also, the Hamaker constant of the Hydrate-Dodecane-Water system(10.35 × 10 -21 J) is positive. It means that, in a water in dodecane emulsion, water droplets will agglomerate with hydrate particles. Inversely, for an dodecane in water emulsion, dodecane droplets will tend to agglomerate to hydrate particles. In this case, we shall notice that the Hamaker constant of the Hydrate-Dodecane-Water system (10.35 × 10 -21 J) is higher than for the Hydrate-Water -Dodecane system (1.989 × 10 -21 J). Consequently, the hydrate attracts more water than dodecane and thus the stable situation is an hydrate particle with a double layer : a water film in contact with the hydrate particle and a dodecane layer on the outside. This situation will happen only if the dodecane surface tension is small enough to allow the droplet-to-shell transition. In the same logic, data from table 8 show that hydrates in suspension in water will aggregate, in order of preference, with metal, sapphire and quartz. Experiment with unconsolidated quartz grains placed in a metallic reactor have to be interpreted with precaution.

Conclusion

In this article, we give an insight of Van der Waals forces acting in systems involving gas hydrates. The microscopic, Hamaker approach to calculate the Hamaker constant of the hydrate-water-hydrate system leads to a value of 4.59 × 10 -21 J which is close to the value of 8.25 × 10 -21 J found with the macroscopic, Lifshitz theory. To apply the latter, we used the Kramers-Kronig relationship and adapted it to use the full information contained in optical experimental data of ice, since water molecules frozen in a rigid structure appear in both ice and gas hydrate solids. Although the Hamaker constant obtained with the Lifshitz theory appears to us more representative of the reality, the difference between the two values is rather low in terms of agglomeration kinetics. In effect, the characteristic agglomeration time is proportional to the inverse agglomeration efficiency α, which is proportional to the Hamaker constant at the power 0.18 [43]. Numerically, this means that the two Hamaker constants lead to a difference of 10 % only between the two characteristic agglomeration times. We then applied this theory to different other systems. Notably, we studied the geometry of two spherical hydrate particles interacting across water and we compared different three-phases systems to predict what configurations would be more stable. In this optic, we applied the theory to systems containing hydrates and water but also ice, dodecane, quartz, sapphire, teflon and metals. These data may appear of interest for industrial and experimental purposes.

1 2 .

 2 The function D (x) is equal to

20r 1 r 2 e 3 i

 3 c = bλ 2π with b=3.1 and λ=100 nm.

Figure 1 Figure 2

 12 Figure 1 Van der Waals constant between hydrate molecules

Figure 3 Figure 4 Figure 5

 345 Figure 3 Residual part R ω to calculate the ice dielectric response function

Figure 6

 6 Figure 6 Interaction potential between two hydrate spheres

Figure

  Figure 7 Attractive force between two hydrate spheres

6 ]Figure 1 : 6 ]Figure 2 :Figure 3 :

 61623 Figure 1: Van der Waals constant between hydrate molecules

Figure 4 :

 4 Figure 4: Ice dielectric response function (Elbaum-Seki [23] and this work)

Figure 5 :

 5 Figure 5: Dielectric response functions of liquid water, ice and hydrate

Figure 6 :

 6 Figure 6: Interaction potential between two hydrate spheres

Figure 7 :

 7 Figure 7: Attractive force between two hydrate spheres

  is about 1.3130 × 10 -19 J with a molecular density of ρ w = 3.342 × 10 28 molecules/m 3 for water. The constant A H ww is equal to A H ww = π 2 .C ww .ρ 2 w or 1.5322 × 10 -19 J. Finally, the Hamaker constant of the hydrate/water/hydrate system is :

7

  Attractive force between two hydrate spheresC ind C orient C disp C total

	H 2 O -H 2 O	10	96	33	139
	CH 4 -CH 4	0	0	102	102
	H 2 O -CH 4	9	0	58	67

Table 1 :

 1 Contributions to the induction, orientation and dispersion interactions (10 -79 J.m 6 )

	Geometry	f (d,geometry)
	Two spheres	2r1r2 d 2 +2r1d+2r2d +	2r1r2 d 2 +2r1d+2r2d+4r1r2 + ln	d 2 +2r1d+2r2d d 2 +2r1d+2r2d+4r1r2	/6
	Two close spheres	r 1 r 2 / [6d(r 1 + r 2 )]
	Sphere -Half space	r d + r d+2r + ln	d d+2r	/6
	Sphere -close Half space r/(6d)
	Two plates	1/d 2 + 1/ (d + 2e) 2 -2/ (d + e) 2 /(12π)
	Two half spaces	1/(12πd 2 )

Table 2 :

 2 Geometric function f (d,geometry)[START_REF] Elimelech | Particle Deposition and Aggregation[END_REF] 25

		C	ρ	A ii
		(10 -79 J.m 6 ) (10 28 m -3 ) (10 -19 J)
	Hydrocarbons	50	3.3	0.5
	Water	139	3.3	1.5

Table 3 :

 3 Hamaker constant for water and hydrocarbons in vacuum 10 -2 6.25 × 10 -4 1.5 × 10 -2 6.9 × 10 -2 3.50 × 10 -3 3.8 × 10 -2 9.2 × 10 -2 1.28 × 10 -3 2.8 × 10 -2 2.0 × 10 -1 5.69 × 10 -3 2.5 × 10 -2 4.2 × 10 -1 1.35 × 10 -2 5.6 × 10 -2 d=74.8 1/τ = 1.05 × 10 11 rad.s -1

	ω j	f j	g j
	(eV )	(eV 2 )	(eV )
	8.21	3.26	0.63
	10.0	3.87	0.84
	11.4	12.0	2.05
	13.6	63.6	3.9
	17.8	114	7.33
	25.2	24.3	5.34
	2.07 ×

Table 4 :

 4 Parameters for the dielectric representation of liquid water

	ω j	f j	g j
	(eV )	(eV 2 )	(eV )
		IR Parameters [24]
	1.716 × 10 -2	2.110 × 10 -3	8.553 × 10 -2
	0.1140	6.754 × 10 -3	5.881 × 10 -2
	0.4003	1.089 × 10 -2	2.594 × 10 -2
		UV Parameters [25]
	8.740	3.081	0.7672
	11.55	11.25	2.666
	14.49	43.12	4.374
	17.54	47.93	5.291
	21.81	41.25	8.066
	32.83	144.8	22.80
		UV Parameters [26]
	8.603	5.027	0.7308
	10.43	5.124	1.342
	12.79	8.885	1.968
	15.06	48.44	3.701
	18.07	54.57	4.688
	22.82	39.92	6.876
	32.97	85.35	15.93

Table 5 :

 5 Parameters for the dielectric representation of ice [23] Quartz 1.93 2.093 × 10 14 1.359 2.032 × 10 16 4.29 Dodecane 0.026 5.540 × 10 14 0.965 0.873 × 10 16 1.991

	C IR	ω IR	C UV	ω UV	ǫ 0

Table 6 :

 6 Parameters for the dielectric representation of quartz and dodecane

	a 0	1.8277	a 5	-0.20779
	a 1 -0.09573 a 6	0.11237
	a 2 0.07746	a 7	0.15226
	a 3	-0.1089	a 8 -0.070614
	a 4 -0.16035

Table 7 :

 7 Coefficients for the polynomial interpolation of the hydrate dielectric

	response function
	M	H-W-M H-V-M
	Hydrate	8.252	9.460
	Dodecane	1.989	0.899
	Quartz	2.890	2.780
	Teflon	1.991	0.874
	Sapphire	4.506	5.329
	Gold	5.018	9.711
	Copper	5.029	9.722
	Silver	5.031	9.724
	Ice	6.975	9.923

Table 8 :

 8 Hamaker constants for different systems involving hydrates (10 -21 J)

The revolution time of the electron of the Bohr atom is equal to the inverse of the first ionisation frequency which is about 3 × 10 15 s -1 .

The molecular density, expressed in molecules/m

, is equal to ρ = N A .ρ V Mwhere ρ V is the density (kg.m -3 ), N A = 6.02283 × 10 23 the Avogadro number and M the molar mass of the component (kg.mole -1 )

These distances are to be compared with the hydrate lattice parameter which is 1.203 nm for sI hydrates.

The center of an hydrate molecule is the position of a given water molecule that constitute it. Refer to the description of the hydrate crystalline structure by Stackelberg[START_REF] Stackelberg | Feste Gashydrate II. Struktur und Raumchemie[END_REF].

Hamaker constants are sometimes expressed in multiples of kT. At 0 • C, the kT product is 3.7708 × 10 -21 J.

The pulsation ω, the frequency ν and the wavelength λ are connected by the relations ω = 2πν and λ = c ν .

In our geometrical configuration, the function f is 1/(12πd 2 )

It corresponds to a pulsation lower than 2.248×10 14 rad.s -1 or a wavelength greater than 8.385 µm or a photon energy below 0.1479 eV. The constant h 1eV is equal to 6.5789 × 10 -16 .

The result of the formula is correct when the parameters (ω j , f j , g j ) are expressed in Joule (1eV = 1.60252 ×

-19 J) and the pulsation in rad.s -1 (h = 1.0543 × 10 -34 J.s.rad -1 )

Polarisability means total polarisability (electronic and orientation)

The method and the expression is from Parsegian[START_REF] Parsegian | Van der Waals Interactions[END_REF]

Footnotes

Note 1 The revolution time of the electron of the Bohr atom is equal to the inverse of the first ionisation frequency which is about 3 × 10 15 s -1 .

Note 2

The molecular density, expressed in molecules/m 3 , is equal to ρ =

NA.ρV M

where ρ V is the density (kg.m -3 ), N A = 6.02283 × 10 23 the Avogadro number and M the molar mass of the component (kg.mole -1 ) Note 3 These distances are to be compared with the hydrate lattice parameter which is 1.203 nm for sI hydrates.

Note 4

The center of an hydrate molecule is the position of a given water molecule that constitute it. Refer to the descripon of the hydrate crystalline structure by Stackelberg [START_REF] Stackelberg | Feste Gashydrate II. Struktur und Raumchemie[END_REF].

Note 5 Hamaker constants are sometimes expressed in multiples of kT. At 0 • C, the kT product is 3.7708 × 10 -21 J.

Note 6

The pulsation ω, the frequency ν and the wavelength λ are connected by the relations ω = 2πν and λ = c ν .

Note 7 In our geometrical configuration, the function f is 1/(12πd 2 ) Note 8 It corresponds to a pulsation lower than 2.248 × 10 14 rad.s -1 or a wavelength greater than 8.385 µm or a photon energy below 0.1479 eV. The constant h 1eV is equal to 6.5789 × 10 -16 .

Note 9

The result of the formula is correct when the parameters (ω j , f j , g j ) are expressed in Joule (1eV = 1.60252 × 10 -19 J) and the pulsation in rad.s -1 (h = 1.0543 × 10 -34 J.s.rad -1 )

Note 10 Polarisability means total polarisability (electronic and orientation)

Note 11

The method and the expression is from Parsegian [START_REF] Parsegian | Van der Waals Interactions[END_REF] List of Tables Ice dielectric response function (Elbaum-Seki [START_REF] Elbaum | [END_REF] and this work) 31 5 Dielectric response functions of liquid water, ice and hydrate . . 32 6 Interaction potential between two hydrate spheres . . . . . . . . 32 7 Attractive force between two hydrate spheres . . . . . . . . . . . 33


Ghassan Jomaa

Patrick Goblet
email: patrick.goblet@mines-paristech.fr

Christophe Coquelet

Vincent Morlot

Kinetic modeling of polyurethane pyrolysis using non-isothermal thermogravimetric analysis

Keywords: kinetic parameters, polyurethane, pyrolysis, thermogravimetry, kinetic modeling

The pyrolysis of polyurethane was studied by dynamic thermogravimetry analysis (TGA). The studied polyurethane is used as organic binder in casting process to make sand cores and molds. A semi-empirical model is presented that can be used to describe polyurethane pyrolysis occurring during TGA experiments. This model assumes that the polyurethane is pyrolysed by several parallel independent reactions. The kinetic parameters of polyurethane pyrolysis were evaluated by fitting the model to the experimental data obtained by TGA over a wide variety of heating rates. A nonlinear least-squares optimization method is employed in the fitting procedure. A hybrid objectives based simultaneously on the mass (TG) and mass loss rate (DTG) curves has been used in the least-squares method. The values of the activation energy obtained by the nonlinear fitting were then recalculated by the methods of Kissinger and Friedmand. Furthermore, the parameters obtained in the present paper were then compared with those reported in the literature.

Introduction

Polyurethane is widely used in foundry industry as an organic binder to harden sand cores.

The latter are inserted into a metallic mold to obtain internal shapes of casting parts. Pouring of a molten metal into the mold causes the thermal decomposition of polyurethane and gas emissions which can represent a severe problem for the quality of the casting parts. Indeed, under the effect of the thermal decomposition of the polyurethane, the pressure of the cores increases (due to the produced gas). If the local gas pressure in the sand cores exceeds the local metallostatic pressure of the solidifying liquid metal at core-casting part interface, gas bubbles can grow into the metal. Depending on whether the gas bubbles escapes through the metal or not, two possible scenarios can occur. In the first one the gas can escape through the casting part, so its effect might A c c e p t e d M a n u s c r i p t be limited to the formation of bubble trails defect [START_REF] Divandari | Mechanisms of bubble damage in castings[END_REF][START_REF] Campbell | Complete Casting Handbook : Metal Casting Processes, Techniques and Design[END_REF]. In the second one, the gas bubbles stay entrapped within the solidified metal. The entrapped gas bubbles contribute thus to the formation of blowholes defects in the casting. Both bubble trails and blowholes defects affect the quality of the casting parts by deteriorating their mechanical properties and thus these parts are generally rejected. Nowadays, in order to prevent casting defects caused by the gas emissions and to assure the soundness of castings, one can use numerical simulation which allows the determination of the velocity and the pressure of the gas in the sand cores. The models used in the numerical simulation of gas emissions are based on differential equations which describe the transport of mass, momentum and energy within the sand cores. These equations contain source terms including the thermal decomposition rate of the organic binder which are calculated using an appropriate kinetic model. The reaction scheme and the parameters of the kinetic model have to be determined experimentally. Therefore, the study of the thermal decomposition of the polyurethane represents a necessary step for the numerical simulation of the gas emissions occurring during the pouring and the solidification of the casting parts.

Generally after pouring the molten metal into the mold, the sand cores are totally covered by the metal. As there is no oxygen within the sand cores, the thermal decomposition of the polyurethane binder can be assumed to be done only by pyrolysis [START_REF] Campbell | Complete Casting Handbook : Metal Casting Processes, Techniques and Design[END_REF]. One of the most common used thermal analysis techniques to study pyrolysis kinetics of organic solids is the thermogravimetry analysis (TGA) [START_REF] Flynn | General treatment of the thermogravimetry of polymers[END_REF]. TGA technique consists in measuring the mass of a substance as a function of temperature (or time) while the substance is subjected to a controlled temperature programme.

The thermogravimetry analysis of a solid matter can be conducted either at isothermal conditions (constant temperature) or non-isothermal/dynamics conditions (generally the temperature varies linearly with time). The investigation of the pyrolysis by TGA is carried out in inert atmosphere such as nitrogen, argon or Helium. However, one can study combustion by using a reactive atmosphere such as oxygen or hydrogen. The results obtained by TGA are principally the curves of mass loss (noted TG) and of rate of mass loss (note DTG) as functions of temperature (or time).

Even though several works have been carried out to determine the kinetic parameters of the pyrolysis of polyurethanes used in applications other than the foundry [START_REF] Font | Pyrolysis study of polyurethane[END_REF], only experimental results have been reported about the pyrolysis of polyurethanes used in foundry [START_REF] Lytle | Determination of thermal decomposition products from a phenolic urethane resin by pyrolysis-gas chromatography-mass spectrometry[END_REF][START_REF] Dungan | Pyrolysis of foundry sand resins : A determination of organic products by mass spectrometry[END_REF][START_REF] Wang | Characterization of hydrocarbon emissions from green sand foundry core binders by analytical pyrolysis[END_REF]. The objective of this work is to propose a kinetic model of the polyurethane pyrolysis and to calculate the kinetic parameters using TG and DTG curves obtained by thermogravimetry analysis of polyurethane. The determined kinetic parameters could be used later within a computational fluid dynamics model for gaseous emissions occurring in sand core during casting process [START_REF] Jomaa | Étude des dégagements gazeux survenant pendant la coulée de pièces d'aluminium[END_REF].

Experimental results

The polyurethane binder studied in this work was prepared by the reaction of phenolic resin and polyisocyanate resin. Both resins are in a liquid form and in combination with organic solvents. A blend of an equal weight of phenolic and polyisocyanate resins were mixed to form a reaction mixture. The formation of polyurethane occurred in the presence of a gaseous amine catalyst (dimethylethanolamine). The obtained solid polyurethane was ground using a mortar and pestle. Three samples were prepared with initial masses of 16. )

β=20 • C min -1 β=60 • C min -1 β=80 • C min -1 (
Temperature ( • C) β=20 • C min -1 β=60 • C min -1 β=80 • C min -1 (b) DTG curves

Kinetic Models of organic solid pyrolysis

There are abundance of published research works dealing with the modeling of the pyrolysis of solid matter using TGA. In these works, several kinetic models were employed to describe the pyrolysis of organic solids. We can mention among the used kinetic models : the single global reaction model [START_REF] Caballero | Kinetics of the thermal decomposition of tannery waste[END_REF][START_REF] Park | Kinetic analysis of thermal decomposition of polymer using a dynamic model[END_REF][START_REF] Conesa | Comments on the validity and utility of the different methods for kinetic analysis of thermogravimetric data[END_REF][START_REF] Font | Kinetic law for solids decomposition. application to thermal degradation of heterogeneous materials[END_REF], the parallel independent reactions model [START_REF] Font | Comparison between the pyrolysis products obtained from different organic wastes at high temperatures[END_REF][START_REF] Várhegyi | Kinetic modeling of biomass pyrolysis[END_REF][START_REF] Rodríguez | Kinetic models for pyrolysis and combustion of sewage sludge[END_REF][START_REF] Hu | Kinetic study of chinese biomass slow pyrolysis : Comparison of different kinetic models[END_REF][START_REF] Li | Kinetic study of corn straw pyrolysis : Comparison of two different three-pseudocomponent models[END_REF], and the distributed activation energy model [START_REF] Miura | A simple method for estimating f(e) and k0(e) in the distributed activation energy model[END_REF][START_REF] Scott | An algorithm for determining the kinetics of devolatilisation of complex solid fuels from thermogravimetric experiments[END_REF][START_REF] Sonobe | Kinetic analyses of biomass pyrolysis using the distributed activation energy model[END_REF][START_REF] Cai | New distributed activation energy model : Numerical solution and application to pyrolysis kinetics of some types of biomass[END_REF][START_REF] Várhegyi | Thermal decomposition of wheat, oat, barley, and brassica carinata straws. a kinetic study[END_REF]. In the single global reaction model, the pyrolysis of organic solid is described by only one reaction. This model is the simplest kinetic model that one can use. It is generally employed to describe thermal decomposition of a pure homogeneous solid.

The parallel reaction model assumes that the organic solid is constituted of many solid fractions (pseudo-components) each of which is decomposed independently by one reaction. Therefore, the interaction of pyrolysis reactions of different pseudo-components is neglected. This model is essentially used to describe the pyrolysis of organic solid blends and mixtures. In the distributed activation energy model, the organic matter is assumed to be decomposed by an infinite number of parallel independent reactions. Each reaction has its own activation energy and the sum of all reactions is given by a distributed activation energy. This model is used generally to take into account the heterogeneity of the studied solid matter.

Kinetic Model of polyurethane pyrolysis

In this section, the kinetic model used in our study to describe polyurethane pyrolysis will be presented. This model is based on a parallel independent reactions model. In this model, the organic solid, noted M, will be considered as a mixture of np pseudo-components, noted M j .

The pseudo-components are supposed to decompose independently from one another, and not to A c c e p t e d M a n u s c r i p t influence each other. The pyrolysis of each pseudo-component produces a gas noted G j and a solid noted C j . Therefore the pyrolysis of M could be described by the reaction scheme R. 1.

M j -→ s j C j + (1 -s j ) G j j = 1, ..., np (R. 1)
Where s j is the mass-based stoichiometric coefficient of pyrolysis reaction of solid pseudocomponent M j . np is the number of parallel reactions. The conversion rate of each reaction j is given by the following equation :

dα M j dt = k j 1 -α M j n j j = 1, np (1)
Where n j is the order of the reaction. α M j is the conversion ratio of pyrolysable part of pseudocomponent M j . It is defined as :

α M j = m 0 M j -m M j m 0 M j -m ∞ M j j = 1, ..., np (2)
Where m 0 M j is the initial mass of M j . m M j the mass of M j during the pyrolysis reaction. m ∞ M j is the final mass of M j when the reaction is complete.

The kinetic coefficients k j are function only of temperature. k j are given in an Arrhenius form :

k j = Z j exp -Ea j RT j = 1, np (3)
Where Z j , Ea j and R are the pre-exponential factors, activation energies, and the universal gas constant, respectively. The overall conversion rate of the solid matter M is given by a linear combination of conversion rates of M j and, is expressed by the equation [START_REF] Font | Pyrolysis study of polyurethane[END_REF].

dα M dt = np j=1 c j Z j exp -Ea j RT 1 -α M j n j (4)
Where the coefficient c j represents the fraction of the overall mass loss due to the pseudocomponent M j . The sum of these coefficients has to be equal to unity :

np j=1 c j = 1 (5)
In the same way, the conversion ratio of the solid matter M is given by a linear combination of conversion ratio of different pseudo-components M j :

α M = np j=1 c j α M j (6)
α M can be calculated using TGA results by the following equation :

α M = m 0 M -m M m 0 M -m ∞ M (7)
A c c e p t e d M a n u s c r i p t m 0 M , m ∞ M j and m M j are the initial, final and current weights of M. If we consider that the TGA experiments have been carried out in non-isothermal conditions, the conversion rate of the organic matter M might be written as :

dα M dt = dT dt dα M dT = β dα M dT (8)
Therefore, we have :

dα M dT = 1 β np j=1 c j Z j exp -Ea j RT 1 -α M j n j (9)
where β = dT/dt is the heating rate of sample during the TGA experiments. If the sample is subjected to a linear temperature ramp, the temperature can be expressed by the following equation :

T = β t + T 0 ( 10
)
Where t is time and T 0 is the initial temperature of the sample. For each reaction, there are four kinetic parameters (c j , Z j , Ea j and n j ).

Using the definition of ratio conversion α M given by the equation ( 7), the rate conversion of M can be evaluated by equation [START_REF] Conesa | Comments on the validity and utility of the different methods for kinetic analysis of thermogravimetric data[END_REF].

dα M dT = -1 m 0 M -m ∞ M dm M dT (11)
If we consider also np as an unknown, the model would have 4 × np + 1 kinetic parameters which have to be determined by fitting kinetic model equations to experimental data obtained by TGA.

Determination of kinetic parameters

There are several methods for determining the kinetic parameters of the pyrolysis of solid materials. We may classify them as linear and nonlinear methods. It is worthy to note that this classification is based upon the mathematical analysis of the experimental results. An classification according to experimental conditions of pyrolysis analysis such as isothermal and nonisothermal may be also used [START_REF] Vyazovkin | Isothermal and non-isothermal kinetics of thermally stimulated reactions of solids[END_REF].

Linear Methods

These methods consist of finding a linear relationship between the kinetic parameters using the reaction rate given by the considered kinetic model. Then using the experimental results the coefficients of the linear relation are determined by a linear regression. Those method are generally applied for a simple kinetic model in which the solid material is supposed to be decomposed in a single step. It is so described by one elementary reaction. This is the case where we have one reaction in the pyrolysis model presented above. In this case the reaction rate is given by the following equation :

dα M dT = Z β exp -Ea RT (1 -α M ) n (12)
A c c e p t e d

M a n u s c r i p t

This expression of reaction rate is used when the reaction mechanism is characterized by a homogeneous chemical kinetic [START_REF]Analysis of non-isothermal reaction kinetics : Part 1. simple reactions[END_REF]. For more general case of other reaction mechanisms, the reaction rate can be expressed as shown below :

dα M dT = Z β exp -Ea RT f (α M ) (13)
Where f (α M ) is the reaction function which can take many mathematical forms depending on the controlling reaction mechanism. For example, some forms of this function can be found in [START_REF] Vyazovkin | Isothermal and non-isothermal kinetics of thermally stimulated reactions of solids[END_REF][START_REF] Popescu | Integral method to analyze the kinetics of heterogeneous reactions under non-isothermal conditions a variant on the ozawa-flynn-wall method[END_REF][START_REF] Galwey | Chapter 3 kinetic background to thermal analysis and calorimetry[END_REF][START_REF] Zhang | Thermal decomposition kinetics of the Zn(II) complex with norfloxacin in static air atmosphere[END_REF][START_REF] Órfão | Kinetic analysis of thermogravimetric data obtained under linear temperature programming-a method based on calculations of the temperature integral by interpolation[END_REF][START_REF] Vlaev | Non-isothermal kinetics of pyrolysis of rice husk[END_REF][START_REF] Chrissafis | Kinetics of thermal degradation of polymers[END_REF].

Many linear methods have been reported in the literature [START_REF] Vyazovkin | Isothermal and non-isothermal kinetics of thermally stimulated reactions of solids[END_REF][START_REF] Chrissafis | Kinetics of thermal degradation of polymers[END_REF][START_REF] Khawam | Application of solid-state kinetics to desolvation reactions[END_REF][START_REF] Yao | Thermal decomposition kinetics of natural fibers : Activation energy with dynamic thermogravimetric analysis[END_REF][START_REF] Sbirrazzuoli | Integral, differential and advanced isoconversional methods : Complex mechanisms and isothermal predicted conversion-time curves[END_REF][START_REF] Slopiecka | Thermogravimetric analysis and kinetic study of poplar wood pyrolysis[END_REF]. The Friedman and Kissinger Methods are presented above.

Friedman Method

This method was introduced by Friedman [START_REF] Friedman | Kinetics of thermal degradation of char-forming plastics from thermogravimetry. application to a phenolic plastic[END_REF]. The linear relationship employed by this method is obtained by taking the logarithm of the equation ( 13) which leads to the expression [START_REF] Várhegyi | Kinetic modeling of biomass pyrolysis[END_REF].

ln β dα M dT = lnZ + ln f (α M ) - Ea RT (14)
For a given conversion degree α M , plotting ln β dα M dT as a function of 1 T for several heating rates yields a straight line whose slope is equal to -Ea R . Therefore the activation energy of the reaction Ea is obtained from the slope. We note that the estimation of the activation energy by

Friedman method does not require to know the reaction function f (α M ). This is why this method is considered as a free-model method. In order to estimate the pre-exponential factor Z from the y-intercept of the straight line, the function f (α M )) must be known.

Kissinger Method

This method was firstly introduced by Kissinger [START_REF] Kissinger | Reaction kinetics in differential thermal analysis[END_REF]. It is based on the temperature of the maximum rate of conversion (temperature of peaks in DTG). Indeed, the derivative of the conversion rate is equal to zero at the peak temperature of the DTG curve. The derivative of the conversion rate is given by the equation [START_REF] Rodríguez | Kinetic models for pyrolysis and combustion of sewage sludge[END_REF].

d 2 α M dT 2 = Ea RT 2 + Z β exp -Ea RT d f (α M ) dα M dα M dT (15)
If we assume that the maximum reaction rate occurs at the peak temperature T max , as the derivative of the conversion rate at T max is equal to zero, we obtain the following equation :

Ea RT 2 max - m Z β exp -Ea RT max = 0 ( 16
)
Which may also be rearranged in the following form

β T 2 max = m Z R Ea exp -Ea RT max (17
ln β T 2 max = ln m Z R Ea - Ea RT max (18)
For a set of DTG curves with different heating rates β, the plot of ln β T 2 max as a function of 1 T max would lead to a straight line whose slope -Ea R gives the activation energy. Kissinger method is also considered as model-free method. If the form of the function f (α M ) is known, the preexponential factor Z constant can be further determined from the y-intercept of the obtained straight line.

Non-linear Methods

This type of methods is based on nonlinear least-square regression [START_REF] Anderson | Comparison of different non-linear evaluation methods in thermal analysis[END_REF][START_REF] Anderson | Effective search of starting values for kinetic parameters estimation[END_REF][START_REF] Opfermann | Kinetic analysis using multivariate non-linear regression. i. basic concepts[END_REF][START_REF] Várhegyi | Kinetics of the oxidation of bismuth tellurite, bi2teo5 : Empirical model and least squares evaluation strategies to obtain reliable kinetic information[END_REF][START_REF] Reverte | Kinetic model identification and parameters estimation from {TGA} experiments[END_REF][START_REF] Zhang | Research on pyrolysis characteristics and kinetics of super fine and conventional pulverized coal[END_REF]. The kinetic parameters are calculated in this method by minimizing an objective function based on the sum of weighted square of the error between the experimental data of thermogravimetry analysis and model equations. These methods allow a direct fitting of the model equations to the experimental data without any rearrangement or approximation. Moreover, they are robust methods to estimate kinetic parameters of models with complex reaction schemes .

Depending on whether the conversion ratio or the conversion rate of the pseudo-compounds is used, the objective function may be defined in two ways. In the first one, the objective function is defined by the conversion ratio :

O f I (a) = 100 × N TGA l=1 N l exp i=1 α exp,l M (t i ) -α l M (t i , T i , a) 2 N l exp ( 19
)
a is the vector of the kinetic parameters to be estimated. It is defined as follows :

a =
O f D (a) = 100 × N TGA l=1 N exp i=1 dα exp,l M (t i ) dt - dα l M (t i ,T i ,a,α α α) dt 2 N l exp dα l M (t i ) dt | 2 max ( 20
)
This form of the objective function is called the differential form. The division by

dα exp,l M (t i ) dt | 2 max
serves to normalize the conversion rate of different samples. In fact, since the TGA experiments are carried out at different heating rates, the difference between the conversion rate of different samples could be very important. Therefore, in order to consider the effect of all experiments in the objective function, the conversion rate of each sample is divided by its maximum value [START_REF] Conesa | Thermogravimetric studies on the thermal decomposition of polyethylene[END_REF][START_REF] Várhegyi | Kinetics of charcoal devolatilization[END_REF].

A c c e p t e d M a n u s c r i p t

It is important to note that both forms of the objective functions have been used in the literature to determine the kinetic parameters. Studies such as the ones of [START_REF] Rodríguez | Kinetic models for pyrolysis and combustion of sewage sludge[END_REF][START_REF] Várhegyi | Kinetics of the oxidation of bismuth tellurite, bi2teo5 : Empirical model and least squares evaluation strategies to obtain reliable kinetic information[END_REF][START_REF] Caballero | Kinetic study of the pyrolysis of neoprene[END_REF] and [START_REF] Yan | Analysis of volatile species kinetics during typical medical waste materials pyrolysis using a distributed activation energy model[END_REF] have employed the integral form. As examples of studies that have used the differential form, we mention [START_REF] Hu | Kinetic study of chinese biomass slow pyrolysis : Comparison of different kinetic models[END_REF][START_REF] Grønli | Thermogravimetric analysis and devolatilization kinetics of wood[END_REF][START_REF] Vamvuka | Devolatilization and combustion kinetics of low-rank coal blends from dynamic measurements[END_REF] and [START_REF] Várhegyi | Combustion kinetics of biomass materials in the kinetic regime[END_REF]. The difference between the estimation of kinetic parameters by the integral and differential forms of objective function have been rarely discussed in various publications. A discussion of the difference between these two forms was conducted by [START_REF] Várhegyi | Thermal decomposition of wheat, oat, barley, and brassica carinata straws. a kinetic study[END_REF]. According to Várhegyi et al. [START_REF] Várhegyi | Thermal decomposition of wheat, oat, barley, and brassica carinata straws. a kinetic study[END_REF], it is somewhat difficult to determine the most appropriate form to determine the kinetic parameters. A good choice for the estimation of kinetic parameters of reactions with slow conversion rates would be the integral form. In other cases, the minimization of the differential form is the most sensible to estimate the kinetic parameters choice.

In this paper, we use a hybrid form of the objective function that simultaneously takes into account the minimization of both integral and differential forms [START_REF] Scott | An algorithm for determining the kinetics of devolatilisation of complex solid fuels from thermogravimetric experiments[END_REF] and [START_REF] Sonobe | Kinetic analyses of biomass pyrolysis using the distributed activation energy model[END_REF]. This hybrid form is given by the following function :

O f = (1 -λ)O f I 2 + λO f D 2 ( 21
)
where λ is a number varying between 0 and 1. When λ equals to 0, we get the integral form of objective function. When λ equals to 1, the form of the objective function is the differential one. In order to determine the kinetic parameters, the minimization of the objective function [START_REF] Cai | New distributed activation energy model : Numerical solution and application to pyrolysis kinetics of some types of biomass[END_REF] has been achieved using the pattern search method of Hooke-Jeeves [START_REF] Hooke | direct search" solution of numerical and statistical problems[END_REF][START_REF] Torczon | On the convergence of pattern search algorithms[END_REF] method. It is an iterative method of minimization which does not require the calculation of derivatives of the objective function. This method has been used by Várhegyi et al. [START_REF] Várhegyi | Combustion kinetics of biomass materials in the kinetic regime[END_REF] to determine the kinetic parameters. [START_REF] Várhegyi | Combustion kinetics of biomass materials in the kinetic regime[END_REF] has shown that even tough this method converges slowly toward the optimal parameters, it is an efficient method.

The evaluation of the objective function requires the calculation of the conversion ratio and conversion rate of polyurethane pyrolysis. The calculation of these values was obtained by the numerical integration of the differential equations given by the equations of rate conversion of the different pseudo-components (Eq. 4). The integration is conducted by the fourth-order Runge-Kutta method.

Results and discussion

Figures 1 shows that the decomposition of the different samples starts at 100 • C. DTG curve illustrates clearly that the pyrolysis of polyurethane is occurring in at least two steps which correspond to the successive reactions of thermal decomposition of the polyol resins and the polyisocyanate resin. The first one takes place between 100 and 300 • C. About 50 % of the total decomposed mass of polyurethane is lost in the first step. The second one starts approximately at 300 and extends to 800 • C. Therefore, the polyurethane decomposes in a broader range in the second step. Otherwise, in this step about 50 % of total decomposed polyurethane is lost. It is important to note that the final char residue represents about 26.75 % of initial mass of samples.

Decomposition of polyurethane by two parallel independent reactions

In order to determine the kinetic parameters of polyurethane pyrolysis, we supposed firstly that the decomposition of the polyurethane is done by two parallel independent reactions. In other words, the decomposition of polyurethane in each step will be described by one reaction.

As initial values for the kinetic parameters we have taken : 1 × 10 5 kJ mol -1 for activation energy, 1.0 × 10 10 s -1 for the pre-exponential factor, 0.5 for the coefficient c j and 1 for the reaction order.

8

A c c e p t e d

M a n u s c r i p t

The parameter λ is set to 0.5. The optimized kinetic parameters for the polyurethane pyrolysis are

shown in table 1. The final value of the objective function is 10.81 %. From the optimal coefficient c j , we find that about 0.45 of the polyurethane is decomposed in the first stage while about 0.55 is decomposed during the second one. It is worth noting that as the sum of the c j coefficient has to be equal to 1, the obtained values c j have been normalized to satisfy this condition. These results may be explained by the fact that each stage represents the decomposition of one of the organic resins from which the polyurethane is synthesized. From the results, it can be seen also that the polyurethane pyrolysis has a reaction order close to two in the 2 step and to 3 in the second one.

- Moreover, the figures indicate that the model did not reproduce the shoulders appearing on the first and second stages of the DTG curves. We note that a shoulder represent generally a peak temperature which overlap another neighboring peak temperature [START_REF] Grønli | Thermogravimetric analysis and devolatilization kinetics of wood[END_REF].

Decomposition of polyurethane by four parallel reactions

It is somewhat difficult to describe each step by one reaction. That is why we have supposed that each step occurs by two parallel independent reactions highly overlapping each other. Therefore, the polyurethane pyrolysis is supposed to be done by four parallel independent reactions.

In order to determine the kinetic parameters of these reactions by least square method, we chose equal initial values of the kinetics parameters : 0.25, 1.0 × 10 10 s -1 , 1 × 10 5 kJ mol -1 , 1 for c j , Z j , Ea j and n j respectively. We take λ equal to 0.5. The initial value of the objective function equals 103.47 %. The best value of the objective function obtained by the minimization process is about of 5.56 %. Table 2 shows the estimated kinetic parameters obtained from the pyrolysis of polyurethane. Figure 3 shows the curves of the conversion ratio and conversion rate obtained by TGA analysis and the used kinetic model at the slowest and fastest heating rates : 20 and 80

• C min -1 .
We can observe a good agreement between experimental and calculated result. We note that the optimal objective function decreased from 10.81 % in the case of two reactions model to 5.56 % for the four reactions model. The decrease of the objective function proves that the latter one fits better the experimental data. The use of four reactions model allow to describe successfully the shoulders appearing in the DTG curves for both slowest and fastest heating rates. Furthermore, the four reactions model determined more precisely the position of the peaks along the temperature axis. However, the fit still slightly inaccurate to describe the peak height especially for the fastest heating rate.

Page 10 of 18

Effect of thermal-lag on the kinetic analysis

In the above sections, we have taken the temperature of samples uniform and equal to the TGA furnace temperature. The latter one increases linearly as a function of time at the constant heating rate β (cf. eq.( 10)). In reality, there is a thermal lag (deviation) between the true sample and furnace temperatures. Indeed, they could be equal only on the surface of the sample. Many papers have investigated thermal lag occurring during TGA analysis [START_REF] Antal | Cellulose pyrolysis kinetics : revisited[END_REF][START_REF] Lin | Kinetics and mechanism of cellulose pyrolysis[END_REF][START_REF] Lyon | Thermokinetic model of sample response in nonisothermal analysis[END_REF][START_REF] Vyazovkin | {IC-TAC} kinetics committee recommendations for collecting experimental thermal analysis data for kinetic computations[END_REF][56]. These papers reported that the thermal lag could be affected by many factors such as heating rate, initial size/mass of the sample, thermal effect of the reactions (heat generation/absorption by exothermic/endothermic reactions), the composition and the thermo-physical properties of the sample and the carrier gas of TGA, etc [START_REF] Antal | Cellulose pyrolysis kinetics : revisited[END_REF][START_REF] Lin | Kinetics and mechanism of cellulose pyrolysis[END_REF][START_REF] Lyon | Thermokinetic model of sample response in nonisothermal analysis[END_REF][START_REF] Vyazovkin | {IC-TAC} kinetics committee recommendations for collecting experimental thermal analysis data for kinetic computations[END_REF].

The temperature of the polyurethane samples was measured during the TGA measurements

11

A c c e p t e d M a n u s c r i p t conducted in the current study. Figure 4 shows that there is an important deviation between the measured and furnace temperatures (henceforth noted T m and T f respectively) for the three chosen heating rates 20, 60 and 80 • C min -1 . The maximum thermal lags (T m -T f ) are equal to 67.99, 86.81 and 87.70 • C for these three heating rates respectively. The variation of the sample temperature as a function of time can be re-described by a linear relation with new heating rates which are 21.28, 65.46 and 88.15 • C min -1 . One can explain these higher thermal lags essentially by three factors which are the heating rates, initial mass of samples and the self-heating of the samples. Indeed as explained by many studies [START_REF] Antal | Cellulose pyrolysis kinetics : revisited[END_REF][START_REF] Vyazovkin | {IC-TAC} kinetics committee recommendations for collecting experimental thermal analysis data for kinetic computations[END_REF], the thermal lag can dramatically increase with heating rate. Moreover, the use of an important initial sample mass can lead to the creation of temperature gradient and therefore to the increasing of thermal-lag (particularly at higher heating rates). It was recommended that the heating rate and the initial mass samples should be less than 10 • C min -1 and 10 mg respectively in order to minimize the thermal-lag [START_REF] Lyon | Thermokinetic model of sample response in nonisothermal analysis[END_REF]. The violation of these recommendations in the measurements conducted in the current study (the initial mass of the samples are 16.7 , 16.0 and 16.4 mg and the heating rates are 20, 60 and 80 • C min -1 respectively) could explain the important thermal lag values occurring during TGA analyses. The third factor which can be noted from the positive deviation between the sample and the furnace temperatures is the self-heating. Indeed, as the temperature of the sample increases faster than that of the furnace, one can conclude that there is a heat generation occurring during polyurethane pyrolysis reaction (i.e. exothermic reaction). The self-heating of the sample can lead to a gradient temperature inside the sample (especially with higher initial mass) [START_REF] Lyon | Thermokinetic model of sample response in nonisothermal analysis[END_REF][START_REF] Vyazovkin | {IC-TAC} kinetics committee recommendations for collecting experimental thermal analysis data for kinetic computations[END_REF]. To minimize the effect of the self-heating, one can use lower initial sample mass [START_REF] Vyazovkin | {IC-TAC} kinetics committee recommendations for collecting experimental thermal analysis data for kinetic computations[END_REF]. In order to illustrate the effect of the thermal lag on the kinetic analysis of polyurethane pyrolysis, we will recalculate the kinetic parameters of the two reactions model using the sample temperature obtained at the three heating rates 20, 60 and 80 The recalculation of the kinetic parameters using the sample temperature shows that the accuracy of the estimated parameters was severely affected by the thermal lag. In order to improve the accuracy of the kinetic analysis, new TGA experiments with lower heating rates and initial mass of samples should be done. Furthermore, the development of a more detailed heat transfer model within the sample will allow a more accurate kinetic analysis of polyurethane pyrolysis. To determine the activation energy of each reaction by Kissinger method, T 2 max is plotted as a function of -1 T max for the peak of each step at different heating rates. T max is the temperature corresponding to the maximum rate of decomposition. As explained above, the resulting plot is a straight line whose slope is Ea R . The estimated activation energy for the first and second steps are equal to 66.60 and 210.50 kJ mol -1 , respectively. The value of the square of correlation coefficient (r 2 ) is equal to 0.9667 for the first step and 0.8939 for the second one, which can indicate that the experimental data are well fitted by linear plot for the two steps.

The determination of activation energy by the Friedman method was carried out for two conversion degrees of 0. The presented results show that the obtained kinetics parameters (activation energies and preexponential factors) in the present work are the lowest compared to the the published ones. The difference between the kinetic parameters estimated in this study and those calculated by other workers may probably be attributed the fact that types of studied polyurethane are different.

A c c e p t e d M a n u s c r i p t

The kinetic models used to describe polyurethane pyrolysis may also influence the obtained values kinetic parameters. Other factors include the experimental conditions of thermogravimetric analysis and the calculation procedure used to determine the kinetic parameters can cause the difference between the values of the estimated parameters.

Conclusion

A kinetic analysis for the polyurethane pyrolysis was conducted in this paper. Thermal analysis of the polyurethane was carried out with the thermogravimetry technique which allows a non-isothermal determination of conversion ratio and conversion rate curves as a function of heating rate. Different experiments for polyurethane pyrolysis were conducted at different heating rates. The experimental results showed that the polyurethane pyrolysis is mainly completed in two steps, the first one takes place between 100 and 300 • C, and the second occurs between 300 and extends to 800 • C.

A kinetic model has been presented to describe solid matter pyrolysis during the TGA analysis. This model is based on the parallel independent reactions models. In order to determine the kinetic parameters of the different reactions of polyurethane pyrolysis, a nonlinear least-square procedure has been used. The proposed model has been employed to describe the polyurethane pyrolysis using two and four reactions. The results shows that the use of four reactions model allowed to described more accurately the experimental data, particularly to reproduce the highly overlapped reactions (shoulders). A comparison of the kinetic parameters calculated using the sample and the furnace temperatures showed that thermal-lag can affect considerably the kinetic analysis. In order to improve the reliability of the kinetic analysis, one should minimize the deviation between the sample and furnace temperature (thermal-lag). The kinetics parameters of polyurethane pyrolysis estimated in the paper would be used later in CFD codes to simulate gas emission occurring during casting process.

7 ,

 7 16.0 and 16.4 mg. The analyses of the samples were conducted by a Perkin-Elmer TGA 7 thermobalance. The samples were heated in nitrogen flow at different heating rates : 20, 60 and 80 • C min -1 . The mass loss and the mass loss rate curves of the polyurethane are given in figure 1. The results shows that the average value of the residual mass of polyurethane equals around 26.75 % of initial mass of the various samples.

Figure 1 :

 1 Figure 1: Experimental data of polyurethane pyrolysis obtained by TGA analysis conducted in nitrogen atmosphere and at various heating rates : 20, 60 and 80 • C min -1

) Page 7 of 18 A

 18 c c e p t e d M a n u s c r i p t Where m = -d f (α M ) dα M | T =T max . Taking the logarithm of this equation leads to the linear relationship of Kissinger :

2 .

 2 The figures show that the model can globally follow the shape of the peak temperatures in the different stages of polyurethane decomposition. However, the two reactions model could not determine accurately the position of the peaks along the temperature axis and their height.

Figure 2 :

 2 Figure 2: Experimental and calculated pyrolysis TG and DTG curves at the slowest and fastest heating rates, 20 and 80 • C min -1 , for 2-parallel independent reactions model

Figure 3 :

 3 Figure 3: Comparison between experimental and calculated data at the slowest and fastest heating rates, 20 and 80 • C min -1 , for 4-parallel independent reactions model

β 1 β

 1 = 80.0 • C min -= 60.0 • C min -1 β = 20.0 • C min -1

Figure 4 :

 4 Figure 4: Variation of furnace and measured temperatures with respect to time at various heating rates : 20, 60 and 80 • C min -1

  15 and 0.85 which are chosen to be in the first and second stage respectively. Plotting of ln β dα M dT as a function of 1 T for each conversion degree leads to the estimated values of activation energy of 111.5 and 202.9 kJ mol -1 for the first and second stage respectively. The square of the correlation coefficient of Friedman plot of the first and second stage are equals to 0.9999 and 0.9909 respectively. The results show that the values of the activation energy obtained by the different methods are globally quite different. The values of the activation energies obtained by non-linear fitting are generally lower than those obtained by Kissinger and Friedman methods. However, the value of activation energy obtained for the first stage using the non-linear fitting method (59.80 kJ mol -1 ) A c c e p t e d M a n u s c r i p t

Table 1 :

 1 Kinetic parameters of polyurethane pyrolysis obtained using the Least-Squares Evaluation of three experiments

	Reactions	c j	c norm j	Z j ( s -1 )	Ea j (kJ mol -1 )	n j	Of (%)
	1 2	0.454 0.450 2.45 × 10 4 0.556 0.550 1.20 × 10 6	59.8 107.2	1.905 10.81 2.882

Examples of comparison between the calculated and experimental results of polyurethane pyrolysis at the slowest and fastest heating rates, 20 and 80 • C min -1 , are presented on the figure

Table 2 :

 2 Kinetic parameters of polyurethane pyrolysis obtained using the Least-Squares Evaluation of three experiments

	Reactions	c j	c norm j	Z j ( s -1 )	Ea j (kJ mol -1 )	n j	Of (%)
	1	0.334 0.331 9.95 × 10 4	59.650	4.455
	2 3	0.160 0.159 1.50 × 10 10 0.235 0.233 3.187 × 10 7	118.12 117.09	0.857 2.336	5.57
	4	0.280 0.277 4.26 × 10 8	148.07	2.575
				10

  • C min -1 . Even though we can use the new heating rates obtained by the linear regression of the sample temperature (67.99, 86.81 and 87.70 • C) in the new kinetic analysis of polyurethane pyrolysis, we will use the raw values of the sample temperature measured by the thermocouple during the TGA experiments. The same initial values of the kinetic parameters used in the previous kinetic analysis with two reactions model were used in the current one. The new optimized kinetic parameters were obtained for a objective function of 10.84 %. There parameters are shown in table3. The new results of kinetic analysis of polyurethane pyrolysis illustrate that taking into account the thermal lag leads to different kinetic parameters. The new activation energies of the first and second reactions are greater than that obtained for the furnace temperature ( 60.5 and 110.6 kJ mol -1 instead of 59.8 and 107.2 kJ mol -1 respectively). The exponential factor of the first reaction was decreased from 2.24 × 10 4 to 2.45 × 10 4 s -1 , however that of second reaction was increased from 1.20 × 10 6 to 1.22 × 10 6 s -1 . The orders of first and second reactions were increased to become 2.038 and 3.022 instead of 1.905 and 2.882 respectively. One can also note that the objective function Of of the current kinetic analysis (has a value of 10.84 %) is just slightly higher than the one obtained in the case of use the furnace temperature (equals to 10.81 %).

	M a n u s c r i p t
	A c c e p t e d

Table 3 :

 3 Kinetic parameters of polyurethane pyrolysis calculated using the sample temperature

	Reactions	c j	c norm j	Z j ( s -1 )	Ea j (kJ mol -1 )	n j	Of (%)
	1 2	0.463 0.458 2.24 × 10 4 0.547 0.542 1.22 × 10 6	60.5 110.6	2.038 10.84 3.022
	6.4. Estimation of activation energy by Kissinger's and Friedman's methods
	A comparison was conducted between the kinetic parameters estimated by the non-linear and
	linear fitting methods. The values of activation energy of each step was determined by two linear
	methods, namely Kissinger and Friedman methods. Table 4 presents the energy of activation of
	polyurethane decomposition calculated by the different methods.

Table 4 :

 4 Activation energies of polyurethane decomposition estimated by the different methods that obtained by Kissinger method (66.60 kJ mol -1 ). For the second stage, the values obtained by the Kissinger and Friedman methods are quite similar, 210.5 and 202.9 kJ mol -1 respectively. The difference between the results obtained by the different methods may be attributed to the dissimilar mathematical treatment in each method.6.5. Comparison of current results with published onesMany studies have investigated the pyrolysis of polyurethane. A comparison between of the estimated kinetic parameters published in the literature and those obtained by the current study is given in table5. Different types of polyurethane have been analyzed in this studies. Font et al

	Methods	Ea 1 (kJ mol -1 ) Ea 2 (kJ mol -1 )
	Non-linear	59.80	107.2
	Kissinger	66.60	210.5
	Friedman	115.5	202.9
	is close to

[START_REF] Font | Kinetic law for solids decomposition. application to thermal degradation of heterogeneous materials[END_REF]

studied the pyrolysis of commercial polyurethane obtained by the reaction bewteen the polyadipate of 1,4-butanediol or 1,6-hexanediol with diphenylmethane p,p%-diisocyanate. Rein et al [57] studied the flexible polyurethane foam. Prasad et al [58] studied polyurethane foam formed out of a reaction between toluene diisocynate (TDI)and a polyol. Pau et al [59] studied polyurethane foams made from the reaction of mainly toluene diisocyanate and polyalkoxy polyether polyol.

Table 5 :

 5 Comparison between the kinetic parameters of polyurethane pyrolysis obtained in this work with other published resultsEa j (kJ mol -1 ) Z j ( s -1 ) n j

	Font et al [12]	133.6 190.4	2.55 × 10 12 0.951 9.76 × 10 15 0.668
	Rein et al [57]	124.0 148.0	1.58 × 10 8 2.00 × 10 11 0.210 1.140
	Prasad et al [58]	135.0 175.0	1.69 × 10 8 8.75 × 10 9	1.000 1.160
	Pau et al [59]	179.0 231.0	3.61 × 10 14 9.510 4.72 × 10 16 1.230
	Present Work	59.8 107.2	2.45 × 10 4 1.20 × 10 6	1.905 2.882

Acknowledgement

This study is fully-funded by Montupet SA. The authors would like to express their deep acknowledgment to Montupet. We thank INERIS (Institut National de l'Environnement Industriel et des RISques) for their technical support in TGA experiments. We are grateful to Dr. Vahid Ebrahimian and M. Aghiles Garah for their help in this work. Finally, we thank the reviewers whose thoughtful guidance helped to improve the manuscript.