Francis Billaud
email: francis.billaud@ensic.inpl-nancy.fr

Julien Gornay
email: julien.gornay@ensic.inpl-nancy.fr

Lucie Coniglio
email: lucie.coniglio@ensic.inpl-nancy.fr

Pyrolysis of secondary raw material from used frying oils

Keywords: Used frying oils, Pyrolysis, Biodiesel, Bioenergy

HAL is

Introduction

Increasing global population is depleting the world's supply of fossil fuels. Pyrolysis of used frying oils could be an attractive approach for the production of light hydrocarbons useful as petrochemical raw material, or the production of dihydrogen utilizable in fuel cells, or the production of synthetic gas (mixture of dihydrogen (H 2 ), carbon monoxide (CO), and carbon dioxide (CO 2 )) in favourable ratios to obtain biodiesel by Fischer-Tropsch catalysis (H 2 /CO = 2) or to obtain bio-methanol ((H 2 -CO 2 )/(CO+CO 2 ) = 2) [START_REF] Song | Technologies for direct production of flexible H 2 /CO synthesis gas[END_REF]. In addition, these used frying oils are wastes, then not expensive (only the price of collecting should be taken into account), and are mainly used as biofuel or as feedstocks in cement works. Used frying oils have, however, the advantage to contain valuable chemical species, mainly triglycerides and free fatty acids.

1-8

Literature review

Contrary to the pyrolysis of triglycerides for which many studies were published in the literature, there is few work concerning pyrolysis of used frying oils. [START_REF] Dandik | Converting used oil to fuel and chemical feedstock through the use of a fractionating pyrolysis reactor[END_REF] studied the conversion of used frying oil to hydrocarbon fuels and chemicals, by using a special pyrolysis reactor equipped with fractionating columns. In Dandik and Aksoy work, pyrolysis experiments were conducted at different temperatures (400°C et 420°C) and with different fractionating column lengths (180, 360, and 540 mm). The main products observed by the authors were condensable hydrocarbons (C 5 -C 17 paraffins, olefins, and in fewer amounts, aromatics, cycloparaffins, and cycloolefins) and gases (H 2 , CO, and CO 2 plus C 1 -C 6 paraffins and olefins, consisting mostly of C 1 -C 3 species). The used frying oils tested by the authors were composed by linoleic, oleic, palmitic and stearic acids (mainly as triglycerides). The authors observed an increase in reactant mixture conversion when increasing temperature and decreasing length of columns. Furthermore, an increase in temperature increased olefin fraction in the liquid phase. According to [START_REF] Dandik | Converting used oil to fuel and chemical feedstock through the use of a fractionating pyrolysis reactor[END_REF], judiciously selected operating conditions of pyrolysis could direct the reaction to products that may be used either as fuel or as raw material in chemical industry.

Experimental section

Feedstock composition

The feedstock used (brand name VEGETAMIXOIL ® from Ecogras company, France) is a sample obtained after filtration and decantation of used frying oils (secondary raw material). The main physicochemical properties of VEGETAMIXOIL ® and mass composition are presented in table 1. While the detailed molecular composition of VEGETAMIXOIL ® is unknown, weight percents in carbon (C), oxygen (O), and hydrogen (H) elements are given. As the melting points of VEGETAMIXOIL ® are above room temperature, the feedstock must be heated until 70°C to be injected in the laboratory reactor system.

Overall setup

The pyrolysis pilot plant can be split into five parts (Figure 1): • Injection of the feedstock (secondary raw material from used frying oils) and of the diluents (nitrogen and/or water); compressed air was used to eliminate the coke produced during pyrolysis.

• Preheating of the mixture, feedstock plus diluents.

• Reaction unit with the oven (Thermolyne F 79300 type) and the reactor placed inside.

• Product trapping by condensation and gas/liquid phase separation.

• Liquid and gaseous effluent analysis, and measurement of the coke formed.

• Feedstock and Gas injection. The feedstock was injected by a piston pump (Hamilton Roy). The flow rates of the pump were between 2.3 10 -3 and 1.7 10 -1 L/min. The gases (nitrogen or air) were injected separately by mass flow regulators. Nitrogen (with addition of water at several tests) was used as diluent during pyrolysis reaction and air was used at the end of each experiment to oxidize the coke formed. A gas meter situated after the reaction and condensation units gave the flow rates of the nitrogen and gaseous products formed during cracking (or the flow rates of the air and coke oxidation products, CO and CO 2 ).

• Preheating. The feedstock (or the air after cracking reaction) and the diluent (nitrogen and/or water) were preheated separately in stainless steel tubes surrounded by a roll of heating resistors. The diluent and the feedstock were preheated to 300°C. At these temperatures, the feedstock was prevented from cracking before entering the reactor.

• Oven and Reactor. An electric oven (brand name Thermolyne, type F 79300) was used. The temperature remained steady around the middle of the oven in a portion of 150 mm. The oven had its own regulator system. The reactor is composed of a stainless-steel tube (550 mm long and 14.3 mm in inside diameter). A thermocouple K type (chromel/alumel) placed in the centre of the oven gave the effective temperature inside the isotherm part of the reactor. A pressure controller situated at the outlet of the reactor measure the pressure.

• Trapping of products. At the reactor outlet, a set of flasks maintained at 20°C (cryostat Julabo) condensed all the condensable cracking products. The obtained liquid fraction was weighted and analyzed at the end of the experiment. The non-condensable cracking products (C 1 -C 4 gaseous hydrocarbons, H 2 , CO, and CO 2 ) together with the diluent (nitrogen) were directed to the analysis circuit, and then to the online gas meter. Hence, the total volume of the gaseous products could be measured during the time of experiment (which starts when steady state pyrolysis conditions were achieved). • Product analysis (liquid and gaseous effluents, and coke). The liquid fraction, containing no oxygenated hydrocarbon, was analyzed by gas chromatography. A Schimadzu GC 17A chromatograph was used, equipped with a PONA capillary column (apolar stationary phase = grafted methyl silicone, film thickness = 0.5 mm, length = 50 m, intern diameter = 0.21 mm). Slow temperature programming of the furnace is required to separate all the products (80°C for 20 min, then +2°C/min up to 300°C). The injector (with 1/100 split) and the flame-ionization detector (FID) were at 300°C, and the carrier gas was nitrogen. Primarily, a qualitative analysis by gas-chromatography/mass-spectrometry (GC-MS) coupling was also realized in order to verify the product molecular structure of the sample.

The cracking gaseous fraction contained dihydrogen (determined by catharometer chromatography equipped with a silica-gel filled column, length = 5 m, intern diameter = 6 mm, carrier gas = nitrogen), CO and CO 2 (determined by a Cosma Cristal 300 analyzer with no dispersive infrared ray absorption), oxygenated and non-oxygenated hydrocarbons (determined by FID-detection GC, programming of the furnace: -60°C for 4 min, +10°C/min up to 180°C, 180°C for 10 min, and then +10°C/min up to 300°C).

At the end of each experiment, the amount of coke formed during the reaction was determined by oxidation with air at 850°C.

The analyzed products can be classified as follows: • H 2 • CO and CO 2 • C 1 -C 4 cut (methane, ethylene, ethane, propylene, propane, 1-butene, n-butane, and traces of acetylene) • C 4 + cut (linear 1-olefins and n-paraffins) • Coke

Description of the parametric study

The effect of various parameters on the VEGETAMIXOIL ® pyrolysis were studied:

• Reactor temperature. An increase in temperature leads to high conversions of the reactant mixture and to the production of small molecular species (H 2 , CH 4 , C 2 H 4 …).

• Residence time. An increase in residence time leads to increasing conversions of the reactant mixture and has a significant effect upon the chemical nature and yields of the formed products.

• Diluent nature. Water or equimolar mixture of water and nitrogen were used.

• Reaction initiator. Hydrogen peroxide (H 2 O 2 ) added in small amounts in the reactant mixture should increase the conversion at low temperature (650°C).

• Reaction inhibitor. Thiophene (C 4 H 4 S) added in very small amounts in the reactant mixture (300 ppm), as often made in hydrocarbon pyrolysis process [START_REF] Billaud | n-Decane pyrolysis at high temperature in a flow reactor[END_REF], should reduce catalytic wall effects that might occur at high temperature (800 °C).

During pyrolysis of VEGETAMIXOIL ® , various products were observed. Among them, key products obtained with interesting yields were focused on: (i) dihydrogen for its applications in fuel cells or as one of the components of synthetic gas, (ii) carbon monoxide (CO) and carbon dioxide (CO 2 ) for their applications in biodiesel production by Fischer-Tropsch catalysis (H 2 /CO = 2) or in bio-methanol production ((H 2 -CO 2 )/(CO+CO 2 ) = 2) [START_REF] Song | Technologies for direct production of flexible H 2 /CO synthesis gas[END_REF], (iii) methane (CH 4 ), ethylene (C 2 H 4 ), and propylene (C 3 H 6 ) for their applications as fuels or petrochemicals.

Experimental results

As mentioned previously (Table 1), the detailed molecular composition of VEGETAMIXOIL ® in not given by Ecogras company. Analytical methods required to determine detailed molecular compositions of used frying oils are very complex, expensive, time consuming, and lead finally to partial results because some components of these oils can not be successfully identified, yet [START_REF] Dobarganes | Formation and analysis of high molecular-weight compounds in frying fats and oils[END_REF][START_REF] Graille | Réactions chimiques induites par la friture[END_REF]. Furthermore, composition of the VEGETAMIXOIL ® samples may change from one collecting source to another. Therefore, the VEGETAMIXOIL ® sample used as reactant mixture has been characterized in this work in terms of C, H, and O element weight fractions rather than in terms of chemical component weight fractions.

As the goal of this work is to produce essentially small molecular species (H 2 , light hydrocarbons, mixtures of H 2 and CO), operating conditions of pyrolysis leading to a total or quasi-total conversion of the reactant mixture were adopted. To obtain a total conversion of the feedstock, the reactor was heated at high temperature (above 700°C). The operating conditions of the pyrolysis experiments realized in this work are shown in table 2. 0.9 0.9 1 1.8 1.1 1.1 1 0.9 a T A and P A are respectively ambient temperature and ambient pressure.

Catalytic wall effects

