Vincent Repellin

René Guyonnet
email: guyonnet@emse.fr

Evaluation of heat treated wood swelling by differential scanning calorimetry in relation with chemical composition

Keywords: differential scanning calorimetry, fiber saturation point, heat treatment, swelling, wood

Retification® is a heat treatment which decreases the swelling of wood and increases its resistance to fungal attack. In this study, differential scanning calorimetry (DSC) was applied in order to determine the fiber saturation point (FSP) of natural and retified® wood. FSP values were used to determine the total swelling of natural and heat-treated wood. The DSC method was compared to the volumetric shrinkage approach. The influence of the heat treatment temperature and duration on the swelling of wood was investigated. Relationships between chemical changes and the reduction of swelling were analysed thoroughly.

The equivalence of the DSC method and the volumetric shrinkage method is shown. FSP in association with anhydrous density is a good indicator for the evaluation of the overall swelling of heat-treated wood. Reduction of wood swelling with increasing temperature and duration of thermal treatment is often attributed to hemicelluloses destruction. This study shows that the reduction of beech wood swelling can not only be attributed to the disappearing of adsorption sites that goes with the hemicelluloses destruction. It is suggested that other phenomena such as structural modifications and chemical changes of lignin also play an important part.

Introduction

Heat treatments of wood become widely used for wood protection against fungal attack, and for reduction of wood swelling [START_REF] Kamden | Durability of heat-treated wood[END_REF][START_REF] Mouras | Propriétés physiques de bois peu durables soumis à un traitement de pyrolyse ménagée[END_REF][START_REF] Niemz | Effect of high temperature on the change in color, dimensional stability and mechanical properties of spruce wood[END_REF]. Retification® is a mild pyrolysis of wood that takes place under nitrogen atmosphere (Guyonnet and Bourgois 1986). Oxidation phenomena and cellulose hydrolyse due to the presence of oxygen and water steam are then reduced. Coloration and chemical modifications take place in the mass of wood, and allow further machining of the material [START_REF] Bourgois | La mesure de la couleur : une méthode d'étude et d'optimisation des transformations chimiques du bois thermolysé[END_REF]. Moreover, retification® avoids the use of chemical compounds that may cause pollutions. The industrial development of this type of heat treatment requires rapid methods for the control of end use properties. This study focussed on the swelling of heat-treated wood.

Wood cell wall is thought to be a composite material made of cellulose microfibrils embedded in a water reactive matrix of hemicelluloses and lignin. The ability of the matrix to adsorb water is thus of critical importance in the swelling of wood [START_REF] Barber | A theoretical model of shrinking wood[END_REF][START_REF] Cave | A theory of the shrinkage of wood[END_REF]. It is known that heat treatments reduce wood hygroscopicity by modifying this water reactive matrix [START_REF] Obataya | Hygroscopicity of heat-treated wood I -effects of after treatments on the hygroscopicity of heat-treated wood[END_REF]Obataya and Tomita 2002) in different ways.

In order to quantify the overall swelling, the fiber saturation point (FSP) was used. The FSP of wood may be defined as the moisture content (MC) at which cell walls are saturated with water and cell cavities are free of liquid water. At the FSP all the adsorption sites of water are busy. Any increase in MC above the FSP results in negligible dimensional increase in the material, whereas any decrease in MC below the FSP leads to a shrinkage of the material. As a consequence, the FSP is traditionally measured by the volumetric shrinkage method. This method requires slow drying, and takes much time to meet the needs of industry.

Nevertheless, the state of water in wood can be quantitatively and qualitatively measured by other methods. For instance 1 H MNR was used [START_REF] Labbé | Moisture content and extractive materials in maritime pine wood by low Field 1 H MNR[END_REF] to determine the water state and distribution in maritime pine wood. A somewhat more simple method, using differential scanning calorimetry (DSC), can be employed to determine the state of water and the FSP of various cellulosic materials [START_REF] Nakamura | Studies on Bound Water of Cellulose by differential scanning calitymetry[END_REF][START_REF] Kabeya | Effect of grinding on bound water in cellulosic materials[END_REF]. The FSP of green wood has already been successfully measured by DSC [START_REF] Simpson | Determination of the fibre saturation point in whole wood using DSC[END_REF]. In this study, FSP of heat-treated wood was determined by DSC.

The FSP values measured by DSC were compared to FSP values measured by volumetric shrinkage. The relevance of the FSP as an indicator of the total swelling of wood was then evaluated. The influence of the heat treatment temperature and duration on the swelling reduction of wood was examined in relation with changes of chemical composition.

Materials and Methods

Heat treatments and samples

Nine batches of common beech (fagus sylvatica), and four batches of maritime pine (pinus pinaster) were carried out in a pilot kiln. The corresponding values of the parameters (duration and temperature) are mentioned in Table 1. Each batch consisted in 16 boards of size 25x85x850 mm. Before treatment, the boards had been equilibrated in a climatic chamber at 20°C and 65% relative humidity during three months.

One board was randomly chosen among each treated batch. In this board, a 2 cm edge square cross section beam was cut across the natural directions of wood. In each beam 10 cubes and a thin lamella of thickness inferior to 1 mm were cut. The cubes were used for volumetric shrinkage measurement of FSP. Each thin lamella was cut again in parallelepipedic samples, that were used for FSP measurements by DSC. Their dry weight ranged from 2 mg to 5 mg. A part of the retified® wood of each batch was ground and sieved for the determination of the chemical composition of natural and heat-treated beech wood.

FSP measurement by Volumetric shrinkage

The cubic samples were immersed in water until constant weight and volume were reached. Then, they were slowly dried. Their weight and length along the three directions were measured. In this way, the volumetric swellings of the cubes versus their moisture content were plotted (Fig. 1). The FSP was then given by the intersection of the two straight lines obtained by a linear regression of the four first points and the four end points of these curves.

FSP measurement by Differential Scanning Calorimetry Water in wood:

The state of water in cellulosic materials was accurately studied by DSC [START_REF] Nelson | The determination of moisture transition in cellulosic materials using differential scanning calorimetry[END_REF][START_REF] Nakamura | Studies on Bound Water of Cellulose by differential scanning calitymetry[END_REF][START_REF] Kabeya | Effect of grinding on bound water in cellulosic materials[END_REF]. In wood as well as in cellulosic materials, water can be present under three main distinct states: free water, non-freezing bound water and freezing bound water. Free water is held in the wood large pores at a liquid state. On the DSC curves, free water presents a melting endothermic peak at a temperature equivalent to pure water (0°C). Non-freezing bound water consists of water molecules adsorbed on the cell wall and presents no melting peak. The quantity of non-freezing bound water present in wood corresponds to the FSP. Freezing bound water is mainly located in small pores and its physical properties are affected by the vicinity of the cell walls. The endothermic 'melting' peak of freezing bound water can be observed at temperatures inferior to 0°C, according to the Gibbs-Thompson relationship.

DSC device and calibration

The device was a SETARAM 131 DSC. Samples were sealed in hermetic stainless steel crucibles. A cooling system using liquid nitrogen allowed to reach temperatures inferior to 0°C. Heat flow and Temperature have been calibrated with three standards: mercury, deionised water and gallium. Afterwards, experimental conditions were set as follows: increase in temperature from -40°C to +40°C, with a heating rate of 2°Cmin -1 . A preliminary study verified the linearity domain of melting enthalpy of pure deionised water droplets weighed from 0.8 mg to 11.0 mg.

Method for FSP measurements

The "direct method" of Simpson and Barton [START_REF] Simpson | Determination of the fibre saturation point in whole wood using DSC[END_REF] was applied to natural and retified® beech wood and maritime pine wood. At first, the samples were oven dried at 103°C and weighed in their anhydrous state (m 0 ). Then they were immersed in deionised water for 48 hours. The water saturated samples were weighed (m s ). Integration of the endothermic peak of free water gave its total heat of melting (ΔH obs , in J).

Free water is assumed to behave like pure water, its melting enthalpy (Δh f ) is taken equal to 334 Jg -1 . It is then possible to calculate the mass of non freezing bound water and the FSP, in percent, thanks to equation (1):

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Δ Δ - - = f obs 0 S 0 h H m m m 100 FSP
(1).

