N Popoff

D Horezan

Boivin

V Marcadon
email: vincent.marcadon@onera.fr

C Davoine

D Lévêque

A Rafray

F Popoff

N Horezan

D Boivin

High temperature mechanical behavior of tube stackings -Part I: Microstructural and mechanical characterization of Inconel R 600 constitutive material

Keywords: Metallic tube stackings, Cellular materials, High temperature behavior, Mechanical characterization, Microstructural characterization

    High temperature mechanical behavior of tube stackings -Part I: Microstructural and mechanical characterization of Inconel 600 constitutive material

Introduction

Metallic cellular structures, such as metal foams, honeycombs, truss lattices, and hollow-sphere or tube stackings, for instance, have been abundantly studied for several years because of their potential for the development of new multi-functional and lightweight structures [START_REF] Evans | Multifunctionality of cellular metal systems[END_REF][START_REF] Ochsner | Multifunctional metallic hollow sphere structures: manufacturing, properties and application[END_REF]. In order to address this issue, tube stackings have been considered here as a 'model' cellular structure. Indeed, the case of tube stackings is particular, since it is a quite new cellular architecture that has been barely studied, despite presenting an extruded character that is very convenient to develop 2D modeling. Moreover, the elementary tubes can be bought easily and both their geometrical and mechanical properties are very reproducible. Thus, regular stackings can be manufactured with only few defects, which enables a very rich dialogue between experiment and modeling involving large deformations and contact issues. The main mechanism that governs the mechanical behavior of tube stackings is the localized plasticity in the walls of the constitutive cells [START_REF] Marcadon | About the use of microindentation to determine the constitutive material's mechanical behavior of cellular structures[END_REF][START_REF] Marcadon | Mechanical behavior of hollow-tube stackings: Experimental characterization and modeling of the role of their constitutive material behavior[END_REF]. Therefore, this is significantly different from the case of honeycombs loaded transversally, for which buckling and cell wall instabilities are observed [START_REF] Papka | Biaxial crushing of honeycombs -Part I: Experiments[END_REF][START_REF] Papka | In-plane biaxial crushing of honeycombs -Part II: Analysis[END_REF]. It is mainly a question of ratio between the cell wall thickness and the cell diameter. The reader can refer to the companion paper of this one (see Part II [START_REF] V Marcadon | High temperature mechanical behavior of tube stackings -Part II: Comparison between experiment and modeling[END_REF]) for a more complete review on the mechanical behavior of cellular structures. Here, focus is on both the mechanical and microstructural properties of Inconel R 600 as constitutive material of the tube stacking structure studied.

Inconel R 600 is a commercially available Ni-Cr-Fe alloy that exhibits a good combination of strength, formability and environmental resistance at high temperatures. This combination of properties makes this material ideally suited for structural components in chemical processing equipment, nuclear reactors or aerospace applications, for instance. In order to provide a good understanding of the stress corrosion cracking resistance of Inconel R 600, microstructural changes that occur when this alloy is heat-treated in the temperature range of 600 • C to 900 • C have been widely studied [START_REF] Airey | Microstructural Aspects of the Thermal Treatment of Inconel Alloy 600[END_REF][START_REF] Park | Effect of heat treatment on fatigue crack growth rate of Inconel 690 and Inconel 600[END_REF]. Bruemmer et al [START_REF] Bruemmer | High voltage electron microscopy observations of microdefformation in alloys 600 tubing[END_REF] suggested that the distribution of the grain boundary carbides may induce a mechanical effect, since carbides reduce the crack tip stress state by crack tip blunting. The effect of heat treatments on the mechanical properties of Inconel R 600 has been studied also, in particular the role of the carbides located at grain boundaries [START_REF] Wang | Effects of grain boundary carbides on the mechanical properties of Inconel 600[END_REF]. Extensive grain boundary cracks are observed on the fracture surfaces of tensile and impact specimens with continuous grain boundary carbides, decreasing the area reduction observed at the rupture. The shortest delay time for precipitation is observed between 800 • C and 900 • C. At higher temperatures, grain growth occurs when the alloy is heated to about 980 • C. At this temperature, the carbide particles finely dispersed into the microstructure, which inhibit grain growth, begin to coalesce. Carbide solution begins at about 1040 • C, whereas solution annealing, observed above 1190 • C, results in the dissolution of all of the carbides [START_REF] Hall | The Microstructural Response of Mill-Annealed and Solution-Annealed INCONEL 600 to Heat Treatment[END_REF], hence an increased grain size. Moreover, it has been shown that the mechanical properties of Inconel R 600 are also affected by the amount of carbon in solid solution [START_REF] Ivanchenko | Internal friction studies of carbon and its redistribution kinetics in Inconel 600 and 690 alloys under dynamic strain aging conditions[END_REF]. Although most of the carbon resides in carbides, and only a small fraction exists as interstitial solution because of the high chromium content (about 16 wt.%) in Inconel R 600, the interaction of carbon with mobile dislocations results in dynamic strain ageing. In order to have a better understanding of the mechanisms involved in ultra-fine grain processing techniques, grain boundary engineering, or welding, the effects of both the deformation level and temperature on the microstructure of Inconel R 600 have been investigated recently [START_REF] Li | The role of deformation temperature and strain on grain boundary engineering of Inconel600[END_REF]. Depending on the deformation level and the temperature, various mechanisms, such as static and dynamic recrystallizations or dynamic recovery, may dominate and affect the development of microstructures. In particular, these mechanisms contribute to reducing the strain hardening during the deformation. Dynamic recrystallization was also observed in Inconel R 600 friction stir welded samples [START_REF] Song | Effect of welding speed on microstructural and mechanical properties of friction stir welded Inconel 600[END_REF]. Blaizot et al [START_REF] Blaizot | Constitutive model for nickel alloy 690 (Inconel 690) at various strain rates and temperatures[END_REF] recently proposed a model accounting for such complex interactions, but for Inconel R 690.

They showed that, for this alloy, the dislocation density is the main microstructural parameter governing the mechanical behavior for the experimental conditions that they considered. This work is aimed at studying the elasto-viscoplastic properties of Inconel R 600, in order to be able to propose a rich dialogue between experiment and modeling when further characterizing the mechanical behavior of tube stacking structures at high temperatures (see Part II [START_REF] V Marcadon | High temperature mechanical behavior of tube stackings -Part II: Comparison between experiment and modeling[END_REF]).

Although this characterization is presented in Part II only, Section 2 describes the tube stacking studied and its process route to precise the brazing heat treatment applied to the constitutive material. Then Section 3 addresses the characterization of the elasto-viscoplastic properties of Inconel R 600, as constitutive material of the tubes, by means of tensile tests performed on single tubes at different strain rates. Various temperatures were investigated, from room temperature to 800 • C in order to study the influence of both the brazing heat treatment and the test temperature. In order to investigate microstructural evolutions induced by both the heat treatment and the test conditions, electron backscattered diffraction (EBSD) analyses were carried out post-mortem on the samples. They are described in Section 4 in terms of deformation mechanisms, twinning and microstructural changes in the grains.

Studied material

The method for fabricating a 'model' cellular structure consisting of regularly stacked circular tubes is briefly reported here. The small sandwich structure was processed by using a brazing heat treatment. It consisted of two Inconel R 600 skins with a thickness of 1 mm and a 5×5-tube stacking core. The tubes were made of Inconel R 600 also, and their outer diameter and thickness were equal to 4 mm and 250 µm, respectively. These seamless tubes were manufactured by Alloyshop using extrusion and their chemical composition is given in Table 1. This nickelchromium-iron alloy is a standard engineering material for applications that require resistance to corrosion at high temperatures.

The tubes were bonded to each other and to the skins by using Ni-1.5B-19.0Cr-7.3Si-0.08C (wt.%) braze foils. They were stacked in a graphite die lined with Inconel R 600 skins, positioned at the top and at the bottom of the stacking in order to provide the sandwich structures. The internal faces of the graphite die were coated with a boron nitride layer to avoid the chemical interactions between carbon and the metallic structure. The main advantage of this braze ensues from its good oxidation resistance properties, compared to those of the Ni-11.0P wt.% braze that we have used previously for room temperature applications, for instance [START_REF] Marcadon | About the use of microindentation to determine the constitutive material's mechanical behavior of cellular structures[END_REF][START_REF] Marcadon | Mechanical behavior of hollow-tube stackings: Experimental characterization and modeling of the role of their constitutive material behavior[END_REF]. A stacking following a square pattern in its transversal plane was produced by this technique (Fig. 1(a)), since this particular braze was available under the shape of 50 µm thickness foils only. The 'square stacking' was achieved by alternatively stacking a layer of tubes and three braze foils. The bonding was achieved by applying a heat treatment in a vacuum furnace. The Ni-1.5B-19.0Cr-7.3Si-0.08C wt.% alloy is liquid above 1150 • C. Thus, the heat treatment chosen here consisted in an increase in the temperature at 100 • C min -1 up to 1170 • C, then followed by a dwell of 15 min at 1170 • C and, finally, followed by the natural cooling of the furnace. At 1170 • C the liquid concentrated at the contact lines by capillarity, hence metallurgical joints were created between the tubes after cooling.

Energy dispersive spectroscopy (EDS) analyses showed that the braze joint consisted of three distinct phases, brittle intermetallic chromium borides and nickel silicides in a ductile Nirich solid solution. These phases were visible in scanning electron microscopy (SEM) chemical contrast images (Fig. 1(b)). Moreover, the mean grain size in the tube walls was rather large. It varied between 100 and 150 µm. This could result from the brazing heat treatment, despite its short duration. SEM analyses revealed that the braze joint width was around 1.43 mm, although there was a slight difference between the horizontal joints and the vertical ones resulting from the fact that the braze foils were initially placed horizontally between the tube layers themselves and between the tube layers and the skins; the mean width for the horizontal braze joints was equal to 1.63 mm, whereas it was equal to 1.23 mm for the vertical ones. The length of the tubes was equal to 20 mm after the electro-discharge cutting of the samples. The density of the cellular structure was about 2.29 g cm -3 compared to that of the constitutive material, which was equal to 8.25 g cm -3 .

Mechanical behavior of the constitutive material

This section is aimed at characterizing the mechanical properties of the tube material Inconel R 600 in-situ. Indeed, the effective work-hardening of cellular structures results from both architectural effects (e.g., stress localization) and the constitutive material hardening properties [START_REF] Marcadon | About the use of microindentation to determine the constitutive material's mechanical behavior of cellular structures[END_REF][START_REF] Marcadon | Mechanical behavior of hollow-tube stackings: Experimental characterization and modeling of the role of their constitutive material behavior[END_REF][START_REF] Amsterdam | The influence of straininduced damage on the mechanical response of open-cell aluminum foam[END_REF][START_REF] Amsterdam | On the plastic collapse stress of open-cell aluminum foam[END_REF][START_REF] Mangipudi | The microstructural origin of strain hardening in two-dimensional open-cell metal foams[END_REF]. This first step is thus necessary to be able to propose a predictive model of the mechanical behavior of the tube stacking structure and to enrich the comparison with the experiment detailed in Part II of this work [START_REF] V Marcadon | High temperature mechanical behavior of tube stackings -Part II: Comparison between experiment and modeling[END_REF]. In order to address the influence of the test temperature on the elasto-viscoplastic behavior of Inconel R 600, uni-axial tensile tests on single tubes have been performed at various temperatures and strain rates.