To minimize coke formation due to the chemical composition of the reactor material, pure water and mixture of water and nitrogen were used. With addition of water, consumption of coke takes place according to the reaction (1). Thus, coke formed is rapidly consumed by water, and a steady state is rapidly achieved.

C + H 2 O = CO + H 2 (1)

Effects of temperature and diluent

To study temperature and diluent effects, pyrolysis of VEGETAMIXOIL ® was carried out at two temperatures (700°C and 800°C), with either pure water or equimolar mixture of water and nitrogen as diluent. Other operating parameters (residence time and dilution) are unchanged. Figures 2 and3 show that pyrolysis at high temperature promotes the formation of small species. Thus, at 800°C, CO, CO 2 and H 2 are formed in higher proportions than at 700°C to the detriment of light hydrocarbons (methane, ethylene, and propylene). Water has the advantage to promote light hydrocarbons at the two temperatures investigated here, and to reduce the formation of CO and CO 2 . Formation of H 2 is, however, also reduced in that case (decreasing of 5 % at 800°C).

N 2 + H 2 O H 2 O N 2 +

Effect of residence time

To determine the effect of residence time, an experiment at 700°C was carried out with water as diluent. Residence time was reduced by half by increasing VEGETAMIXOIL ® flow rate. Results are shown on figure 4 in terms of molar fractions of main products. As it can be observed, variations are not significant. This result can be explained by the fact that the conversion of the feedstock is complete, so in that case residence time has poor influence.

Effect of reaction initiator

In an attempt to improve the product yields with lower pyrolysis temperatures, addition of an initiator of reaction was tested. As the feedstock pyrolysis is a gas phase process involving elementary reactions with radical species, hydrogen peroxide (H 2 O 2 ) was selected as initiator agent. Indeed, at the temperature of pyrolysis investigated in this test (650°C), hydrogen peroxide decomposes to form the following species:

H 2 O 2 → 2 HO • (2) H 2 O 2 + HO • → H 2 O + HOO • (3)
Results of the two experiments carried out at 650°C, the first one with pure water and the second one with hydrogen peroxide stabilized in water (concentration of 2.68 mol/l) show that with H 2 O 2 the molar fraction of dihydrogen (H 2 ) is multiplied by two, whereas the yields in light hydrocarbons are divided by four. Nevertheless, the molar fraction of dihydrogen remains low (10.5%). Moreover, yields in CO and CO 2 are also multiplied by two and by four, respectively. The benefit expected with this initiator agent is not positive.

Effect of reaction inhibitor

Pyrolysis at high temperature can induce heterogeneous reactions at reactor wall that may have negative effects on the formation of the desired products. Therefore, as currently encountered in hydrocarbon pyrolysis studies [START_REF] Billaud | n-Decane pyrolysis at high temperature in a flow reactor[END_REF], a reaction inhibitor (thiophene) was added in small quantities (300 ppm) into VEGETAMIXOIL ® . Results of the pyrolysis at 800 °C, with water as diluent, are shown in figure 6.

With regards to the obtained results, thiophene has an important role in VEGETAMIXOIL ® pyrolysis by directing the reaction to the production of the focused products. Indeed, dihydrogen production is in that case divided by four, whereas molar fractions of light hydrocarbons, like methane or ethylene, are doubled. Moreover, thiophene addition leads to a decrease in CO (6.3% instead of 18.8%) and in CO 2 (5.4% instead of 8.1%).

Conclusions

This study points out a promising alternative for upgrading used frying oils after treatment (secondary raw material). It is demonstrated that, under specific conditions of pyrolysis, large amounts of dihydrogen and of light hydrocarbons, like methane and ethylene, can be poduced. The best pyrolysis operating conditions were observed at 800°C with pure water as diluent. Under theses operating conditions, dihydrogen molar fraction exceeds 40%, while light hydrocarbon molar fractions (methane, ethylene, and propylene) exceed 27 %, with very low production of CO and CO 2 . Furthermore, dihydrogen and carbon monoxide were obtained in a ratio H 2 /CO = 2, favourable for low temperature FT process. Using pure water as diluent offers the triple advantage to be as yet an economical feedstock that can easily be supplied and easily separated from diluent and gaseous products by quenching. Also, it has been observed that addition of thiophene directs the pyrolysis reaction to the production of light hydrocarbons (particularly ethylene) to the detriment of dihydrogen and CO (whose production is divided by 3).

Figure 1 .

 1 Figure 1. Schematic description of the laboratory experimental set up.

Figure 4 .

 4 Figure 4. Effect of residence time on molar fractions of VEGETAMIXOIL ® pyrolysis main products (700°C, pure water as diluent).

Figure 6 .

 6 Figure 6. Effect of inhibitor additions (thiophene) on the product molar fractions obtained from VEGETAMIXOIL ® pyrolysis (800°C and H 2 O as diluent).

Table 1 .

 1 Weight composition and physicochemical properties of VEGETAMIXOIL ® .

		Mean value	Uncertainty
			(in the same units as the mean value)
	Weight composition (%)
	Saturated fatty acids	35	± 3
	Mono-unsaturated fatty acids	47	± 3
	Poly-unsaturated fatty acids	15	± 1
	C content	73.6	± 3
	O content	9.7	± 1
	H content	12.2	± 0.5
	Physicochemical properties
	Density (90°C)	0.88	± 0.1
	Viscosity (mm 2 /s) (90°C)	11.1	± 3
	Minimum melting point (°C)	30
	Maximum melting point (°C)	45
	Lower heating value (MJ/kg)	36.4	± 0.8

Table 2 .

 2 Operating conditions of VEGETAMIXOIL® pyrolysis.

	Experiments	1	2	3	4	5	6	7	8
	Temperature of the reactor isothermal part (°C)	800 800	700	700	650	650 7700 800
	Pressure at the reactor inlet (Torr)	881 852	863	1060	807	801 8815 874
	Residence time	τ	τ	τ	τ/2	τ	τ	τ	τ
	Inhibitor (Thiophene) (ppm)	-	-	-	-	-	-	-	Yes
	Initiator (H 2 O 2 )	-	-	-	-	-	Yes -	-
	Inlet volumic flows
	N 2 (l/min) (TPN) H 2 O (ml/min) (T A , P A ) a VEGETAMIXOIL ® flow (ml/min) (T A , P A ) a	0.5 0.4	-0.8	0.45 0.36	-1.45	-0.93 0.93 00.88 0.9 ---

  Effect temperature and nature of diluent in the production of methane (CH 4 ), ethylene (C 2 H 4 ), and propylene (C 3 H 6 ) during VEGETAMIXOIL ® pyrolysis. Effect of temperature and nature of diluent in the production of dihydrogen (H 2 ), carbon monoxide (CO), and carbon dioxide (CO 2 ) during VEGETAMIXOIL ® pyrolysis.

		40	H 2 CO
	Molar Fraction (%)	20 30	CO 2
		10
		0
	H 2 O CH 4 800°C C 2 H 4 C 3 H 6 Diluent and Pyrolysis temperature (°C) 700°C H 2 O N 2 + H 2 O 800°C 700°C Diluent and Pyrolysis temperature (°C) H 2 O H 2 O Figure 2. N 2 + H 2 O 0 5 10 15 20 25 30 Molar Fraction (%) Figure 3.

Acknowledgments

Ecogras S.A. is gratefully acknowledged for the supplying of VEGETAMIXOIL ® .


M Filippi

E Carrera

Flutter analysis of fixed and rotary wings through a one-dimensional unified formulation

Keywords: Rotating beams, Finite element method, Higher-order theories, Composite structures

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

According to Collar's definition, aeroelasticity is ''the study of the mutual interaction that takes place within the triangle of the inertial, elastic, and aerodynamic forces acting on structural members exposed to an air stream, and the influence of this study on design''. The fluid structure interaction (FSI) has brought to catastrophic events due to sudden failures of bridges, airplanes, helicopter blades, etc. A correct and safe design must require, therefore, an accurate prediction of aeroelastic phenomena. Unfortunately, FSI analyses are, in most cases, too computationally expensive; thus a tradeoff between accuracy and cost is often necessary.

For fixed wing configurations, the simplest aerodynamic theories based on strip approaches were proposed in the first half of the twentieth century [1 5]. These models were properly com bined with simplified structural theories to develop reliable aeroe lastic tools. For instance, aeroelastic analyses were performed on wings [6 8] and civil structures [START_REF] Banerjee | A simplified method for the free vibration and flutter analysis of bridge decks[END_REF] by using the beam plate approach. On the basis of this methodology, the flexural, torsional and secondary stiffness of the cross section are reduced to the equivalent beam quantities (EI, GJ, K, etc.). To overcome the intrin sic limitations of this technique (no coupling due to the material distribution), Librescu et al. proposed a refined 1D theory including non classical effects, such as warping and transverse shear defor mations and the non uniform torsion [10 12]. Moreover, Patil et al. developed a non linear 1D formulation for laminated box beams by adopting an asymptotic technique to compute the cross sectional mechanical properties [START_REF] Cesnik | Aeroelastic analysis of composite wings[END_REF][START_REF] Patil | Aeroelastic tailoring of composite box beams[END_REF]. Later, Palacios et al. used the same formulation to evaluate the aeroelastic behav ior of highly flexible wings [START_REF] Palacios | An intrisic description of the nonlinear aeroelasticity of very flexible wings[END_REF]. Nowadays, the rapid increase of the computational power has motivated the development in the area of computational fluid dynamics (CFD). Over the last 30 years, the CFD based aeroelasticity progressed from full potential theo ries (strip theory, panel methods) to problems governed by the Navier Stokes equations. Recent works offer interesting overviews of CFD aeroelastic tools [START_REF] Bennet | An overview of recent developments in computational aeroelasticity[END_REF][START_REF] Huttsell | Evaluation of computational aeroelasticity codes for loads and flutter[END_REF] furthermore providing useful information about the emerging trends in the aeroelasticity field [START_REF] Livne | Future of airplane aeroelasticity[END_REF]. From a structural point of view, the CFD simulations can be coupled with finite element models, in which non linear elements may be used to include either geometrical effects or large deforma tions. A detailed paper about non linear structural models can be found in [START_REF] Demasi | Dynamic aeroelsticity of structurally nonlinear configurations using linear modally reduced aerodynamic generalized forces[END_REF].

