L Gallimard
email: laurent.gallimard@u-paris10.fr

A constitutive relation error estimator based on traction-free recovery of the equilibrated stress

Keywords: finite element method, equilibrated stress recovery, constitutive relation error, error bounds, goal-oriented error estimation

A new methodology for recovering equilibrated stress fields is presented, which is based on traction-free subdomains' computations. It allows a rather simple implementation in a standard finite element code compared with the standard technique for recovering equilibrated tractions. These equilibrated stresses are used to compute a constitutive relation error estimator for a finite element model in 2D linear elasticity. A lower bound and an upper bound for the discretization error are derived from the error in the constitutive relation. These bounds in the discretization error are used to build lower and upper bounds for local quantities of interest.

INTRODUCTION

The finite element (FE) discretization of a continuous mechanical model leads to a partial loss of the information contained in the continuous model and, thus, to the introduction of discretization errors. Methods have been developed over many years to evaluate the global quality of FE analyses [START_REF] Babuska | A posteriori estimates for the finite element method[END_REF][START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF][START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF]. For linear problems, all these methods provide a global energy-based estimate of the discretization error. Most of the time, such global information is insufficient for dimensioning purposes in mechanical design. In many common situations, the dimensioning criteria involve local quantities (stresses, displacements, intensity factors, etc.). The development of error estimators for such quantities was initiated in the 80s [START_REF] Babuska | The post-processing approach of the finite element method. Part II: the calculation of the stress intensity factor[END_REF][START_REF] Kelly | Procedures for residual equilibration and local error estimation in the finite element method[END_REF]. Recently, numerous works have been published which provide error estimates and bounds for several local quantities of interest [START_REF] Rannacher | A feedback approach to error control in finite element methods: application to linear elasticity[END_REF][START_REF] Peraire | Bounds for linear-functional outputs of coercive partial differential equations: local indicators and adaptive refinement[END_REF][START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors[END_REF][START_REF] Ladevèze | Local error estimators for finite element analysis[END_REF][START_REF] Ohnimus | Local error estimates of fem for displacements and stresses in linear elasticity by solving local Neumann problems[END_REF][START_REF] Prudhomme | Practical methods for a posteriori error estimation in engineering applications[END_REF][START_REF] Gallimard | Error estimation of stress intensity factors for mixed-mode cracks[END_REF]. In these works, the bounds for the quantities of interest are evaluated starting from the computation of an estimation of a global energy error.

As far as we know the error estimators based on the recovery of equilibrated stress (or equilibrated residuals) are the only type of estimators ensuring easily computable bounds for the error in linear analysis [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF][START_REF] Ohnimus | Local error estimates of fem for displacements and stresses in linear elasticity by solving local Neumann problems[END_REF][START_REF] Gallimard | Error estimation of stress intensity factors for mixed-mode cracks[END_REF][START_REF] Ladevèze | Error estimation and mesh optimization for classical finite elements[END_REF][START_REF] Machiels | A 'flux-free' nodal Neumann subproblem approach to output bounds for partial differential equations[END_REF][START_REF] Ainsworth | A Posteriori Error Estimation in Finite Element Analysis[END_REF][START_REF] Luce | A local a posteriori error estimator based on equilibrated fluxes[END_REF][START_REF] Parés | Subdomain-based flux-free a posteriori error estimators[END_REF][START_REF] Nicaise | An a posteriori error estimator for the Lamé equation based on H (div)conforming stress approximations[END_REF]. In this paper we consider the approach based on the concept of constitutive relation error (CRE) introduced in [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF]. One difficulty of this approach is the construction of equilibrated stress fields from the FE solution and the data. The classical technique of construction proposed [START_REF] Ladevèze | Error estimation and mesh optimization for classical finite elements[END_REF], which requires equilibrated traction recovery, is difficult to introduce in a standard FE code because the tractions are expressed on the edges of the elements of the mesh.

This paper proposes a new technique to recover equilibrated stress fields, which is easier to introduce in an FE code. This technique is an extension of the CRE estimators in the work of Machiels et al. [START_REF] Machiels | A 'flux-free' nodal Neumann subproblem approach to output bounds for partial differential equations[END_REF] on residual error estimators. We propose a new lower bound for the error in the constitutive relation, as well.

The paper is organized as follows: In Section 2, we briefly recall the basics of constitutive relation error estimators in linear elasticity. In Section 3 we introduce a new technique to build equilibrated stress fields. Then, in Section 4, we introduce the approach for obtaining a lower bound of the energy norm. In Section 5, a procedure to build the bounds for a quantity of interest is presented. Finally, numerical examples illustrating the behavior of the global and local bounds are presented in Section 6.

THE ERROR IN CONSTITUTIVE RELATION

The problem to be solved

Let us consider a 2D elastic structure defined in a domain bounded by * . The external actions on the structure are represented by a prescribed displacement u d over a subset * 1 of the boundary, a surface force density T defined over * 2 = * -* 1 , and a body force density b defined in . Hooke's operator of the material is denoted by K. Thus, the problem can be formulated as: Find a displacement field u and a stress field r defined in which verify

• the kinematic constraints:

u ∈ U and u | * 1 = u d (1)
• the equilibrium equations:

r ∈ S and ∀u * ∈ U 0 -r : ε(u * ) d + b.u * d + * 2 T.u * d = 0 (2)
• the constitutive relation:

r = Kε(u) (3)
U is the space in which the displacement field is being sought, S = L 2 [ ] 3 the space of the stresses, U 0 the space of the fields in U which are zero on * 1 , and ε(u) denotes the linearized deformation associated with the displacement

[ε(u)] i j = 1 2 (u i, j +u j,i ) (4)
The kinematic's constraints are taken homogeneous in the developments for the sake of simplicity (i.e. u d = 0). The weak form of the problem reads: find u ∈ U 0 such that

-Kε(u) : ε(u * ) d + b.u * d + * 2 T.u * d = 0 ∀u * ∈ U 0 ( 5
)
The Galerkin FE method provides an approximation u h to u in an FE space U h ⊂ U. The finite-dimension space U h is associated with an FE mesh of characteristic size h. Let P h denote a partition of into elements E k . This partition is assumed to verify = k E k (i.e. the mesh cover the whole domain) with E i ∩ E j = ∅ for any i different from j. The discretized problem is:

find u h ∈ U h 0 such that -Kε(u h ) : ε(u * h ) d + b.u * h d + * 2 T.u * h d = 0 ∀u * h ∈ U h 0 ( 6
)
where

U h 0 = {u h ∈ U h ; u h | * 1 = 0}.
The corresponding stress field is calculated using the constitutive relation

r h = Kε(u h ) (7)
The discretization error e h is the difference between the FE solution for the displacement and the actual solution of the problem defined by Equation ( 5)

e h = u-u h (8)
Traditionally, the energy norm is used as the measure of the error

e h = u-u h u = r-r h ( 9
)
where

u u = ε(u) : Kε(u) d 1/2 and r = r : K -1 r d 1/2 (10)

Definition of the CRE

The approach based on the CRE relies on a partition of the equations of the problem to be solved into two groups [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF]. In linear elasticity, the first group consists of the kinematic constraints (1) and the equilibrium equations (2); the constitutive relation (3) constitutes the second group. Let us consider an approximate solution of the problem, denoted by ( û, r), which verifies the first group of equations:

• the field û verifies (1) (i.e. û is kinematically admissible),

• the field r verifies (2) (i.e. r is statically admissible).

The fields ( û, r) are said to be an admissible solution. If ( û, r) verifies the constitutive relation (3) in , then ( û, r) = (u, r) (i.e. the exact solution is found). If, however, ( û, r) does not verify the constitutive relation, the quality of this admissible solution is measured by the CRE error ê( û, r), which is defined with respect to the constitutive relation

ê( û, r) = ( r-Kε( û)) : K -1 ( r-Kε( û)) d 1/2 = r-Kε( û)
A relative error ε is defined by

ε( û, r) = ê( û, r) [ r * K -1 r * d ] 1/2 with r * = 1 2 ( r+Kε( û))

Application of the CRE to FE computations

A key point to develop a CRE estimator is the construction of an admissible solution ( ûh , rh ) from the FE solution (u h , r h ) and the data.

• Since the FE displacement field verifies the kinematic constraints, one takes

ûh = u h in
• However, the stress field r h does not verify the equilibrium equations [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF]. A method to recover a statically admissible stress field rh from r h and the data have been under development for several years [START_REF] Ladevèze | Local error estimators for finite element analysis[END_REF][START_REF] Ladevèze | Error estimation and mesh optimization for classical finite elements[END_REF][START_REF] Gallimard | An enhanced error estimator on constitutive relation for plasticity problems[END_REF][START_REF] Florentin | Evaluation of the local quality of stresses in 3d finite element analysis[END_REF]. The main features of the recovery method are described in Section 3.1.

The CRE associated with the admissible solution (u h , rh ) is denoted as ê(u h , rh )

ê(u h , rh ) = ( rh -Kε(u h )) : K -1 ( rh -Kε(u h )) d 1/2 = rh -Kε(u h )

Upper bound estimate of the discretization error in the energy norm

An essential property is that the CRE estimator ê(u h , rh ) is an upper bound of the discretization error measure e h . The Prager-Synge theorem [START_REF] Prager | Approximation in elasticity based on the concept of functions space[END_REF] leads to the following inequality:

e h ê(u h , rh ) (11)
The proof of this property is easily obtained by introducing the exact solution r = Kε(u) in the expression of the CRE estimator

ê(u h , rh ) 2 = rh -Kε(u h ) 2 = rh -r+Kε(u)-Kε(u h ) 2 = rh -r 2 + Kε(u)-Kε(u h ) 2 +2 ( rh -r) : ε(u-u h ) d = rh -r 2 + u-u h 2 u +2 ( rh -r) : ε(e h ) d ( 12
)
As rh is an equilibrated stress, and e h belongs to U 0 , the following relation holds:

rh : ε(e h ) d = r : ε(e h ) d
and hence

( rh -r) : ε(e h ) d = 0 ( 13
)
The proof is completed by introducing Equation [START_REF] Ladevèze | Error estimation and mesh optimization for classical finite elements[END_REF] and the definition of discretization error in the energy norm in Equation ( 12)

ê(u h , rh ) 2 = rh -r 2 +e 2 h (14)

RECOVERY OF AN EQUILIBRATED STRESS FIELD

Standard method

This method has been first introduced in [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF], and further developed in several papers [START_REF] Ladevèze | Local error estimators for finite element analysis[END_REF][START_REF] Ladevèze | Error estimation and mesh optimization for classical finite elements[END_REF][START_REF] Gallimard | An enhanced error estimator on constitutive relation for plasticity problems[END_REF][START_REF] Florentin | Evaluation of the local quality of stresses in 3d finite element analysis[END_REF]. It is based on a prolongation condition used to link rh to the FE stress r h

∀E ∈ P h ∀i ∈ I E ( rh -r h ) grad ( i ) dE = 0
where E denotes an arbitrary element, I is the set of the mesh nodes, and i the FE scalar shape function associated with the node i. It involves two steps:

• The first step consists of constructing, on the element edges, surface force densities Th that represent the vector fields rh n E as

[ rh n E ] | = E Th with ∈ *E
where E is a function that is constant on each edge and whose value is either 1 or -1, so that on the common edge of two adjacent elements E and E : E + E = 0. Moreover, these force densities are generated in such a way that the volume load b and the surface loads E Th are in equilibrium on each element E of the mesh. More details of this procedure can be found in [START_REF] Ladevèze | Error estimation and mesh optimization for classical finite elements[END_REF]. • The second step consists of constructing on each element E a stress field rE = rh |E solution of the equilibrium equations:

div rE +b = 0 in E rE n E = E Th on *E (15)
For a given set of surface force densities, among all the solutions of ( 15) the best field rE is the one that is the solution to the minimization problem:

min rE verifying(15) 1 2 E ( rE -r h ) : K -1 ( rE -r h ) dV
By duality, this is equivalent to seeking a displacement field u E ∈ U(E) such that

∀u * ∈ U(E) E Kε(u E ) : ε(u * ) dE - E bu * dE - *E E Th u * d = 0 ( 16
)
where U(E) designates the space of the restriction to E of the fields in U. The stress field rE is given by

rE = Kε(u E )
Thus, one can obtain an approximation of rE by solving the problem ( 16) by a classical FE method on E. In practice, it is sufficient, in order to obtain a good approximation, to consider either a discretization of E with a single element but an interpolation of degree p +k, where p is the degree of the interpolation used in the FE method and k a positive integer, or a subdivision of the element E along with an interpolation of degree p. A complete study can be found in [START_REF] Strouboulis | The Finite Element Methods and Its Reliability[END_REF]. However, provided the loading is sufficiently smooth, it is possible to recover rigorously equilibrated stresses [START_REF] Ladevèze | Error estimation and mesh optimization for classical finite elements[END_REF][START_REF] Nicaise | An a posteriori error estimator for the Lamé equation based on H (div)conforming stress approximations[END_REF].

The main drawback of this technique is to require, for the construction of the surface force densities Th , a data structure that is not natural in a classical FE code (i.e. nodes and integration points defined on the edges of the elements). In the next subsection we present a method to build an equilibrated stress field that requires only the knowledge of the classical integration points and the nodes of the mesh.

Method based on traction-free recovery

This method is based on the partition of unity and uses local subdomains different from the mesh elements. It has been used for the residual error estimators [START_REF] Machiels | A 'flux-free' nodal Neumann subproblem approach to output bounds for partial differential equations[END_REF][START_REF] Parés | Subdomain-based flux-free a posteriori error estimators[END_REF][START_REF] Morin | Local problems on stars: a posteriori error estimators, convergence and performance[END_REF][START_REF] Choi | Adaptive computations of a posteriori finite element output bounds: a comparison of the 'hybrid-flux' approach and the 'flux-free' approach[END_REF]. The principles of this method, associated with a new prolongation condition, can be used to recover statically admissible stress fields without computing the surface loads E Th on the edges of the elements of the mesh. Let I v be the set of vertices of P h and i the corresponding piecewise linear (or bilinear for square elements) shape functions. The support of i is the set of the mesh elements that contains the vertice i and is denoted by i (see Figure 1). The sum of the functions k i is a partition of unity

i∈I v i = 1 (17)
Now, the statical stress field rh is sought as a sum of stress fields ri h = i rh computed in each subdomain i such that ri h vanishes ini . Moreover, rh is linked to the FE stress r h by the

(( rh -r h ) grad i ) u * d = 0 ∀u * ∈ U 0 ( i ) (18)
where U 0 ( i ) is the local restriction of the set U 0 on the support i . Using Equation ( 17), the following equality holds:

rh = i∈I v i rh = i∈I v ( i rh ) = i∈I v ri h
and ri h satisfies the following equilibrium equation in i :

div ri h = div ( i rh ) = i div rh + rh grad i = -i b+ rh grad i ( 19
)
as well as the following boundary conditions on * i :

ri h n = i T on * i 2 ri h n = 0 on * i -(* i 1 ∪* i 2 ) ( 20
)
where

* i 2 = * i ∩* 2 , * i 1 = * i ∩* 1
The weak solutions of the problem defined by Equations ( 19) and ( 20) are the fields ri

h verifying i ri h ε(u * ) d = i ( i b-rh grad i ) u * d + * i 2 i Tu * d ∀u * ∈ U 0 ( i ) (21)
Note that if * i 2 = ∅ the corresponding term disappears from the equation. Introducing the prolongation condition [START_REF] Nicaise | An a posteriori error estimator for the Lamé equation based on H (div)conforming stress approximations[END_REF] in Equation ( 21), one obtains

∀u * ∈ U 0 ( i ) i ri h ε(u * ) d = i ( i b-r h grad i )u * d + * i 2 i Tu * d = i ( i r h ε(u * )-r h : ε( i u * )) d + i b( i u * ) d + * i 2
Equation ( 22) is solvable on i if, for any rigid body motion u rb of i , the right-hand term is equal to zero

- i r h : ε( i u rb ) d + i b( i u rb ) d + * i 2 T( i u rb ) d = 0 ( 23
)
As r h satisfies the equilibrium equations in the FE sense, Equation ( 23) is satisfied if i u rb belongs to the FE space U h 0 . If i u rb does not belong to the FE space U h 0 (for instance, for linear triangular elements), it is necessary to add forces F i (M) to the loading such that the right-hand term is equal to zero for any rigid body motion of i and that i∈I v F i (M) = 0. This point will be developed in the Appendix.

Practical implementation

The construction of an admissible stress rh requires the resolution of local problems [START_REF] Strouboulis | The Finite Element Methods and Its Reliability[END_REF]. The stress field ri h is sought as the sum of the FE stress field and a correction stress field ri h = i r h +Kε(d ûi )

Equation ( 22) becomes: find

d ûi in U 0 ( i ) such that ∀u * ∈ U 0 ( i ) i Kε(d ûi ) : ε(u * ) d = i ( i b-r h grad i )u * d - i i r h ε(u * ) d + * i 2 i Tu * d (24)
In practice the resolution is performed in the same way as problem [START_REF] Luce | A local a posteriori error estimator based on equilibrated fluxes[END_REF] in the second step of the standard method (see Section 3.1). An approximation of d ûi is computed by solving the problem (24) by a classical FE method on i . A subdivision of the initial FE mesh on i is introduced: P h,S ( i ) where S is the subdivision per edge in all directions. Let U h,S ( i ) be the associated FE space, U h,S 0 ( i ) the space of the fields in U h,S ( i ) which are equal to 0 on * 1 . The problem [START_REF] Choi | Adaptive computations of a posteriori finite element output bounds: a comparison of the 'hybrid-flux' approach and the 'flux-free' approach[END_REF] becomes: find

d ûi ∈ U h,S 0 ( i ) such that ∀u * ∈ U h,S 0 ( i ) i Kε(d ûi ) : ε(u * ) d = i ( i b-r h grad i )u * d - i i r h ε(u * )+ * i 2 i Tu * d ( 25
)
The main advantage of the traction-free recovery method over the standard method is the simplicity of its implementation. The traction-free error estimator requires only one additional information to the standard data structure of an FE code. This information is the set of the elements that surround a given vertex node, and is easily determined from the connection table. Moreover, it does not require the computation of the flux along the edges of the elements. By comparison the standard method needs a data structure to express the flux on the edges of the elements and a connection between the edges, the elements and the nodes of the mesh, as well as integration points and weight on these edges. Finally, the boundary conditions of the local problems are easy to introduce in the standard structure of an FE code.

LOWER BOUND

Definition of the lower bound

Let us consider a displacement field ũh that satisfy the kinematic constraints (Equation ( 1)). Then, as the admissible stress field built to compute the upper bound rh satisfies the equilibrium equation ( 2)

rh : ε( ũh -u h ) d = r : ε( ũh -u h ) d
By subtracting on each side of the equation the quantity r h : ε( ũh -u h ) d and by using the Cauchy-Schwarz inequality, we obtain the following relation:

( rh -r h ) : ε( ũh -u h ) d r-r h ũh -u h u = e h ũh -u h u ( 26
)
where • and • u are the classical energy norms defined in Equation [START_REF] Ohnimus | Local error estimates of fem for displacements and stresses in linear elasticity by solving local Neumann problems[END_REF]. From Equation ( 26) a lower bound ẽ( ũh -u h , rh -r h ) of the discretization error e h can be derived

ẽ( ũh -u h , rh -r h ) e h ( 27
)
with

ẽ( ũh -u h , rh -r h ) = | ( rh -r h ) : ε( ũh -u h ) d | ũh -u h u

Practical construction of the lower bound

To compute the lower bound, it is sufficient to construct an admissible displacement field ũh . This displacement field can be decomposed as the sum of the FE displacement field and a correction displacement field ũh , where ũh belongs to U 0 :

ũh = u h + ũh
Following the method proposed in Section 3.2 to compute the admissible stress rh , ũh is sought as a sum of displacement fields ũi h computed on each subdomain i ũh =

i∈I v ũi h
To ensure the continuity of ũh , a simple way is to search ũi h such that ũi h is equal to zero on * i . Then ũi h is built in the same FE space as ûi , and ũi h is the solution of the following FE problem defined on i : find ũi

h ∈ U h,S 0,* i ( i ) such that ∀u * ∈ U h,S 0,* i ( i ) i Kε( ũi h ) : ε(u * ) d = i ( i b-r h grad i ).u * d - i i r h : ε(u * ) d ( 28
)
where U h,S 0,* i ( i ) is the space fields in U h,S 0 ( i ), which are equal to zero on * i . This FE problem is similar to the FE problem solved to compute the stress field ri h [START_REF] Strouboulis | A posteriori estimation and adaptative control of the error in the quantity of interest. Part 1: a posteriori estimation of the error in the von Mises stress and the stress intensity factor[END_REF], but the Neumann boundary conditions have been replaced by Dirichlet boundary conditions.

APPLICATION TO THE COMPUTATION OF A GOAL-ORIENTED ERROR ESTIMATOR

In this section the global upper and lower bounds in the energy norm (defined, respectively, in Sections 2 and 4) are used to obtain bounds for quantities of interest. We briefly recall the techniques developed in the literature for the case where the quantity of interest is a linear functional L of the displacement. In this case, the objective of the calculation is to assess the quality of I h = L(u h ) by estimating |I -I h |, where I = L(u) (an example of quantity of interest is proposed in Section 6.2). We refer the reader to References [START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors[END_REF][START_REF] Ladevèze | Local error estimators for finite element analysis[END_REF][START_REF] Ohnimus | Local error estimates of fem for displacements and stresses in linear elasticity by solving local Neumann problems[END_REF][START_REF] Prudhomme | Practical methods for a posteriori error estimation in engineering applications[END_REF][START_REF] Strouboulis | A posteriori estimation and adaptative control of the error in the quantity of interest. Part 1: a posteriori estimation of the error in the von Mises stress and the stress intensity factor[END_REF] for a detailed description of this approach and to [START_REF] Beckers | An optimal control approach to a posteriori error estimation in finite element method[END_REF][START_REF] Heintz | On adaptive strategies and error control in fracture mechanics[END_REF][START_REF] Ruter | Goal-oriented a posteriori error in linear elastic fracture mechanics[END_REF][START_REF] Gallimard | Evaluation of the local quality of the von mises's stress and l2-norm of the stress[END_REF] for the extension to non-linear quantities of interest.

Owing to the linearity assumption, one has

I -I h = L(u)-L(u h ) = L(u-u h ) = L(e h )
thus, the estimate of |I -I h | is equivalent to the estimate of |L(e h )|.

Definition of the auxiliary problem

Let us consider the following auxiliary problem: find u aux ∈ U 0 and r aux = Kε(u aux ) such that

Kε(u * ) : ε(u aux ) d = L(u * ) ∀u * ∈ U 0 (29)
Replacing u * by e h , one gets

L(e h ) = Kε(e h ) : ε(u aux ) d (30)
The function u aux indicates how the discretization error affects the quantity L(e h ). If u aux could be exactly computed, one could determine L(u) directly from the input data, as from (2), [START_REF] Gallimard | Evaluation of the local quality of the von mises's stress and l2-norm of the stress[END_REF] and because u aux belongs to U 0 , one would have

L(u) = b.u aux d + * 2

T.u aux d

Unfortunately, the problem for the function u aux is as complicated to solve as the reference problem for the solution u. However, we can compute an approximate value for u aux .

The relation [START_REF] Peraire | Bounds for linear-functional outputs of coercive partial differential equations: local indicators and adaptive refinement[END_REF] is the starting point of the goal-oriented error estimators developed in [START_REF] Kelly | Procedures for residual equilibration and local error estimation in the finite element method[END_REF][START_REF] Rannacher | A feedback approach to error control in finite element methods: application to linear elasticity[END_REF][START_REF] Peraire | Bounds for linear-functional outputs of coercive partial differential equations: local indicators and adaptive refinement[END_REF][START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors[END_REF][START_REF] Ladevèze | Local error estimators for finite element analysis[END_REF][START_REF] Ohnimus | Local error estimates of fem for displacements and stresses in linear elasticity by solving local Neumann problems[END_REF]. The approaches proposed in these papers differ on three points:

• the technique developed to approximate u aux , • the technique used to obtain lower and upper bounds of |L(e h )|,

• the error estimators used to calculate the bounds.

Approximate solution of the auxiliary problem

Let u aux h be the FE approximation of u aux defined in the FE space U h 0 :

Kε(u * h ) : ε(u aux h ) d = L(u * h ) ∀u * h ∈ U h
Because of the orthogonality property, one has

Kε(e h ) : ε(u * h ) d = 0 ∀u * h ∈ U h 0
Replacing u * h by u aux h the above relation transforms to

Kε(e h ) : ε(u aux h ) d = 0 (32)
By combining ( 30) and (32), one obtains the following relation for the error estimate:

I -I h = L(e h ) = Kε(e h ) : ε(e aux h ) d ( 33
)
where e aux h = u aux -u aux h is the discretization error of the auxiliary problem.

Bounds of the quantity of interest

Following [START_REF] Prudhomme | Practical methods for a posteriori error estimation in engineering applications[END_REF], the error in the quantity of interest can be written as

I -I h = Kε(e h ) : ε(e aux h ) d = 1 4 se h +s -1 e aux h 2 u - 1 4 se h -s -1 e aux h 2 u ( 34
)
where s is a scaling factor. The value of s is chosen so that se h u and s -1 e aux h u have the same amplitude. This implies that s = e aux h u / e h u However, in practice, the discretization errors e aux h u and e h u will be replaced by the computed error estimators in the constitutive relation.

Replacing e h by (u-u h ) and e aux

h by (u aux -u aux h ) in Equation (34) one easily obtains

I -I h = 1 4 u + -u + h 2 u -1 4 u --u - h 2 u ( 35
)
where

u ± = su±s -1 u aux and u ± h = su h ±s -1 u aux h
In order to determine bounds for the error in the quantity of interest, the error estimates of both the initial and the auxiliary problem are necessary. The methodology proposed in Sections 3 and 4 can be used to build a statically admissible stress field raux h and a kinematically admissible displacement field ũaux h for the auxiliary problem. Taking advantage of the linearity of the problems, and introducing the global lower and upper bounds defined, respectively, in Equations ( 27) and [START_REF] Prudhomme | Practical methods for a posteriori error estimation in engineering applications[END_REF], bounds for u + -u + h u and u --u - h u can be built ẽ( ũ+

h -u + h , r+ h -r + h ) u + -u + h u ê(u + h , r+ h ) ẽ( ũ- h -u - h , r- h -r - h ) u --u - h u ê(u - h , r- h ) ( 36
)
where

r+ h = s rh +s -1 raux h and r- h = s rh -s -1 raux h ũ+ h = s ũh +s -1 ũaux h and ũ- h = s ũh -s -1 ũaux h
From Equations ( 34) and ( 36), the classical upper and lower bounds are derived [START_REF] Prudhomme | Practical methods for a posteriori error estimation in engineering applications[END_REF] low I -I h upp (37

)
where

low = 1 4 ẽ( ũ+ h -u + h , r+ h -r + h ) 2 -1 4 ê(u - h , r- h ) 2 upp = 1 4 ê(u + h , r+ h ) 2 -1 4 ẽ( ũ- h -u - h , r- h -r - h ) 2
The bounds for the quantity of interest I are easily derived from (37)

I upp h = I h + upp and I low h = I h + low
The introduction of lower bounds in the energy ẽ( ũ+

h -u + h , r+ h -r + h ) and ẽ( ũ- h -u - h , r- h -r - h
) leads to much sharper bounds than the classical approach based on the error estimator in the constitutive relation, where the 'local' lower and upper bounds are constructed with only the CRE estimator computed on the initial and auxiliary problem.

CRE low = -1 4 ê(u - h , r- h ) 2 CRE upp = 1 4 ê(u + h , r+ h ) 2
In this case, the bounds for the quantity of interest I are

I upp,CRE h = I h + CRE upp and I low,CRE h = I h + CRE low 6. NUMERICAL EXAMPLES
This section is devoted to the study of the error estimates and bounds presented above. In the first example, we analyze the behavior of the global lower and upper bounds in the energy norm on a plate with two holes. The next examples deal with the study of the goal-oriented error estimation for a crack opening problem. In all the cases studied the exact solutions are unknown. The estimates and the error bounds are then verified by comparing results with solutions obtained on fine meshes, which are believed to provide accurate solutions. The quality of the global CRE estimator ê(u h , rh ) and of the global lower bound ẽ( ũh -u h , rh -r h ) is estimated by an effectivity index:

u = ê(u h , rh ) e ref and L = ẽ( ũh -u h , rh -r h ) e ref
where e ref = e h when the exact error is available and e ref is computed from a refined mesh when the exact error is not available.

Behavior of the global energy error

In this section, the study of the effectivity indexes and the convergence rate of the global errors estimators ẽ and ê shows the quality of the global lower and upper bounds presented in Sections 2 and 4. Two structures are studied. The first structure is a thin elastic plate with two holes proposed in [START_REF] Peraire | Bounds for linear-functional outputs of coercive partial differential equations: local indicators and adaptive refinement[END_REF], this is an in-plane stress linear elastic problem loaded with an unit tension along the vertical edges (see Figure 2). The symmetry of the problem allows us to study only one-fourth of the plate. The structure is meshed with 6-nodes triangular elements (as shown in Figure 3). The second structure is a holed plate submitted to bending. The plate clamped on the left edge and loaded with an uniform load is depicted in Figure 4. The initial coarse mesh (6-nodes triangular elements) is shown in Figure 5. The material properties are E = 100 GPa and = 0.3. Owing to the lack of analytical solution for these problems, we use a reference solution computed on a fine mesh. The reference solution is then used to compute the reference error e ref .

The behavior of the global CRE estimator (the upper bound) and global lower bound is analyzed by refining the initial mesh. The local problems [START_REF] Strouboulis | A posteriori estimation and adaptative control of the error in the quantity of interest. Part 1: a posteriori estimation of the error in the von Mises stress and the stress intensity factor[END_REF] defined on i are solved on a subdivision of the initial FE mesh P h,3 ( i ). Figure 6 presents the evolution of the effectivity indexes for the CRE estimator and for the global lower bound. In this example, the mean effectivity index of the CRE estimator is U = 1.08. The mean effectivity index of the lower bound is not as sharp: L = 0.73. Figure 7 presents the convergence rates of ê, e ref and ẽ. The three quantities have similar convergence rates, which are equal to 1. We do not obtain the theoretical convergence rate p = 2 of bilinear elements because of the presence of the corner singularities. Table I presents the numerical results for the second example, a function of the number of dof: the relative error estimator ε, the global upper bound ê, the global lower bound ẽ, the reference error e ref , the effectivity indices L and U . The evolutions of the effectivity indices are represented in Figure 8. The results obtained for this problem are similar to those obtained in the first example, with an upper bound sharper than the lower bound (mean effectivity indexe: U = 1.11 for the upper bound, L = 0.66 for the lower bound).

Application to goal-oriented error estimation in a crack opening problem

This section illustrates the behavior of the goal-oriented error estimator presented in Section 5 compared with the simpler version proposed in [START_REF] Gallimard | Error estimation of stress intensity factors for mixed-mode cracks[END_REF]. The first example shows that for a simple tension state we obtain an improvement of the lower bound, but not for the upper bound. The second example illustrates that for a more complex state of tension, the upper bound is also greatly improved.

For a crack opening problem, a classical quantity of interest is the mean vertical displacement along the crack boundaries The computed quantity is:

I = 1 a a Γ 1
I h = 1 a 1 u h •(-y) d + 2 u h .(y) d
The material properties are E = 100 GPa and = 0. The first example is a thin plate with a single edge horizontal crack subjected to a tension state represented in Figure 9. Four meshes are considered; the coarsest mesh is represented in Figure 10. On each mesh, the quantities I h , upp and low as well as CRE low and CRE upp have been computed and are reported in Table II.

A first observation is that the lower bound low for I -I h is greater than zero. This means that the lower bound for the quantity of interest I h + low is greater than the computed quantity of interest I h . The second observation is that the upper bound obtained directly with the CRE estimator CRE upp is close to the improved local upper bound upp defined in Section 5. However, if the lower bound is considered then the introduction of the improved lower bound low leads to really sharper bounds compared with CRE low (see Table II). This improvement is illustrated in Table III, where we have computed the size of the intervals (I and their ratio. In Figure 11 we have represented the evolution of the bounds for the quantity of interest as a function of the number of DoF.

The second example is a thin plate with a single edge oblique crack subjected to a tension state represented in Figure 12. Four meshes are considered, the coarsest mesh is represented in Figure 13. On each mesh the quantities I h , upp and low as well as CRE low and CRE upp have been computed and are reported in Table IV factor 2 as shown in Table IV. In Figure 14 we have represented the evolution of the bounds for the quantity of interest as a function of the number of DoF.

CONCLUSION

This paper has introduced a construction for the admissible stress fields used in the CRE estimator, which is more simple than the classical techniques proposed in previous papers. This technique leads to very good effectivity indexes. Based on these admissible stress fields, a lower bound for the energy error is defined. This lower bound can be combined with the CRE estimator (which is an upper bound for the energy error) so as to obtain lower and upper bounds for quantities of interest by a classical approach. This combined approach allows to obtain sharper bounds for a quantity of interest in comparison with alternative works based only on the constitutive relation approach.

Figure 2 .

 2 Figure 2. Model problem of the elastic plate with two holes.

Figure 3 .F=1Figure 4 .

 34 Figure 3. Plate with two holes: one-fourth of the plate is meshed.

Figure 5 .

 5 Figure 5. Elastic plate submitted to bending: coarse mesh.

Figure 7 .

 7 Figure

Figure 9 .

 9 Figure 9. Horizontal crack opening problem.

Figure 10 .

 10 Figure 10. Horizontal crack opening problem: coarsest mesh.

Figure 11 .

 11 Figure 11. Horizontal crack opening problem: upper and lower bounds for I .

Figure 12 .

 12 Figure 12. Oblique crack opening problem.

Figure 13 .

 13 Figure 13. Oblique crack opening problem: coarsest mesh.

Figure 14 .

 14 Figure 14. Oblique cracked plate problem: upper and lower bounds for I .

  Figure 6. Plate with two holes: upper and lower bound as a function of DoF.

		1.5
		Upper bound η U Lower bound η L
	Effectivity Index	1
		0.5
				10 3
			Number of DoF
		10 3
		Upper bound e
		Lower bound e
		Reference error e ref
	Global Errors
		10 4
		10 2	10 1	10 0	10 1
			h average
	Copyright q 2008 John Wiley & Sons, Ltd.		Int. J. Numer. Meth. Engng 2009; 78:460-482
					DOI: 10.1002/nme

Table I .

 I Elastic plate submitted to bending.

	Number of Dof	62	122	224	476	784	1572
	ε (%)	43.03	21.27	15.48	11.92	8.618	6.887
	ê	0.3708	0.1988	0.1490	0.1165	8.4690E-02	6.7958E-02
	ẽ	0.1271	0.1241	9.8705E-02	7.1684E-02	5.7886E-02	4.4531E-02
	e ref	0.2952	0.1881	0.1375	0.1010	7.9173E-02	6.2185E-02
	U	1.256	1.056	1.083	1.153	1.069	1.092
	L	0.430	0.659	0.717	0.709	0.731	0.716
			1.4	Upper bound η U Lower bound η L
			1.2
		Effectivity Index	0.8 1
			0.6
			0.4
				10 2	10 3
				Number of DoF
	Figure 8. Elastic plate submitted to bending: upper and lower bound as a function of DoF.

  3. The structure is studied in plane stress and meshed with 6-nodes triangular elements. The reference solution u ref is computed on a very refined mesh. A reference quantity of interest I ref is computed by replacing u by u ref in Equation (38). The local problems[START_REF] Strouboulis | A posteriori estimation and adaptative control of the error in the quantity of interest. Part 1: a posteriori estimation of the error in the von Mises stress and the stress intensity factor[END_REF] defined on i are solved on a subdivision of the initial FE mesh P h,3 ( i ).

Table II .

 II Horizontal crack opening problem: upper and lower bounds for I -I h .

	Number of Dof	154	316	590	1168
	I h	2.9816E-05	3.2939E-05	3.4195E-05	3.4771E-05
	low	2.8459E-06	1.3049E-06	6.8234E-07	3.8664E-07
	upp	6.1811E-06	3.1367E-06	1.7428E-06	1.1824E-06
	CRE low	-4.048E-06	-2.004E-06	-1.108E-06	-0.731E-06
	CRE upp	6.398E-06	3.236E-06	1.796E-06	1.215E-06

Table III .

 III Horizontal crack opening problem: upper and lower bounds for the quantity of interest.

	Number of Dof		154	316	590	1168
	I	upp h -I low h	3.34E-06	1.83E-06	1.06E-06	7.96E-07
	I	upp,CRE h	-I low,CRE h	6.80E-06	3.44E-06	1.91E-06	1.29E-06
	I	upp,CRE h I upp h -I low -I low,CRE h h	2.04	1.88	1.80	1.62
					4	x 10 -5
				bounds for the quantity of interest I	2.6 2.8 3 3.2 3.4 3.6 3.8		Upper bound Lower bound Direct Upper bound Direct Lower bound I ref
						10 2	10 3
							Number of DoF

  . Both the upper bound upp and the lower bound low are sharper than CRE upp and CRE low . The range obtained for the quantity of interest is divided by a Table IV. Oblique crack opening problem: upper and lower bounds for I -I h .

	Number of Dof		230	580	1152	2048
	I h	0.1461E-05	0.1493E-05	0.1541E-05	0.1564E-05
		low	-0.0708E-06	-0.1162E-06	-0.0391E-06	-0.0259E-06
		upp	0.5977E-06	0.3164E-06	0.1755E-06	0.0907E-06
		CRE low	-0.5220E-06	-0.4589E-06	-0.1922E-06	-0.0956E-06
		CRE upp	0.8387E-06	0.5510E-06	0.2779E-06	0.1385E-06
	I	upp,CRE h I upp h -I low -I low,CRE h h	2.04	2.33	2.19	2.01
				x 10 -6
			bounds for the quantity of interest I	1 1.5 2 2.5	Upper bound Lower bound Direct Upper bound Direct Lower bound I ref
				0.5
					10 3
					Number of DoF

Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2009; 78:460-482 DOI: 10.1002/nme

(31) Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2009; 78:460-482 DOI:

10.1002/nme

Int. J. Numer. Meth. Engng 2009; 78:460-482 DOI: 10.1002/nme

u•(-y) d +

u.(y) d (38)where a is the crack length and y is the normal unit vector along the crack edge. Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2009; 78:460-482 DOI: 10.1002/nme

Contract/grant sponsor: Publishing Arts Research Council; contract/grant number: 98-1846389

APPENDIX A

When i u rb does not belong to U h it is necessary to add forces F i to the loading on each domain i such that the resultant and the moment vanish on that domain and such that the sum of F i over the whole structure is zero. To simplify the presentation, we will limit ourselves to the case where T is equal to zero on * 2 . The forces are added on the vertices of the initial mesh P h . Each vertex i ∈ I v is associated with a domain i , the position of the vertex j is denoted by M j . In a domain i , the force added on the vertex j is defined by

Let us denote the set of vertices of the FE mesh belonging to i by I v ( i ). Equation ( 22) is modified and the admissible stress ri h in i is sought as a solution of the following modified problem:

Proposition A.1

The nodal forces g i j do not modify the global equilibrium on the structure (i.e. i∈I v g i j = 0).

Proof

Let us consider the sum of the nodal forces over the set of the vertices I v . As i = 0 ini we have the following relations:

The set ( i ) i∈I v is a partition of unity: i∈I v i = 1 and grad i∈I v i = 0. The global force applied to the vertex j is:

Proof

On a domain i the rigid body motion is written on the basis of the linear shape functions

where q j rb is the displacement of the vertice j


Jean Slawinski
email: jeanslawinski@yahoo.fr

Veronique Billat

Slawinski Js

Billat Vl

Facsm

Slawinski Jean

Changes in

Keywords: Energy cost, oxygen slow component, mechanical work and fatigue

published or not. The documents may come

Introduction

Paragraph Number 1

The process of determining the energy cost of running (C r ) during an exercise of moderate intensity is based on two assumptions. The first is that C r is independent of running time for constant-speed running of 3-30 minutes duration. The second is that a steady-state of oxygen uptake ( ) occurs following 3-4 min of exercise ( [START_REF] Morgan | Factors affecting running economy[END_REF].

2 O V &
However, during constant-load exercise of supra-lactic threshold intensity, the steadystate of is delayed after 3-4 min. A number of physiological factors have been postulated as contributing to this delayed steady state of during a supra-lactic threshold exercise ( [START_REF] Billat | The V(O2) slow component for severe exercise depends on type of exercise and is not correlated with time to fatigue[END_REF]13,[START_REF] Gaesser | The slow component of oxygen uptake kinetics in humans[END_REF]. However, biomechanical factors have received less attention as potential causes which could partially explain this increase in C r .

2 O V & 2 O V &

Paragraph Number 2

In running studies, biomechanical factors have often been used to address the differences in C r among individuals (17, [START_REF] Williams | Relationship between distance running mechanics, running economy, and performance[END_REF]. Kram (22), as well as Farley and McMahon [START_REF] Farley | Energetics of walking and running : insights from simulated reduced-gravity experiments[END_REF], had shown that C r could be determined by the force generated by the muscle during running. Others have also shown that an inverse relationship exists between the stiffness of the body and the C r (15). However, the interactions between mechanical and metabolic variables appear to be very complex. A general mechanical approach has been based on the observation of the mechanical energy change of the centre of mass and of the limbs of the body (14, 34). Indeed, the mechanical energy change associated with the movement of the different segments of the body during running, represents the energy used by the active muscle. Poole et al. [START_REF] Poole | Contribution of excising legs to the slow component of oxygen uptake kinetics in humans[END_REF] have demonstrated that 86% of this increase in C r was related to the increase in active muscle . Therefore the increase in C 2 O V & r observed at the end of supra lactic threshold exercise may be related to the variation in mechanical energy.

Paragraph Number 3

Some studies have significantly related the differences in C r to the energy variations of the centre of mass during running (11,12). Above the lactic threshold and under the influence of fatigue, the mechanical cost associated with the movements of an anatomical point taken as equivalent of the body centre of mass of the runner, has

Paragraph Number 4

Therefore, in order to complete the understanding of the relationship between mechanical and metabolic variables, this study aims to determine whether a change in C int could occur during an exhausting over-ground run and whether this change is related to a modification of C r . The hypothesis was that the increase in C r , generally observed after the 3 rd min of a supra-lactic threshold exercise, is associated with a modification of C int .

Methods

Subjects and protocol

Fourteen subjects (mean height 1.73 ± 0.06 m, mean body mass 62 ± 9 kg, and mean age 21.9 ± 2.8 yrs) volunteered to participate in this study. This population was composed of 3 females and 11 males. Five subjects were physical education students (soccer, rugby and tennis players), the others subjects trained regularly in running (they had a national or a regional level of competition). Table 1 summarises this information. Prior to participation, all the subjects were informed of risks and stress associated with the experimental protocol and gave a written voluntary informed consent and approval received by ethics committee in accordance with the guidelines of the hospital of Paris S t Louis.

Specialist of running (n = 9)

Non specialist of running (n = 5) -An incremental test (3-min stages) to exhaustion (the voluntary stop of the exercise). This test was performed on a track in order to determine the maximal oxygen consumption ( ). This parameter was defined as the highest 30 s

Total (n = 14) v (m.s 2max O V & -1 )
max 2 O V & O V &
value attained during the test. The velocity associated with

(v ) ( O V & max 2 O V & 9)
. The velocity at the lactate threshold (vLT) was defined as the speed measured at value which corresponds to the starting point of an accelerated lactate accumulation between 3.5 and 5 mmol.

L 2 O V & -1 (1).
-An exhaustive run at 95% of v until exhaustion (tlim95).

Throughout the tests, the subjects adopted the required velocity using an audio rhythm the last minute of the exercise). The running speed was documented with the video camera, for each time the runner crossed in front of the camera (Table 2). The energetic cost of running (C r ) was also evaluated throughout the exercise, by the oxygen uptake measurement.

rd min (n = 14) End (n = 14)

Velocity (m.s -1 )

4.87 ± 0.40 4.89 ± 0.37 SR (Hz)

2.98 ± 0.12 3.00 ± 0.16 SL (m)

1.62 ± 0.14 1.63 ± 0.11 C ke (j.kg -1 .m -1 )

1.33 ± 0.33 1.38 ± 0.29 C pe (j.kg -1 .m -1 ) 0.47 ± 0.11 0.48 ± 0.10 C int (j.kg -1 .m -1 ) 0.55 ± 0.07 0.60 ± 0.07 † C r (j.kg -1 .m -1 ) 4.10 ± 0.39 4.32 ± 0.42 † † Significantly different from the 3 rd minute of the exercise (p ≤ 0.05).

Table 2: Changes in the stride rate (SR), the stride length (m), the velocity of the exercise, kinetic (C ke ), potential (C pe ), internal mechanical cost (C pe ) and energetic cost of running (C r ) between the 3 rd and the last minute of the exercise.

Paragraph Number 6 Material

Oxygen uptake measurement. Throughout the exercises, the respiratory and pulmonary gas exchange variables were measured using a breath-by-breath portable gas analyzer (Cosmed K4b 2 , Roma, Italy). Before each test, O 2 and CO 2 analysers were calibrated using ambient air and sample gas references of 16% O 2 and 5% CO 2 .

The flowmeter was calibrated with a volume of air of 3-l (Quinton instruments, Seattle, USA). The accuracy of this system has been tested and is acceptable for and measurement during supra-lactic threshold exercise ( [START_REF] Mcfarlane | Automated metabolic gas analysis system[END_REF].

2 O V & 2 CO V &

Paragraph Number 7

Mechanical measurement. The video sequences obtained during the tlim95 were digitised without distortion on a PC as a series of bitmap images by a Perception Video Recorder card from Silicon Co. Then the images recorded were displayed and analysed on a screen of a regular PC. ANIMAN 3D uses the numerical human model MAN3D [START_REF] Slawinski | Use of lumbar point for the estimation of potential and kinetic mechanical power in running[END_REF][START_REF] Slawinski | Difference in mechanical and energy cost between highly, well and non trained runners[END_REF][START_REF] Tavernier | A 3D motion analysis system using a numerical human model[END_REF] which is delimited by an envelope composed of 155 crowns.

These crowns are superimposed on the different segments of the manikin (figure 1).

The position and posture of MAN3D can be adjusted to each runner. The morphological properties of MAN3D are deduced from the size and the weight of the subject. The inertial properties of the segments are from works of Dempster and Gaughran (16). Then MAN3D is projected onto images. When the projections of MAN3D were superimposed onto the images of the runner, the position and posture of the runner were the same as that of MAN3D. Posture and position adjustments were optimized by locating the positions of some particular points of the runner on the images which must coincide with some equivalent points of MAN3D. There are 5

particulars points which help the researcher to project Man3D on the runner. These points are represented in the following picture (figure 2):

-The top of the head.

-The hand in direct view.

-Two points on the back.

-The end of the foot in direct view.

Our method has been compared with the measure of a force plate-form ( ). The results have shown that the difference between both systems on the computation of the mechanical energy was about 11% and remains constant for different running speeds.

When the segments (hand, forearm and arm) are out of direct view, the system allows Man3D to be projected onto a partially hidden segment with accuracy. Indeed, when the elbow is in direct view, it is possible to deduce the position of the forearm and of the arm with accuracy. However, when the elbow and the forearm is out of direct view (as in figures 1 or 2) the system does not allow one to deduce with the same accuracy the position of the hand, the forearm or the arm. But this case (when the arm is totally hidden) does not represent more than one or two images during one step. Moreover, the position of the hidden arm can be deduced from the previous and the following positions of the arm.

Paragraph Number 8 Data analysis

Energy cost (C r ). The breath-by-breath oxygen uptake data were reduced to 5-s stationary averages in order to reduce the noise so as to enhance the oxygen kinetics characteristics. These data were finally fitted to two exponential functions using a least square fit method: a single-exponential function comprising a delayed linear component (eq. 1) and a double-exponential function comprising two exponential terms which start at two distinct time delays from the onset of exercise (eq. 2) (4). The Fisher test was used to choose the model for which the fit was associated with the highest F value (4).

O V & (t) = y 0 + A 1 ×(1-e -(t-TD1)/ τ 1 ) ×U 1 (fast component) [1] + A 2 ×(1-e -(t-TD2)/ τ 2 ) ×U 2 (slow component) [2]
Where y 0 is the baseline (ml.min

2 O V & -1
), A 1 and A 2 are the asymptotic amplitudes for the exponential terms (ml.min -1 ), τ 1 and τ 2 are the time constants (min), and TD1 and TD2 are the time delay from the onset of exercise (s).

U 1 = 0 for t < TD1 and U 1 = 1 for t ≥ TD1 U 2 = 0 for t < TD2 and U 2 = 1 for t ≥ TD2

The energy cost of running (C r ) was evaluated, 3 mins after the beginning of the exercise and 1 min before its end (12):

1 1 95 O 2 r M V E ) t ( O V C 2 - -⋅ ⋅ ⋅ = & Where (21.3 J.ml ) is the energy equivalent of 1 mlO 2 O E 1 2
O - 2 for a respiratory exchange ratio of 1 (12), and V 95 is the velocity during the test tlim95 (measured with the video camera).

Paragraph Number 9

Mechanical cost. All the trajectories of the different centres of mass of the model (segments and body centres of mass) were smoothed using a polynomial method in order to obtain by derivation the speed displacement of the different centres of mass.

The speed of running was measured with the video camera after the trial. In the referential of the camera, the horizontal position of the CM increases linearly according to time (figure 3). The slope of the relationship between time and the horizontal position of the CM allows one to calculate the speed of the CM of the runner. Moreover, the speed was measured with photo cells placed in the field of the camera. The polynomial method used was a polynomial of high multiple order. After plotting the residual according to the order of the polynomial (37), the polynomial degree was determined individually for the trajectory of the centre of mass of each subject. Variations of potential and kinetic energies (ΔE pe and ΔE ke respectively) were then calculated (in joules).

ΔE pe = ( ) min max H H g M - × × [3] Δ ( ) 2 E ke = min max V V M 2 [4]
"M" is the body mass (kg), "g" is the gravitational acceleration (9.81 m.s ( )

∑ = + = 1 i s s s s int ω I V m 2 E [5]
Where "s" represents a segment, "m 17 2 2 1 s " is the mass of the segment (kg) considered and "I s " is the moment of inertia of the segment (kg.m 2 ). " s V " is the velocity of the centre of mass of the segment (m.s -1 ) with respect to the referential linked to the CM. " s ω "

is the angular velocity of the centre of mass of the segment (rad.s -1 ) with respect to the referential linked to the CM of the segment.

are equal to the positive variation of E ke , E pe and E int during one step, divided by the the number of images for a stride and the step length (SL) was computed by dividing the running speed by SR.

The percentages of variation between the 3 rd min and the end of the exercise of C r and C int , were also calculated (%ΔC r and %ΔC int ).

mass and the step amplitude of the subject.

The step rate (SR) was computed by counting the energy cost of running (C ) increased significantly within the 3 rd and the last minute (table 2).

Paragraph Number 13

ithin the 3 rd and the last minute of exercise the velocity of running, the stride rate (SR) and the stri e l ngth (SL) did not vary significantly (table 2). At the same time, C ke and C pe also remained constant (table 2). Indeed, C ke varied from 1.33 ± 0.33 1 .m -1 to 1.38 ± 0.29 J.k 0.48 ± 0.10 J.kg -1 .m -1 (P = 0.67). These results demonstrated that during track running there was no modification of SR, C or C at the end of an exhaustive exercise.

Paragraph Number 14

However int However, the percentage of variation of C int (%ΔC int ) was significantly different between the specialist and the non inally, there was no correlation between the C r and C int (figure 4C) or between %ΔCr and %Δcint (r = 0.06; P = 0.84). The increase in C r observed after the 3 rd min of a supra-lactic threshold exercise is not associated with a modification of C int .

specialist of running (p = 0.03). Indeed, the non specialists increase C int by 20%, within the 3 rd and the last minute of exercise, while the specialists do not increase it (4%). This increase in C r between the 3 rd and the last minute of exercise is not accompanied by an increase of C ke and C pe . Therefore the increase in C r cannot be related to C ke or C pe . Following the idea that smaller vertical oscillation of the centre of mass is associated with a low C r , there was a significant correlation between C r and C pe (figure 4B). As already demonstrated in the literature (11, 12), low but significant correlatio have been established

was a significant correlation between C r and C ke or C pe , the increase in C r could not be related to the mechanical cost associated with the movements of the centre of mass during the step.

Paragraph Number 19

A detailed analysis of mechanical data showed that the mechanical cost associated with the movements of the limbs around the centre of mass (C int ) increases significantly by 9.7% within the 3 rd and the last minute of the exercise. This increase in C int during track running has never been described and depends on mechanical adaptations which have yet to be explained. Numerous factors influence C int , speed [START_REF] Willems | External, internal and total work in human locomotion[END_REF], load, gradient, gait, or stride frequency [START_REF] Minetti | A model equation for the prediction of mechanical internal work of terrestrial locomotion[END_REF]. In order to simplify the calculation of C int (in Watts), this has been modelled by Minetti and Saibene (24) as:

( )( ) [ ]

Paragraph Number 20

Thus, the internal mechanical cost depends on the speed and on the mass of the segments, as well as on the segmental organisation during running. Using this equation, B significantly from 0.90 component period. This decrease is mainly the result of the respective decrease in contact time (-2.6%) and stride frequency (-2%). Contrary to the results of Borrani et al. (10), the present study showed that C int increased within the third minute and the last minute of the exercise. This difference can be explained by the method used.

Indeed, the study of Borrani et al. (10) was carried out on a treadmill whereas the differences observed between both studies may be explained by specific mechanical int could be linked to an increase in the leg speed during the stride or to a change in the consistent with an increase in the leg speed. As C ke or C pe , the stride rate (SR) was not modified during an exhaustive run. Therefore, during over-ground running, the lack of any change did not support a role for SR during fatigue as a mechanism for the increase in C r , as previously proposed by Borrani et al. (10) and Avogadro et al. [START_REF] Avogadro | Changes in mechanical work during severe exhausting running[END_REF]. Furthermore, two studies carried out in order to understand the possible influence of cadence on the increase of C r during strenuous exercise (3, 7) have already demonstrated that in cycling as in running, the increase of C r was not influenced by cadence.

Paragraph Number 21

int after the 3 rd min of exercise could be related to a change in fications. Indeed, the use of a global method of calculation of the not allow one to determine precisely the segments responsible Thus, the increase in C the geometric configuration of the segments of the body. One probable hypothesis would be connected to the effect of fatigue on the biomechanical stride characteristic [START_REF] Gazeau | Biomechanical events in the time to exhaustion at maximum aerobic speed[END_REF][START_REF] Verkerke | Measuring changes in step parameters during an exhausting running exercise[END_REF]. These effects of fatigue on C int , could depend on the level of experience in running. Indeed, the Student test has showed that the great variability in the %ΔC int depends on the skill/experience of the runner. On track, the non specialists increase C int by 20%, within the 3 rd and the last minute of exercise, while the specialists increase it by only 4%. However, it is difficult to identify segments responsible for these modi mechanical work does for the observed modifications. Only the use of a local method of the calculation of the mechanical work such as the method of calculation of powers developed in joints would allow a better understanding of the observed mechanical modifications. One of the perspectives of the present work would be to study the effects of fatigue in the mechanical work calculated by the method of the joints moments. Hence, the observed modifications of the style could be explained. However, the use of the joints method in running is limited to the study of the power developed at the level of the joints of lower limbs [START_REF] Belli | Moment and power of lower limb joints in running[END_REF][START_REF] Winter | Moments of force and mechanical power in jogging[END_REF].

Paragraph Number 22

The present results showed that both C int and the energy cost of running (C r ) increase within the 3 rd minute and the end of the supra-threshold exercise. However, contrarily of the present study, the figure 4C showed that there was no y of individual variations of C r and C int shows that some runners who displayed an increase in C int did not display an increase in C r . This result showed that the increase of C r could not be related to the increase in C int and that two independent mechanisms might explain the increase in C r and in C int .

Paragraph Number 23

The results reported in this study showed that during track running, exhaustion has no effect on C pe or C ke . Only the C r and C int increased significantly at the end of a supralactic threshold exhaustive run. However, the C r is only correlated with C ke and C pe , and not with C int . Thus, the increase in C r is not associated with a modification of the mechanical work produced in running. The increase in C int after the 3 rd minute of exercise could be related to a change in the geometric configuration of the segments of the body. The geometric configuration might depend on the skill/experience of the runner. The modification of C int suggests that the technique of running is a basic component of the training of endurance runners and above all in fatigue status.

However, the increase in C int during over-ground running to exhaustion has yet to be

  been significantly correlated with C r and increased significantly between the beginning and the end of the exercise (12). Even if the interactions between biomechanical and physiological factors are still open to debate, this result suggests that the increase in C r , generally observed after the 3 rd min of a supra-lactic threshold exercise, may be due to a combined action of both physiological and mechanical mechanisms. Recently, Borrani et al. (10) and Avogadro et al. (2) suggested that the increase in C r in running did not result from a change in the external mechanical cost under the effect of fatigue. Nevertheless, the decrease in stride frequency suggested an alternative mechanical explanation such as a modification of the mechanical cost associated with the movements of the limbs around the centre of mass which is called the internal mechanical cost (C int ).

  which gave the time to cover 20-m. Visual marks were set at 20-m intervals along the track with audio signals determining the speed needed to cover 20-m intervals. To estimate the kinetic, potential and internal mechanical cost of running (C 2max O V & ke , C pe and C int ), the mechanical analysis was performed during the tlim95. The runners were filmed by a video camera using a sampling frequency of 25 frames per second (Sony, Beta SP VCRs, resolution of 572 lines). The camera was positioned 9 meters from the left of the runner. Two consecutives steps were analysed at the beginning of the test (during the second lap of the track, around the 3 rd min after the beginning of the exercise) and two more steps at the end of the test (during the last lap of track, around

29

  This procedure allows the trajectory of the body centre of mass and the trajectories of the limb centres of mass to be obtained. The uncertainty in the determination of the different centres of mass trajectories depends on the resolution of the video camera used.

Figure 1 :

 1 Figure 1: MAN3D human model representation and its 17 segments.

Figure 2 :

 2 Figure 2: Example of projection of Man 3D on the runner, with the 5 particulars points (black).

2

 2

  the maximal and minimal horizontal velocities of the .s -1 ). al energy (Δ ) is calculated as the maximal minus the minimal E int duri -2 ), "H max and H min " are maximal and minimal heights of the body centre of mass (CM) during one step (m). "V max and V min " a CM during one step (m Paragraph Number 10 The internal energy (E int ) represents the sum of the mechanical energies associated with the different segments of the body. The present model used in this work (figure 1) is composed of 17 different segments. The variation of intern int E ng the step.

Figure 3 :

 3 Figure 3: Example of horizontal position of the CM in the referential of the camera, and calculus of the speed of the CM of the runner.

  , significant correlations have been established between C 0.43 and P = 0.03; r = 16% variance between measures. Therefore, such as it was calculated in the present work, C or C constituted a global descriptor of the C . ke rParagraph Number 15The internal mechanical cost (C ), contrary to C and C increased significantly between the 3 and the last minute of the exercise (P ≤ 0.05, table2). The average percentage of variation of C (%ΔC ) was 9.7% with a great inter-individual variability. Indeed, eleven runners presented an increase of C , (comprised) of between 4 and 35% and 3 runners presented a decrease of C (comprised) of between -4 and -15%. Moreover, no correlation was observed between C and C , C SR or SL (respectively, r = 0.07 and P = 0.73; r = -0.2 and P = 0.31; r = 0.06 and P = 0.77; r other mechanical parameters.

F the 3 Figure 4 :

 34 Figure 4: Correlation, between C r and C ke (A), C pe (B) and C int (C) measured in the fourteen runners, at

  et al.(10) have recently showed that on the treadmill, C int decreases (0.07) to 0.87 (0.04) J.kg -1 .m -1 (3.5%) throughout the slowpresent results have been obtained during over-ground running. Therefore the adaptations induced by treadmill or over-ground running(18,[START_REF] Wank | Kinematics and electromyography of lower limb muscles in overground and treadmill running[END_REF]. This increase in C segmental organisation. The stability of SR, observed in the present study, was not int mainly depends on the SR, and on the running speed (V). However other parameters influence C int , including the mass of the lower and upper limb (m L and m U ), the fractional distance of the lower limb centre of mass from the proximal joint (a), the length of the upper limb as a fraction of the lower limb one (b) and the average radius of gyration of limbs, as a fraction of the limb length (g).

  explained. time than intense but submaximal runs. Eur. J. Appl. Physiol.81: 188-196, 2000. 10. Borrani, F., R. Candau, S. Perrey, G.Y. Millet, G.P. Millet, and J.D. Rouillon. Does the mechanical work in running change during the slow component? Med. Sci. Sports Exerc. 35: 50-57, 2003. Bourdin, M., A. Belli, L. Arsac, C. Bosco, J.R. Lacour. Effect of vertical loading on energy cost and kinematics of running in trained male subjects. J. Appl. Physiol. 79: 2078-85, 1995. 12. Candau, R., A. Belli, G, Millet, D, Georges, B, Barbier, J.D. Rouillon. Energy cost and running mechanics during a treadmill run to voluntary exhaustion in humans. Eur. J. Appl. Physiol. Occup. Physiol. 77: 479-85, 1998. 13. Carter, H., A. Jones, T. Barstow, M. Burnley, C. Williams, and H. Doust. Effect of endurance training on oxygen uptake kinetics during treadmill running. .J Appl. Physiol. 89: 1744-52, 2000. 14. Cavagna, G.A., and M. Kaneko. Mechanical work and efficiency in level walking and running. J. Physiol. 268: 647-81, 1977. 15. Dalleau G., Belli A., Bourdin M., Lacour J. R. The spring-mass model and the energy cost of treadmill running. Eur. J. Appl. Physiol. Occup. Physiol. 77: 257-63, 1998. 16. Dempster, W., and G., Gaughran. Properties of body segments based on size and weight. Am. J. Anat. 120: 33-54, 1967. 17. Donelan, J., D. Shipman, R. Kram, A. Kuo. Mechanical and metabolic requirements for active lateral stabilization in human walking. J. Biomech. 37: 827-835, 2004.

Acknowledgments

The authors gratefully acknowledge the Fédération Française de Ski and Dr. Michel Tavernier for their assistance during the experiment.

to the hypothesis correlation between C r and C int . Indeed, the stud