Tensile tests on single tubes

The experiments have been performed using an electro-mechanical testing machine (see Fig.

2

). The length of the tubes was 170 mm and their extremities were handled in holding jaws.

Small 20 mm length rods were inserted into the ends of the tubes to allow their holding and to avoid their collapse and their sliding during the tests, which would have altered the measurements. The loading displacement was controlled at various crosshead speeds: 0.4 mm min -1 , 4.5 mm min -1 and 47.0 mm min -1 . Three test temperatures were considered: room temperature (labelled RT in the captions in the following), 600 • C and 800 • C. The longitudinal elongation of the tubes was measured using an extensometer placed in the furnace, whereas the holding jaws were outside the furnace. Such an experimental device was aimed at preventing some error in the measurement of the strain resulting from the heat gradient that existed in the portion of the samples outside of the furnace, hence an improved identification of the plastic properties and Young modulus especially. The real strain rates applied on the portion of the samples inside the furnace, in which the strain was homogeneous, have been computed from the strain measurements provided by the extensometer. Several samples were tested for each test temperature and strain rate, not only after the brazing heat treatment but also before. In the following, the as-received tubes are labelled NHT, for no heat treatment, and the treated ones are labelled BHT, for brazing heat treatment. The stress-strain curves are plotted in Figs. 3(a) and 3(b) and the results are detailed in the next section. They show a quite good test reproducibility. However, the main conclusion from these tests is the fact that the heat treatment drastically decreases the strength of the tubes, so the mechanical behavior of the tube constitutive material is strongly altered. Such a decrease can be explained by the fact that, before the brazing heat treatment, the tubes were greatly hardened because of their cold drawing process. The heat treatment acts as a recovery treatment, hence the drop in the yield stress, but the recovery of higher hardening capabilities. A decrease in the mechanical behavior of Inconel R 600 is also observed when the test temperature is increased. Furthermore, although no viscosity effect is observed at room temperature and 600 • C, a significant one appears at 800 • C. It is also worth noting that, for the tests at 800 • C, the mechanical responses before and after the heat treatment are considerably closer than for those at room temperature and 600 • C.

Therefore, the test temperature acts as a recovery treatment also, if it is high enough. In the literature, the Young modulus of Inconel R 600 is equal to 214 GPa at 22 • C, 180 GPa at 600 • C and 164 GPa at 800 • C (http://www.specialmetals.com/products/inconelalloy600.php), respectively. Thus, the mean values identified experimentally for this quantity are very close to those of the literature despite some experimental scattering. As concerns σ 0.2m , it varies between 170 MPa and 345 MPa at room temperature in the literature, but this quantity depends strongly on the process route and the history of the plastic loading explaining why there is no precise value. We can only conclude that the experimental values measured are consistent with those of the literature.

Identification of an elasto-viscoplastic model from room temperature to 800 • C

The mechanical properties of Inconel R 600 have been identified for the BHT material by fitting the mechanical responses of a homogeneous equivalent medium (HEM) to the experimental stress-strain curves between 0 and 8% of tensile strain, and by using the Levenberg-Marquardt optimization method. A first identification has been attempted, using a classical mono-potential elasto-viscoplastic model, but it failed to capture the saturation of the viscous effect observed at 800 • C when the strain rate was increased (see Fig. 3(b)). Consequently, we have proposed a new ticity to be dissociated [START_REF] Besson | Non-Linear Mechanics of Materials (Solid Mechanics and Its Applications[END_REF]. It allows all of the features of the behavior of Inconel R 600 observed from room temperature to 800 • C to be obtained.

The model equations are the following. The constitutive material was assumed to be elastoviscoplastic isotropic. The elasticity was supposed to be governed by Hooke law:

σ = E m 1 + ν m ε e + ν m 1 -2ν m tr(ε e )i (1)
where σ, ε e and i denote the stress tensor, the elastic part of the strain tensor ε and the secondorder identity tensor, respectively. E m and ν m are the Young modulus and the Poisson ratio for the constitutive material. The plastic and viscous parts of the strain tensor are noted ε p and ε v , respectively, with the classical decomposition of the total strain ε = ε e + ε p + ε v .

A Von Mises criterion was chosen for the yield surface:

F (σ) = J(σ) -R(p) (2)
where J(σ) = (3/2)s : s and s denotes the deviatoric stress tensor:

s = σ -(1/3)tr(σ)i. R(p)
is the instantaneous yield stress, which results from the constitutive material hardening, with

p = t 0
(2/3) εp : εp dt denoting the cumulative plastic strain. Note that this is a purely plastic potential.

Since the stress-strain curves showed an inflection of the hardening when increasing the strain level, on the basis of our previous works [START_REF] Marcadon | Mechanical behavior of hollow-tube stackings: Experimental characterization and modeling of the role of their constitutive material behavior[END_REF] a sum of non-linear exponential hardening terms has been used to account for such a plasticity saturation phenomenon:

R(p) = σ ym + 3 i=1 Q im (1 -exp(-b im p)) (3)
where σ ym denotes the yield stress and the Q im and the b im are the coefficients of the non-linear hardening terms. Accounting for such a saturating hardening was critical in view of modeling aspects, because the stress concentration around the braze joints can result in considerably high localized cumulative plastic strain values, inducing significant overestimation of the plastic contribution.

The viscous contribution was governed by Norton law:

εv = J(σ) D m nm-1 s D m (4)
where n m and D m denote the Norton exponent and the drag stress of the constitutive material, respectively. εv is the time derivative of the cumulative viscous strain. Note that this power law contains no plastic contribution and has no activation threshold. Other models exist in the literature to account for behaviors exhibiting different viscosity regimes, by using hyperbolic sine based equations [START_REF] Perrin | Creep constitutive equations for a 0.5Cr-0.5Mo-0.25V ferritic steel in the temperature range 600-675 • C[END_REF][START_REF] Vakili-Tahami | High-temperature creep rupture of low alloy ferritic steel butt-welded pipes subjected to combined internal pressure and end loadings[END_REF]. These models are commonly used in creep-fatigue applications and they provide very good results [START_REF] Chaboche | Viscoplastic constitutive equations of combustion chambers materials including cyclic hardening and strain aging[END_REF][START_REF] Barrett | An improved unified viscoplastic constitutive model for strain-rate sensitivity in high temperature fatigue[END_REF][START_REF] Barrett | A unified viscoplastic model for high temperature low cycle fatigue of service-aged P91 steel[END_REF]. However, they assume coupled plastic and viscous contributions because they are only mono-potential models. Although the particular case of hyperbolic sine models has not been evaluated here, when trying to use a mono-potential model with two Norton terms to account for two distinct viscous regimes, before opting for a bi-potential formulation, this was not efficient enough to capture the behavior features observed. For these two temperatures, all of the stress-strain curves were almost superimposed, whatever the strain rate. For each test temperature, the fitted curves have been optimized simultaneously taking into account all of the samples tested at this given temperature. Due to the lack of experimental data on the tangential strain, only the longitudinal strain was measured by the extensometer, the Poisson ratio was assumed to be equal to 0.31 and to not depend on the test temperature. Such an assumption was rather relevant according to the literature, given that the Poisson ratio for Inconel R 600 varies in the range of 0.30-0.32 between room temperature and 800 • C (http://www.specialmetals.com/products/inconelalloy600.php). Similarly, since the mean values measured for the Young modulus were in a very good agreement with those of the literature (see Table 2), we set the Young modulus to be equal to that of the literature for the fitting procedure (in order to improve the identification of the yield stress and the beginning of the plasticity). Coefficients b 1m , b 2m and b 3m were also fixed, since some first estimations of their values showed that these parameters were quite independent from the test conditions (temperature and obviously the strain rate). They were set to 500, 50 and 5, respectively, to account for various non-linearity regimes associated with different strain ranges. The aforementioned fitting procedure was then applied at 800 • C. Similarly to the fits at room temperature and 600 • C, the Young modulus was set to its literature value, whereas the Poisson ratio and the b im were set to the same values as before (i.e. 0.31, 500, 50, and 5, respectively). Obviously, at 800 • C the stress-strain curves were fitted simultaneously, but this time distinguishing each strain rate to account for viscosity. The constitutive properties determined for the three different test temperatures are provided in Table 3. The fitted stressstrain curves are plotted with the experimental ones in Figs. 4(a) and 4(b).

It is worth noting that the constitutive bi-potential elasto-viscoplatic model proposed ob-tains the various features of the mechanical behavior observed for Inconel R 600 very well. At room temperature and 600 • C, the values identified for n m are high, thus the viscosity contribution can be considered as negligible. On the contrary, at 800 • C this value is equal to 3.9, which is classical for Ni-base superalloys. As concerns the hardening capabilities, the yield stress decreases significantly when the test temperature is increased. Moreover, the total hardening capability value (i.e.

3 i=1 Q im ) decreases when the test temperature is also increased and is more than twice as high at room temperature and 600 • C as at 800 • C. At room temperature and 600 • C, the main contribution to the total hardening is that related to high strain levels (Q 3m ), whereas at 800 • C it is that related to medium strain levels (Q 2m ). Therefore, the hardening saturation is more important at 800 • C, but it is offset by an additional viscous strain. It can also be noticed that the experimental values measured for σ 0.2m cannot correspond to σ ym by definition and because of the fast hardening contribution Q 1m . They are rather comprised between σ ym + Q 1m and

σ ym + Q 1m + Q 2m .

Post-mortem microstructural characterization

In order to investigate microstructural changes resulting from the tensile tests, post-mortem EBSD analyses have been carried out on tested samples.

Sample preparation

One of the two parts of the broken tubes was cut in a longitudinal median plane. The cut tubes were cold-mounted in an epoxy resin and polished flat in a plane parallel to the tube axis (Fig. 5). The analyzed surface was chosen to set the tube axis parallel to the transversal direction (TD) of the conventional crystallographic macroscopic frame. As a result, TD was parallel to the load direction and the normal direction (ND) of the conventional crystallographic macroscopic frame corresponded to a direction transverse to that of the load. The scanning was performed near the fracture area and the EBSD patterns were indexed in the face-centered cubic lattice corresponding to the crystallographic system of nickel (space group: Fm3m). A grain was defined as a group of at least five neighboring indexed points having the same crystal orientation within a tolerated angular deviation of 5 • . In order to highlight the crystallographic evolution resulting from the test conditions, two microstructures were compared: Sample 1 (NHT, 0.4 mm min -1 at RT) and Sample 21 (BHT, 47.0 mm min -1 at 800 • C). For the all scans, outer surface of the tubes is at the top, whereas inner surface is at the bottom. Moreover, height of scans is not the same because the area reduction at the final stage of the rupture was considerably higher for sample 21 (tested at 800 • C) than fo sample 1 (tested at RT), hence thinner tube walls.

EBSD analyses

EBSD analyses revealed a grain morphology typical of highly-strained materials (see Fig. 6).

The inverse pole figures referring to TD and ND (IPF TD and IPF ND, respectively) are given in Fig. 6 too. The grains are elongated toward the load direction, i.e., the tubular specimen axis. The deformation spreads over the entire surface analyzed and textures are observed both in TD and ND resulting from the rotation of the slip systems in the grains during the plastic deformation of the tubes. IPF TD shows that there were two preferential orientations of the cubic unit cell: a major <111> direction in blue, and a minor <100> direction in red, as given by the corresponding key color. It means that either the diagonal or the edge of the facecentered cubic lattice was parallel to TD. These two preferential orientations were also revealed by the direct pole figures. Indeed, (001) and (111) plane pole figures show intensity peaks in the directions close to TD. The orientation mapping referring to ND also revealed a texture.

It appeared that the grains are mostly green, meaning that <101> direction was preferentially parallel to ND. This particular texture resulted from the tubular shape of the specimens and the ductility of their constitutive material: both longitudinal and transverse deformations occurred at the final stage of the rupture, leading to a large reduction in the tube diameter locally.

EBSD analyses also revealed a heterogeneous grain size in the thickness of the tube walls.

The finer grain microstructure observed at the inner and outer surfaces of the tubes can be explained by a higher plastic deformation locally, resulting from the initial cold drawing process of the tubes. The tubes are drawn through a mould, using a guide bar to keep the tube straight and to give the inner diameter after drawing. Thus, during this process the outer diameter and the wall-thickness of the tube are reduced resulting in a significant plastic deformation in the material at the end of processing. Grain size is quickly reduced by the successive cold drawing and annealing sequences.

Residual strain analysis

In addition of the crystallographic information previously described (texture and grain size), EBSD data post-treatment can also provide information on the residual plastic strains into the grains. The technique enables analyses with high spatial resolution and sensitive examination of the plastic strain, which are otherwise difficult to obtain by X-ray diffraction or transmission electron microscopy because of either spatial / angular resolution or sampling limitations, respectively. For instance, interesting qualitative maps describing the organization of the defect density can be established through the kernel average misorientation (KAM) analysis. KAM has been used to represent the average misorientation between a given point of measure and its nearest neighbors that belongs to the same grain. Another quantity, the grain orientation spread (GOS, i.e., the average orientation deviation between each point in a grain and the average orientation of the grain), has also been determined in each grain. Thus, GOS maps clearly show the most strained grains, whereas KAM maps rather provide the local density of dislocations. The blue color in all of the maps indicates very slight misorientations. The green, yellow, orange and red colors represent progressively increasing levels of internal or local misorientations. Grains are outlined with black boundaries. The same color scales, misorientations from 0 to 18 • for the GOS and from 0 to 5 • for the KAM, are used for the two samples in Fig. 7. Note that for the GOS analysis, twin boundaries have been considered as grain boundaries.

Fig. 7 clearly shows that the residual strains into the grains were lower for Sample 21, which was deformed at 800 • C, than for Sample 1, which was deformed at room temperature.

In particular, both GOS and KAM maps obtained for Sample 1 indicate that a large amount of strain energy has been stored in the material during the test by forming dislocations. The analysis of the GOS map obtained for Sample 21 also shows a majority of grains colored in green and a family of small grains colored in blue at the bottom of the image, i.e., the inner surface of the tube. They result from dynamic recrystallization that occurred during the tests at 800 • C. Recrystallization occurrence has been confirmed by the fact that no texture was observed in these grains on both the IPF TD and ND given in Fig. 6. In the GOS maps of both samples, small grains colored in blue and dispersed into the tube wall were also observed, revealing the presence of precipitate phases. The largest precipitates were titanium nitrides, which present a crystalline structure similar to that of Inconel R 600, titanium atoms in titanium nitrides forming a cubic lattice with centered faces. Due to the fact that they are brittle, they broke during the tests, and small levels of residual strains were observed.

The comparison between the KAM maps obtained for the two samples reveals a lower density of dislocations in the grains in the middle area for the tube tested at 800 For the lower strain rates (0.4 mm min -1 and 4.5 mm min -1 ), the kinetics of dynamic recovery at 800 • C suffice to allow the annihilation of the dislocations and, as a result, only a slight strain hardening is observed. On the contrary, for the highest strain rate (47.0 mm min -1 ), the test elapses too fast to allow complete dynamic recrystallization and recovery, hence work hardening and plasticity predominate. At lower temperatures, it can be expected that dynamic recrystallization and recovery cannot be thermally activated, thus the absence of any observed viscous effect whatever the strain rate.

Conclusions

This work has concerned both the mechanical and microstructural characterizations of Inconel R 600, as constitutive material of a cellular structure consisting of stacked and brazed tubes, from room temperature to 800 • C.

For that purpose, tensile tests at various strain rates and temperatures have been performed.

They have shown the strong influence of the brazing heat treatment on the elasto-viscoplastic behavior of Inconel R 600. The tubes having been strongly hardened during their cold drawing process, the brazing heat treatment acts as a recovery one, hence a lower yield stress and higher hardening capabilities. The test temperature has a significant influence too. By increasing, it results in a decrease of the strength of Inconel R 600. Moreover, at the highest test temperature, the difference between the mechanical response of the as-received tubes and that of the heattreated ones tends to vanish. Viscosity has been revealed only at the highest test temperature.

Post-mortem EBSD observations carried out have shown that both dynamic recrystallization and recovery occur during the tensile tests performed at 800 • C. The competition between workhardening on one side and dynamic recrystallization and recovery on the other side at 800 • C results in a strong viscosity. For the lower strain rates, the thermal energy supplied during the tests allowed the system to use the strain energy stored during the processing of the tube to limit the strain hardening. This was no more true in the case of the higher strain rate at 800 • C, for which dynamic recrystallization and recovery had no longer enough time to be complete.

In that case, work-hardening and plasticity were predominant. Obviously, for the two lower temperatures considered in this study, work-hardening and plasticity were also predominant whatever the strain rate. From these tensile tests, elasto-viscoplastic properties have then been identified for each test temperature (room temperature, 600 • C and 800 • C) assuming a bi-potential to distinguish the respective contributions of the plasticity and the viscosity, and to capture the viscous effect saturation observed when increasing the strain rate. A non-linear hardening has been used to account for a hardening saturation in the modeling.

In the following, the mechanical properties of Inconel R 600 characterized here have been used as input data for the modeling of the compression behavior of the aforementioned tube stacking structure. These results are detailed in a companion paper addressing both experimental and numerical characterizations of the high temperature behavior of tube stackings (see Part II [START_REF] V Marcadon | High temperature mechanical behavior of tube stackings -Part II: Comparison between experiment and modeling[END_REF]).

Figure 1 :

 1 Figure 1: a) Tube stacking performed b) EDS image of a braze joint between two neighboring tubes (chromium borides are in black, nickel silicides are in light grey and the Ni-rich solid solution is in medium grey) .

Figure 2 :

 2 Figure 2: Experimental setup for the tensile tests on single tubes.

Figure 3 :

 3 Figure 3: Stress-strain curves resulting from the tensile tests on both NHT and BHT single tubes a) at room temperature and 600 • C b) at 800 • C.

  min -1 7.5×10 -3 s -1 142.0 143.2 Sample 22 800 • C 47.0 mm min -1 7.5×10 -3 s -1 125.0 145.6 Sample 23 800 • C 4.5 mm min -1 8.8×10 -4 s -1 --Sample 24 800 • C 47.0 mm min -1 7.5×10 -3 s -1 -identification of the behavior of Inconel R 600 by considering a bi-potential elasto-viscoplastic model. Such a more complex model enables the respective contributions of viscosity and plas-