Mercury Porosimetry

A mercury porosimetry test was carried out using a Micromeritics Autopore IV. The sample was a 0.156 g cube of natural beech wood, cut across wood natural directions. A penetrometer for solid sample analysis was used. It had a bulb volume of 5 cm 3 and a stem volume of 0.392 cm 3 . The pressure was increased from 0.0037 MPa to 185.7661 MPa.

Chemical analysis of heat-treated beech wood

Extractives, lignin, and carbohydrates contents of natural and heat-treated beech wood were determined.

Extractives content was determined according to ASTM method D 1107-96. The initial sample is made of 2 g of wood ground and sieved between 212 µm and 560 µm. These samples were extracted in a soxhlet with a mix of ethanol and toluene. The weight lost during extraction gives the extractives content. The resulting extractive free wood was used to determine the Klason's lignin content according to the ASTM method D 1106-96. This method consists in an acid hydrolysis. The extractive free wood powder is firstly mixed with a 72% sulphuric acid solution. Secondly, water is added. And finally, this solution is filtered. The mass of the solid fraction gives the amount of Klason's lignin in the extractive free wood. The filtrate is kept for further analysis because it contains the soluble fraction of lignin and different simples sugars resulting from the wood polysaccharides hydrolysis. Soluble fraction of lignin is determined in a semi quantitative way by UV spectroscopy [START_REF] Singleton | Colorimetry of total phenolics with phosphomolybdicphosphotungstic acid reagent[END_REF][START_REF] Alberto | Influence des extractibles d'essences du Mozambique sur l'hydratation du ciment[END_REF]. The amount of soluble lignin is given in mg catechine equivalent per g of extractive free anhydrous wood. The amount of sugar resulting from the hydrolysis of wood polysaccharides is measured by high pressure liquid chromatography (HPLC). Glucose, xylose, mannose, cellobiose and glucuronic acid were identified. Their amount in the filtrate (for example the filtrate xylose content) was quantitatively determined in weight percent of extractive free anhydrous wood.

Results and discussion

Shape of the DSC melting peak of water in wood saturated samples.

The DSC curves of pure water had a unique and simply shaped endothermic peak. An endothermic peak matching to the melting of free water was also present for water saturated wood samples. This peak presented a shoulder almost for all specimens (Fig. 2) as it had already been noticed by [START_REF] Simpson | Determination of the fibre saturation point in whole wood using DSC[END_REF]. It is assumed here that this shoulder is due to freezing bound water located in wood pores. Thanks to the Gibbs-Thompson's law equation ( 2) it is possible to calculate the decrease of the water melting temperature present in pores:

( ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ Δ - - + = - = Δ f ls m 0 0 h r σ V exp 1 273 T T T T (2)
where T 0 is the melting temperature of pure water in Celsius degrees. T is the melting temperature of water held in pores of diameter 2r. V m is the molar volume of frozen water (19.6 m 3 mol -1 ). The value of the liquid/solid surface tension (σ ls ) is chosen equal to 20.4 mNm -1 which corresponds to the case of water melting in porous silica media (cited by [START_REF] Maloney | Hydration and Swelling of pulp fibers measured with differential scanning calorimetry[END_REF].

On a DSC curve, the melting temperature is represented by the onset temperature. Since the shoulder is located at the end left of the peak of free water (to the low temperatures), the onset temperature measured on water saturated wood samples corresponds to the melting temperature of freezing bound water. Consequently, the difference in onset temperatures between pure water and water in saturated wood samples should be linked to the size of the pores containing the freezing bound water. From the experimental point of view the difference in average onset temperatures of pure water and free water in beech wood is 0.06°C.

The porometry test showed three populations of pore diameters around 20 µm, 0.44 µm and 0.015 µm (Fig. 3). The corresponding temperature shifts were calculated according to the Gibbs-Thompson's law: 0.9 10 -3 °C, 0.04°C and 1.18°C respectively. As regards the first population (20 µm), the temperature shift is too little be detected. As the difference in average onset temperatures of pure water and water in wood is of 0.06°C, the second population (pore diameter around 0.44 µm and temperature shift of 0.04°C) can match to the shoulder of Fig. 2.

Moreover, pores of diameter inferior to 9 µm are attributed to fibres lumen (Trenard 1980).

Thus the shoulder can match to freezing bound water present in fiber lumen. The third pore population (0.015 µm) should result in a temperature shift of about 1.2°C. However, no corresponding peak was visible on the DSC curves. On the porometry curves obtained by

Trenard in 1980, one can notice that this population is present on 10 mm thick samples and not on 240 µm thick samples (thickness inferior to the fibre length). It is suggested that the distribution around 0.015 µm visible on Fig. 3 is linked to the sampling and does not correspond to an existing pore population. That would explain why no peak is visible at the corresponding temperature on the DSC curves.

Effect of the thermal treatment on wood swelling and precisions about FSP Obataya showed that dry heat treatments result in reversible and irreversible effect on wood hygroscopicity (Obataya et al. 2002). Both reversible and irreversible effects may affect wood swelling. It is thus necessary to indicate which effect the FSP accounts for.

Reversible effects are due to softening and relaxation of the cell wall matrix at high temperature. This phenomenon induces the tightening of the matrix structure during cooling.

The mobility of the polymer chains is then reduced. At low relative humidity, the tightened structure prevents water adsorption. However, if the heat-treated material is moistened at high relative humidity, polymer chains recover their mobility, and water adsorption occurs again.

Irreversible effects are due to the chemical modifications of the material. The sample preparation for DSC as well as for volumetric shrinkage measurement requires water saturation. By this way reversible effects are cancelled. Consequently, the FSP measurement accounts only for the irreversible effects induced by heat treatment.

Comparison of the two methods for FSP measurement and evaluation of overall swelling

Thanks to a suitable choice of parameters, a large range of FSP values could be reached: from 15% to 33%. Fig. 4 shows the linear regression between the FSP measured by volumetric shrinkage and the FSP measured by DSC. The reliability of the DSC method is demonstrated by a slope equal to 1.03 and a correlation coefficient equal to 0.93.

Swelling of wood (G in %) is defined by equation (3):

0 0 S V V V 100 G - × = (3)
where V 0 is the anhydrous volume, and V s the water saturated volume of the specimen.

Experimental swelling was measured on the cubic samples. By assuming that the swelling of wood corresponds to the adsorbed water a calculated swelling (G c ) can be assessed by equation ( 4) [START_REF] Skaar | Wood-water relations[END_REF]:

FSP d G anh C = (4).
Where d anh is the apparent anhydrous density of wood. This relation does not take into account the pore size variations during drying that was studied in 1962 by Boutelje [START_REF] Boutelje | On shrinkage and change in microscopic void volume during drying as calculated from measurements on microtome cross section of Swedish pine (pinus sylvestry L.)[END_REF]). On Fig. 5, the linear regression between G and G c exhibits a slope of 0.94 and a correlation coefficient of 0.92. The FSP associated with anhydrous density can be thus used as a good indicator for the evaluation of wood swelling.

Influence of temperature and duration on the FSP.

Table 2 shows the FSP values measured by DSC on natural and heat treated wood. An increase in retification® temperature leads to a decrease in FSP of common beech and Maritime pine. An increase in retification® duration for beech wood also leads to a decrease in wood FSP. The effect of a long time exposure at low retification® temperatures seems to be equivalent to short times at high temperature. The same average swelling reduction level could be reached with treatments E, I and J that correspond respectively to 5 min at 260°C, 180 min at 220°C and 600 min at 200°C. This observation is consistent with the idea that the reduction of wood hygroscopicity follows an Arrhenius relation as a function of the timetemperature combination, as mentioned by [START_REF] Skaar | Wood-water relations[END_REF].

Reduction of FSP in relation to chemical composition of beech wood.

Tables 2 and3 present the chemical composition of natural and heat treated beech wood. In this part, it is proposed to analyse the relationship between reduction of FSP and wood chemical changes due to thermal treatment.

Variations of extractives content with retification® parameters.

Extractives content increases with increasing temperature and length of treatment. However, we assume in the following discussion, that extractives content has negligible influence on the FSP value. Consequently, FSP decrease by heat treatment is discussed in function of the chemical composition of extractive free anhydrous wood.

Variation of lignin content with retification® parameters.

The soluble fraction of lignin is almost stable whatever parameter values (Table 2). The amount of Klason's lignin increases with increasing temperatures. At 220°C, the quantity of lignin increases significantly for treatments longer than 60 min. At 200°C, it decreases for treatments of 60 min and 180 min length, and increases again with a 600 min treatment.

More precisely, the Klason's lignin content of treatments B and J (5 min at 200°C and 600 min at 200°C respectively) are very close (24.3±0.1% and 24.6±0.2% respectively). In the same time the FSP is seriously lowered (30.6%±2.2% and 15.3%±1.8% respectively). This observation shows that lignin content has a minor influence on the FSP.

Variations of polysaccharides content with retification® parameters

Two types of sugar polymers are present in hardwood: hemicellulose (glucuronoxylans, glucomannans) and cellulose [START_REF] Sjöström | Wood chemistry, fundamentals and applications[END_REF]. During acid hydrolysis, these molecules are depolymerised into simple sugar dissolved in the filtrate. Consequently, cellobiose and glucose content in the filtrate match to the cellulose proportion in wood, as well as glucuronic acid and xylose filtrate content match to glucuronoxylane proportion in wood, and mannose filtrate content matches to mannanes proportion in wood.

The main variation (Table 3) regarding polysaccharides is a decrease of the xylose content in the filtrate with increasing temperatures and durations of retification®. Glucuronic acid presents the same tendency. This observation shows that glucuronoxylan proportion in heat treated beech wood is lower than in natural beech wood. Hemicelluloses are the most thermally sensitive constituents of wood. They are the first components to be destroyed, as it is reported by numerous authors [START_REF] Beall | Thermal degradation of wood components, a review of the literature[END_REF]DeGroot et al.1988;Avat 1993;Weilland 2000). Mannose content is low and almost stable whatever the parameter values.

