Y Bury

M Saudreau

J Borée

G Charnay

Experimental and numerical study of a con¿ned variable density turbulent jet submitted to a strongly pulsed coÀow

An experimental and numerical work dedicated to the study of a con¿ned variable density jet subjected to a time-varying coÀow with high acceleration/deceleration levels is presented in this paper. A global analysis and a hyperbolic model based on momentum equation written for the centerline excess velocity are ¿rst developped to obtain physical understanding and to estimate relevant effects governing such a complex unsteady Àow. In order to observe the jet as it interacts with the main Àow and to quantify those effects, two-component laser Doppler velocimetry and optical concentration measurements based on Mie scattering are used and adapted to unsteady conditions. Light, heavy and air jets are successively studied. During acceleration phases, convective phenomena govern part of the Àow and a front formation appears which strength increases as initial jet density is lower. This front propagates at a constant velocity and deeply modi¿es mixing processes and mass transfer in comparison with steady coÀowing jet. Comparisons with RANS numerical simulations using n model are presented. Agreement with measurements of longitudinal velocity ¿eld is acceptable.

Introduction

Turbulent mixing of a jet exhausting in a different density coÀow is met in numerous industrial situations. The present work is directly linked with development of natural gas vehicles (NGV) which have a high potential to reduce urban air pollution. In such spark ignition engine, the injected gaseous fuel in the intake ports ( ' fDD) is submitted to a pulsed air Àow where acceleration can reach fff} values (} is the gravity acceleration ¿eld). Knowledge of the mixing processes in those particular unsteady conditions is useful to optimize combustion [START_REF] Fulachier | InÀuence of density variations on the structure of low-speed turbulent Àows: a report on euromech 237[END_REF]. Previous study concerning light jet in a pulsed crossÀow [START_REF] Raud | Behaviour of a variable density heated jet subjected to a pulsated crossÀow[END_REF][3] [START_REF] Raud | Experimental study of the behaviour of con¿ned variable density jets in a time varying crossÀow[END_REF] has shown the complex 3D inÀuence of this unsteadiness on the global behavior of the jet. In the present work and owing to the complexity of measurements in real engines, a model experiment setup, considering a variable density turbulent jet exhausting in a strongly pulsed coÀowing duct Àow, has been developped to isolate and analyse the main physical mechanisms governing such unsteady Àows. The experimental setup and the measurements techniques, quoted below, will ¿rst be described.

A physical analysis based on a global integral approach and a hyperbolic model derived from momentum equation written for the centerline excess velocity will show that the Àow is mainly governed by the balance between unsteady and convective effects. This shows a front formation during acceleration phase, convected downstream at a ¿nite velocity, depending on the initial jet density. Buoyancy effects due to the Àow acceleration ¿eld are found to be globally weak.

Two-component laser Doppler velocimetry and laser planar tomography visualizations adapted to unsteady conditions, as well as & " simulations con¿rm this behaviour. Measurements also bring relevant informations concerning strong modi¿cation of turbulence production and anisotropy during acceleration and deceleration phases. Mixing processes and mass transfer are deeply modi¿ed by the time-varying coÀow, in comparison with corresponding steady coÀowing jet.

Experimental Setup

Flows Generation and Main Characteristics

The experimental setup is sketched in ¿gure 2. It consists of a ( ' e66 diameter cylindrical tube jet con¿ned in a Sf Sf66 transparent square duct. The use of a sonic nozzle ensures a time-constant mass Àux jet. Nozzle location is optimized in order to obtain a fully developped turbulent pipe Àow with a L ' Sf6*r axial velocity at the jet exit. The total length of the duct generating the coÀow is 2ef6. The channel Àow is generated by a vacuum pump and pulsed thanks to a motorized butterÀy valve at a s ' .DM5 frequency. Pulsation is tuned with the channel acoustics and corresponds to a half-wave mode [START_REF] Charnay | Periodic Àow in a wind-tunnel produced by rotating shutters[END_REF]. Duct entrance is a pressure node. As a result, for about thirty jet diameters downstream the duct inlet, including test section, coÀow velocity L is quasi-uniform and varies from D to f6*r with time (see ¿g.4). The jet therefore develops in a spatially uniform, time varying coÀow. Resulting values of acceleration/deceleration L rise up to eff}. To bring to the fore main effects governing such a Àow, light, homogeneous and heavy jets are studied. Density variations are obtained considering an air channel Àow and successively air, CO2 ( ' S) and an air-helium mixture ( ' fDD) jet Àow. The corresponding Reynolds numbers, based on outlet jet diameter, vary from bc fff to fc fff, depending on the initial density. The study focuses from f to 2D diameters downstream the jet outlet. Particular care has been taken considering experimental boundary conditions, in view of numerical simulation comparison. To do so, a convergent type collector and a grid, placed at the duct entrance (¿g.2), just upstream the studying zone, generate uniform coÀow velocity pro¿le and isotropic homogeneous turbulence, and minimize the inÀuence of the jet tube's wake.

Duct acoustics was determined by pressure and velocity measurements. In order to observe the jet as it interacts with the main Àow and to quantify those effects, two-component laser Doppler velocimetry (Dantec BSA) and optical concentration measurements based on Mie scattering are used and adapted to unsteady conditions. An original pressurized and regulated oil seeding system was specially designed in this way (see ¿g 1). Considering LDV measurements, coÀow and jet Àow seeding rates are set equal to limit resulting measurements bias. Considering tomography measurements, only the coaxial jet is seeded. An encoder is connected to the motorized butterÀy valve in order to phase average unsteady results. More than Dff samples per encoder degree are used to average LDV data and 2ff images for tomography.

Axial velocity measurements from f to 2D jet diameters and radial pro¿les respectively at f, D, f, D, 2f diameters downstream jet outlet were achieved. For each injected Àuid (air, CO2 and air-helium mixture jets), measurements were performed under the four following coÀow conditions: Unsteady coÀow, steady coÀow at respectively L ' D6*r, L ' .D6*r, L ' f6*r (corresponding to unsteady phases A, B, C, D in ¿g.4). A global analysis is based on momentum excess equation integrated over the cross-section of the duct. In an unsteady situation this balance equation writes as follows ( [START_REF] Bury | Etude expérimentale modèle du mélange de carburants gazeux (GNV) dans les moteurs alternatifs du secteur automobile[END_REF]) :

Results and Discussion

Y Y| ]] 4EL L _r ' Y Y% ]] 4L EL L _r YL Y| ] ] E4 4 _r (1)

E| E|2 E|