Figure 4 :

 4 Figure 4: Stress-strain curves resulting from the tensile tests on the BHT single tubes and their fitted mechanical responses a) at room temperature and 600 • C b) at 800 • C.

Figure 5 :

 5 Figure 5: Schematic illustration of the EBSD analyzed section.

Figure 6 :

 6 Figure 6: EBSD analysis of samples 1 and 21: IPF referring to TD (top figures) and ND (middle figures) of the crystallographic macroscopic frame, the corresponding key color and direct pole figures (bottom figures). Sample 1 is on the left whereas sample 21 is on the right. 17

Figure 7 :

 7 Figure 7: GOS and KAM maps for the two analyzed samples (1 and 21) and their corresponding color scale. GOS maps are top figures whereas KAM maps are bottom ones. Sample 1 is on the left whereas sample 21 is on the right. (d) Large vertical straight lines observed in the KAM map obtained for sample 21 only are artefacts resulting from the polishing process.

Table 1 :

 1 Chemical composition of Inconel R 600.

	Element	Cr Fe Si Mn	C	Cu	S	Ni
	wt.% (weight content) 15.5 8.0 0.5 1.0 0.15 0.5 0.01 bal.

Table 2 :

 2 Mechanical properties of Inconel R 600 identified experimentally from the tensile tests on single tubes. Samples 23 and 24 are mentioned in the table since their stress-strain curves are plotted in Figs.3(b) and 4(b), but E m and σ 0.2m could not be measured in those cases because of too noisy curves.

		Temperature Loading rate	Strain rate	E m (GPa) σ 0.2m (MPa)
	NHT Sample 1	RT	0.4 mm min -1	5.0×10 -5 s -1 205.0	321.0
	Sample 2	RT	47.0 mm min -1 5.7×10 -3 s -1 225.0	325.0
	Sample 3	600 • C	0.4 mm min -1	5.0×10 -5 s -1 202.0	218.0
	Sample 4	600 • C	47.0 mm min -1 5.9×10 -3 s -1 162.0	217.0
	Sample 5	800 • C	0.4 mm min -1	1.2×10 -4 s -1 112.0	123.0
	Sample 6	800 • C	4.5 mm min -1	1.1×10 -3 s -1 141.0	176.0
	Sample 7	800 • C	47.0 mm min -1 8.5×10 -3 s -1 155.0	167.0
	BHT Sample 8	RT	4.5 mm min -1	5.1×10 -4 s -1 270.0	239.3
	Sample 9	RT	4.5 mm min -1	5.0×10 -4 s -1 259.0	240.7
	Sample 10 600 • C	0.4 mm min -1	6.7×10 -5 s -1 180.0	172.4
	Sample 11 600 • C	0.4 mm min -1	6.7×10 -5 s -1 145.0	166.7
	Sample 12 600 • C	4.5 mm min -1	6.0×10 -4 s -1 186.0	157.1
	Sample 13 600 • C	4.5 mm min -1	6.0×10 -4 s -1 160.0	171.0
	Sample 14 600 • C	4.5 mm min -1	6.0×10 -4 s -1 204.0	167.1
	Sample 15 600 • C	47.0 mm min -1 6.2×10 -3 s -1 175.0	153.3
	Sample 16 800 • C	0.4 mm min -1	1.3×10 -4 s -1 119.0	116.1
	Sample 17 800 • C	0.4 mm min -1	1.3×10 -4 s -1 114.0	115.2
	Sample 18 800 • C	0.4 mm min -1	1.3×10 -4 s -1 115.0	103.9
	Sample 19 800 •