Glucose and Cellobiose content increase slightly and monotonously with increasing parameter values. Cellulose is resistant in the temperature range of retification®. Consequently, its proportion in extractive free anhydrous wood should increase proportionally to the decrease of the other constituents.

Relationship between polysaccharides content and FSP

Both glucuronoxylan content and FSP decrease. It is often assumed that the destruction of hemicelluloses goes together with a vanishing of water adsorption sites (mainly hydroxyl groups) that induces a reduction of wood swelling. In order to verify if hemicelluloses destruction is the only reason for the FSP reduction, we propose to calculate a FSP value related only to the polysaccharides content of natural and retified® wood. The filtrate xylose content (FXC in %) and the filtrate glucose content (FGC in %) stand for the proportion of glucuronoxylane and cellulose in wood. We assume that they allow to calculate the FSP, thanks to equation ( 5): [START_REF] Bourgois | The color measurement : A fast method to study and to optimize the chemical transformations undergone in the thermically treated wood[END_REF].

( ) ( ) C gl M FGC gl n wa M xy M FXC xy n wa M FSP + + = (
Where n xy , (resp. n gl ,) is the number of water molecules that can be adsorbed by a xylose (resp. glucose) unit of glucuronoxylane (resp. cellulose). M wa , (resp. M xy , M gl ,) is the molar mass of water (resp. xylose, glucose). FXC and FGC values are issued from experimental data (Table 3), they are given in weight percent of extractive free anhydrous wood. The term C accounts for the contribution to the FSP of the other constituents. C is chosen constant because the calculated FSP is chosen to be only related to wood polysaccharides content.

Equation ( 5) should be understood as follow. For 100g of free extractive anhydrous wood, one have (FXC)g of xylose unit in the filtrate, i.e. Values of n xy and n gl are issued from literature [START_REF] Nakamura | Studies on Bound Water of Cellulose by differential scanning calitymetry[END_REF][START_REF] Kabeya | Effect of grinding on bound water in cellulosic materials[END_REF]. It was shown that 4 molecules of water can adsorb on one xylose unit of larch xylan and that 2 water molecules can adsorb on one glucose unit of wood cellulose of cristallinity 44%.

Once the FSP value is calculated thanks to FXC and FGC, it is possible to compare the experimental and calculated reduction of FSP for each treatment. We define FSP variations from natural to retified® wood by:

(6). wood) ted (heat trea FSP - wood) l FSP(natura FSP = Δ
Assuming that hemicelluloses destruction is the only reason for the FSP reduction, the calculated ΔFSP should be equal to the experimental ΔFSP. Fig. 6 compares experimental and calculated ΔFSP values. The first observation is that experimental and calculated ΔFSP show the same tendency in spite of a rather large scattering. It is consequently expected that the FSP value depends on the hemicellulose content. Nevertheless the slope of the linear regression is 3.3, which means that the experimental variations of the FPS are to 3.3 times higher than the expected calculated variations. We admit that the number of water molecules that can adsorb on cellulose and glucuronoxylane in natural and retified® beech wood is close to the values found by Kabeya and Nakamura. As a consequence the decrease of wood swelling can not be explained simply in terms glucuronoxylane destruction. Some other reasons for the reduction of wood swelling have to be found.

Modifications of the lignin carbohydrate complex by heat treatments

It has been shown that some cross linking and condensation phenomena take place in lignin during retification® [START_REF] Weiland | Etude physico-chimique du traitement thermique du bois. Optimisation de paramètres du procédé de rétification®[END_REF]. As a consequence, the ability of lignin to adsorb water may be decreased. Though the lignin content does not vary significantly (treatments B, F, H, J), it is not excluded that chemical modifications of lignin induce significant decrease in the FSP. In addition, lignin is closely linked to the other constituents, and forms the lignin carbohydrate complex (LCC). Working on milled wood lignin, Obataya (Obataya et al. 2002) suggested that the modification of lignin and of the LCC during the heat treatment is mainly responsible for the irreversible reduction of wood hygroscopicity.

In this study, Klason's lignin from natural and heat treated wood (treatment E) were equilibrated at 70% relative humidity and 20°C after drying at room temperature. The MCs were low due to the drastic procedure used for the extraction of Klason's lignin destroying the LCC. Lignin from natural wood presented a MC of 7.05% and lignin from treated wood had a MC of 3.72%. Reduction of lignin hygroscopicity during heat treatment should thus take part in the reduction of wood swelling.

Moreover, lignin and hemicelluloses are closely linked to each other in cell wall. The lignin and hemicellulosic matrix of wood cell wall is generally thought to be isotropic and amorphous. Some authors however [START_REF] Vincent | From cellulose to cell[END_REF][START_REF] References Yan | Direct observation of the main cell wall components of straw by AFM[END_REF] report that hemicelluloses are regularly distributed in the cell wall structure: hemicelluloses are located in a lattice between the cellulose fibrils. The remove of hemicelluloses should thus modify this particular structure, and the ability of the wood cell wall to adsorb water. It would explain why in our study the FSP is so dependent on the FXC. Thus, structural modifications due to hemicelluloses destruction could be also involved in the reduction of wood swelling by heat treatment.

Conclusion

This study deals with the use of DSC for evaluation of the swelling of heat-treated wood. The method is based on the measurement of the FSP. The correlation of the DSC method with the volumetric shrinkage method is shown with linear regression. Moreover, the FSP associated with anhydrous density is a good indicator for the evaluation of the overall swelling of heattreated wood. This study shows that the reduction of wood swelling of beech wood can not be explained only in terms of the disappearing of hemicelluloses adsorption sites that goes with the hemicelluloses destruction. It is suggested that other phenomena such as structural modifications and chemical changes of lignin also play an important part. Extractive content in % of anhydrous wood (ASTM D1107-96), Klason's lignin content in % of extracted anhydrous wood (ASTM D1106-96), and soluble fraction of lignin in mg equivalent cathechine per g of extracted anhydrous wood (U.V. spectroscopy). The values presented are averages on four samples and the standard deviations are presented in brackets.

Table 3. Chemical composition of natural and heat-treated beech wood: carbohydrates and glucuronic acid content in weight percent of extracted anhydrous wood (HPLC).

  extractive free anhydrous wood.

Fig. 1 .FSPFig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .

 123456 Fig. 1. Measurement of the FSP by the volumetric shrinkage method ; the FSP is the abscissa

Fig. 1 .

 1 Fig. 1. Measurement of the FSP by the volumetric shrinkage method ; the FSP is the abscissa

Fig. 2 .

 2 Fig. 2. DSC melting peaks of pure deionised water and free water present in saturated wood

Fig. 3 .

 3 Fig. 3. Pore size repartition of natural beech wood (carried out with mercury porosimetry), on

Fig. 4 .

 4 Fig. 4. Comparison of two different methods for the measurement of the FSP of heat-treated

Fig. 5 .

 5 Fig. 5. Comparison of experimental swelling and swelling calculated from FSP and

Fig. 6 .

 6 Fig. 6. Comparison of the experimental and calculated variations of the fibre saturation point

Table 1 .

 1 Parameters of the heat treatments of retification® (time and temperature).

	Treatment number	species	Temperature (°C)	Time (min)
	A	Common beech	Untreated	Untreated
	B	Common beech	200	5
	C	Common beech	220	5
	D	Common beech	240	5
	E	Common beech	260	5
	F	Common beech	200	60
	G	Common beech	220	60
	H	Common beech	200	180
	I	Common beech	220	180
	J	Common beech	200	600
	K	Maritime Pine	Untreated	Untreated
	L	Maritime Pine	180	5
	M	Maritime Pine	230	5
	N	Maritime Pine	240	5
	O	Maritime Pine	260

Table 2 .

 2 FSP in %H 2 O and chemical composition of natural and heat-treated wood:

	Extractive content in % of anhydrous wood (ASTM D1107-96), klason's lignin content in %
	of extracted anhydrous wood (ASTM D1106-96), and soluble fraction of lignin in mg
	equivalent cathechine per g of extracted anhydrous wood (U.V. spectroscopy). The values
	presented are averages on four samples and the standard deviations are presented in brackets.
					Soluble fraction of
			Extractive	Klason's Lignin	lignins
	Treatment	FSP (%)	content (% of anhydrous	content (% of extractive free	(mg catechine equivalent per g of
			wood)	anhydrous wood)	extractive free
					anhydrous wood)
	A	32.9 (2.2)	0.7 (0.2)	24.3 (0.1)	13.5 (2.7)
	B	30.6 (2.1)	1.1 (0.2)	24.3 (0.6)	11.3 (0.8)
	C	28.0 (1.6)	2.8 (0.1)	23.8 (0.4)	11.5 (1.3)
	D	18.1 (0.8)	4.1 (0.3)	26.5 (0.5)	7.9 (0.3)
	E	16.6 (0.9)	5.0 (0.5)	35.4 (1.8)	9.1 (0.5)
	F	29.6 (1.8)	1.5 (0.1)	19.9 (0.1)	12.2 (1.5)
	G	23.2 (1.3)	4.2 (0.2)	26.6 (0.3)	14.8 (2.4)
	H	22.0 (1.0)	3.3 (0.2)	20.3 (0.1)	12.7 (1.6)
	I	15.8 (0.9)	5.0 (0.3)	29.7 (0.2)	12.2 (1.0)
	J	15.3 (1.8)	4.1 (0.2)	24.6 (0.2)	9.9 (0.1)
	K	31.7 (2.1)	-	-	-
	L	32.4 (1.5)	-	-	-
	M	24.1 (1.2)	-	-	-
	N	24.3 (1.2)	-	-	-
	O	13.8 (2.2)	-	-	-

Table 3 .

 3 Chemical composition of natural and heat-treated beech wood: carbohydrates and glucuronic acid content in weight percent of extractive free anhydrous wood (measured by HPLC).

	Treatment	Glucose (FGC) in (weight % of extractive free anhydrous wood)	Cellobiose in (weight % of extractive free anhydrous wood)	Xylose (FXC) in (weight % of extractive free anhydrous wood)	Glucuronic acid in (weight % of extractive free anhydrous wood)	Mannose in (weight % of extractive free anhydrous wood)
	A	33.4 (0.8)	2.3 (0.3)	20.6 (0.9)	12.1 (1.7)	3.3 (0.6)
	B	34.3 (1.3)	2.3 (0.4)	18.0 (1.6)	10.9 (2.2)	3.8 (0.7)
	C	35.1 (1.0)	2.7 (0.4)	17.4 (0.9)	9.4 (2.1)	3.5 (0.6)
	D	38.8 (0.6)	3.4 (0.2)	13.5 (0.2)	8.9 (1.4)	3.4 (0.5)
	E	39.7 (1.1)	4.5 (0.8)	6.7 (0.5)	3.9 (2.1)	3.8 (0.3)
	F	34.1 (0.4)	3.5 (0.2)	18.5 (0.1)	11.3 (0.5)	4.8 (0.6)
	G	34.6 (0.7)	3.9 (0.3)	14.6 (0.6)	10.0 (1.1)	6.5 (3.8 )
	H	35.0 (0.5)	4.7 (0.3)	11.2 (0.7)	7.9 (0.6)	4.7 (1.0)
	I	38.8 (1.4)	4.4 (0.1)	7.2 (02.)	10.1 (0.3)	3.1 (0.6)
	J	35.8 (0.7)	4.4 (0.1)	10.1 (0.4)	4.8 (1.3)	4.2 (0.4)

Table 1 .

 1 Parameters of the heat treatments of retification® (time and temperature).

Table 2 .

 2 Table 2. FSP in %H 2 O and chemical composition of natural and heat-treated wood:


Jérôme Lohéac
email: jerome.loheac@irccyn.ec-nantes.fr

Jean-Francois Scheid
email: jean-francois.scheid@univ-lorraine.fr

Jean-François Scheid

Time optimal control for a perturbed Brockett integrator Application to time optimal axi-symmetric micro-swimmers

Keywords: Time optimal controllability, State constraints, Numerical resolution, Brockett integrator, Micro-swimmers AMS classification. 34H05, 49J15, 49K15, 49K20, 93C15

published or not. The documents may come

Time optimal control for a perturbed Brockett integrator

Introduction

In this paper, we consider the system:

ḣ = V (α), λ (1.1a) α = λ , (1.1b)
where (h, α) ∈ R × R n are the state variables, λ ∈ R n is the control variable and V : R n → R n is a given map. In (1.1a), •, • denotes the Euclidean product in R n . The system (1.1a), (1.1b) is completed with the initial conditions: h(0) = 0 and α(0) = 0 .

(1.1c)

The first aim of this work is the following. Given T > 0, h f ∈ R * and c > 0, find a control λ : [0, T ] → R n such that the final target h f is reached at time T :

h(T ) = h f and α(T ) = 0 , (1.2)
together with the state constraint:

|α(t)| c (t ∈ [0, T ]) , (1.3) 1
where | • | denotes the Euclidean norm in R n . Providing that V ∈ C ∞ (R n , R n ) and the Jacobian matrix ∇V (0) ∈ M n (R) is not symmetric, we prove that this control problem can be solved. More precisely, there exists a minimal time T ⋆ = T ⋆ (h f , c) > 0 such that the control problem (1.1)-(1.3) admits a solution λ ∈ BV ([0, T ], R n ) together with the control constraint:

|λ(t)| 1
(for a.e. t ∈ [0, T ]) .

(1.4)

The second aim of this work which actually constitutes the main issue of this note, is to present an efficient numerical method for the computation of time optimal controls for the system (1.1)- (1.4). The method relies on explicit optimal solutions given in [START_REF] Lohéac | Time optimal control for a nonholonomic system with state constraint[END_REF] for the time optimal control problem (1.1)-(1.4) with V (α) = M α where M ∈ M n (R) is a non-symmetric matrix. The numerical method is also based on a full direct discretization as depicted in [START_REF] Trélat | Contrôle optimal[END_REF]§ 9.II.1]. Let us mention that (1.1) with V (α) = M α is a natural extension of the Brockett integrator:

ḣ = α 2 λ 1 -α 1 λ 2 , α1 = λ 1 , α2 = λ 2 .
Roughly speaking, up to some adaptations, our numerical method is a full direct discretization method (see [30, § 9.II.1]) initialized with the analytic solution obtained in [START_REF] Lohéac | Time optimal control for a nonholonomic system with state constraint[END_REF] for M = ∇V (0). Numerical examples will show that this initialisation procedure provides better results than arbitrary initial guesses. It improves the robustness and reduces the computational time.

This strategy is applied to the time optimal control of micro-swimmers. We will see that the efficiency of the micro-swimmer we consider is low and the minimal time is high. Our strategy also provides a preliminary estimate of the minimal control time and hence we also get an estimate of the time discretization step required for suitable results. In this note, we only focus on numerical results. No theoretical proof of robustness nor algorithmic complexity are given. We emphasize that even if direct discretization methods are robust for control problems, knowing an approximation of the optimal solution is even better.

Moreover, in this note, for the sake of clarity, we only present the dimension case n = 2, but a similar methodology can be obtained for any n 2. This note is organized as follows. We first give controllability and time optimal controllability results for the control problem (1.1)-(1.4) in section 2. Analytic results obtained in [START_REF] Lohéac | Time optimal control for a nonholonomic system with state constraint[END_REF] are recall in section 3. In section 4, we present the numerical strategy and its application to time optimal micro-swimmers is studied in section 5.