Regarding the aeroelasticity of rotary wings, the problem becomes more complex due to the geometric nonlinearities that must be taken into account in both elastic and aerodynamic terms. This need has driven the development of suitable structural and aerodynamic theories as confirmed by the considerable number of articles available in the literature. Almost 50 years ago, Loewy provided a thorough overview on a broad range of topics related to the dynamics and aeroelasticity of rotary wings such as flap lag flutter, pitch lag flutter and ground and air resonance [START_REF] Loewy | Review of rotary-wing v/stol dynamic and aeroelastic problems[END_REF]. Although being more limited in scope, other contributions to the study of classical flutter conditions and unsteady aerody namic theories were proposed in [START_REF] Ham | Helicopter blade flutter. Agard report[END_REF][START_REF] Dat | Aeroelasticity of rotary wing aircraft[END_REF], respectively. In that per iod, the aerodynamic models for the instability predictions were essentially extensions of steady and quasi steady strip theories conceived for fixed wing configurations. The first important unsteady approach for hover, based on Theodorsen's model, was proposed by Loewy [START_REF] Loewy | A two-dimensional approach to the unsteady aerodynamics of rotary wings[END_REF]. Though in an approximate manner, this theory takes into account the effect of the spiral returning wake beneath the rotor. Other valuable models were proposed on the basis of Greenberg's theory [START_REF] Greenberg | Airfoil in sinusoidal motion in a pulsating stream[END_REF], where a pulsating velocity varia tion and a constant pitch angle were included in Theodorsen's model. These modifications were extended to the case of rotary wings in [START_REF] Friedmann | Effect of modified aerodynamic strip theories on rotor blade aeroelastic stability[END_REF], in which Theodorsen's, Loewy's and Poisso's theories were used for studying the flap lag torsional coupling. The authors pointed out that the above modifications should be included to realistically reproduce the aeroelasticity of rotary wings, where the assumptions commonly used in deriving strip theories (1 cross sections are assumed to perform only simple harmonic pitching and plunging oscillations about a zero equilibrium posi tion; 2 the velocity of oncoming airflow is constant; 3 usual potential small disturbance unsteady aerodynamics are assumed to apply) are intrinsically violated. The modified strip theory with the quasi steady approximation was adopted to study the stability of composite helicopter blades with swept tips both in hover and in forward flight [START_REF] Yuan | Aeroelasticity and structural optimization of composite helicopter rotor blades with swept tips[END_REF]. The swept tip blades were modeled using non linear 1D finite element with the inclusion of transverse shear deformations and out of plane warping. These results were recently used to verify the accuracy of a nonlinear structural for mulation valid for slender, homogeneous and twisted blades [START_REF] Piccione | A finite element model for structural and aeroelastic analysis of rotary wings with advanced blade geometry[END_REF]. Other analyses on advanced geometry blades were carried out in [START_REF] Bir | Aeromechanical stability of rotorcraft with advanced geometry blades[END_REF] with the purpose of evaluating the effects of sweep and droop on both rotor aeroelastic stability and rotorcraft aero mechanical stability. The evolution of the mechanics of helicopter blades, focusing on the aeroelastic and aerodynamic issues, was exten sively discussed in [ 29 32]. Moreover, even if the aeroelasticity of large wind turbine blades is inherently different, several related aeroelastic problems were solved using rotary wing theories [33 35]. Furthermore, the response problems of an isolated wind turbine blade and a complete rotor/tower configuration were detailed discussed in [36 38].

Within this work, we propose an advanced 1D formulation to predict the flutter conditions of various fixed and rotary wing con figurations. The higher order beam theories are obtained using the Carrera Unified Formulation (CUF), which enables to derive, at least theoretically, an infinite number of models [39 41]. The Taylor type expansions (referred to as TE) represent a particular class of the 1D CUF theories. The TE models have been adopted for the study of mechanical behaviors of thin walled structures [START_REF] Carrera | Performance of cuf approach to analyze the structural behavior of slender bodies[END_REF][START_REF] Carrera | Free vibration analysis of thin-walled cylinders reinforced with longitudinal and transversal stiffeners[END_REF] and non conventional cross sections made of isotropic, composite [START_REF] Carrera | Laminated beam analysis by polynomial, trigonometric, exponential and zig-zag theories[END_REF][START_REF] Carrera | Free vibration analysis of laminated beam by polynomial, trigonometric, exponential and zig-zag theories[END_REF] and functionally graded [START_REF] Daoud | Free vibration of FGM layered beams by various theories and finite elements[END_REF] materials. Other encouraging results have also been obtained in the dynamics of rotors by analyzing both blades [START_REF] Carrera | Free vibration analysis of rotating composite blades via carrera unified formulation[END_REF] and spinning shafts [ 45 47]. Furthermore, aeroelastic analyzes were performed on fixed wing configurations by combining the TE elements with aerodynamic panel methods [START_REF] Varello | Static aeroelastic response of wing-structures accounting for in-plane cross-section deformation[END_REF][START_REF] Carrera | A refined structural model for static aeroelastic response and divergence of metallic and composite wings[END_REF] and unsteady strip theories [START_REF] Filippi | Aerodynamic and mechanical hierarchical aeroelastic analysis of composite wings[END_REF].

The structural model: Carrera Unified Formulation

CUF states that the displacement field, u T ðx; y; z; tÞ ðu x ; u y ; u z Þ T , is an expansion of generic functions, F s ðx; zÞ for the displacement vector, u s ðyÞ: uðx; y; z; tÞ F s ðx; zÞu s ðy; tÞ s 1; 2; . . . ; T ð1Þ where T is the number of terms in the expansion and, according to Einstein's generalized notation, s indicates summation. In this work, Eq. (1) consists of Taylor like expansions denoted as 'TE'. For example, the third order displacement field (TE3) is:

u x u x 1 þ x u x 2 þ z u x 3 þ x 2 u x 4 þ xz u x 5 þ z 2 u x 6 þ x 3 u x 7 þ x 2 z u x 8 þ xz 2 u x 9 þ z 3 u x 10 u y u y 1 þ x u y 2 þ z u y 3 þ x 2 u y 4 þ xz u y 5 þ z 2 u y 6 þ x 3 u y 7 þ x 2 z u y 8 þ xz 2 u y 9 þ z 3 u y 10 u z u z 1 þ x u z 2 þ z u z 3 þ x 2 u z 4 þ xz u z 5 þ z 2 u z 6 þ x 3 u z 7 þ x 2 z u z 8 þ xz 2 u z 9 þ z 3 u z 10
The remarkable feature of these models is that classical beam theories can be obtained as particular cases of Taylor expansions. Although, the capabilities of the TE approaches have been evaluated in several works, we suggest referring to [START_REF] Carrera | Beam structures. Classical and advanced theories[END_REF] for a detailed description. The stresses and the strains are grouped as follows:

p f zz xx xz g T r p fr zz r xx r xz g T n f zy xy yy g T r n fr zy r xy r yy g T ð2Þ
where the subscript ''p'' stands for the terms lying on the cross section, while ''n'' stands for the terms lying on the other planes, which are orthogonal to the cross section. The strain displacement relations and Hooke's law are, respectively:

p D p u n D n u ðD np þ D ny Þu ð3Þ r p Cpp p þ Cpn n r n Cnp p þ Cnn n ð4Þ
where D p and D n are linear differential operators. Both laminated box beams and cylinders can be considered constituted by a certain number of either straight or curved plates of orthotropic material, whose material coordinate systems generally do not coincide with the physical coordinate system (x, y, z) of the structure. The matrices of material coefficients of the generic material k are

For the sake of the brevity, the explicit forms of the coefficients of the C matrices are shown in [START_REF] Carrera | Variable kinematic one-dimensional finite elements for the analysis of rotors made of composite materials[END_REF]. The finite element method is used in order to consider arbitrary cross section profiles. The generalized displacement vector is u s ðy; tÞ N i ðyÞq si ðtÞ ð 6Þ

where N i ðyÞ are the shape functions and q si ðtÞ is the nodal displace ment vector q si ðtÞ q ux si q uy si q uz si n o T ð7Þ

2-D steady and unsteady aerodynamic theories for fixed-wing

Despite steady and unsteady strip theories were already com bined with the 1D CUF elements in [START_REF] Filippi | Aerodynamic and mechanical hierarchical aeroelastic analysis of composite wings[END_REF], for the sake of complete ness, this section provides a brief description of these aerodynamic approaches. According to the assumptions of the small disturbance theory, Theodorsen considered a flat plate with a control surface, which moves in vertical translation h(t) and rotate about an axis at x ba through an angle a(t). As derived in [START_REF] Bisplinghoff | Aeroelasticity[END_REF], the final expres sion of the lift for unit of span is

L pq a b 2 ½ € h þ U _ a ba€ a þ 2pq a UbCðkÞ _ h þ Ua þ b 1 2 a _ a ! ð8Þ
where 'U' is the free stream velocity, 'b' the semichord, 'q a ' the air density and 'CðkÞ' the deficiency lift function, which depends on the so called reduced frequency k xU b . The term 'a' defines the posi tion of the rotation axis with respect to the center of the section and it is dependent on the support condition, the lamination scheme and the applied load. Due to the difficulty involved in defining its correct value for swept and asymmetric laminated structures, the elastic offset 'a' is assumed to be null. Theodorsen identified the term 'CðkÞ' as a Hankel function of the second kind, which is, in turn, comprised of Bessel functions of the first and second kinds [START_REF] Bisplinghoff | Aeroelasticity[END_REF][START_REF] Rodden | Theoretical and computational aeroelasticity[END_REF] CðkÞ FðkÞ þ iGðkÞ H

1 ðkÞ H ð2Þ 1 ðkÞ þ iH ð2Þ 0 ðkÞ ð2Þ
The first contribution of Eq. ( 8), contains the virtual mass (or non circulatory terms) that acts as a form of inertia of the fluid sur rounding the airfoil and, since the single and doubly differentiated terms relative to the mass properties of the structure are small, they are usually neglected. The remaining terms determine the value of the lift and represent the circulatory part. The steady state simpli fication can be easily obtained by replacing the function CðkÞ with the unitary value and omitting the term containing _ a.

Lu2pq a Ub½ _ h þ Ua ð 10Þ
Diederich's approximation is used to correct the sectional lift coef ficient 'C L ' relating to the influence of the aspect ratio ' AR 2Lw cm ' and the sweep angle ('K')