Table 3 :

 3 Mechanical properties of Inconel R 600 identified from the tensile tests on single tubes.

		RT		600 • C		800 • C
	E m (fixed) 214.0 GPa		180.0 GPa		164 GPa
	ν m (fixed) 0.31		0.31		0.31
	σ ym	223.1 MPa		145.7 MPa		114.6 MPa
	Q 1m	7.6 MPa		11.9 MPa		12.4 MPa
	b 1m (fixed) 500.0		500.0		500.0
	Q 2m	37.8 MPa		12.5 MPa		133.2 MPa
	b 2m (fixed) 50.0		50.0		50.0
	Q 3m	482.5 MPa		495.6 MPa		13.3 MPa
	b 3m (fixed) 5.0		5.0		5.0
	n m	56.6		56.6		3.9
	D m	606.2 MPa.s	1 nm	606.2 MPa.s	1 nm	1228.3 MPa.s	1 nm

Acknowledgements

This work was performed within the 'Aero-Thermodynamic Loads on Lightweight Advanced Structures II' project investigating high-speed transport. ATLLAS II, coordinated by ESA-ESTEC, was supported by the European Union within the 7th Framework Program Theme 7

Transport, Contract No.: ACP0-GA-2010-263913. Further information on ATLLAS II can be found at http://www..esa.int/techresources/atllas II. We are also indebted to our colleagues M. Thomas, P. Kanouté and S. Kruch for many fruitful discussions.


François François Boulier
email: boulier@lifl.fr

François Lemaire
email: lemaire@lifl.fr

Marc Moreno Maza

Franc ¸ois Boulier

¸ois Franc

Lemaire

PARDI !

Keywords: differential algebra, PDE, characteristic sets, change of rankings, gcd

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

INTRODUCTION

In this paper, we propose an algorithm which solves the following problem: given a characteristic set C of a prime differential ideal a w.r.t some ranking R and another ranking R = R, compute a characteristic set C of a w.r.t. R. The algorithm that we present, called 1 PARDI applies for systems of partial differential polynomial equations. It specializes to systems of ordinary differential polynomial equations and is then called 2 PODI. It specializes to nondifferential polynomial equations where it is called 3 PALGIE.

As far as we know, Ollivier was the first to solve the problem addressed in this paper. Let's quote [19, page 95]: "one can [design] a method for constructing a characteristic set of a finitely generated prime differential ideal as soon as one can effectively test membership to this ideal". An al- 1 PARDI is an acronym for Prime pARtial Differential Ideal. In French, "pardi" is an oldfashioned swearword such as, say, "egad" in English. 2 PODI is an acronym for Prime Ordinary Differential Ideal. 3 PALGIE is an acronym for Prime ALGebraic IdEal. However, since "algie" means "suffering" in French, one might also understand PALGIE as "polynomial suffering" say. gorithm is given in SCRATCHPAD in [19, page 97]. In all approaches, a known characteristic set provides the membership test algorithm. Such a problem was also considered in [START_REF] Boulier | Efficient computation of regular differential systems by change of rankings using Kähler differentials[END_REF]. However, the algorithms presented in [START_REF] Boulier | Efficient computation of regular differential systems by change of rankings using Kähler differentials[END_REF] compute differential polynomials which are not necessarily part of the desired characteristic set but only help computing it. They are complementary to PARDI.

The restriction to prime ideals is realistic. Indeed most differential systems coming from real problems generate differential prime ideals. Quite often, nondifferential polynomial systems in positive dimension either generate prime ideals or can be decomposed into prime ideals. Assuming that prime ideals are given by characteristic sets is realistic too. In the differential case, it happens quite often (e.g. dynamical systems in nonlinear control theory) that the input equations already form characteristic sets w.r.t. some rankings. Let us observe that one can decide primality for any ideal presented as the set of the polynomials pseudoreduced to zero by a characteristic set. The algorithm we propose generalizes to ideals which are not necessarily prime. However, for the reasons explained above and the legibility of the paper, we prefer to restrict ourselves to the prime case. Our algorithm applies to perform invertible changes of coordinates on the dependent and independent variables. Such changes realize ring isomorphisms between two differential rings φ : R → R, and one-to-one correspondences between the differential ideals of R and the ones of R. However the image C of a characteristic set C of a is usually not a characteristic set of the ideal a = φa and there is usually no ranking w.r.t. which a characteristic set of a could be easily deduced from C. The idea is then to apply PARDI over C but to test membership in a by performing the inverse changes of coordinates and testing membership in a using C.

Our approach offers several advantages. It identifies the algebraic subproblems which occur in the differential computations and solves them by a purely algebraic method. This improves the control of the coefficients growth and avoids many useless computations only due to differential considerations. This very important advantage w.r.t. all other approaches permits us to handle some unsolved problems, even using a preliminary implementation. The three variants were implemented: PARDI in MAPLE, PODI in C and PALGIE in MAPLE, C and Aldor. The application to the change of variables was implemented in MAPLE. A last contribution is the conceptual simplicity of our algorithm, which contrats with the high technicity of its implementation. Ev-erybody knows that the common roots of two univariate polynomials over a field are given by their gcd. Our algorithm applies this very simple idea and replaces any two univariate polynomials by one of their gcd over the fraction field of some quotient ring which makes much more sense than speaking of full remainders as in the previous approaches.

DEFINITIONS AND NOTATIONS

Let X be an ordered alphabet (possibly infinite). Let R = K[X] be a polynomial ring where K is a field. Let p ∈ R \ K be a polynomial. The leader of p, denoted ld p, is the greatest indeterminate x which occurs in p. The polynomial p can be written as p = a d x d + • • • + a1 x + a0 where d = deg(p, x) and the polynomials ai are free of x. The polynomial ip = a d is the initial of p. The rank of p is the monomial x d . The reductum of p is the polynomial pip x d . If x d and y e are two ranks then x d < y e if x < y or x = y and d < e. The separant of p is the polynomial sp = ∂p/∂x. Let A ⊂ R \ K be a set of polynomials. Then IA (resp. SA) denotes the set of the initials (resp. the separants) of its elements. We denote HA = IA ∪ SA. The set A is said to be triangular if its elements have distinct leaders. Let q be a polynomial. We denote prem(q, p) the pseudoremainder of q by p, both polynomials being viewed as univariate polynomials in the leader of p. We denote prem(q, A) "the" pseudoremainder r of q by all the elements of A i.e. any polynomial r obtained from q and the elements of A by performing successive pseudoreductions and such that prem(r, p) = r for every p ∈ A. Without further precisions, r is not uniquely defined. Fix any precise algorithm.

If A is a subset of a ring R then (A) denotes the ideal generated by A. Let a be an ideal of R. If S = {s1, . . . , st} then the saturation a : S ∞ of a by S is the ideal a :

S ∞ = {p ∈ R | ∃a1, . . . , at ∈ N s.t. s a 1 1 • • • s at t p ∈ a}.
Reference books for differential algebra are [START_REF] Ritt | Differential Algebra[END_REF] and [START_REF] Kolchin | Differential Algebra and Algebraic Groups[END_REF]. We also refer to the MAPLE VR5 and VI diffalg package and to [START_REF] Boulier | Representation for the radical of a finitely generated differential ideal[END_REF][START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF][START_REF] Petitot | Quelques méthodes de Calcul Formel appliquéees à l'étude des équations différentielles[END_REF][START_REF] Hubert | Factorization free decomposition algorithms in differential algebra[END_REF]. A derivation over a ring R is a map δ : R → R such that δ(a + b) = δa + δb and δ(a b) = (δa)b + a(δb) for every a, b ∈ R. A differential ring is a ring endowed with finitely many derivations which commute pairwise. The commutative monoid generated by the derivations is denoted by Θ. Its elements are the derivation operators θ = δ a 1 1 • • • δ am m where the ai are nonnegative integer numbers. The sum of the exponents ai, called the order of the operator θ, is denoted by ord θ. The identity operator is the unique operator with order 0. The other ones are called proper.

If φ = δ b 1 1 • • • δ bm m then θφ = δ a 1 +b 1 1 • • • δ am+bm m . If ai > bi for each 1 ≤ i ≤ m then θ/φ = δ a 1 -b 1 1 • • • δ am-bm m .
A differential ideal a of R is an ideal of R closed under derivation i.e. such that a ∈ a ⇒ δa ∈ a. Let A be a nonempty subset of R. We denote [A] the differential ideal generated by A which is the smallest differential ideal which contains A.

Let U = {u1, . . . , un} be a set of differential indeterminates. Derivation operators apply over differential indeterminates giving derivatives θu. We denote ΘU the set of all the derivatives. Let K be a differential field. The differential ring of the differential polynomials built over the alphabet ΘU with coefficients in K is denoted R = K{U }. A ranking is a total ordering over the set of the derivatives [12, page 75] satisfy-ing the following axioms 1. δv > v for each derivative v and derivation δ, 2. v > w ⇒ δv > δw for all v, w ∈ ΘU and derivation δ. Let us fix a ranking. The infinite alphabet ΘU gets ordered. Consider a polynomial p ∈ R \ K. Then the leader, initial, . . . of p are well defined. Axioms of rankings imply that the separant of p is the initial of every proper derivative of p. Let rankp = v d . A differential polynomial q is said to be partially reduced w.r.t. p if no proper derivative of v occurs in q. It is said to be reduced w.r.t. p if it is partially reduced w.r.t. p and deg(q, v) < d. A set A of differential polynomials is said to be differentially triangular if it is triangular and if its elements are pairwise partially reduced. It is said to be autoreduced if its elements are pairwise reduced. Autoreduced implies differentially triangular. If A is a set of differential polynomials and v is a derivative then Av = {p ∈ ΘA | ld p ≤ v}. Thus Rv denotes the set of the differential polynomials with leader less than or equal to v.