Controllability and time optimal controllability

Using Chow Theorem, the following controllability result holds.

Proposition 2.1. Let c > 0 and V ∈ C ∞ (R 2 , R 2
) be given and assume that the Jacobian matrix ∇V (0) is not symmetric i.e. :

∇V (0) = (∇V (0)) ⊤ .
Then, for every h f ∈ R * and every T > 0, there exists λ ∈ C 0 ([0, T ], R 2 ) with |λ(t)| 1 for all t ∈ [0, T ] such that the solution of (1.1) satisfies the final condition (1.2) together with the state constraint (1.3).

Proof. Let us write f

1 (h, α) = V 1 (α), 1, 0 ⊤ and f 2 (h, α) = V 2 (α), 0, 1 ⊤ for V (α) = (V 1 (α), V 2 (α)) ⊤ .
The system (1.1) becomes: d dt

h α = λ 1 f 1 (h, α) + λ 2 f 2 (h, α) .
The Lie bracket of f 1 and f 2 at the point 0 is given by:

[f 1 , f 2 ] 0 =   ∂ α1 V 2 (0) -∂ α2 V 1 (0) 0 0   ,
which does not vanish if ∇V (0) is not a symmetric matrix. Thus, under this assumption, the Lie algebra generated by {f 1 , f 2 } and evaluated at the point 0 is of dimension 3. In addition, this Lie algebra is independent of h and the fist order Lie bracket is continuous with respect to α. Therefore, there exists ε > 0 such that for every (h, α) ∈ R × B 0 (ε), the Lie algebra generated by {f 1 , f 2 } evaluated at the point (h, α) is of dimension 3. The result follows from Chow's Theorem (see for instance [30, chap. 5, Proposition 5.14], [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF] or [START_REF] Jurdjevic | Geometric control theory[END_REF]).

This result combined with the Filippov Theorem (see for instance [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF][START_REF] Cesari | Optimization-theory and applications[END_REF][START_REF] Hartl | A survey of the maximum principles for optimal control problems with state constraints[END_REF]) leads to the following optimal control result:

Proposition 2.2. Let c > 0 and V ∈ C ∞ (R 2 , R 2 ) such that ∇V (0) = (∇V (0)) ⊤ .
Then, the set of times T > 0 such that there exists λ ∈ BV (0, T )

(T n ) n ∈ (R * + )
N and associated controls λ n ∈ BV (0, T n ) 2 solving the control problem such that T n → 0. Thus, the corresponding time optimal control does not make sense. [START_REF] Lohéac | Time optimal control for a nonholonomic system with state constraint[END_REF]Proposition 3.3], it can be proved that if λ is a time optimal control, then |λ(t)| = 1

As in

for almost every t ∈ [0, T ⋆ ]. Consequently, the control variable can be written as λ(t) = R(Θ(t)) for almost every t, with Θ(t) ∈ R and

R(θ) = cos θ -sin θ sin θ cos θ (θ ∈ R) , (2.5)
denotes the rotation matrix of angle θ.

An approximate linearized time optimal control problem

For a small parameter c > 0, we have V (α) = V (0) + ∇V (0)α + o(c) for every α ∈ B 0 (c). Instead of considering the time optimal controllability of the system (1.1), we first consider the approximated linear system:

ḣ = V (0) + ∇V (0)α, λ , (3.6a)
α = λ , (3.6b)
with the initial condition (1.1c). We assume that the Jacobian matrix of V at α = 0 is not symmetric. Thus, Proposition 2.2 ensures the existence of a minimal time T ⋆ such that the solution (h, α, λ) of the control problem (3.6) with initial condition (1.1c) satisfies (1.2) together with (1.3) and (1.4). In addition, according to [START_REF] Lohéac | Time optimal control for a nonholonomic system with state constraint[END_REF]Proposition 7], this optimal time T ⋆ and time optimal controls can be explicitly computed.

Proposition 3.1. Let h f = 0, c > 0 and V ∈ C ∞ (R 2 , R 2 ). Let γ ∈ R such that γJ = 1 2 ∇V (0) -(∇V (0)) ⊤ with J = 0 -1 1 0
. Assume γ = 0 and define:

d ⋆ = 2|h f | π|γ| and τ = πc 2 .
Then, the minimal time T ⋆ of the control problem (3.6), (1.1c), (1.2), (1.3) and (1.4) is given by:

T ⋆ =    πd ⋆ if c d ⋆ , π(d ⋆ ) 2 2c + τ otherwise.
Moreover, the time optimal control λ ⋆ is continuous and given by:

λ ⋆ (t) = R (Θ(t) + Θ 0 ) 1 0 (t ∈ [0, T ⋆ ]) ;
where, for every θ ∈ R, R(θ) is the rotation matrix of angle θ defined by (2.5), Θ 0 is any real constant and Θ(t) is given by:

• If c d ⋆ , Θ(t) = sign(γh f )2π t T ⋆ (t ∈ [0, T ⋆ ]) ; • if c < d ⋆ , Θ(t) =                sign(γh f )π t τ if t ∈ [0, τ ) , sign(γh f )π t + τ 2τ if t ∈ [τ, T ⋆ -τ ] , sign(γh f )π t + τ -T * 2 τ if t ∈ (T ⋆ -τ, T ⋆ ] .
Proof. Since α(0) = α(T ), it can be easily proved that λ is a control on [0, T ] for the system (3.6), (1.1c)-(1.4) if and only if λ is a control on [0, T ] for the following system: Remark 3.2. Let us explain the formulae of the optimal solution given in Proposition 3.1. The optimal deformation α ⋆ (t) = t 0 λ ⋆ (s) ds possesses the following characteristics (see Figure 1):