c la dC L da p AR p AR þ c la0 cosðKÞ c la0 cosðKÞ ð 11Þ
where 'L w ' is the length of the wing, 'c' is the mean chord and, 'c la0 ' is the lift curve slope equal to 2p. Finally, since the quantity 'bp' may be approximated with the integral R b b b x bþx dx, which reproduces the distribution of pressure over a slightly inclined, thin and uncam bered airfoil, Eqs. ( 8) and [START_REF] Librescu | Aeroelastic divergence of swept-forward composite wings including warping restraint effect[END_REF] become

L 2p AR cosðKÞ p AR þ 2p cosðKÞ Z b b  b x b þ x dx q a UCðkÞ _ h þ Ua þ b 1 2 a _ a ! ð12Þ
L u 2p AR cosðKÞ p AR þ 2p cosðKÞ Z b b b x b þ x dx q a U½ _ h þ Ua ð 13Þ

Unsteady aerodynamic theory for blade

According to [START_REF] Bielawa | Rotary wing structural dynamics and aeroelasticity[END_REF], the general expression for the lift distribution on a rotary wing yields dL dy

pq a b 2 X 2 € h X 2 þ _ a X ba € a X " # þ 2pq a bX 2 CðkÞ _ h X þ ya þ b 1 2 a _ a X " # ð14Þ
where the reduced frequency 'k' is 2xc 3XR and, 'R' is rotor radius. Although this formulation is limited as it relies on Theodorsen prob lem, it can be easily generalized since the unsteadiness is described by the lift deficiency function.

In the general case, the rotor blade develops oscillating loads, which generate a radial distribution of vorticity blown downward at velocity 'u'. The wake can be described via models with increas ing levels of complexity but, historically, the first successful attempt to incorporate the ''returning wake'' into the unsteady aerodynamics problem was proposed by Loewy [START_REF] Loewy | A two-dimensional approach to the unsteady aerodynamics of rotary wings[END_REF]. According to Loewy's theory, the returning wake is modeled as an infinite number of layers, which are equally spaced by a fixed distance, 'h l ' semi chords. Here, the factor 'h l ' quantifies the contribution of the underlying vortices to the lift distribution on the reference blade. The main difference with respect to Theodorsen's model lies in the form of the lift deficiency function

C 0 ðk; m; h l Þ H ð2Þ
1 ðkÞ þ 2J 1 ðkÞ Wðkm; h l Þ H ð2Þ
1 ðkÞ þ iH ð2Þ
0 ðkÞ þ 2½J 1 ðkÞ þ iJ 0 ðkÞWðkm; h l Þ ð15Þ
where 'J n (k)' and 'H (2) n (k)' are the first kind Bessel and second kind Hankel functions of order n, respectively. The wake weighting func tion 'W(k,m,h l )' for a rotor with 'Q' blades is defined as Wðk; m; hÞ

1 þ P Q 1 q 1 ðe kh l Q e i2p m Þ Q q e iUq e khQ e i2p m 1 ð16Þ
where 'm' and 'U q ' are the ratio of oscillatory frequency to rotational frequency ( x=X) and the phase angle between the qth blade and the reference blade. Due to periodicity, Eq. ( 16) can be expressed following the formulation for a single blade rotor with the modified values ' m' and ' ĥ'

Wðk; m; ĥÞ 1

e k ĥQ e i2p m 1 m x QX ĥ 2pk 0 Q ðb=rÞ ð17Þ
where 'k 0 ' C T =2 p ' (where C T is the thrust coefficient) is the inflow ratio and, 'r' is the generic radial distance.

Equations of motion

The general form of Hamilton's Principle states that the varia tion of the kinetic and potential energy and the variation of the work exerted by non conservative forces must equal zero during any time interval. It may be expressed as

Z t 1 t 0 d T ðU þ U r0 Þ ð Þ dt þ Z t 1 t 0 dL ext dt 0 ð18Þ
where d refers to kinematically admissible perturbations during the indicated time interval and T: total kinetic energy of the system U + U r0 : total potential energy of the system, including both strain energy ('U') and potential of any conservative external forces ('U r0 ') L ext : work done by non conservative forces acting on the system, including any arbitrary external loads

The application of this variational concept leads directly to the equations of motion (EoM) for any given system. The energy expressions for a rotary wing are T 1 2

Z V qð _ u T _ u þ 2u T X T _ u þ u T X T Xu þ 2 _ u T Xr þ 2u T X T XrÞdV ð19Þ
U 1 2 Z V ðu T D T CDuÞdV ð20Þ
U r 0 Z V ðr 0 T nl ÞdV r 0 X 2 q r h R þ 1 2 R 2 r h r 1 2 r 2 ! ð21Þ
where 'q' is the material density, 'r' is the distance from the neutral axis, 'r h ' is the hub radius and the constant rotational speed matrix is

X 0 X 0 X 0 0 0 0 0 2 6 4 3 7 5
The external work generated by the lift can be written as

L ext Z y Z x u z ðx; y; z top Þ Lðx; y; z top Þ dxdy ð22Þ
where 'z top ' is the upper z coordinate of the cross section. According to Eqs. ( 8) and ( 14), Eq. ( 22) for fixed and rotary wings becomes, respectively Introducing CUF (Eq. ( 1)) and FE (Eq. ( 6)) approximations and substituting Eqs. ( 19) [START_REF] Dat | Aeroelasticity of rotary wing aircraft[END_REF] in Eq. ( 18), we obtain

L ext cost Z y Z x u T b x b þ x s I L _ u þ u T b x b þ x s I L u ;x 2 4 þ u T b 1 2 a b x b þ x s I L _ u ;x 3 5 dx dy L ext cost Z y Z x u T X y b x b þ x s I L _ u þ u T X 2 y 2 b x b þ x s I L u ;x dx dy 2 4 3 5 þ cost Z y Z x u T X y b x b þ x s I L b 1 2 a _ u ;x dx dy
Z t 1 t 0 dq T si M ijss € q sj þ dq T si D ijss TOT _ q sj þ dq T si K ijss TOT q sj þ dq T si F is X r dt 0 ð23Þ
where the matrices are written in terms of CUF fundamental nuclei.

For fixed wings, we have

D ijss TOT D ijss L þ D ijss Lm K ijss TOT K ijss þ K ijss L
whereas, for rotary wings

D ijss TOT D ijss LX þ D ijss LmX þ D ijss X K ijss TOT K ijss þ K ijss LX þ K ijss X þ K ijss r0
The subscripts 'L' and 'X' denote the contributions due to the aero dynamic load and rotation, respectively. The compact expressions of the fundamental nuclei yield ; I i; y j l ; I i; y j; y l Z l N i N j ; N i N j; y ; N i; y N j ; N i; y N j; y dy I i j ly ; I i j l y 2 I i; y j; y lr 0 Z l y N i N j ; y 2 N i N j ; r 0 N i; y N j; y dy

M ijss I ij l / ðF s qIF s Þ. D ijss X I ij l / ðF s qIF s Þ . 2X D ijss L cost I ij l / b x b þ x s F s ðx; z top ÞI L F s ðx; z top Þ 0 @ 1 A . x D ijss Lm cost I ij l b 1 2 a / b x b þ x s F s;x ðx; z top ÞI L F s ðx; z top Þ 0 @ 1 A . x D ijss LX cost X I ij ly / b x b þ x s F s ðx; z top ÞI L F s ðx; z top Þ 0 @ 1 A . x D ijss LmX cost X I ij ly b 1 2 a / b x b þ x s F s;x ðx; z top ÞI L F s ðx; z top Þ 0 @ 1 A . x K ijss X I ij l / ðF s qIF s Þ . X T X K ijss r0 I i;yj;y lr0 / ðF s qIF s Þ . X T X K ijss I ij l / D T np ðF s IÞ Ck np D p ðF s IÞ þ Ck nn D np ðF s IÞ h i þ D T p ðF s IÞ Ck pp D p ðF s IÞ þ Ck pn D np ðF s IÞ h i . þ I ij;y l / D T np ðF s IÞ þ D T p ðF s IÞ Ck pn h i F s . þI Ay þ I i;yj l I T Ay / F s Ck np D p ðF s IÞ þ Ck nn D np ðF s IÞ h i . þ I i;yj;y l I T Ay I Ay / F s Ck nn F s . K ijss L cost I ij l / b x b þ x s F s;x ðx; z top ÞI L F s ðx; z top Þ 0 @ 1 A . x ð24Þ
K ijss LX cost X 2 I ij l y 2 / b x b þ x s F s;x ðx; z top ÞI L F s ðx; z top Þ 0 @ 1 A . x
To obtain the natural frequencies and the normal modes of the rotor, we solve the homogeneous equation:

M € q þ D TOT _ q þ K TOT q 0 ð25Þ
Assuming a periodic solution q qe ixt ð26Þ

substituting Eq. ( 26) and its derivative in Eq. ( 25), we obtain

qe ixt K TOT þ G TOT ix M x 2 Â Ã 0 ð27Þ
The quadratic eigenvalues problem (QEP) of generic order R in Eq. ( 27) is written in a classical linear system of order 2ÂR:

M € q þ D TOT _ q þ K TOT q 0 _ q þ _ q 0 & ð28Þ
Now, introducing:

a q _ q & ' _ a _ q € q & ' ð29Þ
the equations of motion assume the form: R T

1 ix I 0 ð30Þ
where:

T 1 R ðK TOT Þ 1 ðD TOT Þ ðK TOT Þ 1 ðMÞ I 0 " # ð31Þ
Since the aerodynamic matrices depend on the reduced frequency 'k', the solutions are found through an iterative procedure [START_REF] Nguyen | Flight dynamics of flexible aircraft with aeroelastic and inertial force interactions[END_REF].

The frequencies are introduced by trial and error for any speed in the considered range until Eq. ( 30) is obtained. The flutter velocity is determined when the real part of one eigenvalue is null.

Numerical results

The illustrated approach enables to study the dynamics of wings and blades subjected to aerodynamic loads. Therefore, the following section presents results related to vibrational and flutter analysis performed on fixed and rotary configurations made of iso tropic and composite materials.

Rectangular cross section beams