One distinguishes the partial reduction algorithm, which is denoted partial rem from the full reduction algorithm, denoted full rem. See [12, page 77]. Let q and p be two differential polynomials. The partial remainder partial rem (q, p) is the pseudoremainder of q by the (infinite) set of all the proper derivatives of p. The full remainder full rem (q, p) is the pseudoremainder of q by the set of all the derivatives of p (including p). Let A be a set of differential polynomials. We denote partial rem (q, A) and full rem (q, A) respectively the partial remainder and the full remainder of q by all the elements of A. Let v = ld q and A = A ∩ Rv. The partial remainder q of q by A is partially reduced w.r.t. all the elements of A and there exists a power product h of elements of S A such that h q ≡ q mod (Av). The full remainder q of q by A is reduced w.r.t. all the elements of A and there exists a power product h of elements of H A such that h q ≡ q mod (Av).

A pair {p1, p2} of differential polynomials is said to be a critical pair if the leaders of p1 and p2 are derivatives of some same differential indeterminate u (say ld p1 = θ1u and ld p2 = θ2u). Denote θ12 the least common multiple between θ1 and θ2. One distinguishes the triangular situation which arises when θ12 = θ1 and θ12 = θ2 from the nontriangular one which arises when θ12 = θ2 (say). In the last case, the critical pair is said to be a reduction critical pair. In the triangular situation, the ∆-polynomial ∆(p1, p2) is

∆(p1, p2) = s2 θ12 θ1 p1 -s1 θ12 θ2 p2.
In the nontriangular one,

∆(p1, p2) = prem(p2, θ2 θ1 p1).
If A is a set of differential polynomials then critical pairs(A) denotes all the critical pairs that can be formed with any two different elements of A and ∆(A) = {∆(p1, p2) | {p1, p2} ∈ critical pairs(A)}. A critical pair {p1, p2} is said to be solved by a system A = 0, S = 0 of differential polynomial equations and inequations if there exists a derivative v < θ12u such that ∆(p1, p2) ∈ (Av) : (S ∩ Rv) ∞ . Here is a criterion: if HA ⊂ S and full rem (∆(p1, p2), A) = 0 then the critical pair {p1, p2} is solved by A = 0, S = 0.

A differential system A = 0, S = 0 of a differential polynomial ring R is said to be a regular differential system w.r.t. a ranking R if (1) A is differentially triangular, (2) S contains the separants of the elements of A and is partially reduced w.r.t. A, (3) A is coherent 4 i.e. every critical pair in critical pairs(A) is solved by A = 0, S = 0.

A subset C of a differential ideal a is said to be a characteristic set of a if C is autoreduced and a contains no nonzero polynomial reduced w.r.t. C. If C is a characteristic set of a and HC contains no zero divisor in the quotient ring R/a then a = [C] : H ∞ C and p ∈ a ⇔full rem (p, C) = 0. This is the case when a is prime.

For algorithmic notations, we use pseudocode notation and an imperative programming style: we allow functions to modify some of their arguments. For this purpose, we put (as in ADA) the keywords in, out and in out in front of the formal parameters of the functions. The keyword in means that the formal parameter is readable only ; the keyword out means that the formal parameter is writable only ; the keyword in out means that the formal parameter is both readable and writable.

In the different algorithms, we often write the test p ∈ a where p is some differential polynomial. Since a characteristic set C of a is known, the test can translate as full rem (p, C) = 0. These full reductions can be very expensive in practice and can be replaced by more subtle tests using first some inexpensive criteria (e.g. full reductions after evaluating derivatives which are not derivatives of leaders over Z/nZ for some prime integer n).

THE OLD ALGORITHM

The Rosenfeld-Gröbner algorithm, which is implemented in the diffalg package of MAPLE VR5, solves the problem under consideration in this paper. See [START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF] for proofs and [START_REF] Boulier | Triangularisation de systèmes différentiels[END_REF] for the pseudocode of Rosenfeld-Gröbner. In the general case, Rosenfeld-Gröbner needs to split the solutions of the current system which annihilate some differential polynomial p from the solutions which do not annihilate p. The differential polynomial p is usually the initial or the separant of some differential polynomial used with Ritt's reduction. When the characteristic set C of a (prime) is known, only one branch needs being considered: the first one if p ∈ a ; the second one if p / ∈ a. The specialized Rosenfeld Gröbner function is a specialized pseudocode of Rosenfeld-Gröbner. It handles quadruples A, D, P, S where A is the set of the differential polynomial equations already processed, D is the set of the critical pairs to be processed, P is the set of the differential polynomial equations to be processed, S is the set of the differential polynomial inequations. Before stating the loop invariants, we need to give the definition of a critical pair nearly solved by a quadruple A, D, P, S . Two axioms are enough if the basic complete subfunction is used ; a third one is necessary if the advanced one is used.

1. every critical pair solved by A ∪ ∆(D) = 0, S = 0 is nearly solved by A, D, P, S , 2. every critical pair in D is nearly solved by A, D, P, S . 4 Same definition as that of Rosenfeld [START_REF] Rosenfeld | Specializations in differential algebra[END_REF] if HA ⊂ S.

The algebraic subproblems