ḣ = 1 2 ∇V (0) -(∇V (0)) ⊤ α, λ , (3.7a
) α = λ , (3.7b
• If c d ⋆ , the optimal trajectory t → α ⋆ (t) is a circle of diameter d ⋆ , starting from 0;
• if c < d ⋆ , the optimal trajectory t → α ⋆ (t) is composed of three arcs of circle. The first arc of circle is an half-circle of diameter c, starting from 0 . The second one lies on the circle of diameter 2c, centred at the origin 0. Finally, the third arc of circle is a half-circle of diameter c, reaching the origin 0 at the final time t = T ⋆ .

Numerical computation of a time optimal trajectory

In this section, we will explain our numerical approach for computing optimal trajectories, which is mainly inspired from [30, § 9.II.1].

Finite dimensional minimisation problem

Let us rewrite our time optimal control problem as a general optimisation problem. First of all, we consider a time control T (i.e. a time that solves the control problem (1.1)-(1.4)), so that the

0 α 2 α 1 |α| = c λ(0) d ⋆ 0 α 2 α 1 |α| = c λ(0) α(T ⋆ -τ ) α(τ ) (a) Case c d ⋆ . (b) Case c < d ⋆ .
Figure 1: Typical time optimal trajectories of α = (α 1 α 2 ) ⊤ .

minimal time T ⋆ is lower than T and we define the cost function

J 0 : (T, λ) ∈ [0, T ] × BV (0, T ) 2 → T ∈ R.
The optimisation problem reads as follows:

min T = J 0 (T, λ) (T, λ) ∈ [0, T ] × BV (0, T ) 2 , |λ(t)| 1 (for a.e. t ∈ [0, T ]) , |α(t)| c (t ∈ [0, T ]) , α(0) = α(T ) = 0 , h(0) = 0 , and h(T ) = h f ,
where α and h are solutions of ḣ = V (α), λ , α = λ .

According to item 2 of Remark 2.3, any optimal solution λ satisfies |λ(t)| = 1 for almost every t and hence it can be expressed as λ(t) = R(Θ(t))e 1 for almost every t, with e 1 = (1, 0) ⊤ . The state variable α is also given by α(t) = t 0 R(Θ(s))e 1 ds. Then, we introduce the new cost function J 1 : (T, Θ) ∈ [0, T ] × BV (0, T ) → T ∈ R and the above optimisation problem can be simplified to:

min T = J 1 (T, Θ) , (T, Θ) ∈ [0, T ] × BV (0, T ) , t 0 R(Θ(s))e 1 ds c (t ∈ [0, T ]) , T 0 R(Θ(t))e 1 dt = 0 , T 0 V t 0 R(Θ(s))e 1 ds , R(Θ(t))e 1 dt = h f . (4.8)
In order to numerically tackle this problem, we set n ∈ N * an integer and

t 0 < • • • < t n n + 1 time discretization points. For every i ∈ {0, • • • , n}, we consider Θ i ∈ R an approximation of Θ(t i ).
Let us also define a quadrature formula

I(t 0 , • • • , t k ; f 0 , • • • , f k ) approximating t k t0 f (s) ds with f i = f (t i ). Thus, defining the cost function J n : (T, Θ 0 , • • • , Θ n ) ∈ R + ×R n+1 → T ∈ R, the discretized version of (4.8) is: min T = J n (T, Θ 0 , • • • , Θ n ) , T, Θ 0 , • • • , Θ n ∈ R + × R n+1 , |α k | c (k ∈ {1, • • • , n -1}) , I 0, T n , • • • , (n-1)T n , T ; R(Θ 0 )e 1 , • • • , R(Θ n )e 1 = 0 , I 0, T n , • • • , (n-1)T n , T ; V (α 0 ) , R(Θ 0 )e 1 , • • • , V (α n ) , R(Θ n )e 1 = h f , with α k = I 0, T n , • • • , kT n ; R(Θ 0 )e 1 , • • • , R(Θ k )e 1 , for k ∈ {0, • • • , n} . (4.9)
This nonlinear minimisation problem can be solved for instance with an interior-point method, SQP (sequential quadratic programming) or an active-set algorithm (see e.g. [START_REF] Bonnans | Numerical optimization. Universitext[END_REF][START_REF] Byrd | An interior point algorithm for large-scale nonlinear programming[END_REF][START_REF] Gill | Practical optimization[END_REF]) as implemented is the Optimization toolbox of MATLAB through the fmincon function. We will see on practical examples that even if the interior point algorithm is robust with respect to the initialisation, the better the initialisation guess is, the better the result is.

A natural initialisation procedure is to use the time optimal control obtained analytically for the linearized problem with ∇V (0)α in place of V (α).

Numerical example

Let us consider:

V (α) = -α 2 + α 2 2 + 10α 3 2 α 1 + α 2 1 . (4.10)
For the computation of the integrals in the minimisation problem, we use the trapeze method, i.e., for f ∈ C([τ 0 , τ 1 ], R), we approximate τ1 τ0 f (t) dt by:

I(t 0 , • • • , t n ; f 0 , • • • , f n ) = n k=1 t k -t k-1 2 (f k-1 + f k ) , (4.11)
with t k = τ 0 + k(τ1-τ0) n and f k = f (t k ). In order to numerically solve the optimisation problem (4.9) with (4.10), we use the fmincon function of the Optimization toolbox of MATLAB with a SQP (sequential quadratic programing) algorithm.A first observation is that if we choose the initial value (T, Θ 0 , • • • , Θ n ) = (1, 0, • • • , 0), the algorithm is not converging. On the other hand, when initializing with the optimal solution of the approximate linearized problem (3.6),(1.1c),(1.2)-(1.4) (see Section 3), we obtain a reasonable solution to the nonlinear optimisation problem (4.9) (see Figure 2).

The SQP algorithm with the trivial initialisation point (T, Θ 0 , A way to avoid this is to use the initialisation point:

• • • , Θ n ) = (1, 0, • • • , 0)
(T, Θ 0 , • • • , Θ n ) = (1, 0, π, • • • , nπ) , (4.12)
so that the discrete state variables associated to this control are given by α k = 0 (for every k ∈ {0, • • • , n}) and thus, the state constraints are satisfied. We compare in Table 1 the solutions obtained for these to different types of initialization values. These results show that the algorithm is more efficient when the optimal solution of the linearized problem is chosen as an initial guess.

In Table 1, the comparison is made with V given by (4.10), the target h f = 1 and the state constraint given by (1.3) with c = 0.3. The number of discretization points is given by n + 1, T ⋆ is the optimal time obtained, N SQP is the number of iterations inside the SQP algorithm and T CPU returns the CPU time used for computations.

In addition, the optimal solution for the nonlinear optimization problem (4.9) obtained with the initial value (4.12) does not look as good as the optimal solution obtained from the optimal solution of the linearized problem (3.6) (see in comparison Figure 2 and Figure 3).

Θ(t) 0 π 2π 0 1 2 3 4 5 Θ t (c) solution t → Θ(t) mod 2π
Figure 3: Solution of the time optimal problem (1.1)-(1.4) with V given by (4.10). The objective is h f = 1, the constraint parameter c is 0.3 and the number of discretization points n is 300. The SQP algorithm has been initialized with the initial value (4.12).

Application to time optimal micro-swimmers

Understanding the motion of micro-organisms is a challenging issue since at their size the fluid forces are only viscous forces and micro-organisms live in a world where inertia does not exist. Despite the pioneer works modeling and analysing the motions of micro-swimmers (see for instance [START_REF] Taylor | Analysis of the swimming of microscopic organisms[END_REF][START_REF] Lighthill | On the squirming motion of nearly spherical deformable bodies through liquids at very small Reynolds numbers[END_REF][START_REF] Lighthill | Mathematical biofluiddynamics[END_REF][START_REF] Purcell | Life at low Reynolds number[END_REF][START_REF] Childress | Mechanics of swimming and flying[END_REF][START_REF] Shapere | Geometry of self-propulsion at low Reynolds number[END_REF][START_REF] Shapere | Efficiencies of self-propulsion at low Reynolds number[END_REF]), the swimming of micro-organisms has only been recently tackled as a control problem. A lot of controllability results for various swimmers has been obtained (see for instance [START_REF] Alouges | Optimal strokes for low Reynolds number swimmers: An example[END_REF][START_REF] Alouges | Optimal strokes for axisymmetric microswimmers[END_REF][START_REF] Michelin | Efficiency optimization and symmetry-breaking in a model of ciliary locomotion[END_REF][START_REF] San Martin | A control theoretic approach to the swimming of microscopic organisms[END_REF][START_REF] Lohéac | Controllability and time optimal control for low reynolds numbers swimmers[END_REF] for axi-symmetric swimmers, [START_REF] Alouges | Optimally swimming Stokesian robots[END_REF][START_REF] Lohéac | Controllability of 3D low Reynolds number swimmers[END_REF] for general swimmers or [START_REF] Alouges | Enhanced controllability of low Reynolds number swimmers in the presence of a wall[END_REF][START_REF] Gérard-Varet | Rough wall effect on micro-swimmers[END_REF] when the fluid domain is not the whole space R 3 ). Let us also point out that numerical strategy for axi-symmetric swimmers have already been presented in [START_REF] Alouges | Numerical strategies for stroke optimization of axisymmetric microswimmers[END_REF]. The work explained in this note can be seen as an improvement of this numerical strategy when the cost function is reduced to the final time.

In this section, we consider a micro-swimmer performing axi-symmetric deformations. We denote by (e 1 , e 2 , e 3 ) the canonical basis of R 3 and we assume that e 1 represents the symmetry axis. At any time, the swimmer will be diffeomorphic to the unit sphere of R 3 and its shape at rest is the unit sphere S 0 . The control problem is the following: starting from the initial location 0, reach the final position h f e 1 by shape changes such that at the initial and final positions, the swimmer is the unit sphere of R 3 .

For those swimmers, we will study the time optimal controllability of the dynamical system associated to the swimming problem.

Modeling and problem formulation

Axi-symmetric coordinates

Since we are considering axi-symmetric swimmers, we introduce the spherical coordinate system (r, θ, φ) ∈ R + × [0, π] × [0, 2π). For every x = (x 1 x 2 x 3 ) ⊤ ∈ R 3 , the spherical coordinates (r, θ, φ) = r(x), θ(x), φ(x) are such that:

x = r   cos θ sin θ cos φ sin θ sin φ   (x ∈ R 3 , (r, θ, φ) ∈ R + × [0, π] × [0, 2π)) ,
with the associated local system of unit vectors (e r , e θ , e φ ) given by (see Figure 4):

Swimmer's deformations

The shape of the swimmer at rest is the unit ball of R 3 which forms the reference shape denoted by S 0 . We assume that the deformation of the swimmer is axi-symmetric with respect to the symmetry axis e 1 . More precisely, we assume that the deformation X is built from two elementary deformations D 1 and D 2 , that is

X(t, x) = x + α 1 (t)D 1 (x) + α 2 (t)D 2 (x) , (t 0 , x ∈ R 3 ) , (5.13)
where D 1 and D 2 are axi-symmetric and radial deformations, i.e.

D i (x) = r(x)δ i (θ(x))e r (x) (t 0 , x ∈ R 3 , i ∈ {1, 2}) , (5.14) with δ i ∈ C 1 ([0, π], R) and α i ∈ L ∞ (R + , R) such that: α 1 (t)δ 1 (θ) + α 2 (t)δ 2 (θ) > -1 (t 0 , θ ∈ [0, π]) , (5.15)
so that X(t, •) is a C 1 -diffeomorphism on S 0 . In practice, we will consider c > 0 small enough such that (1.3) implies (5.15). We define the domain S(t) occupied by the deformed swimmer at time t in the reference frame attached to the swimmer:

S(t) = X(t, S 0 ) ⊂ R 3 (t 0) .
We also assume that the deformation X does not produce any translation. To this end, we introduce the mass density ρ 0 (x) = 1 in the shape of the swimmer S 0 at rest and we assume that the mass is locally preserved during the deformation, that is to say that the density of the swimmer at any time t 0 is given by:

ρ(t, X(t, x)) = 1 |JacX(t, x)| (t 0 , x ∈ S 0 ) ,
where JacX(t, •) denotes the Jacobian of the mapping X(t, •). According to (5.13) and (5.14) we have:

ρ(t, X(t, x)) = 1 (1 + α 1 (t)δ 1 (θ(x)) + α 2 (t)δ 2 (θ(x))) 3 (t 0 , x ∈ S 0 ) .
With this mass density, we have, for all t 0:

S(t) ρ(t, x) dx = S0 dx := m 0 (5.16)
and the mass center of the swimmer is given by:

S(t) xρ(t, x) dx = S0 (x + α 1 (t)D 1 (x) + α 2 (t)D 2 (x)) dx = π 0 1 0 2π 0 r (1 + α 1 (t)δ 1 (θ) + α 2 (t)δ 2 (θ)) e r (θ, φ) r 2 sin θdφdrdθ = α 1 (t) π 2 π 0 δ 1 (θ) cos θ sin θ dθ + α 2 (t) π 2 π 0 δ 2 (θ) cos θ sin θ dθ e 1 ,
Consequently, we assume:

π 0 δ i (θ) sin(2θ) dθ = 0 (i ∈ {1, 2}) , (5.17)
so that the mass center of the swimmer does not move with the deformation X.

Finally, let us consider the domain S † (t) occupied by the swimmer in the fluid at time t which is given by S † (t) = S(t) + h(t)e 1 .

(

Since we have assumed that the deformation X does not introduce any translation, h(t)e 1 is the mass center position of the swimmer in the fluid domain at time t. The domains S 0 , S(t) and S † (t) are depicted on Figure 5.

In terms of the control theory, α = (α 1 , α 2 ) ⊤ is the system's input and h is its output we aim to control.

Micro-swimmer and fluid flow

It is well known that the system governing the motion of a micro-swimmer is inertia less (see for instance [20, § 5.3] or [START_REF] Childress | Mechanics of swimming and flying[END_REF]) and reduced to the following problem coupling,

• the Stokes equation: where F (t) = R 3 \ S(t) and with lim |x|→∞ u(t, x) = 0 ;

-∆u + ∇p = 0 in F (t) , (5.19a
) div u = 0 in F (t) , (5.19b)
S(t) S 0 0 0 0 h(t) X(t, •) Id + h(t)e 1 S † (t)
(5.19c)

• the boundary condition (continuity of the velocity at the boundary of the micro-swimmer):

u(t, x) = ḣ(t)e 1 + α1 (t)D 1 X(t, •) -1 (x) + α2 (t)D 2 X(t, •) -1 (x) e r X(t, •) -1 (x) ; (5.20)
• quasi-static Newton law:

0 = F(t) • e 1 , (5.21)
where

F(t) = ∂S(t)
σ(u, p)n dΓ with σ(u, p) = ∇u + (∇u) ⊤ -pI 3 is the Cauchy stress tensor.

The well posedness of (5.19) follows from [START_REF] Girault | A well-posed problem for the exterior Stokes equations in two and three dimensions[END_REF], where the limit (5.19c) has to be understood in a weak sense (that is 1 + |x| 2 u(x) ∈ L 2 (F ) 3 ). More precisely, we have (u, p) ∈ W 1 0 (F (t)) 3 × L 2 (F (t)) where we have set:

W 1 0 (F (t)) = ϕ ∈ L 2 loc (F (t)) , ∇ϕ ∈ L 2 (F (t)) 3 , 1 + |x| 2 ϕ ∈ L 2 (F (t)) .
This ensures that u| ∂F (t) ∈ H 2 (∂S(t)) 3 and 1 ∈ H 1 2 (∂S(t))).

A geometric control problem

Let us first notice that in the full system (5.19)-(5.21) the time does not appear directly but only through the parameter α(t). Consequently, we define S(α) as the image of S 0 by the map X(α) : (5.15) holds. For convenience, we also define the corresponding fluid domain F (α) = R 3 \ S(α).

x ∈ R 3 → x + α 1 D 1 (x) + α 2 D 2 (x) ∈ R 3 , for α ∈ R 2 such that
For every α ∈ R 2 satisfying (5.15), we define (u α 0 , p α 0 ) the solution of: and for i ∈ {1, 2}, (u α i , p α i ) is the solution of:

-∆u α 0 + ∇p α 0 = 0 in F (α) , (5.22a) div u α 0 = 0 in F (α) , (5.22b
-∆u α i + ∇p α i = 0 in F (α) , (5.23a) div u α i = 0 in F (α) , (5.23b) lim |x|→∞ u α i (x) = 0 , (5.23c)
u α i (x) = D i X(α) -1 (x) , for all x ∈ ∂S(α) (5.23d)
Then, by decomposing the solution (u(t, •), p(t, •)) of (5.19)- (5.20) in terms of {(u α i , p α i )} i∈{0,1,2} and due to the linearity of the Cauchy stress tensor, relation (5.21) becomes:

ḣ(t)F 0 (α(t)) • e 1 = -( α1 (t)F 1 (α(t)) + α2 (t)F 2 (α(t))) • e 1 , (5.24)
where for every i ∈ {0, 1, 2} and every α ∈ R 2 satisfying (5.15), we have set:

F i (α) = ∂S(α)
σ(u α i , p α i )n dΓ (5.25) with σ(u α i , p α i ) = ∇u α i + (∇u α i ) ⊤p α i I 3 . By Green formula, we have:

e 1 • ∂F σ(u 0 , p 0 )n dΓ = 2 F D(u 0 ) : D(u 0 ) dx > 0 .
Consequently, F 0 (α) • e 1 = 0 and the control problem can be recast as the generalized Brockett system (1.1) with

V (α) = -1 F 0 (α) • e 1 F 1 (α) • e 1 F 2 (α) • e 1 .
(5.26)

This system form has been already establish in the pioneer work of F. Alouges, A. DeSimone and A. Lefebvre [START_REF] Alouges | Optimal strokes for axisymmetric microswimmers[END_REF].

Notice that even if α represents the physical control of the swimming system, it is more convenience for analysis to consider λ = α as the control variable. This will also allow us to control both the shape and the position of the swimmer. We also set the initial conditions (1.1c) for this system and the target position to be reached in a time T > 0 is given by (1.2). These initial and final conditions mean that the micro-swimmer is the unit sphere located at the origin at initial time and a unit sphere located in h f e 1 at final time T .

We point out that a state constraint on α is needed to be able to well define F i (α(t)) for every time t 0. This constraint is given by (5.15).

In addition, according to [21, Lemma 1] or [START_REF] Lohéac | Controllability and time optimal control for low reynolds numbers swimmers[END_REF]Theorem 2.6], there exists c > 0 small enough, such that the mapping α ∈ R 2 → V (α) ∈ R 2 with V defined by (5.26) is of regularity C ∞ on B 0 (2c). Consequently, we chose c > 0 small enough such that this condition holds together with (5.15). Remark 5.1. From [START_REF] Lohéac | Controllability and time optimal control for low reynolds numbers swimmers[END_REF], there exists δ 1 , δ 2 ∈ C 1 ([0, π], R) in (5.15) such that ∇V (0) is not symmetric. Moreover, using the arguments of [START_REF] Lohéac | Controllability of 3D low Reynolds number swimmers[END_REF], this situation is generic.

Numerical computation of a time optimal trajectory

Let us consider the deformations D 1 and D 2 similar to the one given in [START_REF] Lohéac | Controllability and time optimal control for low reynolds numbers swimmers[END_REF], that is to say that D 1 and D 2 are given by (5.14) with

δ 1 (θ) = P 2 (cos θ) and δ 2 (θ) = P 3 (cos θ) , (5.27)
where P 2 (resp. P 3 ) is the Legendre polynomial of order 2 (resp. 3). With these two elementary deformations, we obtain from [START_REF] Shapere | Efficiencies of self-propulsion at low Reynolds number[END_REF],

V These expressions allow us to compute the explicit form of the time optimal controls for the approximate linearized system (3.6).

In order to numerically compute the fluid forces F i (α) • e 1 for i ∈ {0, 1, 2}, we use the spherical harmonics expansion of the exterior Stokes solution given in [START_REF] Brenner | The stokes resistance of a slightly deformed sphere[END_REF] (see also [START_REF] Happel | Low Reynolds number hydrodynamics with special applications to particulate media[END_REF] or [START_REF] Shapere | Efficiencies of self-propulsion at low Reynolds number[END_REF]). Then, we compute the time optimal controls by using the direct discretization method explained in section 4. The initial guess for the discrete nonlinear optimization problem is chosen as the explicit optimal solution of the approximate linear control problem (3.6) given in Proposition 3.1. We expect that the approximate optimal deformation for the linearized problem is close to the optimal solution for the nonlinear problem, so that the direct method will converge quickly.

We apply this computational method with the deformation D 1 and D 2 given above through Legendre polynomials and we choose c = 0.3 and h f = 1 2 . The optimal trajectories for h and α are depicted on Figure 6. The optimal time is T ⋆ ≃ 48.9. We numerically observe on Figure 6 that the optimal trajectory for α is mainly periodic. On Figure 7, we also give the optimal trajectory of h during a period in time and finally we plot on Figure 8 the different shapes of the swimmer under the optimal deformation for different instants in the time period .

Conclusion

In this note, we have presented a numerical strategy for solving a time minimal control problem for a nonlinear system which generalizes the Brockett integrator. This method is based on the explicit computation of time optimal controls for a state constraint Brockett system which approximates the original nonlinear control problem. These explicit solutions are used as initial guesses in the nonlinear problem. Numerical tests have shown the good behavior of this strategy. This numerical method is then applied to a time optimal control problem for micro-swimmer with shape changes. The shapes of the swimmer are small deformations of the unit sphere. Despite two elementary deformations are used for the swimmer's shapes, it is easy to extend this work to the case of a finite number of elementary deformations.

Finally, we observe numerically that the optimal trajectory for α is mainly periodic. It would be interesting to prove rigorously this property. For instance, one could expect the following result : for |h f | large enough and for α ∈ A with A a compact set of R 2 , the trajectory of α is composed by one starting curve followed by a periodic curve and ending with a final curve. Such a result would be a real improvement for numerical computations since we would only have to solve three simpler (but coupled) time optimal control problems.

  ) together with (3.6), (1.1c)-(1.4). The result follows from [22, Proposition 7].

  fails to converge because the discrete state trajectory α 0 , • • • , α n associated to this control does not fill the state constraint |α k | c in (4.9). Indeed, for the initial value (T, Θ 0 ,• • • , Θ n ) = (1, 0, • • • , 0),due to (4.11) we have α k = t k e 1 for all k ∈ {1, • • • , n} and in particular we get |α n | = |T e 1 | = T = 1.