We firstly consider three beams with rectangular cross sections with the geometrical and mechanical properties listed in Table 1.

Cases 1 and 2 are isotropic beams with different material prop erties (aluminum and foam), whereas the case 3 is a laminated plate with six layers of the equal thickness. Fig. 1 shows the varia tions of the natural frequencies with the rotational speed for all three cases.

In Fig. 1(a) and (b), the results are compared with the 3D finite element solutions reported in [START_REF] Truong | Structural dynamics modeling of rectangular rotor blades[END_REF]. For the metallic beam (case 1), the TE3 results perfectly match the reference solution, whereas, for the foam beam (case 2), the order of the theory has been increased for describing the significant 3D effects due to the soft material. Fig. 1(c) and (d) show the frequency variations for the case 3 for two symmetric stacking sequences [45 2 /0] S and [45/ 45/0] S , respectively. It should be noted that, for the first lamination scheme, the TE3 theory predicts higher values of the 5th and 7th frequencies with respect to the TE4 solution. The corresponding modal shapes are dominated by torsional deformations. Conversely, regarding the [45/ 45/0] S scheme, the two Taylor like expansions provide almost equivalent results for all modal shapes. 27 (1.11%) : results not provided.

(%) : percentage difference with respect to the present results ðRef : À Pres:Þ Ref : . Fig. 2. Flutter velocity (X f ) and frequency (f f ) as functions of the inflow ratio. '-Ã-': Case 1; '-Ã-': Case 3 [0/0/90] S . TE3.

Table 3

Flutter velocities [rpm] and frequencies [Hz] for the case 3 for various laminations. k0=0.05. The flutter analyses have been performed to determine the instability conditions of the three fixed wing configurations. Table 2 lists the flutter velocities computed with the unsteady Theodorsen theory within the CUF framework. For cases 1 and 2, the results are compared with those obtained using the DLM com bined with the 1D CUF elements [START_REF] Petrolo | Flutter analysis of composite lifting surfaces by the 1d carrera unified formulation and the doublet lattice method[END_REF]. For case 3, further reference solutions have been taken from [START_REF] Hollowell | Aeroelastic flutter and divergence of stiffness coupled, graphite/epoxy cantilevered plates[END_REF], where Hollowell et al. studied the aeroelastic behavior of orthotropic plates through analytical and experimental simulations.

Case Lamination TE3 TE4 X f f f X F f F 3 ½0 2 =90
The comparisons reveal that the unsteady strip theory predicts higher flutter velocities than those of DLM for the first two cases with a maximum percentage difference equal to about 7%. For case 3, the maximum discrepancy increases up to about 10% for the lamination scheme [45/ 45/0] S . Except for this stacking sequence, it should be noted that the CUF results are close to the experimen tal data.

In addition, further flutter analyses have been performed with the contributions of the rotational speed. As pointed out in Section 3.1, the distribution of the aerodynamic loads depends on the value of the spacing factor, ' ĥ'. Therefore, the first analysis has been devoted to evaluating the variations of the instability condi tions by varying the inflow ratio, 'k 0 ', for the metallic and lami nated blade (with Q 4). The results are shown in Fig. 2.

The x axis reports the logarithmic value of the inflow ratios starting from 0.0242, which corresponds to ĥ 1:0 for the metallic structure. For the case 1, it is interesting to note that the flutter speed suddenly increases, whereas the frequencies decreases, when k 0 overcomes the values 0.03 ( ĥ ' 1:238) and 0.07 ( ĥ ' 2:89). Before and after these thresholds, the curves are nearly flat. Conversely, the flutter conditions for the case 2 show a quasi linear dependence with the inflow ratio before the reaching of the asymptotic values. For both structures, for high values of k 0 , the Loewy approach coincides with the Theodorsen theory. Table 3 reports the instability conditions for the different lamination sequences of the laminated plate assuming the inflow ratio k 0 0:05.

Rotating multi section blade: NACA 0012

The airfoil of the blade is the NACA 0012, and the chord and length are 0.3454 and 5.186 m, respectively. Fig. 3 shows the blade profile that is constituted by different sections whose material properties are listed in Table 4. The structure is cantilevered and it is modeled with 7 4 node beams. Several problem parameters are evaluated including the lamination angle of the skin ply, the aerodynamic theories and the order of TE models. Firstly, for the ply angle h equal to zero, Fig. 4 shows the natural frequencies as a function of the rotational speed by adopting differ ent displacement theories.

The comparison reveals that the three expansions predict simi lar results for the mode shapes dominated by the bending whereas, slight discrepancies are evident for the torsional deformations (4th and 6th).

Secondly, the flutter analyses have been carried out by using the aerodynamic model proposed in Section 3.1 with different expres sions of the lift deficiency function. Assuming that the four blade rotor is operating at a thrust coefficient C T 0:014, the corre sponding values of inflow ratio 'k 0 ' and spacing factors ' ĥ' are 0.084 and 3. Table 5 shows the flutter conditions as functions of the lamination angle. The used displacement theories are TE3 and TE4 expansions.

It is observed that Loewy and Theodorsen theories lead to dif ferent results when the angle is 15 and 30. These discrepancies are due to the possibility to have a negative pitch damping, adopt ing the Loewy approach, which may imply the single degree of freedom flutter in blade torsion. In fact, from Fig. 5, where the damping coefficients are plotted as functions of m for h 15, it is noted that the pitch damping reaches high neg ative values when m is approaching 0. Conversely, the coefficient related to the flap motion is positive for all values of m.

Moreover, it is observed significant differences between the flutter velocities obtained through the TE3 and TE4 theories due to the bending torsion coupling that affects the mode shapes. The percentage differences between the two sets of results vary from 2% to 16% depending on the lamination sequence and the considered aerodynamic theory.

Conclusion

This work aims to present a simple tool for the aeroelastic anal ysis of rotating and fixed wings. Unsteady aerodynamic theories are coupled with refined beam theories obtained through Carrera Unified Formulation. Numerical simulations have been performed considering rectangular and multi section profiles constituted by isotropic and composite materials. The effects of the aerodynamic theories, the lamination sequences and the order of structural models on the flutter conditions have been evaluated. The results have been compared, when possible, with theoretical, experimen tal and numerical solutions. To conclude, the proposed formulation enables to study the dynamics of several rotating blade configura tions with a significant accuracy and a lower number of degrees of freedom than the 3D solution; to predict, if combined with the unsteady aerodynamic the ory, the flutter conditions of isotropic and laminated fixed wings; to predict the flutter instability of blades in hover with sim plified and complex shaped profiles. Within this work, Carrera Unified Formulation has shown excel lent performances in the study of the dynamics and the aeroelas ticity of rotors, therefore paving the way for future research that will include: extension to nonlinear structural analysis (geometrical as well as material nonlinearities); analysis of complex shaped rotating structures (swept tip blades, twisted blades etc.); extension to more sophisticated aerodynamic theories (panel methods, CFD)

Fig. 1 .

 1 Fig. 1. 8 frequencies vs rotational speed for cases 1, 2 and 3.

Fig. 4 .

 4 Fig. 4. Frequency vs rotational speed for the multi-section NACA 0012 blade. h 0.

Table 1

 1 Dimensions [mm] and material properties [GPa], [kg/m 3 ].

	Case	c	t	L w	r h	E 11	E 22	E 33	G 12	G 31	G 23	m 12	m 31	m 23	q
	1	25.4	1.6002	444.5	63.5	68.9	-	-	-	-	-	0.3	-	-	2713
	2	86.36	21.59	863.6	0.0	0.0827	-	-	-	-	-	0.202	-	-	190.99
	3	76.2	0.804 a	305.0	0.0	98.0	7.9	7.9	5.6	5.6	5.6	0.28	0.28	0.28	1520

a total thickness of the six-layer beam.

Table 2

 2 Flutter velocities for the three configurations of Table1.

	Case	Theory	Lamination	Present	DLM + 1D-CUF	Anl. [58]	EXP. [58]
	1	TE3	Isotropic	134.9	131.5 ( 2.5%)	-	-
	2	TE5	Isotropic	64.4	60.2 ( 6.98%)	-	-
	3	TE4	½0 2 =90 S	23.5	23.2 ( 1.29%)	21.0 ( 11.9%)	25 (6.0%)
		TE4	½45= À 45=0 S	44.2	40.3 ( 9.68%)	39.0 ( 13.3%)	>32
		TE4	½45 2 =0 S	27.8	26.2 ( 6.11%)	27.8 ( 0.0%)	28 (0.71%)
		TE4	½30 2 =0 S	26.7	25.9 ( 3.1%)	27.8 (3.96%)

Table 4

 4 Material properties of airfoil NACA 0012. Elastic moduli in [GPa] and densities in [kg/m 3 ].

	Material	E 11	E 11	E 11	G 12	G 31	G 23	m 12	m 31	m 23	q
	1	139.0	11.0	11.0	6.05	6.05	3.78	0.313	0.313	0.313	1578
	2	35.756	15.427	15.427	3.902	3.902	3.902	0.27	0.27	0.27	1826.9
	3	1.446	-	-	-	-	-	0.49	-	-	13683.
	4	0.0827	-	-	-	-	-	0.2	-	-	190.99
	5	0.03105	-	-	-	-	-	0.45	-	-	95.48
	-: the material is isotropic.

Table 5

 5 Flutter velocities [rpm] and frequencies[Hz] for the NACA 0012.

				h
				[À45]	[À30]	[À15]	[0]	[15]	[30]	[45]
	C(k)=1	TE3	X f	424.37	501.81	428.66	212.37	232.90	307.86	305.95
			f f	16.03	16.15	13.10	13.50	17.76	21.05	19.14
		TE4	X f	378.53	470.20	409.85	212.47	224.21	292.59	289.53
			f f	13.69	15.28	12.64	13.00	17.09	19.96	18.30
	Theodorsen	TE3	X f	446.14	518.62	455.69	289.15	262.31	347.59	351.41
			f f	15.26	16.07	13.48	12.07	16.17	19.38	18.00
		TE4	X f	385.21	490.92	437.16	278.64	252.38	330.59	332.50
			f f	13.05	15.27	13.00	11.59	15.59	18.39	17.21
	Loewy	TE3	X f	431.05	252.10	105.04	281.51	247.13	331.36	346.73
			f f	14.95	0.283	0.151	11.57	14.36	17.64	17.16
		TE4	X f	379.01	239.01	103.13	253.72	239.01	315.89	327.82
			f f	13.00	0.270	0.150	9.871	13.94	16.80	16.38