In fact a third term, writing UU 4EL L _r , thouroughly appears on the right-hand side of equation ( 1). This term can be interpreted as the inÀuence of con¿nement on the Àow. However, measurements easily show that this term is negligible for all the cases studied, so we will not take it into account. Equation (1) states that the temporal change of the momentum excess E| in a jet slice is due to the longitudinal variation of the excess momentum Àux E|2 plus a buoyancy term E| if gas density of the jet is different from air density. Phased variations of those terms, for the air-helium jet case at f( and D( downstream the jet outlet, are plotted in ¿gure 3. It may be noted that only terms E| and E| are directly computed from mean velocity and density ¿elds measurements. The momentum excess Àux E|2 is then deduced from balance equation ( 1) : E|2 ' E| E|.

For distances corresponding to f( and D( downstream the jet outlet, ¿gure 3 shows that during deceleration (from phase f to phase bf ) the buoyancy term E| is negative and contributes to decrease the Àow momentum excess. The situation is just the opposite during acceleration (from phase bf to phase Hf). However we can notice that buoyancy effects are weak when compared to convective one. In this way, pulsed Àow is globally governed by the balance between convective and temporal variation terms.

During the deceleration phase, the excess momentum Àux at the jet exit a E| ' L EL L E| considerably increases since the jet Àux is time-independant. This excess momentum Àux therefore "feeds" excess momentum in the coÀowing jet, downstream the jet outlet, and term E| is positive. The reverse situation is obtained during acceleration phase. a E| decreases and the longitudinal variation of excess Àux "evacuates" excess momentum. Figure 3 shows that those two sequences do not occur at the same phase for ' f and ' D, where % is the downstream distance from the jet outlet. Moreover, comparison with time evolution of coÀow velocity L E| (¿g 4) shows that for ' D, a sharp increase of the local excess momentum still occurs at the beginning of the acceleration phase. In fact, considering the axial centerline velocity L and quasi-steady situations, one can see that L is maximum when the coÀow velocity L E| is the highest, and minimum when L E| is the lowest (see ¿g.7 and 5). As perturbations are convected at velocity L , the excess momentum evacuated downstream focuses during acceleration phase in a starting jet fashion and is responsible for the formation of a velocity front. This behaviour is more formally addressed by studying the local evolution of the excess centerline velocity. A hyperbolic model based on momentum equation for the centerline excess velocity L L using a simple hypothesis for turbulence closure was obtained. This model will not be detailed here [START_REF] Bury | Etude expérimentale modèle du mélange de carburants gazeux (GNV) dans les moteurs alternatifs du secteur automobile[END_REF] but it shows that the convective nature of the Àow is responsible for a focusing on characteristic curves (i.e. the formation of a "shock wave") during acceleration phase only.

The unsteady behavior of the pulsed coÀowing jet is seen to be complex but can thus be analysed with fairly general arguments. Time scale comparison can be useful to estimate a priori the unsteadiness "level" of the jet. For a convective phenomenon, | ' is an appropriate time scale for the jet. If L , | sharply increases when one moves downstream. In the present situation, if | is evaluated using the steady jet at mid coÀow velocity L ' ' .D6*r and compared with | '

, the time scale of the periodic pulsation ( s ' .DM5 is the pulsation frequency), for ' D. From the jet exit to f( the jet is therefore quasi-static. On the contrary, the region of interest is fully unsteady.

Local Analysis

The axial evolution of the mean longitudinal velocity at phases A, B, C, D (¿g.4) is displayed in ¿gure 5 for the light jet. Each result is then compared with steady coÀow velocity measurements in ¿g.6 and 7. A very large amplitude variation of the axial velocity is detected in ¿gure 5. The axial velocity at phase C, corresponding to the end of the deceleration phase, is much lower than the velocity for a corresponding steady coÀow L ' D6*r (¿g.7). On the contrary, the mean axial velocity at the end of the acceleration phase (A) is higher than the velocity for a corresponding steady coÀow L ' f6*r. At phases B and D, the instantaneous coÀow velocity seen by the jet is the same ( L =17.5m/s). The acceleration (respectively deceleration) value are however maximum. ¿gure 6 shows that the longitudinal evolutions of axial velocity for phases B and D lie on both sides of the steady coÀow case L ' .D6*r. Similar behaviour is observed for homogeneous and heavy jets. One noticeable feature is the formation of the velocity front during acceleration phase. Signature of front propagation can be obtained by measuring longitudinal velocity L (¿g 14) and its gradient along jet axis (¿g 9). Measurements show that this front appears during acceleration phase (from phase C to phase D) for : D, which was predicted by dimensional analysis, and propagates at a constant velocity T ' b6*r which corresponds to the slope on the E%c | diagram in ¿gure 9.

As long as we are interested in mixing we may wonder if this behaviour modi¿es mixing properties of the jet. In ¿gure 10 and 8 are shown phase averaged density ¿elds and centerline This phenomenon looks like starting jet ( [START_REF] Zhang | Effects on acceleration on turbulent jets[END_REF]) and the head of the perturbation is expected to be a kind of spherical vortex. However in our case, during front formation, ambiant Àuid velocity still strongly varies. Thus the front growth is alterated by the outer pressure gradient.

Another important aspect of this front formation is that generated longitudinal gradient velocity (¿g 9) dramatically increases : (¿g 11) component via turbulent energy production ( [START_REF] Borée | Measurements and image analysis of the turbulent ¿eld in an axisymmetric jet subject to a sudden velocity decrease[END_REF]). In ¿gure 13 structural parameter introduced by Towsend & '

, where (resp. ) is the longitudinal (resp radial) rms velocity, is plotted versus Eoc |. Anisotropy is enhanced as the front propagates and decreases behind the front when energy produced on :component is transferred to : and :components. This strong alternate anisotropy may induce serious problems concerning & " modelisation since it is based on isotropic hypothesis. Computations are presently performed with ESTET computer code. Turbulence modelisation, in a Favre averaged formulation, is based on either & " or -" closure schemes [START_REF] Simonin | Numerical implementation of second-moment closures and application to turbulent jets[END_REF] and boundary conditions are given by experimental measurements. We present in ¿gures 5,8,15 and 16 ¿rst simulation results with standard & " model for axisymetric Àow ( ' S).

Conclusion

Experimental and numerical investigations based on two-component laser Doppler velocimetry, laser sheet tomography visualizations and & " simulations have been performed to study the dynamic behaviour of a variable density jet submitted to a strongly pulsed con¿ned coÀow. This con¿guration modelizes natural gas injection in intake ports of vehicles (NGV). This gaseous fuel has a high potential to reduce urban air pollution. A large part of this work was devoted to isolate and analyse the inÀuence of the unsteady effects on mixing processes and mass transfer between the jet and the coÀow.

A ¿rst global integral approach and a hyperbolic model based on momentum equation written for centerline excess velocity has been developped to predict the behaviour of the Àow. Convective effects globally governe the Àow. This shows a front formation during acceleration phase, convected downstream at a ¿nite velocity, depending on the initial jet density. Velocity and density measurements as well as & " calculations con¿rm this. Therefore, mixing is strongly affected by this front. Turbulence production and anisotropy are either deeply modi¿ed during acceleration and deceleration phases.

In real engine con¿gurations, as the jet behaviour highly depends on the phase of the pulsed air Àow, the instant of injection will be a crucial parameter for optimization of intake mixing process.

Figure 1 :

 1 Figure 1: Experimental setup

Figure 2 :

 2 Figure 2: Presentation of the test section

3. 1 .

 1 Global AnalysisWe consider the situation displayed in ¿gure 2. The jet Àuid has an ejection velocity L and an initial density 4 possibly very different from the ambient density 4 . The jet Àows parallel to the time varying coÀow L E|. The effect of the gravity ¿eld } is neglected in all the cases presented in this paper. First, it can easily be checked that the jet Froude number de¿ned by 8 '[6] is signi¿cantly large. The initial development of the jet is therefore inertial and is not affected by gravity[7][8]. Second, we deal with contrasted situations with L }and | | where | ' is the time scale of the pulsation and | is the characteristic time scale of the movement induced by gravity along the studying zone length u. The large amplitude variation of the coÀow is expected to have a signi¿cant effect on the jet behaviour.

Figure 3 :Figure 4 :

 34 Figure 3: Phased averaged evolution of integral terms from equation (1)

  phase averaged U component velocity U = 30 m/s (A) U = 17.5 m/s (B) U = 5 m/s (C) U = 17.5 m/s (D)

Figure 5 :

 5 Figure 5: Light jet case comparison of L evolution between measurements (symbol) and & " simulation (lines)

Figure 6 :

 6 Figure 6: Light jet case comparison of L evolution between steady and unsteady coÀow for L ' .D m/s

Figure 7 :

 7 Figure 7: Light jet case comparison of L evolution between steady and unsteady coÀow for L ' D m/s and L ' f m/s

Figure 8 :

 8 Figure 8: Light jet case comparison of centerline density evolution between measurements (dashed line) and & " simulation (solid line) At a given station downstream jet exit during acceleration or deceleration phase two cumulative effects are present: change in excess momentum Àux at the jet outlet and change in the outer pressure gradient needed to accelerate the coÀow.One noticeable feature is the formation of the velocity front during acceleration phase. Signature of front propagation can be obtained by measuring longitudinal velocity L (¿g 14) and its gradient along jet axis (¿g 9). Measurements show that this front appears during acceleration phase (from phase C to phase D) for : D, which was predicted by dimensional analysis, and propagates at a constant velocity T ' b6*r which corresponds to the slope on the E%c | diagram in ¿gure 9.As long as we are interested in mixing we may wonder if this behaviour modi¿es mixing properties of the jet. In ¿gure 10 and 8 are shown phase averaged density ¿elds and centerline

Figure 9 :Figure 11 :

 911 Figure 9: Phased-averaged evolution of long. gradient of U component centerline velocity

  Figures 15 and 16 show that the phase averaged radial evolution computed with & " model are qualitatively correct. More work is presently performed to compare the turbulent ¿eld prediction with experiments.

Figure 13 :

 13 Figure 13: Phased-averaged evolution of radial pro¿le for structural parameter & '

Int. Symp. on Turbulence, Heat and Mass Transfer


Yong He

Zhongyong Fan
email: zyfan@fudan.edu.cn

Jia Wei

Suming Li
email: lisuming@fudan.ac.cn

Morphology and Melt Crystallization of Poly(L-lactide) Obtained by Ring Opening Polymerization of L-lactide With Zinc Catalyst

The morphology and melt crystallization of zinc catalyzed poly(L-lactide) (PLLA) were investigated by using differential scanning calorimetry (DSC), polarized optical microscopy, and scanning electron microscopy. Isothermal melt crystallization performed at 95-135°C showed that the morphology depends on the degree of supercooling, as illustrated by crystallite perfection and lamellar thickening behaviors. Double melting peak was observed on DSC thermograms and attributed to the meltrecrystallization mechanism, small and imperfect crystals becoming gradually more stable ones. Circumferential and hexagonal cracks were detected in PLLA spherulites, which were formed during melt-crystallization at 135°C and quenching in liquid nitrogen. Rhythmic growth and thermal shrinkage are suggested to be the two main factors accounting for the formation of periodic cracks. Spherulite growth rates of PLLA were evaluated by using combined isothermal and nonisothermal procedures, and were analyzed by the secondary nucleation theory. The maximum growth rate reached 9.1 m/ min at 130°C. The temperature range investigated (120 -155°C) belongs to the Regime II of crystallization. The value of U* was found to be 1890 cal/mol, instead of 1500 cal/mol commonly used in literature, and K g and were estimated to be 3.03 ؋ 10 5 K 2 and 1.537 ؋ 10 ؊4 J/m 2 , respectively. As a result, no distinct difference between PLLA catalyzed by zinc metal and those prepared with stannous octoate catalyst exists in this work. POLYM.

INTRODUCTION

Degradable polymers have been attracting much attention for potential biomedical and environmental applications. As a typical degradable polyester, polylactide (PLA) can be used to fabricate various osteosynthetic devices, drug delivery systems, and tissue engineering scaffolds, etc. The degradation characteristics of PLA, which are of major importance for various applications, mainly depend on its morphology and crystallinity [1][2][3].

The crystallization and morphology of poly(L-lactide) (PLLA) have been extensively investigated [4 -7]. Pennings et al. examined the morphology and crystal growth of PLLA as a function of molecular weight (MW) and degree of supercooling by using polarized optical microscopy (POM) equipped with a hot stage. The Hoffman and Lauritzen equation was applied to calculate the nucleation constants. Morphological changes from spherulites to axialites were observed at 163°C because of lower degree of supercooling [4]. Miyata et al. obtained the nucleation parameters of PLLA of different MW under various isothermal conditions. The maximal overall isothermal crystallization rate and spherulite growth rate were obtained at 105°C and at 120°C, respectively [5]. On studying the cold crystallization of amorphous PLLA, Kawaguchi et al. proposed the mechanism of phase transformation. The two exothermic peaks detected on differential scanning calorimetry (DSC) were attributed to the cold crystallization and phase transition from ␤ form to ␣ form. A multiple endothermic peak was observed by annealing samples between 105 and 120°C [6]. On the other hand, Iannace reported that the Regime II-III transition occurs below 115°C [7]. In this work, isothermal melt crystallization of PLLA was investigated in detail as a function of isothermal time, with particular attention paid to double melting behaviors.

The spherulite growth rate (G) is generally determined under isothermal conditions. At a specific temperature, the plot of spherulite radius (r) vs. isothermal time (t) is linear and the slope G can be expressed as:

G ϭ dr dt .
Nonisothermal procedure is also employed to obtain G by monitoring the growth of r as a function of temperature (T).

G ϭ dr dT dT dt

where dr/dT is estimated from the first derivative of the r vs.

T plot and dT/dt is the cooling rate as the spherulite growth front advances [8].

Experimental data of spherulite growth rate have been analyzed according to the secondary nucleation theory by Hoffman. The Hoffman and Lauritzen equation for the growth rate of a linear polymer crystal with folded chains is given by:

G ϭ G 0 exp ͫ Ϫ U* R ͑T c Ϫ T ϱ ͒ ͬ exp ͩ Ϫ K g T c ⌬T f ͪ
where K g is the nucleation constant; ⌬T is the degree of supercooling defined by ͑T m 0 Ϫ T c ͒; T m 0 is the equilibrium melting point; T c is crystallization temperature; f is a factor expressed as 2T c /͑T m 0 ϩT c ͒; U* is the activation energy for segment diffusion to the crystallization site; R is the gas constant; T ϱ is the hypothetical temperature expressed as (T g Ϫ C) where all motion associated with viscous flow ceases, C is a constant (30 -50 K); and G 0 is the front factor [9].

It is well known that the activation energy of linear PLLA is a specific parameter depending on the repetitional diffusion of chain segments. However, most authors used the assumed value of U* ϭ 1500 cal/mol, which gave good fit of the growth rate data for a number of polymers [4,9,10]. In this work, the method of the nonlinear least squares curve fitting is used to apply the Hoffman and Lauritzen equation to spherulite growth rate data of PLLA, and the values of U* and K g are determined.

Isothermal spherulite growth of PLLA has been investigated by many authors [11][12][13][14]. Tsuji et al. found that PLLA crystallizes in Regime II or Regime III kinetics depending on T c if the number-average molecular weight (M n ) exceeds the critical values, and solely in Regime III kinetics for all the T c with the M n below critical values [11,12]. They also investigated the spherulite growth of low MW PLLA [13]. On the other hand, Abe et al. reported that the transitions from Regime II to Regime I and from Regime III to Regime II occur at around 147°C and 120°C, respectively [14].

Stannous octoate and zinc metal are two commonly used catalysts or initiators in ring opening polymerization of lactides. The morphology and crystallization of PLLA initiated by stannous octoate have been investigated by many authors [4 -7, 11-13]. Nevertheless, the presence of residues was suspected and indirectly shown, such as octanoyl moieties combined to some of the alcoholic chain ends, octanoic acid, and hydroxy tin octoate. Stannous octoate can be hardly removed from the chain ends, which could affect the structure and morphology of polymers [15]. In contrast, there are no chain end modifications in the case of zinc catalyzed PLLA, which is an ideal material for model studies of polymer crystallization and morphology [15,16]. In this article, we report on melt crystallization behaviors, morphology, and spherulite growth rates of zinc catalyzed PLLA.

EXPERIMENTAL

Poly(L-lactide) (PLLA) was synthesized by ring opening polymerization of L-lactide using zinc metal as catalyst, according to the standard procedure reported by Vert et al. [16]. The obtained polymer was purified by the dissolution/ precipitation method using dichloromethane as solvent and ethanol as nonsolvent to remove low molecular weight species and residual monomer. The weight-average molecular weight (M w ) of PLLA as determined by GPC is 1.36 ϫ 10 5 g/mol, with a polydispersity index of 2.2. The rotatory power ͓␣͔ 20 D measured in chloroform solution is Ϫ153°, which is consistent with literature data [17][18][19].

Differential scanning calorimetry (DSC) experiments were conducted by using a DuPont DSC 910 instrument calibrated with indium standard. Each sample of 5(Ϯ0.1) mg was scanned at a heating rate of 10°C/min under N 2 atmosphere after various thermal treatments. The summit of endothermic peak was taken as the melting temperature (T m ). In the process of isothermal melt crystallization, the samples were melted at 200°C for 3 min, rapidly cooled to the crystallization temperature, and then allowed to crystallize for 30 min.

Spherulite morphologies were examined with an Olympus model BH-2 polarized optical microscope (POM) and a Philips XL30 scanning electron microscope (SEM). The spherulite growth rate was measured by using a Leica DMLP POM equipped with a Linkam THMSE-600 hot stage. Thin films of PLLA were obtained by cold pressing under a pressure of 20 bars with an oil press of 769YP-15Z Model. Each sample was melted at 200°C between two slides, kept for 3 min to erase any thermal history, and then subjected to different thermal treatments after cooling from the melt at 50°C/min. The radius of growing crystals was monitored during crystallization with a video camera system taking photomicrographs automatically at appropriate time intervals. Dry nitrogen gas was applied throughout the hot stage.

RESULTS AND DISCUSSION

Melting Behaviors

Figure 1 shows the DSC thermograms of PLLA samples after various thermal treatments. The glass-transition tem-perature (T g ) of quenched PLLA is detected at 59.5°C with an enthalpy change of 34.5 J/kmol as estimated from heat flow change during the transition (curve a). The cold crystallization temperature (T c,c ) appears at 112.2°C and T m at 170.8°C, in agreement with literature data [5,6,19,20]. The melting endothermic enthalpy is about 48 J/g, which is equal to the exothermic enthalpy of cold crystallization. Various values of PLLA perfect crystal ͑⌬H m 0 ͒ have been reported, ranging from 91 to 148 J/g [5,[21][22][23][24][25]. The value of 91 J/g reported by Wunderlich was used to obtain the degree of crystallinity (X c ), which is equal to 53% [22].

After melt crystallization at various temperatures, the glass transition and cold crystallization are hardly detectable (curves b-h). A small exothermic peak is detected between the two endothermic melting peaks for samples crystallized at 95 and 100°C (curves b, c), indicating that recrystallization occurs upon heating. A double melting peak is observed after melt crystallization at 105, 110, and 115°C (curves d, e, f). In the case of higher melt crystallization temperatures (120°C, 125°C), single melting peak appears (curves g, h).

Figure 2 shows the melting behavior of PLLA after melt crystallization at 110°C for various time periods up to 24 h. Two endothermic peaks are detected in each case, which are denoted as T m,1 and T m,2 , respectively. With increasing isothermal time, T m,1 shifts to higher temperature, while T m,2 slightly decreases. The melting temperature and corresponding enthalpy data are presented in Table 1. It is noteworthy that the total melting enthalpy (⌬H m ) of PLLA, including both endothermic peaks, remains the same (c.a. 45 J/g) despite the isothermal time. In fact, with increasing isothermal crystallization time, initially formed small and imperfect crystals are subjected to crystal perfection and lamellar thickening. Therefore, T m,1 corresponding to thickened lamellae increases, and the value of ⌬H m,1 becomes higher. On the other hand, both T m,2 and ⌬H m,2 diminish due to less melt-recrystallization. Double-melting behavior has been observed for many semicrystalline polymers, such as poly(ether ether ketone), poly(ethylene terephthalate), and iso-polystyrene (iPS) [26 -29]. Melt-recrystallization model and dual lamellar thickness model have been proposed in the case of iPS [28,29]. Yasuniwa et al. investigated the T m,1 and T m,2 dependence on the heating/cooling rates, suggesting that the melt-recrystallization mechanism be taken into account for the double melting behavior of PLLA [30].

When the degree of supercooling is relatively high (Fig. 1, curves b, c), the initially formed lamellae can hardly evolve because of the high transportation energy, resulting in a small endothermal peak of T m,1 followed by meltrecrystallization. On the contrary, with lower degree of supercooling, almost perfect crystals are formed during melt crystallization, thus suppressing melt-recrystallization. As a consequence, a single melting endotherm is observed.

Cracks in PLLA Spherulites

An interesting phenomenon was observed by using POM and SEM: PLLA spherulites isothermally crystallized at 135°C and quenched in liquid nitrogen exhibit circumferential and hexagonal cracks (Figs. 3 and4). No extinction rings are visible on the spherulites, indicating that lamellae twisting and scrolling are not involved in the spherulite growth process, and consequently, are not related to the formation of the periodic cracks here [31]. In literature, optical micrographs exhibiting such cracks were obtained; however, no closer attention was paid to them [32,33]. Similar patterns were reported for other polymers such as poly(3-hydroxybutyrate) (PHB) and poly(ethylene oxide) (PEO). Zhang et al. observed fine circles on banded spherulites of PHB, which were assigned to terraces resulting from the instability of the moving boundary of the melt film [34]. Terraces were also observed by Ding et al. in the case of PEO, which were formed at the front of the growing spherulites just before or exactly when two spherulites impinge on each other [35]. In other words, these fine and concentric circles were not related to lamellae organization in spheru-lites. In the case of our PLLA spherulites, however, the periodic cracks detected seem to be different from the fine circles or terraces of PHB and PEO.

Spiral and target patterns are common morphologies observed not only in excitable media and nonlinear dissipative systems, but also in polymeric, organic, and inorganic materials. It appears that the periodic cracks in PLLA spherulites present similar characteristics. Kyu et al. employed the time-dependent Ginzburg-Landau-equations model C (TDGL model C) to analyze the pattern forming aspects of polymer crystallization. The formation of spiral and target growth patterns in polymer spherulites was demonstrated theoretically and used to explain the twisted lamellar orientation [36]. In a recent article, Lotz and Cheng excluded the possibility that such spherulite growth features account for the twisted, untwisted, or scrolled lamellae. The authors considered unbalanced surface stresses as the origin of twisting and scrolling of polymer crystals [31]. Interest- ingly, the theoretical simulation seems applicable to our results. In the case of PLLA, the periodic cracks formed in spherulites with no extinction rings are quite similar to the target patterns. It can be assumed that these cracks result from the rhythmic growth of target spherulites, as simulated by the phase-field TDGL model C. As crystalline-rich bands and amorphous-rich bands grow rhythmically, analogous to the theoretical simulation, the difference between those two bands induces formation of cracks on cooling process.

Subsequent investigations revealed that such regular and periodic cracks disappeared when cooling the spherulites at ambient temperature, although a few cracks were occasionally observed under unpolarized conditions. Therefore, the periodic cracks resulted from the shrinkage of spherulites on cooling process, which are different from terraces observed in the form of fine circles and extinction rings formed simultaneously with spherulite growth [34]. Different radial and circumferential thermal expansion coefficients that generate large stresses would be responsible for the occurrence of such cracks on cooling. Thermal shrinkage was previously proposed by Keller et al. to explain the cracks of PHB spherulites [37].

According to the arrangement of the crystallographic axes in the spherulites, the b axis is always radial, while a and c are tangential. While the thermal expansion coefficient along the c axis is negligible, there should be a difference between ␣ r , the thermal expansion coefficient along the b axis, and ␣ c concerning the a axis. Thus, a large tensile stress in the radial direction on cooling spherulites would induce circumferential or hexagonal cracking when ␣ r is higher than ␣ c .

With increasing spherulite radius, the outer cracks were gradually smoothed to become circumferential, the inner ones being hexagonal (see Fig. 4). These inner cracks might result from the (1,1,0) growth plane as hexagonal shaped crystals are formed at the early stage of melt crystallization at 135°C. Therefore, it can be concluded that rhythmic growth features and thermal shrinkage are the two main factors accounting for the formation of periodic cracks. It is noteworthy that these two factors have never been considered together to explain crystallization behaviors in literature. Further investigations are underway to better understand this phenomenon.

Spherulite Growth Rate (G) of PLLA

When a PLLA sample is cooled directly from the melt at 1°C/min, it would be difficult to obtain G values at high temperatures because of the long induction period prior to crystallization. A self-nucleation procedure was applied to enlarge the temperature range in the nonisothermal crystallization measurements. At first, the melted sample was cooled at 50°C/min to 130°C. Once the first spherulite was observed, the sample was heated up to 155°C at 50°C/min and kept for 1 min followed by cooling at 1°C/min. Figure 5 shows the variation of the spherulite radius as a function of temperature. The radius of the spherulite increases with decreasing temperature. G values were calculated by taking the first derivative of the r vs. T plot, according to the following formula:

G i ϭ Ϫ 1 2 ͩ r iϩ1 Ϫ r i T iϩ1 Ϫ T i ϩ r i Ϫ r iϩ1 T i Ϫ T iϪ1 ͪ .
The results were smoothed by adjacent-averaging and shown in Fig. 6. With decreasing crystallization temperature, the spherulite growth rate increased to attain a maximum value, G max ϭ 9.1 m/min, at about 130°C and then decreased. The G max value is larger than 5.1 m/min and 6.7 m/min, which were previously obtained respectively at 128°C in literature [4,8]. Such a difference in growth rates could not be ascribed to the influence of molecular weights since all those samples possessed similar molecular weights. On the other hand, PLLA using zinc as catalyst has hydroxyl and carboxyl groups at both chain ends, while crystallization of stannous octoate catalyzed PLLA might be affected by the presence of octanoyl moieties at the chains ends. The higher crystallization temperature for and smaller chain ends could account for the higher G max value obtained in this work. G values were also obtained under isothermal conditions at specific temperatures (130, 140, 145°C). Figure 7 shows the radius growth as a function of time at 130°C, the slope of the plot giving the growth rate. G values of 9.1, 6.6, and 4.3 m/min were obtained for crystallization at 130, 140, and 145°C, respectively. Therefore, spherulite growth rates determined by nonisothermal procedure are consistent with those by isothermal crystallization.

In the spherulite growth rate equation described by Hoffman and Lauritzen, we could set G 0 , U*, and K g as three unknown parameters, T c as the independent variable, G as the dependent variable. C equals to 30°C as an acceptable value according to previous reports, [4] T m 0 ϭ 207°C [4, 21], and T g ϭ 59.5°C. Nonisothermal procedure provided G values at various temperatures, which allowed us to obtain the defined parameters by using the method of the nonlinear least squares curve fitting. The fitting curve shown in Fig. 6 corresponds to the following parameter values: G 0 ϭ 5.1 ϫ 10 9 m/min, U* ϭ 1890 cal/mol, K g ϭ 3.03 ϫ 10 5 K 2 .

Figure 8 presents the plot of ln G ϩ U*/R ͑T c Ϫ T ϱ ͒ vs. 1/͑T c ⌬T f ͒, in which the U* value was 1890 cal/mol, substituting the assumed value of 1500 cal/mol commonly used in literature [4,7,8]. The plot exhibits very good linearity, giving the same K g value. There is no deflection in the plot, showing the absence of Regime I-II transition or Regime II-III transition. The value of K g is in agreement with literature data, and the temperature range belongs to Regime II [4,8].

In Regime II, K g is defined as ͑2b 0 e T m 0 /⌬H m 0 k͒, where b 0 is the layer thickness; is lateral surface energy;

e is fold surface free energy; ⌬H m 0 is heat of fusion per unit volume, and k is the Boltzmann constant. Using the values reported in literature [38], b 0 ϭ 0.517 nm, ⌬H m 0 ϭ 1.11083 ϫ 10 8 J/m 3 , e ϭ 6.089 ϫ 10 Ϫ4 J/m 2 , the value is obtained ( ϭ 1.537 ϫ 10 Ϫ4 J/m 2 ), which is higher than that reported in literature ( ϭ 1.203 ϫ 10 Ϫ4 J/m 2 ) [4,38].

CONCLUSIONS

Isothermally melt crystallization of PLLA synthesized by using zinc catalyst was studied with different degrees of supercooling at 95-135°C. The results showed that the morphology depends on the degree of supercooling, as illustrated by crystallite perfection and lamellar thickening behaviors. Double melting peak was attributed to the meltrecrystallization mechanism, small and imperfect crystals becoming gradually more stable ones. It was also found that the whole endothermal enthalpy remained the same as the lower melting peak turned to be larger with increasing isothermal time.

Circumferential and hexagonal cracks observed in PLLA spherulites were examined. Rhythmic growth and thermal shrinkage are identified as the two main factors accounting for the formation of periodic cracks.

Isothermal and nonisothermal procedures were combined to evaluate the variation of spherulite growth rate with crystallization temperature. The maximum G value is about 9.1 m/min, which is higher than those reported in literature. The data were analyzed by using the secondary nucleation theory. The method of the nonlinear least squares curve fitting was employed, directly giving the value of U*. And crystallization at the temperature range (120°C Ϫ155°C) belongs to Regime II, with U* ϭ 1890 cal/mol, K g ϭ 3.03 ϫ 10 5 K 2 , and ϭ 1.537 ϫ 10 Ϫ4 J/m 2 .

FIG. 1 .

 1 FIG. 1. DSC curves of various PLLA samples. (a) quenched after 2 min melting at 200°C; (b)-(h) isothermally melt crystallized for 30 min at 95, 100, 105, 110, 115, 120, 125°C, respectively.

FIG. 3 .

 3 FIG. 3. Polarized optical micrographs of PLLA spherulites melt crystallized at 135°C. (a), (b) fine circles, (c), (d) hexagonal cracks.

FIG. 5 .

 5 FIG. 5. Variation of the spherulite radius as a function of temperature during cooling at 1°C/min from 155 to 120°C after self-nucleation at 130°C.

FIG. 7 .

 7 FIG. 7. Variation of the spherulite radius with time during isothermal crystallization at 130°C. FIG. 8. Kinetic analysis of growth rate data of PLLA (U* ϭ 1890 cal/min, T ϱ ϭ T g Ϫ 30°C).

TABLE 1 .

 1 Double melting temperature and enthalpy data of PLLA after melt crystallization at 110°C for various time periods.

	Time (h) T m,1 (°C) T m,2 (°C) ⌬H 1 (J/g)	⌬H 2 (J/g) ⌬H m (J/g)
	0.25	164.0	172.1	21.1	24.5	45.6
	0.5	164.2	171.9	21.2	24.2	45.4
	1	164.8	171.8	24.2	20.5	44.7
	4	165.1	171.8	29.1	17.0	46.1
	10	165.7	171.5	31.3	13.9	45.2
	24	166.1	171.4	36.8	9.2	46.0

POLYMER ENGINEERING AND SCIENCE-2006 DOI 10.1002/pen

POLYMER ENGINEERING AND SCIENCE-2006 DOI 10.1002/pen

POLYMER ENGINEERING AND SCIENCE-2006 DOI 10.1002/pen

ACKNOWLEDGMENTS

The authors thank Professor Rongshi Cheng at Nanjing University for the use of POM equipment and valuable discussions.

Contract grant sponsor: Chinese National Nature Science Foundation; contract grant number: 20474013.


M Saudreau

Analysis of a variable density turbulent jet submitted to a pulsed coÀow

An experimental and numerical work dedicated to the study of a variable density jet submitted to a time-varying coÀow with high acceleration / deceleration levels is presented in this paper. To deal with this complex unsteady Àow a simple model taking into account advective effects has been developped. It permits to bring to the fore that advective effects are responsible for a front creation during acceleration phase. This front will interact with the downstream jet structure and modify turbulence and mixing properties. As this study is linked with industrial application we simulated the entire Àow by a standard n % model. Comparisons with LDV measurements are presented and discussed.

INTRODUCTION.

Turbulent mixing of a jet exhausting in a different density coÀow is met in numerous industrial situations. The present work is directly linked with development of natural gas vehicles (NGV) which have a high potential to reduce urban air pollution. In such spark ignition engine, the injected gaseous fuel in the intake ports ( @ 3=88) is submitted to a pulsed air Àow where acceleration can reach 6333j values (j is the gravity acceleration ¿eld). Knowledge of the mixing processes in those particular unsteady conditions is useful to optimize combustion (Fulachier,1989). Previous study concerning light jet in a pulsed crossÀow (Bury,1997Raud, 1999) has shown the complex 3D inÀuence of this unsteadiness on the global behavior of the jet. In the present work and owing to the complexity of measurements in real engines, a model experimental setup, considering a variable density turbulent jet exhausting in a strongly pulsed coÀowing duct Àow, has been developed to isolate and analyse the main physical mechanisms governing such unsteady Àows.

EXPERIMENTAL AND NUMERICAL SETUP.

Experimental setup.

The experimental setup (Bury,2000) is sketched in ¿gure 1. It consists of a G @ 7pp diameter cylindrical tube jet con¿ned in a 93 93pp transparent square duct. The use of a sonic nozzle ensures a time-constant mass Àux jet. Nozzle location is optimized in order to obtain a fully developed turbulent pipe Àow with a X @ 93p@v axial velocity at the jet exit. The total length of the duct generating the coÀow is 5=73p. The channel Àow is generated by a vacuum pump and pulsed thanks to a motorized butterÀy valve at a i @ :8K} frequency. Pulsation is tuned with the channel acoustics and corresponds to a half-wave mode [START_REF] Charnay | Periodic Flow in a Wind-Tunnel Produced by Rotating Shutters[END_REF]. Duct entrance is a pressure node. As a result, for about thirty jet diameters downstream the duct inlet, including test section, coÀow velocity X is quasi-uniform and varies from 8 to 63p@v with time (see ¿g.2). The jet therefore develops in a spatially uniform, time varying coÀow. Resulting values of acceleration/deceleration b X rise up from .733j to :33j.

To bring to the fore main effects governing such a Àow, light, homogeneous and heavy jets are studied. Density variations are obtained considering an air channel Àow and successively air, CO2 ( @ 4=9) and an air-helium mixture ( @ 3=88) jet Àow. The corresponding Reynolds numbers, based on outlet jet diameter, vary from <> 333 to 63> 333, depending on the initial density. The study focuses from 3 to 58 diameters downstream the jet outlet. Particular care has been taken considering experimental boundary conditions, in view of numerical simulation comparison. To do so, a convergent type collector and a grid, placed at the duct entrance (¿g.1), just upstream the studying zone, generate uniform coÀow velocity pro¿le and isotropic homogeneous turbulence, and minimize the inÀuence of the jet tube's wake.

In order to observe the jet as it interacts with the main Àow and to quantify X +w, @ 63 m/s, B: X +w, @ 4:=8 m/s deceleration, C: X +w, @ 8 m/s and D: X +w, @ 4:=8 m/s acceleration ) those effects, two-component laser Doppler velocimetry (Dantec BSA) and optical concentration measurements based on Mie scattering are used and adapted to unsteady conditions. Axial velocity measurements from 3 to 58 jet diameters and radial pro¿les respectively at 3, 8, 43, 48, 53 diameters downstream jet outlet were achieved. For each injected Àuid (air, CO2 and air-helium mixture jets), measurements were performed under the four following coÀow conditions: unsteady coÀow, steady coÀow at respectively X @ 8p@v, X @ 4:=8p@v, X @ 63p@v (corresponding to unsteady phases A, B, C, D in ¿g 2).

Numerical study.

RANS code Estet-Astrid [START_REF] Laurence | Numerical Implementation of Second-Moment Closures and Application to Turbulent Jets[END_REF] has been used to perform and test turbulence closure scheme for this Àow. Continuity, momentum and scalar transport equations are resolved in a Favre averaging formulation. Turbulence closure used is standard n % model adapted to axisymmetric jet, and U % standard model. For scalar quantities coupled with dynamic ¿eld, turbulence closure is performed using a constant Schmidt turbulence number set to 3=: in our case. Boundary conditions are set using experimental phased averaged quantities for both coÀow and jet inlet. To achieve comparison with experiment, we simulate the entire test section (¿g.1), starting from grid position to 63 diameters downstream. For computational facilities we consider a cylindrical duct which radius is set to obtain the same experimental coÀow mass Àow rate. Thus a 2D axisymmetric computation can be performed.

As experimental duct Àow is driven by pressure gradient variation, pulsation is simulated in the same way. We impose a variable outlet pressure S +w, and a constant coÀow inlet charge Thus at each time step, duct velocity X +w, ¿ts to impose pressure gradient. Taking into account linear frictional losses, we reproduce closely experimental velocity variation show in ¿gure.2.

In practical, the ¿rst step of simulation consists of calculating stationary axisymmetric jet with coÀow velocity equal to 63p=v . In a second step starting from this stationary jet, duct velocity is pulsed during two periods. During the ¿rst period, the initial condition is evacuated and only the second one is compared with experimental results.

ON THE LONGITUDINAL DEVELOPMENT OF THE PERTURBATION.

Experimental results.

The axial evolution of the mean longitudinal velocity at phases A, B, C, D (¿g.2) is displayed in ¿gure 3 for the light jet. Each result is then compared with steady coÀow velocity measurements in ¿g.4 and 5. A very large amplitude variation of the axial velocity is detected in ¿gure 3. The axial velocity at phase C, corresponding to the end of the deceleration phase, is much lower than the velocity for a corresponding steady coÀow X @ 8p@v (¿g.5). On the contrary, the mean axial velocity at the end of the acceleration phase (A) is higher than the velocity for a corresponding steady coÀow X @ 63p@v. At phases B and D, the instantaneous coÀow velocity seen by the jet is the same (X @ 4:=8p@v). The acceleration (respectively deceleration) value are however maximum. ¿gure 4 shows that the longitudinal evolutions of axial velocity for phases B and D lie on both sides of the steady coÀow case X @ 4:=8p@v. Similar behaviour is observed for homogeneous and heavy jets.

Physical analysis -derivation of an hyperbolic model

Momentum equation governing our Àow can be derived from Navier Stokes equations using classical hypothesis for quasi-parallel and turbulent Àows. Since the Àow is axisymmetric we will denote by +X+{> u> w,, Y +{> u> w,, Z +{> u> w,, mean velocity component in cylindrical coordinate system +{> u> , for axial, radial and azimuthal direction. In external region Àow is governed by pressure gradient and since con¿nement is negligible [START_REF] Bury | Experimental and numerical study of a con¿ned variable density turbulent jet submitted to a strongly pulsed coÀow[END_REF], governing equation reads [START_REF] Bury | Experimental and numerical study of a con¿ned variable density turbulent jet submitted to a strongly pulsed coÀow[END_REF] :

:=8p@v +X X , ~} . X +X X , . Y +X X , ~} @ 4 . = (1) t1 t2 t3 t4
Time variation of mean excess momentum along the Àuid particle trajectory is due to turbulent diffusion (w6) and to unsteady buoyancy term (w7). In the following we will consider only situations where buoyancy term is strictly zero (air jet) or small compared to other equation terms. Buoyancy contribution is detailled in [START_REF] Bury | Experimental and numerical study of a con¿ned variable density turbulent jet submitted to a strongly pulsed coÀow[END_REF].

The time variation of mean excess momentum X+3> u> w, X +w, is im- In steady coÀowing jets, mean advection and turbulent diffusion are in equilibrium and have therefore time scales of the same order of magnitude : w where X is the mean centerline velocity and [ the downstream distance from jet exit. Note that w increases strongly versus axial distance (w + , in a free steady jet). w is compared to w at phase C in ¿gure 7 where light, homogeneous and heavy jets are considered. One sees clearly that while the early development of the jet is quasi-steady (w w ), both time scale are of the same order of magnitude for all jet at about @ 48. The state of this turbulent Àow is therefore intermediate between quasi-steady and rapidly distorted Àows.

A simple physical model can be built to describe the spatio-temporal evolution of the centerline mean velocity X +{> w, for homogeneous jet Àows. If one assumes that the relative variation of mean excess momentum obtained at a given location { as a Àuid particle travels over distance { is the same in unsteady and steady situation. Considering Àuid particles, their trajectories in +{> w, diagram correspond to characteristic curves + , of equation : @ X +{> w,. Along these curves, if we note by +=, stationary quantities, the hypothesis made above implies that:

@

(2)

This expression can be integrated and it comes that along line + , expres-sion @ is conserved. A similar model has been used successfully in unsteady free jets (Borée,1996) and is able to reproduce experimental observations. It is important to note that equation 2 supposes an equilibrium between turbulent diffusion and instantaneous excess momentum. In particular this statement is true in the quasi-steady initial jet region as X +{> w, X +w, { and then { g{ whatever the phase [START_REF] Bury | Experimental and numerical study of a con¿ned variable density turbulent jet submitted to a strongly pulsed coÀow[END_REF].

In practical solution is advanced in time along lines + , by setting as a constant (¿gure 8). This methodology is the same as characteristic method used for solving hyperbolic equation (Whitham,1973). Considering a Àuid particle starting from jet exit ({ @ 3) at time w , its velocity is X +3> w , and the current velocity coÀow X +w ,. At time w A w the Àuid particle have been convected at local speed X +{> w, to the station {+w, @ U X +{+v,> v,gv. At this new location its velocity is given by: X +{> w, X +w, @ +X +3> w , X +w ,, ~} = ~}

(3) t1 t2

Term t1 represents the initial excess momentum viewed by Àuid particle at time w and term t2 is related to the role of turbulent diffusion according to stationary jet. Thus relation 3 expresses that mean excess momentum at +{> w, is therefore not associated to the current mean excess momentum at ejection X +3> w, X +w, but to X +3> w , X +w ,. This observation illustrates the phase lag due to convective effects in the jet structure. In particular, the dissymetry observed between acceleration and deceleration phase in ¿gure 4 is clear. Deceleration (resp. acceleration) is associated to lower (resp. higher) initial mean excess momentum X +3> w , X +w ,. We present in ¿gure 6 model and measurements results for the homogeneous jet when initial condition corresponding to stationary coÀowing jet with X @ 4:=8p@v. The pulsation used for X +w, is the experimental pulsation shown in ¿gure 2. As a ¿rst result, the Àapping of the jet centerline velocity is well reproduced. Second important point is that we mimic differences between phases B (Symbol ) and D (Symbol F). Thus hyperbolic model permits to explain jet structure differencies between acceleration and deceleration phases.

Following with this model, we can now consider two simple experiments. A ¿rst one we consider a stationary coÀowing jet at time w . At time w A w we decrease coÀow velocity (deceleration) At a given station {, as centerline velocity X +{> w, decreases in time, characteristic curves diverge and unsteady adaptation of the jet can be viewed in a +{> w, diagram as expansion waves.

On the contrary, the second experiment, starting from the same stationary jet, consists to increase coÀow velocity (acceleration). Thus X will increase and jet adaptation in +{> w, diagram will present merging characteristic curves. In a certain manner as propagation of compressive wave, a shock will eventually appear. For accelerating jets this front creation has already been mentioned and leads to a strong modifycation of turbulence and mixing process [START_REF] Zhang | Effects of acceleration on turbulent jets[END_REF].

EFFECT OF COFLOW ACCELERATION PHASE ON FLOW STRUCTURE -COMPARISON WITH & " MODEL 4.1 Unsteady adaptation of the jet structure.

To deal with accelerating case, we plot in ¿gure 9 for different downstream distances from jet exit, unsteady longitudinal velocity gradient adimensionalized by steady longitudinal velocity gradient for coÀow velocity of 8p@v versus time. We note a strong peack at @ 48 and downstream to @ 58. We note also that intensity increases with downstream distance. This high longitudinal velocity gradient appears during acceleration phase (<3 to 463 ) and corresponds to a break in X decaying slope. As previously discussed this phenomemun can be viewed as a consequence of the convergence of characteristic curves. Thus for ? 43 the longitudinal gradient does not strengthen since the jet is quasi-steady. In fact both three jets present an maximum of the local acceleration but the light jet acceleration is twice larger than the homogeneous and heavy jet acceleration. So only the light jet is presented.

Since we know maximum gradient location at any time, we are able to ¿nd its propagation law and deduce its velocity. It propagates at constant velocity Y but depends on density ratio (Figure 10). For heavy jet Y * 79p@v, for homogeneous jet Y * 5;p@v and for the light jet Y * 4<p@v. Similar phenomenum have been found out by Zhang & al (Zhang,1996) for accelerating homogeneous turbulent axisymmetric jet. In case of an imposed constant acceleration of the jet exit velocity, a front is created and is propagated at constant speed. They ¿nd out that front developpement is mainly controlled by upstream velocity variation. In our case during acceleration phase and at a ¿xed downstream position near the jet nozzle X +{> w, increases in a quasi-linear fashion in time. Referring to Zhang article this could explain our constant front velocities. It should be noted that those velocities differ accord- Physically and referring to previous hyperbolic model this phenomenum is due to increase of mean velocity near jet exit in the beginnning of acceleration phase i.e. merging of characteristic curves. The mean excess momentum stored in the jet structured during deceleration is evacuated rapidly downstream in starting jet type.

Effect on turbulence and mixing -Comparison with & " model.

In ¿gure 11 longitudinal x and radial y rms velocity are plotted versus radial distance at location @ 53 for three phases during coÀow acceleration. One before (<3 ), second during (453 ) and third after (493 ) front passage. As a ¿rst result the longitudinal gradient dramatically increases ? x A component via turbulent energy production (from <3 to 453 ). Secondly anisotropy is enhanced as the front propagates and decreases behind the front when energy produced on ? x Acomponent is transferred to ? y A and ? z Acomponents.

n simulation result is shown in ¿gure 12. Longitudinal rms velocity

x is deduced from turbulent energy n as x @ t n. Of course since n model is based on isotropic turbulence hypothesis, it can obviously not mimic anisotropic but turbulence production due to longitudinal gradient is well simulated but its intensity is lower than measured one.

In ¿gure 13 is plotted mean centerline density pro¿le versus downstream distance during front propagation. As a ¿rst results mixing process is modifying since decaying slope presents a strong longitudinal density gradient located near @ 48. n simulation overestimates longitudinal gradient and as a consequencies mispredicted front intensity. This problem deals with model constants which are adjust for axisymetric turbulent jet case. In our case, in Àow region where the jet dominates (w ? w i.e. near jet exit) simulation works well. Thus front velocities are correctly captured (Figure 10). On the contrary in Àow region where outer pressure gradient effects are stronger (w A w ) Àow behaviour is mispredicted and leads to an overestimation of front intensities (Figures 3 and13). In the view of industrial application we choose to simulate our Àow using a standard n model which is widely used for engine calculation. As a ¿rst result this model works pretty good as it mimics thetime evolution of the jet and front creation (¿gure 3) whatever initial jet density (¿gure 10). Nevertheless to deal with anisotropic turbulent production and unsteady buoyancy effects, higher order models are necessary to improve siluation prediction. To do so U and O=H=V= calculations are actually performed and results are still under investigation.

CONCLUSION.

This work based on experimental coupled with numerical study permits to bring to the fore main physical effects governing a turbulent jet submitted to a pulsed coÀowing stream.

The use of a hyperbolic model derived from momentum equation for centerline mean excess velocity has shown that jet response in a non quasi-static fashion. During acceleration phase the external coÀow pulsation leads to a front formation which deeply modi¿es jet properties. Production of turbulence is enhanced in an anisotropic way and mixing process is modi¿ed by comparison with stationary jet mixing.

Standard n model show good agreement with measurements in the limits of its possibility. Other simulation based on U and Large Eddy Simulation will permit to study turbulence modi¿cation in this Àow during front propagation.

In the view of industrial application, this study shows that unsteady Àows can not be deduced directly from stationary ones. Unsteady Àows imply memory effects, different time adaptation, phase lag effects and need inevitably to be studied. Thus simple experiments and numerical studies have to be done to isolate main unsteady physical effects governing such Àows. Thus by constant comparison between experiment and simulation classical physical models used in numerical tool can be tested and improved.

FigureFigure 2 :

 2 Figure 1: Presentation of the test section

Figure 4 :Figure 6 :

 46 Figure 3: Light jet case comparison of X evolution between measurements (symbol) and n simulation (lines)

Figure 7 :

 7 Figure 7: Time scale comparison at phase C (X @ 8p@v)

  Figure 9: Temporal evolution of mean longitudinal velocity gradient for light jet at different downstream distance from jet nozzle

  Figure 12: Experimental turbulence production (Symbols) during acceleration phase. Comparison with n prediction (Solid lines)

Figure 13 :

 13 Figure 13: Density pro¿l at phase 443 . Dashed line : measurements and solid line : n

saudreau@imft.fr

borée@imft.fr


L Montastruc
email: ludovic.montastruc@ensiacet.fr

C Azzaro-Pantel

B Biscans

M Cabassud

A Davin

S Domenech

L Pibouleau

INFLUENCE DES CONSTANTES D'EQUILIBRE ET DES COEFFICIENTS D'ACTIVITE SUR LA PRECIPITATION DES PHOSPHATES DE CALCIUM

Keywords: Clés. précipitation, constantes d'équilibres, calcium, phosphate, coefficient d'activité

The objective of this study is to quantify the influence of equilibrium constants both in liquid and solid phases and activity coefficients involved in phosphate calcium precipitation. The observed discrepancy in liquid phase equilibrium constants reported in the literature leads to determine an uncertainty zone describing the precipitation domain.

Résumé. Cette étude permet de quantifier l'influence des constantes d'équilibres tant de la phase liquide que de la phase solide et des coefficients d'activités sur les courbes d'équilibre de la précipitation des phosphates de calcium. L'écart entre les valeurs relatives aux constantes d'équilibre en phase liquide obtenues à partir d'une analyse bibliographique a conduit à déterminer une zone d'incertitude décrivant le domaine de précipitation.

Le rejet de phosphore dans le milieu naturel aqueux se traduit par un développement excessif des algues planctoniques et induit généralement une élévation du pH du milieu concerné accompagnée d'une diminution de l'oxygène dissous (phénomène d'eutrophisation). La réduction en apports phosphorés provenant des rejets urbains est ainsi considérée comme le facteur clé de la lutte contre la pollution des rivières et des lacs et se doit d'être particulièrement efficace. Le procédé retenu dans le cadre général de ce travail est basé sur la précipitation du phosphate dans un lit fluidisé de sable. Dans ce contexte, l'optimisation de la précipitation complexe des phosphates de calcium implique une bonne connaissance des mécanismes physico-chimiques mis en jeu et nécessite le développement d'un modèle suffisamment représentatif des phénomènes. Nous présentons dans ce qui suit un modèle basé sur la thermodynamique de la précipitation de deux espèces : l'ACP (Phosphate de Calcium Amorphe) et le DCPD (Phosphate de DiCalcium Dihydraté). L'objectif de ce travail est de quantifier l'importance des constantes d'équilibres et des coefficients d'activité sur la précipitation.

EQUATIONS DU MODELE

Les équations du modèle ont pour but de déterminer des courbes d'équilibres relatives à l'évolution du taux de conversion du phosphate en fonction du pH. Elles impliquent l'équation de sursaturation de l'espèce qui précipite, les bilans de conservation de matière (Phosphate et Calcium) ainsi que la relation d'électroneutralité. Lors de cette précipitation, les espèces considérées dans le milieu aqueux sont, d'une part, pour l'acide phosphorique, H 3 PO 4 , H 2 PO 4 -, HPO 4 2-, PO 4 3-, et d'autre part, la concentration en Ca 2+ et les sels de calcium correspondants. L'équation de la réaction de précipitation de l'ACP est la suivante :

2 3 4 3 4 2 3 2 ( ) Ca PO Ca PO + - + → (1)
La sursaturation de l'espèce ACP est définie par le paramètre β tel que :

2 3 4 2 3 3 2 4 ([ ] ) ([ ] ) 1 ln( ) 5 Ca PO ACP ACP Ca PO Ks λ λ β + - + - = (2)
où λ désigne le coefficient d'activité de l'ion considéré et Ks ACP le produit de solubilité de l'ACP L'équation de la réaction de précipitation de DCPD a lieu selon :

2 2 4 4

Ca HPO CaHPO

+ -

+ →

(3) De même, on définit la sursaturation de l'espèce DCPD :

2 2 4 2 2 4 ([ ] )([ ] ) 1 ln( ) 2 Ca HPO DCPD DCPD Ca HPO Ks λ λ β + - + - = (4) Figure 1.
i A i B i AB i K 1 H + H 2 PO - 4 H 3 PO K 2 H + HPO - 4 H 2 PO - 4 K 3 H + PO - 4 HPO - 4 K 4 Ca 2+ H 2 PO - 4 CaH 2 PO 4 + K 5 Ca 2+ HPO - 4 CaHPO 4 K 6 Ca 2+ PO - 4 CaPO 4 - K 7 Ca 2+ OH - CaOH + K w H + OH - H 2 O Tableau 1. Constantes d'équilibre pour le système Ca-PO 4 -H 2 O ( )( ) ( ) i i i i A B K AB =

INFLUENCE DES CONSTANTES D'EQUILIBRE

Les constantes d'équilibres relatives au système Ca-PO 4 -H 2 O peuvent être classifiées en deux groupes : les constantes concernant la phase liquide et celles relatives à la phase solide. L'influence des 8 paramètres pris en compte (6 en phase liquide et 2 en phase solide) a été évaluée à l'aide d'un plan d'expériences à deux niveaux afin de réduire le nombre de simulations. Le premier niveau concerne uniquement des constantes d'équilibres relatives à la phase liquide. Des résultats ont ainsi montré que certaines constantes d'équilibres à savoir {K1, K4, K6, K7} n'avaient aucune influence sur les courbes d'équilibres. Lors du deuxième niveau de planification, l'étude a porté sur l'ensemble des constantes d'équilibres { K2, K3, K5, pKs ACP , pKs DCPD }.

Niveau 1 de la planification d'expériences : Constantes d'équilibre en phase liquide L'influence des contantes d'équilibre de la phase liquide sur la courbe d'équilibre liquide-solide a été étudiée à l'aide des valeurs trouvées dans la littérature (Vieillard, 1978)

INFLUENCE DES COEFFICIENTS D'ACTIVITE

Plusieurs modèles de coefficients d'activités (Zematis, 1986) selon le domaine de variation de force ionique ont été établis. Dans la gamme de force ionique considérée, deux modèles simples sont couramment employés, respectivement celui de Davies et celui de Debye-Huckel (Zematis, 1986). Le modèle de Davies, souvent utilisé dans le système étudié (Seclker, 1994) ne nous semble pas judicieux car il ne tient compte de la physico-chimie de chaque ion présent dans le milieu. Le modèle de Debye-Huckel a été retenu dans une première approche. Son avantage est de permettre de faire varier les coefficients d'activités suivant les ions à travers leurs rayons hydratés. Il prend en compte le fait que les ions ne peuvent se rapprocher indéfiniment les uns des autres. Dans cette expression, le facteur α représente la distance la plus courte à laquelle l'anion, non hydraté, peut s'approcher du cation central hydraté. De plus, les interactions ion-ion et ion-solvant sont prises en compte par l'introduction du paramètre C DH , qui traduit la réduction de la valeur de la constante diélectrique au fur et à mesure que la solution est plus concentrée.

Le modèle de Debye-Huckel a pour expression : Etant donné que le choix du modèle de coefficient d'activité a un impact sur les courbes d'équilibre, l'influence d'une incertitude de plus ou moins 8% sur les coefficients d'activités calculés à partir du modèle de Debye-Huckel a été prise en compte. Cette étude a été réalisée à l'aide d'un plan d'expériences sur l'ensemble des 9 espèces ioniques présentes dans le milieu réactionnel. Des résultats ont montré que seuls les 4 coefficients relatifs à H + , HPO 4 2-, Ca 2+ et H 2 PO4 -ont une influence notable sur la courbe d'équilibre. Pour visualiser l'importance de ces coefficients d'activités sur les courbes d'équilibre, nous avons recherché les valeurs relatives à l'incertitude sur les 4 coefficients qui permettent d'avoir les mêmes courbes d'erreur obtenues avec les diverses constantes d'équilibre tant en phase liquide qu'en phase solide. Ainsi, une incertitude de 20% sur les coefficients d'activité correspond aux courbes d'incertitude obtenues avec les constantes d'équilibres du tableau 1. Néanmoins, la figure 2 montre qu'une incertitude de 10% influence déjà fortement les courbes d'équilibre.

2 10 1 DH i DH DH Log A z C B µ λ µ α µ = - + + (5) avec 3 1 4 ln10 2 o A DH B N e A k T ρ π ε   =       , 2 2 A o DH B e N B k T ρ ε = et C DH = 0,055
Le recours à des modèles plus complexes de coefficients d'activités sont de plus en plus utilisés comme Pitzer ou de Chen (Zematis, 1986) et ont été envisagés dans cette étude. Leur utilisation soulève une difficulté pratique à savoir la connaissance de nombreux paramètres. Le modèle de Pitzer nécessite ainsi 4 paramètres pour obtenir le coefficient d'activité de chaque espèce aqueuse tandis que le modèle de Chen utilise 2 paramètres. Pour utiliser le modèle de Pitzer sans aucune hypothèse simplificatrice, il faudrait plus de 150 paramètres. Pour le modèle de Chen, les données nécessaires, bien que moins nombreuses, sont rarement accessibles dans la littérature. Il a donc été écarté de ce travail. Le modèle de Pitzer est basé sur le développement de la loi de Debye-Huckel à laquelle sont ajoutés des termes rendant compte de la dépendance des forces d'interaction à courte distance à la force ionique. Les paramètres d'interaction pour des espèces formées à partir de la complexation des métaux sont très difficiles à mesurer et donc inexistants. Le formalisme de la méthode de Pitzer est beaucoup plus lourd que celui de la loi de Debye-Huckel. Par souci de concision, seule l'équation donnant le coefficient d'activité moyen γ ± d'un sel M υ+ X υ- (Zematis, 1986) est rappelée ici :

( ) 1.5 2 2 2 ln 2 z z f m B m C γ γ γ υ υ υ υ γ υ υ + - + - + - ± ± ±     = + +           (6) avec (
)

ln10 2 ln 1 3 1 DH f A b b b γ µ µ µ   = - + +   +     (7) où b=1.2 ( ) ( ) { } ( ) ( ) { } 2 1 0 1 1 1 2 1 2 2 2 2 2 2 2 2 2 1 1 0.5 exp 2 1 1 0.5 exp B γ β β α µ α µ α µ α β α µ α µ α µ α ± = + -+ - - + -+ - - (8) 3 2 C C γ φ ± = (9)
où β 0 , β 1 , β 2 et C Φ sont les paramètres de Pitzer obtenus à partir du lissage de données expérimentales.

Les résultats de la planification d'expériences sur les actions des divers coefficients d'activité ont montré qu'il n'est pas nécessaire de prendre en compte l'ensemble des coefficients d'activité : par conséquent, il n'est pas utile de mettre en place le modèle de Pitzer pour tous les coefficients d'activité puisque seuls, ceux relatifs à H + , HPO 4 2-, Ca 2+ et H 2 PO4 -ont une influence sur les courbes d'équilibre. Par ailleurs le modèle de Debye-Huckel a permis de façon préliminaire de calculer les concentrations des diverses espèces ioniques présentes dans le milieu. Ainsi, les résultats montrent que le système est constitué principalement de 6 espèces ioniques K + , Cl -, HPO 4 2-, Ca 2+ et H 2 PO4 -, les concentrations des autres espèces ioniques étant cent fois plus faibles. Ceci nous conduit à les négliger totalement et à réduire le nombre de paramètres lors de la mise en place du modèle de Pitzer. Les résultats obtenus avec les différents modèles de coefficients d'activité sont représentés sur la figure 3. Ils montrent que le choix du modèle de coefficients n'a pas un impact majeur en terme de taux de conversion en fonction du pH. Cet élément est important dans la mesure où le modèle de Pitzer nécessite un temps de calcul relativement important (le double de celui de Debye-Huckel) et une initialisation délicate (en général, elle est effectuée à partir des résultats obtenus par le modèle de Debye-Huckel). Finalement, le modèle de Debye-Huckel paraît représenter les phénomènes de façon suffisamment précise et pourrait servir de base pour l'optimisation du procédé.

CONCLUSIONS ET PERSPECTIVES

Nous avons identifié les constantes d'équilibre du système Ca-PO 4 -H 2 O tant en phase liquide qu'en phase solide qui agissent de façon notable sur la précipitation. Leur impact est d'autant plus important que l'on travaille à un pH inférieur à 8. Par ailleurs, nous avons quantifié l'importance relative des coefficients d'activités et mis en évidence que seuls les activités relatives à 4 espèces ioniques sont à prendre en compte. Ces deux résultats suggèrent deux voies importantes pour la suite de ce travail consacrée à l'optimisation du fonctionnement du lit fluidisé utilisé pour la précipitation des phosphates de calcium :

-

  Le point d'inflexion sur les courbes de la figure 2 correspond au changement de l'espèce qui précipite. Pour les pH inférieurs au point d'inflexion, les espèces ACP et DCPD précipitent tandis que pour les pH supérieurs seul ACP précipite. Selon le pH utilisé, on constate qu'une incertitude sur la valeur des constantes d'équilibres a un impact majeur sur le taux de conversion. Ainsi les écarts les plus importants sont observés pour des valeurs de pH inférieures à 8.

Figure 2 .

 2 Figure 2. Influence des constantes d'équilibre sur la précipitation

Figure 3 .

 3 Figure 3. Influence des coefficients d'activités sur la précipitation

  D'une part, il faut identifier les jeu des constantes d'équilibre qui traduit de façon la plus représentative les phénomènes observés et opérer si possible à pH élevé. -D'autre part, le modèle de Debye-Huckel semble suffisant pour représenter les coefficients d'activités par rapport au modèle plus sophistiqué de Pitzer qui nécessite un temps de calcul plus lourds. Cet élément est à prendre en compte en vue du couplage du modèle avec une méthode d'optimisation. REFERENCES ZEMATIS J.F., CLARK D.M., MARSHALL R., SCRIVNER N.C. Handbook of aqueous electrolytes thermodynamics : Theory and application, Editions DIPPR, 36-56, (1986). VIEILLARD P., Géochime des phosphates. Etude thermodynamique. Application à la génèse et à l'altération des apatites, Thèse, Université Paul Sabatier Toulouse, (1978). Seckler, M.M. Seckler, M.M., Calcium phosphate precipitation in a fluidized bed, Ph.D. Thesis, Technische Universiteit Delft, 1994 constantes d'équilibre incertitude de -20% sur les coefficients d'activité

Fin de la précipitation Précipitation du DCPD Bilan sur le phosphate et le calcium Electroneutralité 0 DCPD β =

  Organigramme du modèle de précipitation des phosphates de calcium Après précipitation de la phase ACP, la sursaturation de la phase DCPD est calculée. Si cette dernière est supérieure à zéro, cela signifie que potentiellement, la phase DCPD peut précipiter. Les notations suivantes (cf. Tableau1) ont été adoptées pour les constantes d'équilibre du système Ca-PO 4 -H 2 O : K

	l'influence du pH sur le taux de conversion. L'alimentation en calcium se fait par l'intermédiaire de CaCl 2 : par conséquent, nous avons pris une concentration en Cl -égale au double de celle de Ca
	Conditions initiales
	[KOK], [CaCl 2 ], [H 3 PO 4 ]
	Précipitation de l'ACP
	Résolution des équations de :
	Bilan massique sur le phosphate et le calcium
	Electroneutralité 0 ACP β =
	Calcul de la sursaturation de DCPD
	Non	β DCPD >0	Oui
			Calcul de la sursaturation de l'ACP
			Fin de la précipitation
	Les équations du modèle sont finalement résolues selon le schéma proposé dans la figure 1.
	Dans le modèle, toutes les relations d'équilibres de K 1 à K w sont transformées de façon à ce que chaque espèce soit exprimée uniquement en fonction des concentrations [Ca 2+ ], [PO 4 3-] et [H + ]. Le
	système à résoudre lors de la précipitation de l'ACP comporte 12 équations non linéaires avec 12
	inconnues (11 concentrations d'espèces aqueuses et un taux de conversion). Pour simplifier la
	résolution numérique en raison d'une initialisation délicate, le nombre d'équations est ramené aux quatre équations de bilans et les inconnues du système sont uniquement les concentrations [Ca 2+ ], [PO 4 3-], [H + ] et le taux de conversion en phosphate. Ce système d'équations est résolu par une
	méthode classique de Newton-Raphson pour différentes concentrations KOH afin d'analyser

2+

initiale.

de la planification d'expériences : Constantes d'équilibre relatives aux phases solide et liquide

  . Valeurs des constantes d'équilibres de la phase liquide relevées dans la littérature Les valeurs de pKmin et de pKmax correspondent aux valeurs minimales et maximales trouvées dans la littérature tandis que la valeur pK utilisée correspond aux valeurs calculées par Vieillard (1978) au moyen de concepts thermodynamiques. Ce choix a été effectué dans la perspective d'une étude sur l'influence de la température sur les courbes d'équilibre.Les différences entre les valeurs des constantes d'équilibre de la phase solide trouvées dans la littérature étant notables, les valeurs de pKs ACP et de pKs DCPD ont été calculées à partir des expériences de Seckler (1994). Une incertitude de 2% sur la valeur du pKs ACP et de 5% sur la valeur du pKs DCPD a été prise en compte soit,

	Log 10 constante	pK min à 25°C	pK utilisé à 25°C	pK max à 25°C	Ecart (%)
	d'équilibre
	K1	2.03	2.147	2.2	8.37
	K2	7.0	7.205	7.207	2.95
	K3	12.00	12.343	12.343	2.85
	K4	1.4	1.408	1.44	2.85
	K5	2.5	2.739	2.739	9.56
	K6	6.459	6.459	6.5	0.63
	K7	1.23	1.23	1.23	0.00
	Tableau 2. Niveau 2


The strain field of light jets The strain field of light jets L. Joly L. Joly 1 1 , J. Suarez , J. Suarez 1 1 , P. Chassaing , P. Chassaing 1,2 1,2

1 ENSICA -Département de Mécanique des Fluides -Toulouse, France. • The baroclinic torque is Reynolds-number and density-ratio sensitive,

• The strain fields of light jets present a folded layered structure and higher pick values,

• According to the theoretical prediction by CK00 the growth rate of transient « pressureless modes » is bounded above by ,

• Growth rates of transient modes may be up to twice more than their homogeneous counterpart before and during the pairing process,

• This may ensure the rapid growth of streamwise-aligned structures yielding side ejections of light fluid,

• The baroclinic torque mechanism produces thinner braids indicating higher spanwise or azimuthal wavenumbers,

• A need for an (ongoing) direct stability approach « à la » Klaassen and Peltier JFM 227 (1991).

2

  Institut de Mécanique des Fluides de Toulouse -Toulouse, France. Outline -Specific features of light jets -Basics of the 2D development of KH billow in inhomogeneous jets and mixing layers -The strain field of a light jet with pairing -Potential growth rate of oblique modes according to the Caufield and Kerswell's criterion

Fisrt•

  observation by Monkewitz & Bechert (Album of Fluid Motion -1988) Side ejections on helium jets from Hermouche, Imft (1996) Question: why vortex induction naturally produces side jets in round light jets only ? t=0 t=2 t=4 t=6 t=8 t=10 t=12 The baroclinic torque in the inhomogeneous mixing layer Kelvin-Helmoltz Instability in the temporal mixing layer s=3, Re = Uδ /ν = 1500 Vorticity Baroclinic torque 5 The order of magnitude of the baroclinic torque Temporal mixing layer s=3, Re = Ud/n = 500,1500,3000 Increasing Re • Small changes in the growth rate • Vorticity redistribution within one KH wavelength • Positive/negative generated circulation increases with the Reynolds/Schmidt numbers Diffusive damping The structure of the 2D light Jet (with pairing) Temporal jet s=1/3, Re = UD/ν Analysis of the strain field The shear and the strain rate across a contour line The model of the braid region according to Caulfield and Kerswell (PoF vol 12, 2000) γ (stretching) -Blue : negative γ (compression) -Shaded : above same positive level Passive scalar Light jet The strain field according to the Caulfield and Kerswell's analysis Temporal jet s=1/3, Re = UD/ν = 2500 Passive scalar Light jet t = 4,5The strain field according to the Caulfield and Kerswell's analysis Temporal jet s=1/3, Re = UD/ν The structure of the infinite Froude number inhomogeneous jet is biased by the baroclinic torque,

13

  Streamwise baroclinic torque on a pair of counter-rotating vortices vorticity Baroclinic torque vorticity • Offset of vortex centers toward the light side Potenttial • Potential amplification by the strain field via the densitygradient Spatial structure of the baroclinic torque on the pattern of counterrotating vortices Counter-rotating vortices on a density-gradient layer

The nature of the baroclinic torque in infinite Froude number jets vanishes in barotropic flows signature of inertia difference


Sofiane Hocine

Sébastien Elgue

Ludovic Montastruc

Catherine Azzaro-Pantel

Michel Cabassud

Luc Pibouleau

Serge Domenech

A hybrid optimization method (GA/SQP) for methyl acetate production process by reactive batch distillation

The initialization problem of classical optimization techniques (for instance, Sequential Quadratic Programming) is generally crucial for the robustness of the code. This paper presents a contribution to handle this kind of problems by implementation of a hybrid optimization technique based on the successive use of a Genetic Algorithm (GA) and of a classical Sequential Quadratic Programming (SQP) method: the GA is used to perform a first search in the solution space and to locate the neighbourhood of the solution. Then, the SQP method is employed to refine the solution. The optimization of methyl acetate production process by reactive batch distillation serves as an illustration to validate the approach.

Introduction

Batch processes are widely carried out for the production of fine chemical or pharmaceuticals with generally many successive reaction and separation steps involved. Synthesis optimization is often limited to the determination of the optimal operating conditions of each step separately. Therefore, such an approach does not lead to the optimal conditions for global synthesis. For example, conversion optimization of a reaction for which separation between the desired product and the by-products is more difficult than between the reactants will involve an important operating cost, due to further difficulties in the separation scheme. Thus, the simultaneous integration of all the process steps in a single optimization approach is of major importance [START_REF] Wajge | RBDOPT: a general purpose object-oriented module for distributed campaign optimization of reaction batch distillation[END_REF]Reklaitis, 1999, Elgue et al., 2002). The purpose of the study is to optimize the reactive distillation process of methyl acetate production. A simulation model was previously developed and validated for describing this process [START_REF] Elgue | Optimization of a Methyl Acetate Production Process by Reactive Batch Distillation[END_REF]. The model was further used in connection with an optimization loop involving an objective function combining two criteria, i.e. operating time and conversion rate with various weights. This study has shown that the initialization problem of classical optimization techniques which may be used (for instance, Sequential Quadratic Programming) is crucial for the robustness of the code. For this purpose, an optimization strategy combining a twostage approach, i.e. a genetic algorithm for initialization purpose and identification of the search zone followed by an SQP method to refine the solution is proposed in this paper: the GA is used to perform a first search in the solution space and to locate the neighbourhood of the solution. Then, the SQP method is employed to refine the solution. In this paper, the basic principles of the model development will be briefly recalled. The main results of the classical optimization approach by SQP will be presented. Then, the basic operations involved in the design of the GA developed in this study (encoding with binary representation of real values, evaluation function, adaptive plan) are presented. The performances of the method are analyzed and compared with the only SQP procedure.

Model presentation

A general model taking into account the various configurations of reaction separation process was developed [START_REF] Elgue | Optimization of a Methyl Acetate Production Process by Reactive Batch Distillation[END_REF]. The mathematical model involves a rigorous description of a batch reactor connected with an overhead distillation column and is formulated as a set of differential and algebraic equations (DAE system) around each plate. Ordinary differential equations (ODE) are due to energy balances, total and component mass balances. Algebraic equations (AE) consist of vapour-liquid equilibrium relationships, summation equations and physical property estimations. To reduce the complexity of the model, the following typical assumptions have been made on each tray: perfect mixing between vapour bubbles and liquid, equilibrium between liquid and vapour bubbles, introduction of Murphree efficiency, negligible vapour holdup compared to the liquid holdup and constant volume of liquid holdup. The mathematical model is connected to Prophy physical property estimation system (enthalpy model, equilibrium constant, hydrodynamic relationship and bubble point temperature) with associated data bank. DISCo [START_REF] Sargousse | DISCo: un nouvel environnement de simulation orienté objet[END_REF], a general solver of DAE systems based on the Gear method is used to solve the mathematical model. The model is presented in detail in [START_REF] Elgue | Optimization of a Methyl Acetate Production Process by Reactive Batch Distillation[END_REF] for methyl acetate production by reactive distillation according to the following reaction:

O (1)
Production of methyl acetate is carried out by addition of acetic acid to methanol with sulfuric acid as homogeneous catalyst, as it can be shown on the reaction scheme (1). A good agreement was obtained between simulations and experimental runs, thus validating the model assumptions and confirming interest in process optimization. For this purpose, let us recall that operating time, reaction yield, energy consumption, safety and environmental constraints are generally the key points of reactive batch distillation processes. In the concerned process, the column is always operating with a total reflux policy during reaction phases and with a constant reflux policy during distillation. Thus, energy consumption (heat provided by the boiler, condenser cooling water) is only a function of the operating time. In this way, process optimization only involves two criteria, i.e. operating time and conversion of reactants. Therefore, a function combining operating time and conversion rate with various weights has been proposed in order to optimize the production process. Hence, a mono-objective optimization method based on a Sequential Quadratic Programming method (SQP) was used. All the different processes, reaction-distillation or their combination have the same dynamic structure (see Figure 1). In fact, according to the total reflux end time (t R ), before or after the time of reaction equilibrium (t eq ), the process is respectively a coupling or a reaction-distillation process. Thus, the total reflux end time appears to be a main variable of the optimization problems, as well as reflux ratio (R) and operating time (t op ). The total reflux end time is a very sensitive parameter of reactive batch distillation. On the one hand, if total reflux policy ends too early, conversion will be reduced owing to reactant distillation. On the other hand, a too long total reflux policy will lead to

H 2 O CH 3 OH C OH C O CH 3 CH 3 O + + H SO 2 CH 3 Water Methyl acetate Methanol Acetic acid (118°C) (100°C) (57°C) (65°C)
additional operating time and, consequently, to waste of time. An adequate value for total reflux time is very difficult to estimate experimentally and optimization represents a very challenging objective for its automatic determination.

Distillation Reaction

Figure 1: Dynamic structure of batch reactive distillation process

Optimization with SQP package

Two kinds of optimization problems were investigated (see Table 1). In the former, the minimum operating time necessary to obtain the desired reactant conversion is determined. In the latter, the objective functions involve a combination of operating time and conversion with various weights.

Model equations

Several optimization problems were solved with SQP method (SQP package from IMSL library) [START_REF] Schittkowski | NLPQL: A FORTRAN subroutine solving constrained non-linear programming problems[END_REF]. The main results are presented in Table 2. The two first solutions present the optimal operating conditions for respectively 93 percent and 95 percent of acetic acid conversion (C onv ). The three other solutions are relative to operating conditions leading to optimal criteria, for which the conversion weight decreases respectively. Results of the first set of optimization problems show that a significant total reflux time (more than 15 min) is required for a high conversion of reactants. The second set show that if conversion is favoured, total reflux time is almost invariable (around 23 min) and only further operating time allows reaching a higher conversion.

Optimization with a GA/SQP hybrid optimization method

A typical feature of complex optimization problems is that the generation of starting points that are used in the search with conventional optimization methods (for example SQP) is very important to guarantee the success of the optimization procedure.

The preliminary study has shown that the only use of an SQP method is very sensitive to the choice of the initial guess and often leads to a failure. Consequently, this work is motivated by the development of a technique for the automatic generation of good starting points.

For this purpose, a hybrid optimization method is proposed in this paper: a Genetic Algorithm (GA) is used to perform a first search in the space and to locate the neighborhood of the solution. Then, the derivative driven optimization tool (SQP) is used to refine the solution. Let us recall that Genetic Algorithms [START_REF] Goldberg | Genetic algorithms in search, optimization and machine learning[END_REF] differ from most classical optimization methods since no assumption about the problem space is required and yet produce a global search. GAs use a guided random search in which many different solutions to a problem are investigated and refined simultaneously to identify near-optimum solutions. A major interest of such methods is that they lead to reasonable solutions even with a poor initial guess. GA implementation requires the definition of parameters: population generation mode, population size, (i.e., the number of individuals forming a population which must be sufficiently large to create sufficient diversity, in order to cover the possible solution space) crossover probability, mutation rate, survival, crossover and mutation mechanisms. In this paper, the initial population has been selected randomly. Nevertheless, constraint respect at a precision of 5.10 -6 is required. Each potential solution, i.e., a real value, has been coded as an haploid chromosome with the so-called "weight box". The weight box consists to encode each digit of real in 4 bits with respectively the following weights 1, 2, 3, 3 with a precision of 10 -20 . By lack of place (see [START_REF] Montastruc | A global optimization technique for Solid-Liquid Equilibrium: Application to calcium phosphate precipitation[END_REF], the choice of weight is not justified here but has been achieved after a sensitivity analysis. For example, the real 0.18 is encoded 1000 0111, the zero digit and the point are not encoded. Each chromosome represents the 3 unknowns. Since a traditional GA is based on a function maximization, the evaluation function has been adapted to also consider a minimization case. More precisely, the treated problem involves both operating time minimization and conversion rate maximization. Moreover, as initially proposed by [START_REF] Goldberg | Genetic algorithms in search, optimization and machine learning[END_REF] and followed by several authors [START_REF] Costa | Evolutionary algorithms approach to the solution of mixed integer non-linear programming problems[END_REF], the penalty function method has been used to take into account the constraints in the optimization procedure. Two cases have thus been taken into account:

Evaluation function = (t max -t op ) + P X conv , for t op > t R



Evaluation function= (t max -t op )+ P X conv -P 2 (t R -t op ) 2 , for t op < t R In these expressions, t R is the time from which a change in reflux ratio occurs. The constraint t op -t R >0 is thus taken into account. If this constraint is not satisfied during individual selection phase, another penalty term is introduced (coefficient P 2 ). Consequently, the individuals will have lower chances to be selected and to contribute to the next generation.

Results

GA parameter setting

The definition domain of the variables which has been taken into account is as follows: 0.6 < r < 0.9; 1.5 < R < 9 with R= r / 1-r; 0 < t op < 130 min; 5 min < t R < 50 min A preliminary study with the GA leads to the following parameter set, which has been arbitrarily fixed, following classical guidelines: Generation number = 40; Population number= 40; Crossover rate = 75%;Mutation rate =25 %.

GA+SQP typical results

Four different values of the penalty coefficient P were considered : P=100;1000, 1700 and 2500. By lack of place, only the results relative to P=1700 are presented (See Table 3). They show that the results obtained by SQP are very similar to those obtained by coupling GA and SQP, which shows that the GA is relatively efficient. It has been pointed out that an increase in penalty coefficient P may lead to significant variations in results, which mean that several solutions may exist. Another investigation was performed coupling GA and SQP at different stages of the GA evolution, i.e., for the best individuals obtained respectively at the 4 th , 10 th , 12 th and 20 th generations (NG). The results are proposed in Table 4. It can be seen that SQP converges to the same solution for NG=10 and NG=20 whereas initialization is not the same. Yet, different results are obtained for NG=4 and NG=12, which means that the number of potential solutions may be high. It must also be noted that the computational time for GA is significantly higher, which confirms that the use of GA at its early stages of evolution is particularly interesting for determination of initialization values, which are generally difficult to find by trial-and-error procedures.

Conclusions

In this paper, a hybrid optimization tool combining a Genetic Algorithm and an SQP method has been developed and tested on optimization of methyl acetate production process by reactive batch distillation. This two-level strategy leads to very efficient search: the GA has been able to provide good starting points for the subsequent SQP method (even if the user has provided rough starting points), thus favouring the local search. The generic feature of the optimization strategy presented in this study finds a widespread application for hybrid optimization problems.

Table 1 :

 1 Optimization problems

		Objective function	Variables	Constraints
	1.	Operating time	Operating time	Model equations
			Total reflux end time	Acid Conversion
			Reflux ratio
	2.	Operating time,	Operating time
		acetic acid conversion	Total reflux end time
			Reflux ratio

Table 2 :

 2 Results of the optimization problems (with SQP)

	Objective	Constraint	Operating	Total reflux	Reflux	Acetic acid
	function		time (min)	end (min)	ratio	conversion
	t op	C onv  93 %	80	15	1.8	93.0 %
	t op	C onv  95 %	118	20	3.0	95.0 %
	t op + (1-C onv ) x 33	/	141	25	2.3	95.9 %
	t op + (1-C onv ) x 17	/	105	23	2.5	94.5 %
	t op + (1-C onv ) x	/	48	9	1.0	85.8 %
	1.7

Table 3 :

 3 Results of GA and GA+SQP for P=1700

	Run		Reflux ratio Operating time Total reflux time Conversion rate (%)
	1	GA	2,02	120 min	12 min	95,00
		GA+SQP	2,58	90 min	10 min	93,66
	2	GA	1,68	120 min	24 min	95,07
		GA+SQP	1,60	120 min	24 min	94,94
	3	GA	2,12	109 min	17 min	94,65
		GA+SQP	2,14	107 min	18 min	94,65
	4	GA	2,48	106 min	24 min	94,56
		GA+SQP	2,75	107 min	19 min	94,65

Table 4 :

 4 Results of GA and GA+SQP for different stages of A evolution (NG, number of generation)

	NG=4	Reflux	Operating	Total reflux	Conversio	Evaluation	Computation
		rate	time	time	n rate (%)	function	time
	GA	3,46	126 min	15 min	95,18	-676 041,5	20 min
	GA+SQP	3,72	126 min	15 min	95,21	-676 254,4	23 min
	NG=10
	GA	2,38	100 min	15 min	94,23	-676 819,0	45 min
	GA+SQP	2,68	100 min	15 min	94,34	-677 650,1	52 min
	NG=12
	GA	2,38	91 min	15 min	94,11	-677 122,6	55 min
	GA+SQP	2,52	91 min	15 min	94,15	-677 364,4	59 min
	NG=20
	GA	2,56	100 min	15 min	94,32	-677 500,1	90 min
	GA+SQP	2,68	100 min	15 min	94,34	-677 650,8	93 min