6 hFigure 2 :

 62 Figure 2: Solution of the time optimal problem (1.1)-(1.4) with V given by (4.10). The objective is h f = 1, the constraint parameter c is 0.3 and the number of discretization points n is 300. Numerically, we obtain a minimal time T = 4.866577.

1 Figure 4 :

 14 Figure 4: Spherical coordinates.

Figure 5 :

 5 Figure 5: The deformation X and the translation he 1 of the swimmer.

1 2 (

 2 ∂F (t)) and consequently, the expression∂S(t)σ(u, p)n dΓ is seen as a duality product (since σ(u, p)n ∈ H -1

e 1 ,

 1 for all x ∈ ∂S(α)(5.22d)

  Optimal trajectory for h. (b) Optimal trajectory for α.

Figure 6 :

 6 Figure 6: Optimal trajectories for the system (1.1), (1.2)-(1.4) with V given by (5.26) together with the deformations D 1 and D 2 defined by (5.14) with (5.27). In addition, we have set c = 0.3 and h f = 1 2 .

Figure 7 :

 7 Figure 7: Optimal trajectory of h between t 0 = 0.5384 and t 1 = 2.3002 corresponding to one period of α.

Table 1 :

 1 Comparison of the results for the two different types of initialisations with the SQP algorithm.

	n	T ⋆	N SQP T CPU (s)	n	T ⋆	N SQP T CPU (s)
	50	5.14982	716	26.532	50	4.87770	79	3.380
	100 4.87109 1340 102.685	100 4.86974	108	7.381
	150 5.31835 3861 764.928	150 4.86779	135	18.636
	200 4.96607 1425 291.104	200 4.86710	191	32.064
	250 5.17521	242	63.278	250 4.86676	200	41.511
	300 5.12574 1096 546.292	300 4.86658	213	64.223
	350 4.87108	685 363.426	350 4.86646	270 123.074
	400 Not converging		400 4.86721	143 102.644
	450 5.00869	287 287.430	450 4.86661	157 155.655
	500 5.07592	364 580.006	500 4.86636	206 389.241
	(a) Initialisation with the initial value (4.12)	(b) Initialisation with the optimal control obtained for the linearized problem (3.6)
				7

3 and the number of discretization points n is 300. Numerically, we obtain a minimal time T = 4.866577.