Edouard Canot

HAL is

Energy approach: Gibbs free potential

Defining the Gibbs free energy as: (P. Attard, 2000;Elliott & Voitcu, 2007)

G(,) =  ( A LV -A SL cos  -2 V  )
where A LV and A SL are the liquid-vapor and solid-liquid surface area then with  fixed (not V!)

F =  ( 2  r 0 +   r 0 2 ) F = ∂ G ∂ ϵ

Distribution of the gaps: D()

χ (κ)= ∂ Ḡ ∂ β Ḡ(κ)=N g ∫ 0 ϵ max D(ϵ) G(ϵ , κ)d ϵ ∂ D(ϵ) ∂ ϵ ∂ G(ϵ) ∂ ϵ χ (κ)=- 1 3 N g [ ∫ 0 ϵ max (1+ϵ) ∂ D ∂ ϵ G(ϵ , κ)d ϵ+D(0) Ḡ(0) ]
Non constant shift due to the volumetric dilatation of the grains' packing Boundary effect at the origin, comes from the integration by parts

  Liquid menisci equilibrium in a wet granular medium in pendular regime -Cohesive strength estimation for random smooth spheres packings É. CANOT, R. DELANNAY Institut de Physique de Rennes, UMR 6251, Campus de Beaulieu, 35042 Rennes, France 1 -Can a simple model well represent the cohesion? 2 -Two smooth spherical grains' system 3 -Statistics for random spheres' packings 4 -Computing the cohesive strength References Sandcastle: a small water content gives strong capillary forces! In pendular regime, capillary forces come from surface tension Maximum of this cohesive strength should occur at the end of the pendular regime… (see Richefeu et al., 2016, for a realistic 2D simulation) Control parameters: R radius of the spherical grains  half of the gap between them  total curvature of the meniscus surface  surface tension  contact angle liquid-solid Variables obtained after computation: r 0 radius of the neck r B radius of the wet surface of the solid part V volume of the liquid meniscus S surface area of the liquid meniscus F attractive force between the two grains Meniscus shape computed numerically: integration of an ODE system, in an axisymmetric coord. system, specifying a constant total curvature  Numerical results for r 0 , V, S, F: stored in a database Main tendencies for the behaviour of the capillary force F(,V) With respect to the volume V: →force ~ Cst for contact (=0) (even for infinitesimal volume!) With respect to the gap : →force is always decreasing Leads to unstable situation as the following: →humidity aggregates grains as clusters: (Staron et al., 2005)

  →numerically obtained via DEM (L. Oger, IPR) Density of Gibbs free energy: The Cohesive Strength : → average of the normal components of the tensile tensor (Richefeu et al., 2016); (avoid the full tensile tensor, as in Scholtès et al., 2009) → computed via a Uniform Volumetric Expansion where  is the dilatation → doesn't need the direction of the contacts, uses only statistical information (see below) Equilibrium → all liquid menisci have the same curvature   related to Relative Humidity via Kelvin equation Recast in terms of instead of because D() is well approximated by a mathematical function: → needs to work with a shifted distribution of gaps isotropic case increasing radius (no gravity) + monodisperse → compacity: 0.64 anisotropic case falling balls (gravity) + quasi-monodisperse → compacity: 0.61 -low radius STD: 0.005

→

  Good behavior but not realistic packing (because lot of true contacts must occur!) gaps: Log-Cauchy distribution no true contact between spheres density of gaps (per R 3 ): N g ~ 1.03 5 -Results and Conclusion gaps: Dirac + Uniform distribution 46% of true contacts between spheres density of gaps (per R 3 ): N g ~ 0.82 → Easy to parameterize -isotropic case: median of log(gap) + sharpness -anisotropic case: % of true contacts + distribution width Albert et al., 1997, Maximum angle of stability in wet and dry spherical granular media, Phys. Rev. E, 56, 6, pp 6271-6274 Halsey & Levine, 1998, How sandcastles fall, Phys. Rev. Letters,

  already emphasized by Bocquet et al., 2002 → avoids true contacts between smooth spheres → allows the decreasing of the cohesive strength w.r.t. liquid volume (or saturation): see the regime II of Halsey & Levine, 1998 → another possibility is the cone-plane contact model of Albert et al., 1997Edouard.Canot@univ-rennes1.fr Renaud.Delannay@univ-rennes1.fr

	Pendular	Pendular
	Funicular	Funicular
		Behaviour not expected!
		Smooth spheres model not appropriate!
		→ Needs introduction of roughness

80, 14, pp 3141-3144 P. Attard, 2000, Thermodynamics analysis of bridging bubbles…, Langmuir, 16, pp 4455-4466 Bocquet et al., 2002, Physics of humid granular media, C. R. Physique 3, pp 207-215 Staron et al., 2005, Multi-scale analysis of the stress state in a granular slope…, Eur. Phys. J. E, 18, pp 311-320 Elliott & Voitcu, 2007, On the thermodynamics stability of liquid capillary bridges, Canadian J. of Chem. Engng, 85, pp 692-700 Scholtès et al., 2009, On the capillary stress tensor in wet granular materials, Int. J. Num. Anal. Meth. Geomech., 33, pp 1289-1313 Richefeu et al., 2016, Lattice Boltzmann modelling of liquid distribution in unsaturated granular media, Comput. Geotech., in press, DOI:10.1016/j.compgeo.2016.02.017 Three-regime model of Halsey & Levine, 1998 → capillary force decrease as liquid volume tends to zero ? Advantages of roughness introduction: → → Strong dependency with the contact angle, as expected → Plots w.r.t. saturation (strongly linked to Relative Humidity and curvature  via the Kelvin equation)


Audrey Juban

Stephanie Briancon

François Puel
email: francois.puel@centralesupelec.fr

Stéphanie Briancon

Prof F Puel

Impact of compression speeds on phase transition of caffeine

Keywords: Processing-induced-transformations (PITs) during direct compression: Direct compression ; Tableting, Polymorphic transformation, Processing induced transition (PITs), Process parameters, Caffeine

published or not. The documents may come

Processing-induced-transformations (PITs) during direct compression: impact of compression speeds on phase transition of caffeine

Introduction

For pharmaceutical industry, processing-induced-transformations (PITs) in drug products manufacturing are of great importance. Several authors have started to work on these transformations [START_REF] Brittain | Effects of mechanical processing on phase composition[END_REF][START_REF] Zhang | Phase transformation considerations during process development and manufacture of solid oral dosage forms[END_REF][START_REF] Descamps | Perspectives on the amorphization/milling relationship in pharmaceutical materials[END_REF][START_REF] Juban | Processing-induced-transformations (PITs) during direct compression: Impact of tablet composition and compression load on phase transition of caffeine[END_REF] . Solid-state transitions that may occur during manufacturing processes can affect the quality of the final dosage form [START_REF] Bauer | Ritonavir: an extraordinary example of conformational polymorphism[END_REF] .

Even now, manufacturers have no control over these phenomena. A recent theoretical approach of the PITs assumes that each step of a process can generate a stress (thermal, mechanical or chemical) which could cause a modification in the thermodynamic equilibrium of the system and lead to a phase transformation of the product [START_REF] Morris | Theoritical approaches to physical transformations of active pharmaceutical ingredients during manufacture processes[END_REF][START_REF] Petit | Transformations physiques de cristaux moléculaires induites par les procédés industriels de fabrication : approche théorique, conséquences et exemples. Cristal 3[END_REF] .

Nowadays, tablets are the most popular pharmaceutical dosage forms in production [START_REF] Wu | A comparative study of compaction properties of binary and bilayer tablets[END_REF][START_REF] Hennigan | Quantitative polymorph contaminant analysis in tablets using raman and near infra-red spectroscopies[END_REF] . In the past few years, the use of direct compression has rapidly increased [START_REF] Zhang | Physical properties and compact analysis of commonly used direct compression binders[END_REF] due to its economic interest and its process which avoids the steps of wet granulation and drying processes. In the pharmaceutical field, it is well known that active pharmaceutical ingredients (APIs) and excipients may exist under different crystalline forms depending on the treatment that they received [START_REF] Chieng | An overview of recent studies on the analysis of pharmaceutical polymorphs[END_REF][START_REF] Nauha | Predicting" polymorphs of pharmaceuticals using hydrogen bond propensities: probenecid and its two single-crystal-to-single phase transitions[END_REF] . For solid dosage form development, it is recommended to use the stable form of the API, as the starting material [START_REF] Chakravarty | Phase transformation in thiamine hydrochloride tablets: influence on tablet microstructure, physical properties, and performance[END_REF] . However, as each polymorphic form has its own physico-chemical properties (such as solubility, dissolution rate, and bioavailability), metastable forms may be used in pharmaceutical development because of the desired properties for the final product. Many parameters in the manufacturing process itself could have an impact on these phase transitions. In the literature, several authors have observed polymorphic transformations occurring during the compression process [START_REF] Chan | Polymorphic transformation of some drugs under compression[END_REF][START_REF] Lefebvre | Polymorphic transitions of carbamazepine under grinding and compression[END_REF][START_REF] Boldyreva | Effect of pressure up to 5.5 GPa on dry powder samples of chlorpropamide form-A[END_REF][START_REF] Otsuka | Nano-and macro-geometrical structural change of caffeine and theophylline anhydrate tablets during hydration process by using X-ray computed tomography[END_REF][START_REF] Mazel | Polymorphic transformation of anhydrous caffeine under compression and grinding: a re-evaluation[END_REF] .

In this work, caffeine (C8H10N4O2) was used as a model API. Two polymorphic forms of caffeine are known and called Form I and Form II. The commercial anhydrous caffeine Form II is stable at room temperature until 145 °C and the Form I is stable at ambient pressure from the transition temperature (145 °C) until its melting point 236 °C19 . This common molecule, which presents simple enantiotropic behaviour and a reversible transformation at solid state, is well known for polymorphic studies [START_REF] Mazel | Polymorphic transformation of anhydrous caffeine under compression and grinding: a re-evaluation[END_REF][START_REF] Pirttimakki | Effects of grinding and compression on crystal structure of anhydrous caffeine[END_REF][START_REF] Hubert | Process induced transformations during tablet manufacturing: phase transition analysis of caffeine using DSC and low frequency micro-Raman spectroscopy[END_REF][START_REF] Hédoux | Polymorphic transformation of anhydrous caffeine upon grinding and hydrostatic pressuring analyzed by low-frequency raman spectroscopy[END_REF] .