We say that the specialized Rosenfeld Gröbner function encounters an "algebraic subproblem" when the differential polynomial p (let's denote it p = p2) provided to complete has the same leader as some differential polynomial p1 ∈ A. Let's look at complete to see how specialized Rosenfeld Gröbner behaves in that case. Assume rankp1 = v d 1 and rankp2 = v d 2 . Since p1 ∈ A and p2 is reduced w.r.t. A we have d2 < d1. The differential polynomial p2 is stored in A, the differential polynomial p1 is removed from A and the pair {p1, p2} is stored in D. After a few loops, specialized Rosenfeld Gröbner extracts the pair {p1, p2} from D and computes ∆(p1, p2). Since p1 and p2 have the same leader we have

∆(p1, p2) = prem(p1, p2, v).
Denote p3 this pseudoremainder and assume ld p3 = v. Then rankp3 = v d 3 with d3 < d2. Arguing as above we see that p3 is stored in A, that p2 is removed from A and that the pair {p2, p3} is stored in D. Let's summarize: the specialized Rosenfeld Gröbner function starts computing a very basic (bad) pseudoremainder sequence when it encounters an algebraic subproblem. p3 = prem(p1, p2, v), p4 = prem(p2, p3, v), . . .

Moreover, at each step i

• many critical pairs between the current pseudoremainder pi and the other elements of A are generated (not only the reduction pair with pi-1),

• the separant of pi (not only the initial) is stored in S.

The points above come from differential considerations. We shall see they are useless. The analysis done in this section holds with the advanced version of complete too.

THE NEW ALGORITHM

The PARDI algorithm identifies as such the algebraic subproblems arising during the differential treatment. It optimizes the computation of the pseudoremainder sequence and avoids completely the computations due to the pointless differential considerations described above. Assume again that the new differential polynomial p = p2 has the same leader as some p1 ∈ A. PARDI relies on two key ideas.

First, replacing both p1 and p2 by their "gcd" g in (R -/a -) [v] where

R -= K[w ∈ ΘU | w < v] and a -= a ∩ R -.
Actually g is the "lsr" i.e. the "last nonzero subresultant" of p1 and p2 in (R -/a -) [v]. It is also one of the gcds of p1 and p2 in G [v] where G denotes the field of fractions of R -/a -. For legibility however, we shall speak of the gcd of p1 and p2 in (R -/a -)[v] though this is slightly incorrect. We point out that the differential polynomial g always has positive degree in v (property (iii) of lsr1).

Second, applying a "master-student" relationship between C and A: every quantity which is zero in R -/a -(i.e. reduced to zero by C) but not reduced to zero by A is stored in the list P of the equations to be processed. Roughly speaking, when P is empty then the "student" A is able to reproduce all the computations first performed by the "master" C.

Computing the gcd (the lsr, sorry)

The key ideas are the following.

1. Starting from a (good) algorithm for computing a pseudoremainder sequence. We choose the algorithm of Lionel Ducos [START_REF] Ducos | Optimizations of the subresultant algorithm[END_REF] but the [START_REF] Lombardi | New structure theorem for subresultants[END_REF] algorithm would fit as well.

2. Verify at each step i that the leading coefficient of the current subresultant pi is nonzero in R -/a -.

If it is nonzero, just continue the [START_REF] Ducos | Optimizations of the subresultant algorithm[END_REF] algorithm. Do not try to normalize this leading coefficient in any sense (idea explained below).

If it is zero, restart the computation of the pseudoremainder sequence between the former subresultant pi-1 and reductum(pi).

3. Do not store in S the separants of the subresultants but only their initials. Do not generate critical pairs.

The first and the third points show the advantages of the new algorithm w.r.t. the old one: the growth of the coefficients is controlled with the efficient way of [START_REF] Ducos | Optimizations of the subresultant algorithm[END_REF]; the differential treatments are completely avoided.

Let's come back to the second point. This is a very simple but very important idea, already applied in [START_REF] Moreno Maza | On Triangular Decompositions of Algebraic Varieties[END_REF], which permits to perform the [START_REF] Ducos | Optimizations of the subresultant algorithm[END_REF] algorithm in the ring (R -/a -)[v].

The only difficult operation is the exact quotient computation in R -/a -. By not normalizing at all the leading coefficients of the subresultants we can perform it as if we were computing in R -[v], which is easy. For pseudoremainder sequence algorithms, the most convenient choice is to represent the residue classes of the coefficients with representatives which make the exact quotient operation easy to perform. The only precaution to take is to make sure that the leading coefficient of each subresultant pi is nonzero.

When it is, one theoretically could go on the [START_REF] Ducos | Optimizations of the subresultant algorithm[END_REF] algorithm after replacing pi by its reductum but the trick explained above would not work anymore: exact quotient operations would become very tricky to perform in R -/a -. Another much simpler possibility5 consists in restarting the [START_REF] Ducos | Optimizations of the subresultant algorithm[END_REF] algorithm over pi-1 and the reductum of pi. This is what we do. However let us observe that in the context of [START_REF] Moreno Maza | Polynomial gcd computations over towers of algebraic extensions[END_REF], an algorithm is given which performs efficiently exact quotients in the ring R -/a -by normalizing polynomials in the sense of [START_REF] Lazard | A new method for solving algebraic systems of positive dimension[END_REF]. The adaptation of that strategy to the general case does not seem to be efficient in practice [START_REF] Aubry | Triangular sets for solving polynomial systems: A comparative implementation of four methods[END_REF].

Let's consider the implementation. The code is a slight modification of [START_REF] Ducos | source of the axiom package pseudoremaindersequence[END_REF] and [START_REF] Moreno Maza | On Triangular Decompositions of Algebraic Varieties[END_REF]. Functions Lazard2 and nsr can be found in [8, Lazard2 and next sousResultant2 functions].

The parameters p and q satisfy the following properties: they both have the same leader v and belong to a. Their initials and separants do not belong to a. The parameters A, P and S are components of a quadruple as for specialized Rosenfeld Gröbner. In particular, the initials of the elements of A belong to S. Moreover, p and q are partially reduced w.r.t. A. The last nonzero subresultant g returned by the function satisfies the following property:

(i) the leading coefficient of g w.r.t. v does not lie in a.

All the subresultants computed from p and q belong to the ideal (p, q) in the ring (R -/a -) [v]. In other words, they belong to (p, q) + a -in R -[v] and, since p, q ∈ a (ii) g ∈ a (iii) ld g = v (the resultant of p and q belongs to a since it is one of the subresultants ; it belongs to R -; it is thus zero in (R -/a -)[v]).

From (i) and (iii) we see that the initial of g does not lie in a. Moreover, since the separants of p, q do not lie in a and since g is a gcd of p, q in G[v] where G denotes the field of fractions of R -/a -, we see that sg / ∈ a.

(iv) the initial and the separant of g do not lie in a.

The lsr1 function records in P all the coefficients which are zero in R -/a -but not reduced to zero by A. It stores in S the initials i1, . . . , in = ig (after making sure that they do not belong to a) of the successive computed subresultants. Denote j the ideal ((A ∪ P ) ∩ R -) : (S ∩ R -) ∞ where the variables A, P and S are considered at the end of the computation. We claim we have:

(v) (p, q) ⊂ (g) : (i1 • • • in) ∞ in (R -/j)[v]. Indeed (p, q) ⊂ (g) : (i1 • • • in) ∞ in (R -/a -)[v] (classical property).
It suffices thus to prove that every leading coefficient which belongs to a -belongs to j. These leading coefficients are either recorded in P by ensure lcoeff1 or reduced to zero by A. In the former case they belong to (P ∩ R -), in the latter case they belong to (A ∩ R -) : (S ∩ R -) ∞ for can be simplified when g = q. One could also perform full remainders instead of partial ones but this would make the [START_REF] Ducos | Optimizations of the subresultant algorithm[END_REF] algorithm less efficient.

function PARDI(in C, R, R) begin A, D, P, S := ∅, ∅, C, HC (HC is computed w.r.t. R) (from now on, the ranking R is implicitely used) while D = ∅ or P = ∅ do take and remove a new equation p ∈ P or a critical pair {p1, p2} ∈ D. In the latter case let p = ∆(p1, p2) p := partial rem (p, A) p := ensure rank(p, C, A, P ) if p = 0 then if ∃q ∈ A such that ld p = ld q then g := lsr1(p, q, C, A, P, S) A, D, P, S := complete( A \ {q}, D, P, S , g) (replaces p and q by their "gcd") else A, D, P, S := complete( A, D, P, S , p) fi fi od return specialized regCharacteristic(A, S) end

The following function simplifies the differential polynomial p while its initial and its separant lie in a.

Proof of the first invariant

We only give a sketched proof of the most interesting issue: proving that the differential ideal is not changed when the algorithm performs the purely algebraic treatment (call to lsr1). We consider the value of the quadruple A, D, P, S just after lsr1 is run. We claim we have a = [A ∪ {p} ∪ ∆(D) ∪ P ] : S ∞ . This relation held before calling lsr1. This function enlarges P with elements of a and enlarges S with elements which do not belong to a (i.e. non zerodivisors elements modulo a for the ideal is prime). This concludes the proof of the claim. 2 The different versions of the complete function preserve this property for g ∈ a and ig sg / ∈ a (properties (ii) and (iv) of lsr1). At the beginning of the next loop, the differential ideal under consideration is b = [(A \ {q}) ∪ ∆(D) ∪ P ∪ {g}] : S ∞ . We must prove that a = b. The inclusion b ⊂ a holds since g ∈ a. For the opposite inclusion a ⊂ b it suffices to prove that the removed differential polynomials p, q ∈ b. According to property (v) of lsr1 p, q ∈ (g) :

(i1 • • • in) ∞ in (R -/j)[v].
The ideal j ⊂ b for it is generated by polynomials of R -and the differential polynomials p, q do not belong to that ring. Therefore p, q ∈ (g) :

(i1 • • • in) ∞ in (R -/(b ∩ R -))[v].
Since g ∈ b and i1, . . . , in are recorded in S by lsr1 we have

(g) : (i1 • • • in) ∞ ⊂ b whence p, q ∈ b and a ⊂ b. 2

Termination proof

The rank of A decreases at each loop. This rank cannot strictly decrease infinitely many times (classical property [12, proposition 3, page 81] of autoreduced sets). It thus suffices to prove that this rank cannot be constant infinitely many times. This situation only arises if PARDI calls lsr1 and if g = q or if all the coefficients of the differential polynomial p (viewed as a univariate polynomial in its leader) picked from P lie in a. In both cases, the algorithm suppresses a differential polynomial of P (i.e. p) and enlarges P with finitely many (possibly none) differential polynomials having leader strictly less than that of p. This cannot happen infinitely many times (classical argument of graph theory6 [13, Satz 6.6] and the fact that [12, page 75] rankings are well orderings). Thus PARDI stops.

Variant of algorithm

It is interesting in practice to keep A as a regular chain in the sense [START_REF] Kalkbrener | A Generalized Euclidean Algorithm for Computing Triangular Representations of Algebraic Varieties[END_REF][START_REF] Aubry | On the theories of triangular sets[END_REF]. It allows to take an inequation into account as soon as it arises: either it is regular and A is left unchanged or it is not and A gets smaller. Let us observe that inequations need anyway to be kept until the end: this optimization does not avoid calling specialized regCharacteristic. Thus, in the worst case, some computation time is lost. However, memory consumption is the most important issue and this variant improves it.

EXAMPLES

The following system C is a characteristic set w.r. We are looking for a characteristic set C of the same prime differential ideal w.r.t. the elimination ranking R

• • • > ux > uy > u > • • • > vxx > vxy > vyy > vx > vy > v.
In the analysis below, we only give the ranks of the differential polynomials. On the diagrams, black circles give the leaders of the elements of A while integers give the degrees of the elements of A w.r. Now P is empty and rankD = {{u, ux}, {vxxx, v 2 xxy }}. The first critical pair is picked and removed from D. After full reduction, the ∆-polynomial has rank vxxy. This is an algebraic subproblem. The situation is very easy for one of the differential polynomials has degree 1 while the other has degree 2 and there is basically no gain w.r.t. the old algorithm here. The gcd (last nonzero subresultant if you prefer) of these polynomials is the one with degree 1. It replaces the degree 2 polynomial in A and generates a critical pair with the differential polynomial with rank vxxx. At the end of the computation of the gcd, the resultant of the two polynomials, which has rank v 4 xx as a differential polynomial but which is zero in the quotient ring (R -/(a ∩ R -))[vxx], was not reduced to zero by A (where Let us observe specialized Rosenfeld Gröbner could not carry this example out in the diffalg package in MAPLE VR5 (it cannot be carried out with diffalg in MAPLE VI because of a bug). Even if we could, the analysis would have been much more painful. Indeed, we have presented our example in terms of gcd and resultant. This makes much more sense than if we had presented it in terms of full remainders using Ritt's reduction.

R -= K[w ∈ ΘU | w < vxx w.r.t. R]). It

Euler's equations for a perfect fluid

Expressed as differential polynomials, Euler's equations for a perfect fluid in two dimensions are Σ 8 < :

v 1 t + v 1 v 1 x + v 2 v 1 y + px = 0, v 2 t + v 1 v 2 x + v 2 v 2 y + py = 0, v 1
x + v 2 y = 0.

There are three differential indeterminates v 1 , v 2 (components of the speed) and the pressure p. They depend on three independent variables x, y (space variables) and the time t. For some orderly ranking, Rosenfeld-Gröbner gets with nearly no computation the characteristic set

C 8 > > < > > : pxx = -2 v 2 x v 1 y -2 (v 2 y ) 2 -pyy, v 1 t = -v 2 v 1 y -px + v 2 y v 1 , v 1 x = -v 2 y , v 2 t = -v 1 v 2 x -v 2 v 2 y -py.
For the elimination ranking (p, v 1 ) ≫ degrevlex(v 2 ) with t > x > y the implementation of PARDI was able to compute a characteristic set of the differential prime ideal. This characteristic set cannot be written in this paper (the computer file is 600 kilobytes large). There are 7 equations involving more than 50 different derivatives. Intermediate computations took more than 500 megabytes on the computers of the UMS MEDICIS (GAGE Laboratory, of the École Polytechnique). We have rankC = {px, py, v 1 , v This example could not be previously solved by any other algorithm. A remaining challenge (given by Pommaret) consists in computing a characteristic set for a ranking eliminating v 1 and v 2 .

Conclusion

We have described a new algorithm for converting a characteristic set of a prime differential ideal from one ranking into another. As far as we know, it is the first time that algebraic subproblems are clearly identified during differential computations (i.e. before Rosenfeld's lemma applies). It is also the first time that advanced practical algorithmic methods coming from both polynomial algebra (gcd computation over quotient rings) and differential algebra (analogue of Buchberger's criteria) are merged.

Loop invariants 1 .

 1 a = [A ∪ ∆(D) ∪ P ] : S ∞ , 2. the set of ranks of A is autoreduced, 3. every critical pair in critical pairs(A) is nearly solved by A, D, P, S (useful for PDE only), 4. HA ⊂ S. function specialized Rosenfeld Gröbner(in C, R, R) begin A, D, P, S := ∅, ∅, C, HC (HC is computed w.r.t. R) (from now on, the ranking R is implicitely used) while D = ∅ or P = ∅ do take and remove a new equation p ∈ P or a critical pair {p1, p2} ∈ D. In the latter case let p = ∆(p1, p2) p := full rem (p, A) p := old ensure rank(p, C, P, S) if p = 0 then A, D, P, S := complete( A, D, P, S , p) fi od return specialized regCharacteristic(A, S) end At the end of the main loop, A, D, P, S = A, ∅, ∅, S is such that A is coherent (because of the emptyness of D and the third loop invariant). The specialized regCharacteristic function, described in appendix, partially reduces S w.r.t. A and makes the elements of A pairwise partially reduced. It thereby gets a regular differential system A = 0, S = 0 such that a = [A] : S ∞ and then computes the desired characteristic set C by a purely algebraic method. Remark that the specialized regCharacteristic function is a specialized version of [7, regCharacteristic] but that, in the different versions of the diffalg package the computation is performed with different (less efficient) algorithms based on Gröbner bases computations [6, 10]. The following function simplifies p as long as its initial or its separant lies in a. function old ensure rank(in p, C, in out P, S) begin while p / ∈ K and (ip ∈ a or sp ∈ a) do if ip ∈ a then P := P ∪ {ip} p := reductum(p) else P := P ∪ {sp} S := S ∪ {ip} p := d pv sp where v d = rankp fi od return p end The following complete function is basic. It enlarges A with the new equation p and removes from A the equations the leader of which is a derivative of the one of p (so that it preserves the second loop invariant of the calling function). It enlarges D with the critical pairs between p and the equations of A. It enlarges S with the initial and the separant of p. Let us observe that the equations removed from A now belong to reduction critical pairs whence this function does not modify the differential ideal, in the sense that [A ∪ {p} ∪ ∆(D) ∪ P ] : S ∞ = [A ′ ∪ ∆(D ′ ) ∪ P ′ ] : S ′∞ . function complete(in A, D, P, S , p) begin A ′ := {p} ∪ {q ∈ A | ld q is not a derivative of ld p} D ′ := D ∪ {{q, p} | q ∈ A and {q, p} is a critical pair} P ′ := P S ′ := S ∪ {ip, sp} return A ′ , D ′ , P ′ , S ′ end

  function ensure rank(in p, C, A, in out P ) begin while p / ∈ K and (ip ∈ a or sp ∈ a) do if ip ∈ a then if prem(ip, A) = 0 then P := P ∪ {ip} fi p := reductum(p) else if prem(sp, A) = 0 then P := P ∪ {sp} fi p := d pv sp where v d = rankp fi od return p end

  t. the orderly ranking R • • • > uxx > uxy > uyy > vx > vy > ux > uy > v > u. of the prime differential ideal [C] : H ∞ C in the ring Q{u, v} endowed with derivations w.r.t. x and y. Ranks appear on the left-hand side of the equal signs. The denominators of the right hand side are the initials of the differential polynomials. = ux, vy = (ux uy + ux uy u)/(4

P

  is empty, rankD = {{vxy, v 4 yy }, {vxx, vxy}, {vxxy, v4 xx }}. The three ∆-polynomials are reduced to zero by A. After specialized regCharacteristic is performed, one gets

  t. their leaders. Initially, A = D = ∅ and S = {4 u, 2 ux, 2 uy} contains HC w.r.t. R. We have P = C. The set of ranks of P (w.r.t. R from now on) is rankP = {ux, ux, u 2x , u 2 y }. The implemented version of PARDI is a draft version: it picks new equations from P first and critical pairs from D only if P = ∅ (the list D is not sorted). It was written in MAPLE VI by the second author. It slightly differs from the algorithm presented in this paper for it performs full remainders instead of partial remainders in the main loop. We chose this variant because it makes the analysis slightly shorter over that example.First, an equation with rank ux is picked and removed from P and stored in A.

		d/dy	d/dy
			2
		1	1
		u	v
		d/dx	d/dx
	d/dy	d/dy
	1
	u	v
	d/dx	d/dx
	Now rankP = {ux, u 2 x , u 2 y }. A new equation with rank u 2 x
	is picked and removed from P . After full reduction, its rank
	is u. It forms a critical pair with the former equation, which
	is withdrawn from A.
	d/dy	d/dy
	1
	u	v
	d/dx	d/dx
	Now rankP = {ux, u 2 y } and rankD = {{u, ux}}. A new equation with rank u 2 y is picked and removed from P . After full reduction, its rank is v 2 xxy .
	d/dy	d/dy
		2
	1
	u	v
	d/dx	d/dx
	Now rankP = {ux} and rankD = {{u, ux}}. A new equa-
	tion with rank ux is picked and removed from P . After full
	reduction, its rank is vxxx. It forms a critical pair with the
	former equation.

  is recorded in P . Now rankP = {v4 xx }, rankD = {{vxxx, v 2 xxy }, {vxxx, vxxy}}. The resultant which has rank v4 xx is picked and removed from P . Two new critical pairs are generated between this differential polynomial and two of the differential polynomials in A. These two differential polynomials are withdrawn from A. Implemented in the advanced version of complete, the analogue of Buchberger's second criterion cancels the two old pairs of D.Now P is empty and rankD= {{vxxx, v 4 xx }, {vxxy, v 4 xx }}.The first critical pair is picked and removed from D. After full reduction, the ∆-polynomial has rank v 3 xx . This is again an algebraic subproblem. The gcd between this differential polynomial and the element of A with rank v 4 xx has degree 1 in vxx. It replaces the differential polynomial with rank v 4 xx . During its computation no critical pair was generated and this is a very important gain w.r.t. specialized Rosenfeld Gröbner. At the end of its computation, the resultant which has rank v 4 xy as a differential polynomial but is zero in the quotient ring, was not reduced to zero by A. It is recorded in P . No critical pair is generated. Now rankP = {v 4 xy } and rankD = {{vxxy, v 4 xx }}. A new equation with rank v 4 xy is picked and removed from P . It is put in A and generates a new critical pair. The first critical pair is picked and removed from D. After full reduction by A, the ∆-polynomial has rank v 2 xy . This is a purely algebraic subproblem. The gcd between this differential polynomial and the element of A with rank v 4 xy has rank vxy. It replaces the degree 4 polynomial in A. It generates one critical pair. At the end of its computation, the resultant, which has rank v 4 yy is stored in P . The resultant is picked and removed from P . It is recorded in A. By the analogue of Buchberger's criterion only one new critical pair is generated (the pair {v 4 yy , vxx} is avoided).

	d/dy	d/dy
	1	1
	u	v
	d/dx	d/dx
	d/dy	d/dy
		4
	1	1
	u	v
	d/dx	d/dx
	P is empty and rankD = {{vxx, v 4 xy }, {vxxy, v 4 v 1 1 1 xx }}. u d/dy d/dy	u 1 d/dy	d/dx	v d/dy	1	1	d/dx
	d/dx	d/dx
	rankP = {v 4
			d/dy		d/dy
							two old critical
								pairs are cancelled
			1		4
			u		v
				d/dx				d/dx

yy } and rankD = {{vxx, vxy}, {vxxy, v 4 xx }}.

  2 xxxxt , v 2 xxxtt , v 2 xxytt , v 2 xxxyyt }. The diagram of the differential indeterminate v 2 is 7

	t
	y
	v2
	x

According to Lionel Ducos himself (we take this opportunity to thank him very much for his comments).

Every infinite tree which is locally finite involves a branch of infinite length.

The authors would like to thank Marc Giusti for his help.

A polynomial is said to be regular w.r.t. a triangular set {p1 < • • • < pn} if it belongs to none of the associated prime ideals of (p1, . . . , pn) : (ip 1 , . . . , ip n ) ∞ . A triangular set {p1 < • • • < pn} is said to be a regular chain if n = 1 or ip ℓ is regular w.r.t. {p1, . . . , p ℓ-1 } for each 2 ≤ ℓ ≤ n.

C is necessarily coherent at the end of the computation but not necessarily during the computation.

they are partially reduced w.r.t. A and the initials of the elements of A belong to S. They thus belong to j. This concludes the proof of the claim. 2 function lsr1(in p, q, C, A, in out P, S) begin v := ld p if deg(p, v) < deg(q, v) then swap p and q fi found := false while not found do δ := deg(p, v)deg(q, v) s := i δ q S := S ∪ {-iq} (p, q) := (q, prem(p, -q)) z := p rankfall := false while not found and not rankfall do q := ensure lcoeff1(q, C, A, P, rankfall) if q = 0 then found := true elif not rankfall then

(p, q) := (q, nsr(p, q, z, s))

The main function

The PARDI function handles quadruples A, D, P, S in the same way as specialized Rosenfeld Gröbner and maintains the same loop invariants. A difference (which simplifies proofs) with specialized Rosenfeld Gröbner: the elements of every critical pair generated by PARDI always have different leaders (i.e. there are no "algebraic" critical pairs). The code

Appendix

The advanced complete function is described in [START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF]. It applies an analogue of Buchberger's second criterion. Pseudocode can be found in [START_REF] Boulier | Triangularisation de systèmes différentiels[END_REF]. The specialized Rosenfeld Gröbner function is a variant of [START_REF] Boulier | Computing canonical representatives of regular differential ideals[END_REF]regCharacteristic] 3. C is a regular chain 8 in the sense of [START_REF] Kalkbrener | A Generalized Euclidean Algorithm for Computing Triangular Representations of Algebraic Varieties[END_REF][START_REF] Aubry | On the theories of triangular sets[END_REF].

At the end of the computation, C is a characteristic set of (A) : S ∞ whence it is a characteristic set 9 of [A] : S ∞ . Observe that C is not provided to specialized regCharacteristic as a parameter since it is not necessary in order to choose the right branch in the splitting tree. Splittings are avoided in regularize: when deg(g, v) = 0, a "pseudofactorization" of q is discovered i.e. for some nonnegative integer α we have i α g q = g pquo(q, g). Then g is necessarily the factor to throw away (the one which does not belong to a). Last, the loop could be optimized by performing pseudoquotients while the pseudoremainder is zero.

The parameters of lsr0 are two polynomials p, q having the same leader v and with leading coefficients w.r. In particular, the polynomial g is itself regular if it is the resultant of p and q. function lsr0(in p, q, in out C) begin v := ld p if deg(p, v) < deg(q, v) then swap p and q fi found := false while not found do if deg(q, v) = 0 then found := true z := q deg(p,v) else δ := deg(p, v)deg(q, v) s := i δ q (p, q) := (q, prem(p, -q)) z := p fi rankfall := false while not found and not rankfall do q := ensure lcoeff0(q, C, rankfall) if q = 0 then found := true else regularize(lcoeff(q, v), C) if not rankfall then δ := deg(p, v)deg(q, v) z := Lazard2(q, lcoeff(q, v), s, δ) if deg(z, v) = 0 then found := true else (p, q) := (q, nsr(p, q, z, s)


Robert Kalaba

Harriet Natsuyama

Firdaus Udwadia

AN EXTENSION OF GAUSS'S PRINCIPLE OF LEAST CONSTRAINT

Keywords: Mechanical systems, Gauss's principle of least constraint, constraint forces, Lagrange multiplier

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

INTRODUCTION

Lagrange and Gauss were among the first system scientists. They considered mechanical systems as being characterized by potential energy, kinetic energy and the constraint function, all of which refer to the system as a whole. Lagrange also emphasized the use of generalized coordinates to describe the current configuration. He focused attention on the need to consider constraints on the system, and introduced what are today called Lagrange multipliers. Constraints are what make a collection of point masses and rigid bodies into a system.

Lagrange proposed a fundamental principle for dealing with the constraints, the principle of virtual work: the constraint forces do no work on the system in a virtual displacement. This has worked well because practical mechanical systems, through design and the use of lubricants, minimize the effects of constraint forces that do work on a system. Using this principle, Gauss was able, in 1829, to establish his principle of least constraint, which provided an alternative variational principle for the motion of mechanical systems.

The aim of this paper is to produce an extension of Gauss's principle which is applicable in situations where constraint forces such as friction do work on the system in virtual displacements.

A key ingredient in the analysis is the employment of pseudoinverses of matrices, a tool that was not available to those early investigators. In Section 2, we sketch Lagrange's and Gauss's approach to general constrained motion of mechanical systems and state the result that they lead to for the actual acceleration vector (Eq. ( 8)). Section 3 goes on to state an extension of the principle of virtual work to cover cases in which there are constraint forces that do work on the system. Then in Section 4 the extended principle of virtual work is used to derive the new equation of motion (Eq. ( 33)). This makes it possible to provide an extension of Gauss's principle in Section 5, one which covers the case in which constraint forces that do work on the system in a virtual displacement are significant. The validity of the new principle is verified by showing that it leads to the same equation of motion, Eq. ( 44), as was established earlier in Eq. (33).

Equation ( 44) is the most general possible equation of constrained motion. In modeling a mechanical system the modeler may neglect constraint forces that do work on the system by putting c 0: Or they may be included by having c -0 and specified by Eq. ( 12). There are no remaining possibilities, as the derivation reveals.

CLASSICAL MODELS OF CONSTRAINED MOTION

Consider a discrete mechanical system whose generalized state vector is q, generalized velocity vector is _ q; and generalized acceleration is € q: These vectors are of dimension n. Then Lagrange's equations of motion take the form

M € q Q; ð1Þ
where M is an n by n positive definite symmetric matrix, and Q is an n-dimensional force vector. Suppose next that this system is subjected to m consistent equality constraints of the form f i ðt; q; _ qÞ 0; i 1; 2; . . .; m:

Upon differentiation with respect to t these equations assume the form

A€ q b; ð3Þ
where A is an m by n rectangular matrix, and b is an m-dimensional vector. Both may depend upon t, q, _ q: To maintain these constraints a constraint force Q c is required so that Eq. ( 1) becomes

M € q Q þ Q c ; ð4Þ
where the constraint force Q c is n-dimensional.
For the determination of the unknown vectors € q and Q c at each moment, 2n independent relations are required. Equation (4) provides n, and if the rank of A is r, Eq. ( 3) provides r more. A customary way of providing the additional ones is through the principle of virtual work. If we let v be an n-dimensional vector such that

Av 0; ð5Þ
where v is a so-called generalized virtual displacement, then the principle is that the constraint force Q c does no work in such a virtual displacement. This provides the additional n r conditions needed for the determination of both € q and Q c , which we shall actually see below.

In equation form this is: if v is a solution of Av 0; then

v T Q c 0: ð6Þ
Alternatively, Gauss, in 1829, suggested determining € q as the vector which minimizes the constraint function G, G ð€ q aÞ T Mð€ q aÞ; ð7Þ subject to the constraint A€ q b; where a is the acceleration that the system would have if there were no constraints. Both of these characterizations lead to the explicit formula for € q as

€ q a þ M 1=2 AM 1=2 þ ðb AaÞ; ð8Þ
where ðAM 1=2 Þ þ is the usual pseudoinverse of the matrix AM 1=2 : This formula [START_REF] Udwadia | Analytical Dynamics: A New Approach[END_REF] shows that the actual acceleration is made up of the free motion acceleration plus a perturbation term. This term is proportional to the vector b Aa; which indicates the extent to which the free motion acceleration a does not satisfy the constraint equation A€ q b: The proportionality matrix is

M 1=2 ðAM 1=2 Þ þ :
The formula in Eq. ( 8) is useful theoretically and especially computationally, for computing environments such as MATLAB contain commands for obtaining pseudoinverses.

INCORPORATION OF CONSTRAINT FORCES THAT DO WORK

But suppose that still other constraint forces are at work, forces that may do positive or negative work on the system in a virtual displacement. Their existence is beyond dispute. How may the standard framework be modified to account for them?

Let us simply replace Eq. ( 6) with the relation

v T Q c v T c; ð9Þ
where c is an n-dimensional vector which may depend upon t, q, and _ q: The vector c, a force, could, for example, represent the force of friction. The right side of the last equation, which may be positive, zero, or negative, can be chosen to account for the forces that do do work on the system in a virtual displacement. Now let us see what effect this has on both Gauss's principle and the explicit formula (8) for the actual generalized acceleration € q:

DERIVATION OF THE NEW EQUATION OF MOTION

The model for determining € q and Q c is now

M € q Q þ Q c ; ð10Þ
A€ q b; ð11Þ
and, for all v such that Av 0;

v T Q c v T c: ð12Þ
We introduce the notation

B AM 1=2 ; ð13Þ
€ r M 1=2 € q; ð14Þ
so that Eq. ( 11) becomes

B€ r b; ð15Þ
the most general solution of which is € r B þ b þ ðI B þ BÞw: ð16Þ
In the last equation w is an arbitrary n-vector. Equations ( 10) and ( 12) yield

v T ðM € q QÞ v T c; ð17Þ
for all v for which Av 0: This may be rewritten as

v T ðM € q Q cÞ 0; ð18Þ
or v T ðM 1=2 € r Q cÞ 0: ð19Þ
Making use of Eq. ( 16) it is seen that

v T M 1=2 B þ b þ ðI B þ BÞw È É Q c 0: ð20Þ
We may rewrite the equation Av 0 as

AM 1=2 M 1=2 v 0; ð21Þ
and put u M 1=2 v; ð22Þ
so that u T v T M 1=2 ð23Þ
and

Bu 0: ð24Þ
Equation (20) becomes u T B þ b þ ðI B þ BÞw M 1=2 Q M 1=2 c n o 0: ð25Þ
Since the vector u is in the null space of the matrix B, Bu 0; it is known that

u T B þ 0; ð26Þ
which means that u T w M 1=2 Q M 1=2 c n o 0: ð27Þ
Since the term in brackets is orthogonal to every vector u for which Bu 0; it must belong to the range space of the matrix B T , so that

w M 1=2 Q M 1=2 c B T z; ð28Þ
where z is an m-dimensional vector. Equation ( 16) becomes

€ r B þ b þ I B þ B M 1=2 Q þ M 1=2 c þ B T z h i : ð29Þ
But it is known that

B þ BB T B T ; ð30Þ
so that the last equation becomes

€ r B þ b þ I B þ B M 1=2 Q þ M 1=2 c h i : ð31Þ
It follows that

M 1=2 € q B þ b þ I B þ B M 1=2 Q þ M 1=2 c h i ; ð32Þ
which implies that € q a þ M 1=2 B þ ðb AaÞ þ M 1=2 ðI B þ BÞM 1=2 c; ð33Þ
where a M 1 Q: This is the desired equation of motion. It is the most general one possible, as the derivation shows. The first two terms on the right side have already been seen in Eq. ( 8).

The new term, due to the constraint forces that do work on the system in a virtual displacement, is the third one. Thus only a small change is needed to obtain the most general equation of constrained motion. Equation ( 31) is especially informative. The first term on the right side, B þ b, is a vector in the range space of the matrix B T and represents the contribution of the constraint forces that do no work on the system. The second term on the right is a vector that lies in the null space of the matrix B. It depends linearly on the generalized impressed force vector Q and the vector c, which characterizes the constraint forces that do work on the system. The first and second terms on the right side of Eq. ( 31) are, of course, orthogonal to each other, and one cannot compensate for the other.

Equations ( 31) and ( 33) are useful in dealing with Coulomb sliding friction. If there is a single particle, the vector c, assuming Cartesian coordinates, lies in the direction of _ q= q T q À Á 1=2 and has a magnitude that depends upon the vector B þ ðb BM 1=2 QÞ B þ ðb AaÞ and the coefficient of friction.

AN EXTENDED PRINCIPLE OF GAUSS

An extended principle of Gauss is ð€ q a M 1 cÞ T Mð€ q a M 1 cÞ min ð34Þ subject to the restriction

A€ q b: ð35Þ
Let us show that this also leads to Eq. ( 33). Introduce

y M 1=2 € q a M 1 c  à ; ð36Þ
so that the extended Gauss principle becomes

y T y min ð37Þ
subject to A M 1=2 y þ a þ M 1 c h i b ð38Þ
or

AM 1=2 y b Aa AM 1 c: ð39Þ
But the shortest length solution of the consistent linear algebraic equation system in Eq. ( 39) is

y AM 1=2 þ b Aa AM 1 c : ð40Þ
It follows from Eq. ( 36) that

€ q M 1=2 y þ a þ M 1 c: ð41Þ
Thus € q M 1=2 ðAM 1=2 Þ þ ðb Aa AM 1 cÞ n o þ a þ M 1 c: ð42Þ
But this can be written as

€ q a þ M 1=2 ðAM 1=2 Þ þ ðb AaÞ M 1=2 ðAM 1=2 Þ þ AM 1 c þ M 1 c; ð43Þ
€ q a þ M 1=2 ðAM 1=2 Þ þ ðb AaÞ þ M 1=2 I ðAM 1=2 Þ þ ðAM 1=2 Þ n o M 1=2 c; ð44Þ
which is again the equation of motion of the system subject to constraint forces that do no work on the system and to constraint forces that do work on the system in a virtual displacement.

DISCUSSION

In this paper, we have presented an extension of Gauss's principle of least constraint and have exhibited the equation of motion to which it leads, Eq. ( 44). In modeling a mechanical system of mass points and rigid bodies the modeler may put c 0; the classical choice, or choose c -0 and specified by the work it does in a virtual displacement.

Robert

  Kalaba is Professor of Biomedical Engineering, Economics, and Electrical Engineering at the University of Southern California. In addition to his teaching and research in Biomedical Engineering, he works closely with colleagues in the Medical School, the Business School, and the Urban Planning and Economics departments. Outside of the university, Professor Kalaba serves on the editorial boards of eight journals, including Applied Mathematics and Computation, Computers and Mathematics with Appli cations, and Information Sciences, and is a reviewer for many others. At present there are two major ingredients to Professor Kalaba's research program in biomedical engineering. The first involves modeling of constrained motion in living systems, based on a recent mathematical theory developed in collaboration with Professor F. Udwadia (Mechanical Engineering). The second is work on artificial neural networks. Here both Kohonen and Rumelhart networks are being used in studies of chemical equilibrium and the obtaining of initial approximations for system identification problems. Harriet H. Natsuyama (H.H. Kagiwada) was a research associate of Richard Bellman in dynamic programming and optimization at the Rand Corporation. She contributed to systems engineering as senior scientist at Hughes Aircraft Company and as Rockwell International Professor of Systems Engineering, California State University. Her Ph.D. in astrophysics and inverse problems is from Kyoto University, and she has degrees from the University of Hawaii in mathematics and physics. She has authored Terrestrial Radiative Transfer: Modeling, Computation, and Data Analysis [Springer Verlag, 1998] and other books in the systems sciences.

Firdaus

  E. Udwadia is a Professor at Departments of Civil Engineering, Aerospace & Mechanical Engineering (School of Engineering), Mathematics (School of Letters, Arts, and Sciences), and Information and Operations Management (Marshall School of Business), University of Southern California. His research interests include applied mechanics and applied mathematics, earthquake engineering, dynamics and control, engineering management and decision systems, and collaborative engineering design.