In a previous work [START_REF] Juban | Processing-induced-transformations (PITs) during direct compression: Impact of tablet composition and compression load on phase transition of caffeine[END_REF] , three parameters that could impact the phase transition of the caffeine Form I (CFI) into Form II (CFII) have been analysed: nature of the diluent, dilution of the API and compression pressure.

The two first formulation parameters have an impact on the CFI transition degree. By contrast, for the two binary systems studied, the authors have shown that the transition degree induced by the compression was not dependent on f the compression pressure.

This study, complementary to this previous work, focused on the evaluation of the impact of another process parameter on the phase transition of caffeine Form I: the compression speed (i.e. velocity of the mobile punch) during tableting. Indeed, it is well known that compression speed, applied during tableting, can have significant effects on the compaction properties of pharmaceutical powder [START_REF] Wu | A comparative study of compaction properties of binary and bilayer tablets[END_REF][START_REF] Garr | An investigation into the capping of paracetamol at increasing speeds of compression[END_REF][START_REF] Garr | The effect of rate of force application on the properties of microcrystalline cellulose and dibasic calcium phosphate mixtures[END_REF][START_REF] Larhrib | Compressing polyethylene glycols: the effect of compression pressure and speed[END_REF][START_REF] Tye | Evaluation of the effects of tabletting speed on the relationships between compaction pressure, tablet tensile strength, and tablet solid fraction[END_REF] . Garr and Rubinstein [START_REF] Garr | An investigation into the capping of paracetamol at increasing speeds of compression[END_REF] have showed that the capping pressures decrease with increasing speed of compression. Moreover, Larhrib et al. [START_REF] Larhrib | Compressing polyethylene glycols: the effect of compression pressure and speed[END_REF] display an increase in the mean yield pressure with the compression speed.

Consequently, the compression speed during tableting may also have an impact on the polymorphic transition of the caffeine Form I into Form II. Up to our knowledge, this point was not yet investigated in the literature. In this study, tablets were made at three different compression speeds. During tableting, the compression speed was constant until the maximum compression pressure defined was reached. The control of this parameter was possible with the compression technique used (electric press), which was not the case in the previous study (compression driven by force). In addition, for each velocity of the mobile punch studied, tablets made from three compression pressures were analysed. The determination of the solid-solid phase transition of caffeine Form I was performed by thermal analysis.

Material and methods

Materials

The caffeine Form I (CFI) was used as a model API. It was obtained by using the same method as described by Juban et al. [START_REF] Juban | Processing-induced-transformations (PITs) during direct compression: Impact of tablet composition and compression load on phase transition of caffeine[END_REF] . In order to ensure that the particles exhibit a 100 % polymorphic purity (Form I), a Differential Scanning Calorimetry (DSC) analysis was performed for each batch of CFI (results not shown).

Microcrystalline cellulose (MCC), Avicel ® PH-102, was obtained from FMC Biopolymer.

The particle size distribution for these two materials was measured with a Mastersizer® 3000 using an AeroS dry dispersion unit (Malvern Instruments, U.K.). The Fraunhofer model was selected in order to transform the dynamic light scattering measurements in a particle size distribution. The resulting volume density function was averaged. The calculation of the different standard percentile (Dv10, Dv50 and Dv90) was performed following the British standard BS2955.1993. Measurements were performed in triplicate on each powder. Their values with their standard deviations are listed in Table 1.

Table.1: Particle size distribution parameters of CFI and MCC (n=3)

Tablets preparation

Cylindrical tablets were prepared with a binary mixture of anhydrous caffeine Form I (a day after its manufacture) and MCC as a diluent. Diluent was mixed in various proportions with caffeine content of 78 and90 wt% of CFI. The blending was performed with a tridimensional mixer (Turbula ® T2F) at a rotational rate of 49 rpm (for 10 minutes). A series of pure CFI tablets (100 wt% of CFI) was also made. All compacts were made at constant mass (about 150 mg) by using a tensile test system (ElectroPuls TM E10000, Instron ® , Élancourt, France), and its data acquisition software (WaveMatrix TM , Instron ® ). A sensor of 10kN was used.

Compressions were performed with an electric press. The compression is driven by displacement until reaching the required force. This technology allows obtaining a constant compression speed. The acquisition of the position of the punch and the pressure is made with a frequency of 10 Hz. The mean relative deviation of the compression speed during tableting is of 5%. Specific tools with flat-punches 6.50 2) were measured just after compression according to the European Pharmacopoeia method section 2.09.05 [START_REF]European Directorate for the Quality of Medicines & HealthCare (EDQM)[END_REF] .

The storage was made at room temperature during three months. 115

Table.2: Tablet characteristics

(-) (wt% CFI) (mm.min -1 ) (MPa) (mg) (mm) (mm) (g.cm -3 ) (g.cm -3 ) (%)
198.

Transition kinetic study by thermal analysis

The Differential Scanning Calorimetry (DSC) device used for the calorimetric measurements was a Q200 ® from TA Instruments. Calibration was performed with indium. Aluminum crucibles with perforated caps were used and the heating rate was 10°C.min -1 from 20°C to 270°C under a nitrogen purge of 50 mL/min.

MCC releases water when heated above 40 • C. Thus and to avoid an increase of pressure in the crucible, a hole was made in the lids for DSC measurements. This method ensured a good reliability of the experiments. Analyses were made in triplicate. DSC takes into account a macroscopic fraction (about 10 mg) of the tablet and was applied in order to estimate the transition degree in the core of tablet. For the studies of CFI transition degree, a kinetic analysis has been performed for each mixture (uncompressed powders and tablets), at the following specific times: 2, 15, 30, 45, 60 and 90 days after tablets fabrication.

Uncompressed powders were considered as the references, in order to investigate the tableting impact by comparing the transition kinetic profiles between tablets and uncompressed powders.

Concerning crucibles preparation, two methods were used depending on the type of the sample: powder or tablet. For powders, around 10 mg of sample was put in the crucible. For a correct and reproducible DSC analysis, the bottom of each crucible must be completely covered by the sample. Indeed, this method allows ensuring a homogeneous temperature distribution throughout the sample. Consequently, a special tool was designed for these thermal analyses: a flat punch with a diameter slightly lower than the crucible one. This tool allow to pack the powder on the bottom of the cap. For tablet a similar method was used.

As describe by Fig. 1, the core of tablet was cut by using a cutter blade. Then, this part was put in the crucible, without grinding step, and flattened by using the special tool in order to entirely cover the bottom of the cap. These methods of preparation, without grinding step, allow ensuring that no mechanical stress was applied on samples. Consequently, the crucibles preparation did not induce the polymorphic transformation of the caffeine Form I.

As described by Juban et al. [START_REF] Juban | Processing-induced-transformations (PITs) during direct compression: Impact of tablet composition and compression load on phase transition of caffeine[END_REF] , the transition degree (τ) of the caffeine Form I into Form II was obtained with Eq. 1:

𝜏𝜏 = ∆ 𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 𝐻𝐻 (𝐼𝐼𝐼𝐼→𝐼𝐼) 𝑜𝑜𝑓𝑓 𝑡𝑡ℎ𝑒𝑒 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑒𝑒 ∆ 𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 𝐻𝐻 (𝐼𝐼𝐼𝐼→𝐼𝐼) 𝑜𝑜𝑓𝑓 𝑠𝑠𝑝𝑝𝑝𝑝𝑒𝑒 𝑓𝑓𝑜𝑜𝑝𝑝𝑠𝑠 𝐼𝐼𝐼𝐼 • ∆ 𝑓𝑓𝑓𝑓𝑡𝑡 𝐻𝐻 (𝐼𝐼) 𝑜𝑜𝑓𝑓 𝑠𝑠𝑝𝑝𝑝𝑝𝑒𝑒 𝑓𝑓𝑜𝑜𝑝𝑝𝑠𝑠 𝐼𝐼 ∆ 𝑓𝑓𝑓𝑓𝑡𝑡 𝐻𝐻 (𝐼𝐼) 𝑜𝑜𝑓𝑓 𝑡𝑡ℎ𝑒𝑒 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑒𝑒 (1)
The experimental validation of this method (results shown in Juban et al. [START_REF] Juban | Processing-induced-transformations (PITs) during direct compression: Impact of tablet composition and compression load on phase transition of caffeine[END_REF] ) was performed, for each composition studied, in order to confirm the Expression (1). Fig. 1: Sample fraction (black part) of tablet for DSC Analysis.

Statistical analysis

In this work, statistical analyses were performed using the R software (The R foundation) with a threshold of 5% corresponding to the alpha risk. As all data were not normally distributed and variances were not equal, consequently, data could not be normalized, a non-parametric Kruskal-Wallis test was performed.

In order to determine significant differences between two parameters, a multiple comparison post hoc test (kruskalmc) was conducted.

Results

Evolution of the CFI transition degree (τ) in uncompressed powders

During three months, a kinetic study was conducted in order to observe the evolution over time of the CFI transition degree (τ) in uncompressed powders (see Fig. 2). Each symbol represents one of the three compositions studied: pure CFI and binary mixtures of CFI/MCC: 90 and 78 wt% of CFI.

Two transition profiles were observed: one for pure CFI and another one for powder mixtures.

Uncompressed powder of pure caffeine Form I showed a first quick transition stage over the first 15 days until reaching a plateau (about 30 %) from the fifteenth day. However, for the two binary mixtures, this quick transition stage has lasted 60 days. After that, a plateau value of about 55 % was reached. For the transition degree values, the difference observed between the two mixtures was not significant.

Evolution of the CFI transition degree (τ) in tablets

Impact of the velocity of the mobile punch during tableting

As with uncompressed powders, a kinetic study of the CFI transition degree occurring in tablets was also conducted during three months. For each composition, results obtained for tablets were compared to the ones obtained in uncompressed powders. In a first step, the impact of compression load on CFI transition degree in tablets was conducted. Whatever the composition studied and the compression speed tested, no significant difference was observed between the three compression loads (data not shown).

Consequently results obtained for tablets made at different compression pressures were gathered.

In order to evaluate the impact of the velocity of the mobile punch, the evolution over time of the CFI transition degree in tablet made at three different compression speeds (50, 500 and 4500 mm.min -1 ) was measured. For clarity reasons, Fig. 3 presents the kinetic study of the CFI transition degree for one composition (100 wt% of CFI). However, the same results were observed for other compositions.

Uncompressed powder and tablets made at three different compression speeds are represented by the cross and full symbols, respectively.

The transition profile in tablets was similar to the one observed with uncompressed powder: a first rapid growth stage until reaching a plateau. The plateau value was more important in tablets than in uncompressed powders. Moreover, for all compositions, a constant gap of around 25% was observed between uncompressed powders and tablets. Whatever the compression speed, the transition degree was almost identical just after compression and the transition profiles were similar all along these kinetic studies. The velocity of the mobile punch did not alter the transition of the caffeine form I into form II after compression. Nevertheless, the content of the transformed phase (caffeine form II) was initially more important in tablets than in uncompressed powder.

Impact of the API dilution

In order to observe the influence of the API dilution on the transition of caffeine form I into form II, Fig. 4 displays the evolution over time of the gap between the CFI transition degree obtained in uncompressed powders and in tablets for the three compositions studied. Each composition is represented by a symbol.

Considering that compression pressure and compression speed have not an impact on the CFI transition degree, the results obtained for tablets made at different compression loads and speeds were regrouped (n = 27 and n = 3 for tablets and uncompressed powders, respectively). Then, for each composition studied, the difference between uncompressed powder and tablets, for the transition of the caffeine Form I into Form II, was calculated. This method allows ensuring a satisfactory comparison between the different compositions. Two situations were observed depending on the tablet formulation. Indeed, for pure CFI, the calculated gap slightly increased during the first 30 days (between 16 % and 26 %) until reaching a plateau from the thirtieth days (around 30 %). By contrast, for binary mixtures of CFI/MCC, gap was kept constant (around 25 %) all along the kinetic study. The addition of a diluent in the formulation seems to have an impact on the CFI transition degree, but this impact appears to be independent of the diluent content.

Discussion

Evolution of the CFI transition degree (τ) in uncompressed powders

Whatever the composition studied the transition of the metastable caffeine Form I into its stable form (Form II) was rapidly observed in uncompressed powders. These results are consistent with those find in the literature by Mazel et al. [START_REF] Mazel | Polymorphic transformation of anhydrous caffeine under compression and grinding: a re-evaluation[END_REF] , Hédoux et al. [START_REF] Hédoux | Low-and High-Frequency Raman Investigations on Caffeine: Polymorphism, Disorder and Phase Transformation[END_REF] and Juban et al. [START_REF] Juban | Processing-induced-transformations (PITs) during direct compression: Impact of tablet composition and compression load on phase transition of caffeine[END_REF] and confirm that the caffeine Form I which is metastable at room temperature and atmospheric pressure tends to move progressively with time towards its stable Form II.

Evolution of the CFI transition degree (τ) in tablets

Several similarities were observed between results obtained in this study and the ones presented by Juban et al. [START_REF] Juban | Processing-induced-transformations (PITs) during direct compression: Impact of tablet composition and compression load on phase transition of caffeine[END_REF] , even though the tableting presses were different. In fact, for all compositions, CFI transition degree was higher in tablets than in uncompressed powders. As reported by some authors in the literature [START_REF] Mazel | Polymorphic transformation of anhydrous caffeine under compression and grinding: a re-evaluation[END_REF][START_REF] Pirttimakki | Effects of grinding and compression on crystal structure of anhydrous caffeine[END_REF] , these results mean that a part of this transformation was induced by the direct compression process.

Moreover, the compression pressures did not have an impact on the transition of the metastable caffeine Form I into its stable Form II. For CFI, whatever the intensity of the compression pressure applied, the transition degrees just after compression were similar. These results mean that the compression pressure has just a triggering effect when a sufficient compression pressure is applied during tableting. Moreover, the results show that the transition degree of the caffeine form I into form II just after compression was not dependent on the mobile punch velocity whatever the composition (pure caffeine or caffeine mixed with a diluent).. For this polymorphic system, the control of the punch velocity is unnecessary since a very small velocity (long compression time) and a high velocity (short compression time) led to the same transition degree. This process parameter seems to be inefficient for controlling polymorphic transition induced by the compression.

Fig. 5 presents the evolution over time of the gap between the CFI transition degree obtained in uncompressed powders and in tablets for the three compositions studied. A comparison was made between results obtained in this work and the ones presented by Juban et al. [START_REF] Juban | Processing-induced-transformations (PITs) during direct compression: Impact of tablet composition and compression load on phase transition of caffeine[END_REF] . Each composition was represented by a symbol. Full or empty symbols correspond to the gaps calculated in this study or the ones presented in the previous work, respectively. During the kinetic study of pure caffeine Form I, the calculated gap for the CFI transition degree between uncompressed powder and tablets displays a slight increase during the first days before reaching a plateau. By contrast, the ones calculated for binary mixtures were kept constant all along the kinetic analysis (around 20%). Results obtained for mixtures of CFI/MCC are in agreement with the ones presented by Juban et al. [START_REF] Juban | Processing-induced-transformations (PITs) during direct compression: Impact of tablet composition and compression load on phase transition of caffeine[END_REF] . No significant difference was observed between these gaps, meaning that the polymorphic transition of the caffeine Form I is no alter by process parameters applied during tablet manufacturing. Consequently, it could be interesting to express the results for the CFI transition degree in gap between uncompressed powder and tablets than in percentage of the transformed phase. This method allows comparing results obtained from different tableting processes.

𝑋𝑋(𝑡𝑡) = 1 -𝑒𝑒 �-� 𝑡𝑡 𝜏𝜏 𝑁𝑁 � 𝑡𝑡 � (2)
Where t the time, X the volume fraction of the transformed phase, τN the characteristic kinetic time, and n the exponent indicating the mechanism which controls the transition. A value of n lower than 1 means that the transition mechanism is only dependent on the nucleation.

The calculated profiles are shown in Fig. 6 and compared with the experimental points. The fitted parameters n and τN are displayed in Table 3. The mean relative deviation between the experimental and the calculated volume fractions is 5 % both for uncompressed powders and tablets. These results confirm that the transition kinetic rate measured in this study could also be correctly represented with this stretched exponential law. The values of n are in a close range from 0.21 to 0.50. This results are similar than the ones observed by Juban et al. [START_REF] Juban | Processing-induced-transformations (PITs) during direct compression: Impact of tablet composition and compression load on phase transition of caffeine[END_REF] . Whatever the compression technology used, the transition mechanism is nucleation dependent.

Conclusion

In order to improve the understanding of the PITs occurring during the direct compression process, a complementary study to a previous one (Juban et al. [START_REF] Juban | Processing-induced-transformations (PITs) during direct compression: Impact of tablet composition and compression load on phase transition of caffeine[END_REF] ) was performed in order to investigate the impact of the velocity of the punch on the polymorphic transition since this process parameters is known to influence tablet mechanical properties. This was possible by using an electric press instead of a compression simulator.

Several similarities with the previous work were observed. The CFI transition degree was higher in tablet than in uncompressed powder, meaning that the direct tableting process has the ability to trigger the polymorphic transformation of the caffeine form I into form II. Moreover, the value of the CFI transition degree was not altered by the maximum compression load. In the same way, this study shows that the punch velocity had not a decisive influence on the level of the polymorphic transition of caffeine form I into form II just after compression. This observation is corroborated by the analysis of the mechanism of the transition occurring in the tablet at rest which is independent on the press technology used. Controlling the punch velocity seems not relevant for preventing or enhancing the transition of a solid form in a tabletting process

The incorporation of a diluent in the formulation had a small impact on the CFI transition degree.

Nevertheless, unlike the results presented by Juban et al. [START_REF] Juban | Processing-induced-transformations (PITs) during direct compression: Impact of tablet composition and compression load on phase transition of caffeine[END_REF] , the impact of the dilution of the caffeine Form I with MCC as diluent wasn't observed in this study. This difference could be due to the press compression technology used for tableting.

2 a

 2 Average calculated from 3 measurements mm of diameter, and die set, were fitted on the device. Compression was made by the upper punch. The compression forces were measured with the sensor located on the upper punch. In this work, tablets were 110 produced without lubricant. Three different compression speeds were studied: 50, 500, 4500 mm.min-1 for three different compression forces (1.66 kN, 3.32 kN, 6.64 kN) corresponding at 50 MPa, 100 MPa and 200 MPa in compression pressure, respectively. The tablet weight and thickness (shown in Table

Fig. 2 :

 2 Fig.2: Evolution over time of the CFI transition degree in uncompressed powders for the three compositions studied, n=3.

Fig. 3 :

 3 Fig.3: Evolution over time of the CFI transition degree in uncompressed powders and tablets of pure CFI made at three different compression speeds, n=3 for uncompressed powder and n=9 for tablets.

Fig. 4 :

 4 Fig.4: Evolution over time of the gap between the CFI transition degree in uncompressed powders and tablets for the three compositions studied.

Fig. 5 :

 5 Fig.5: Evolution over time of the gap for CFI transition degree between uncompressed powders and tablets for the three compositions studied both in this work and results presented by Juban et al. (2016)

Fig. 6 :

 6 Fig.6: Evolution over time of the volume fraction of the transformed phase. Comparison between experimental data and calculated values obtained from the kinetic laws (2) for uncompressed powder and tablets, respectively.

Table . 3: Comparison of the values of τN and n calculated in this study with the kinetic laws (2) for uncompressed powder and tablets, respectively, with the ones presented by Juban et al., 2016.

 .

Acknowledgements

The support of the Ministère de l'Enseignement Supérieur et de la Recherche is gratefully acknowledged.

We would like to thank Medel'Pharm especially T. Ménard, G. Tardy and B. Villa as well as the technician G. Conod-Nardi B. Ponsard and O. Pollet for their technical support. We also wish to thank the IVTV platform (ANR-10-EQPX-06-01).

Declaration of interest

The authors report no conflicts of interest. Only authors are responsible for the content and writing of this paper.