Laurent Pierre
email: laurent.pierre@bordeaux.ensam.fr

Denis Teissandier
email: denis.teissandier@u-bordeaux1.fr

Jean Pierre Nadeau
email: jean-pierre.nadeau@bordeaux.ensam.fr

Tolerancing analysis taking into account thermomechanical strains

Keywords:

To improve the performance of a helicopter turboshaft engine requires optimising the energy yield of the different components, and more particularly controlling clearance between the tips of the high pressure turbine blades and the stator. Dimension-chain tools take into account the manufacturing dispersion of the parts and assembly defects. This ensures the interchangeability of the different components and guarantees that a turbine can carry out different service functions, as the turbine is modelled in infinitely rigid solids. However, this approach does not take thermomechanical effects into account. And yet, the different operating regimes of a helicopter engine make it indispensable that the effects caused by the thermodynamic cycle should be integrated. The aim of this article is to show how using dimension chain and thermomechanical tools can contribute to controlling clearances at the tip of a high pressure turbine blade.

2-Tolerancing analysis taking into account thermomechanical strains

The geometric models used in 3D dimension-chains are generally based on the following hypotheses: no form default in real surfaces, no local strain in surfaces in contact and no deformable parts.

Two different behaviours will be considered for every system studied: -rigid behaviour: all parts are at 20°C and are considered as being infinitely rigid.

-thermomechanical behaviour: some parts are subjected to thermomechanical stresses. Thermomechanical strains on the parts are then taken into account in the 3D dimension chains.

The problem was restricted to those systems where the topological structure of the contact graphs (or of the joints) in rigid behaviour and thermomechanical behaviour remain the same. This means that there is no supplementary joint and no suppression of a joint between the two behaviours. On the other hand, a joint defined by different parameters (minimal clearance, maximal clearance, nature of contact …) may change.

2.1-Rigid behaviour: dimension-chain of infinitely rigid solids

Modelling the geometrical defects in a part

Real surfaces are modelled by substituted surfaces [START_REF] Bourdet | Geometrical behavior laws for computer aided tolerancing[END_REF]. A substituted surface is an ideal surface (i.e. geometrically perfect) of the same type as the nominal surface of which it is a particular physical representation. This is described in figure 1, which illustrates the three different types of modelling used in tolerancing: the nominal model (CAD model), the skin model and the substituted surfaces model. The nominal model is the one used in a geometric modeller: it is by definition geometrically perfect. The skin model provides a representation of the real surfaces of the part and acts as a graphic support for expressing geometric specifications [START_REF] Dantan | Geometrical product specifications -model for product life cycle[END_REF].

Nominal model (CAD model)

The skin model of the shaft consists of three real surfaces: one is nominally plane and two are nominally cylindrical. Lastly, the substituted surface model consists of three surfaces: one plane and two cylindrical. The geometric defects of a real surface are defined by the relative position between a substituted surface and the corresponding nominal surface. Figure 2 illustrates the defaults of the relative position of the substituted surfaces designated by 1,1; 1,2 and 1,3 of the shaft 1. Each surface is designated by a pair of numbers: 1,2 refers to (substituted) surface 2 of part 1. Surfaces 1,1n; and 1,2n 1,3n refer respectively to the nominal surfaces corresponding to surfaces 1,1; 1,2 and 1,3. Figure 2 shows part 1 (large circle) in diagram form and also shows surfaces 1,1; 1,2 and 1,3 (small circles). This acts as a graphic support to show the structure of a mechanism on which one can visualise the dimension-chains [START_REF] Ballu | Choice of functional specifications using graphs within the frame work of education[END_REF].

Observe 1,1/1,2 d   position deviations between surfaces 1,1 and 1,2. Position deviations between surfaces 1,2 and 2,2, which are still called simply geometric deviations, characterise the geometric defaults of surfaces 1,1 and 1,2. When modelling these deviations using substituted surfaces it is not possible to take into account the form defaults of surfaces 1,1 and 1,2.

1,1/1,2 d   can be formalised mathematically by a small displacement torsor [BB1], by a matrix [T1], etc. In this article, the small displacement torsor is used [START_REF] Clément | A Study of Optimal-Criteria Identification Based on the Small-Displacement Screw Model[END_REF].

The following relations can be written:

1,1/1,2 1,1/1,1n 1,1n /1,2n 1,2n /1,2 d d d d                    (1)
Given that by definition the nominal model of a part has no defect, (2) can be expresses as follows :

1,1/1,2 1,2 /1,2n 1,2n /1,2 d d d              (2)

Modelling 3D dimension-chains in a mechanism

Let us consider figure 3 : two parts, 1 and 2, are in cylindrical pair type contact via their respective surfaces 1,2 and 2,2 and in ball and plane pair contact via their respective interfaces 1,3 and 2,3. In a diagram, this is represented by two edges linking respectively apexes 1,2 ; 2,2 and 1,3 ; 2,3. Let us suppose that a functional condition (FC) limiting the relative position of surfaces 1,1 and 2,1 must be respected. This is a condition of coaxiality such that, for every direction orthogonal to the x-axis, the displacement d of point A on the axis of surface 1,1 in relation to the axis of surface 2,1 respects the following equation (see figure 3):

min max dd d  (3)
In the diagram, FC is within a rectangle on an edge linking apexes 1,1 and 2,1.

The diagram describing the parts and surfaces of a system and also the joints and the functional condition is commonly called a joint graph or contact graph [START_REF] Ballu | Choice of functional specifications using graphs within the frame work of education[END_REF].

Here the specifications for coaxiality of 1,1 in relation to 1,2 and of 2,1 in relation to 2,2 on parts 1 and 2 respectively are specifications that impact on the FC. Similarly, diameter specifications for surfaces 1,2 and 2,2 also influence the FC.

x y z 1 1,2 1,3 1,1 2 2,2 2,3 2,1 A x y z x y z 1 1,2 1,3 1,1 2 2,2 2,3 2,1 A 1 1,2 1,3 1,1 2 2,
d 1,2 t 1,1 A L 1,1 t 2,2 B D 2 2 d 2,2 B L 2, 2 A D 1,2 d 1,2 t 1,1 A L 1,1 A D 1,2 d 1,2 t 1,1 A t 1,1 A L 1,1 t 2,2 B D 2 2 d 2,2 B L 2, 2 t 2,2 B D 2 2 d 2,2 B t 2,2 B t 2,2 B D 2 2 d 2,2 B L 2, 2 transfer transfer D 2,2 D 1,2 FC 1,2 1,1 1 2,1 2,2 2 1,3 2,3 CP BPP D 2,2 D 1,2 FC 1,2 1,1 1 2,1 2,2 2 1,3 2,3 CP BPP FC 1,2 1,1 1 2,1 2,2 2 1,3 2,3 BPP CP FC 1,2 1,1 1 2,1 2,2 2 1,3 2,3 BPP CP x y z 1 1,2 1,3 1,1 2 2,2 2,3 2,1 A x y z x y z 1 1,2 1,3 1,1 2 2,2 2,3 2,1 A 1 1,2 1,3 1,1 2 2,
d 1,2 t 1,1 A L 1,1 t 2,2 B D 2 2 d 2,2 B L 2, 2 A D 1,2 d 1,2 t 1,1 A L 1,1 A D 1,2 d 1,2 t 1,1 A t 1,1 A L 1,1 t 2,2 B D 2 2 d 2,2 B L 2, 2 t 2,2 B D 2 2 d 2,2 B t 2,2 B t 2,2 B D 2 2 d 2,2 B L 2, 2 transfer transfer D 2,2 D 1,2 FC 1,2 1,1 1 2,1 2,2 2 1,3 2,3 CP BPP D 2,2 D 1,2 FC 1,2 1,1 1 2,1 2,2 2 1,3 2,3 CP BPP FC 1,2 1,1 1 2,1 2,2 2 1,3 2,3 BPP CP FC 1,2 1,1 1 2,1 2,2 2 1,3 2,3 BPP CP Figure 3 : Transfer of FC.
Generally speaking, the two location specifications and the two diameter specifications are the result of the transfer of FC to parts 1 and 2.

Each one appears on an arc oriented to the circles of parts 1 and 2. The arrow indicates the toleranced element and the origin of the arc indicates the surface of reference. For the diameter specifications, these are specifications that are intrinsic to surfaces 1,2 and 2,2. The influent specifications, the influent joint (CP) and the FC formalise an influent cycle (see figure 3). Equation (4), justified by the influent cycle, formalises at point A dependences between the FC, the location deviations of parts 1 and 2 and the cylindrical pair joint deviations (between surfaces 1,2 and 2,2) [START_REF] Teissandier | A Computer Aided Tolerancing Model : Proportioned Assemblies Clearance Volume[END_REF] :

1,1/2,1 1,1/1,2 1,2/2,2 2,2/2,1 d = d + d + d                 (4)
One result here is that the ball and plane pair joint between 1,3 and 2,3 is non-influent on the FC. The deviations of surfaces 1,3 and 2,3 are also non-influent on FC. This means that respect for the FC is independent of any geometrical defaults in surfaces 1,3 and 2,3 and of the ball and plane pair joint between 1,3 and 2,3.

Formalisation of relations between functional condition, geometric specifications and contact specifications

In this article, deviation hulls [START_REF] Giordano | Clearance Space and Deviation Space, Application to three-dimensional chains of dimensions and positions[END_REF] to define the admissible limits for geometrical defaults in a part (defined by a specification) and clearance hulls to define the admissible limits for relative displacements between two surfaces in contact (defined by a clearance) [GS1] will be used.   the Small Displacement Torsor, characterising the relative position of 1,1 in relation to 1,2. In the base   ,,

A D 1,2 d 1,2 t 1,1 A L 1,1 Tolerance zone A xis A = A xis of 1,2 t 1,1 A xis of 1,1 A B x y z L 1,1 x y z A D 1,2 d 1,2 t 1,1 A L 1,1 A D 1,2 d 1,2 t 1,1 A t 1,1 A L 1,1 Tolerance zone A xis A = A xis of 1,2
x y z at point A 1,1/1,2 d   is expressed as follows:

1,1,/1,2 1,1/1,2 A A,1,1,/1,2 d       ρ ε (5) with 1,1/1,2 1,1,/ 2,1 1,1/1,2 1,
        2 2 2 1,1 A,1,1/1,2 A,1,1/1,2 1,1/1,2 2 2 2 1,1 B,1,1/1,2 B,1,1/1,2 .. 2 .. 2 t t                          ε y ε z ε y ε z D (6)
To simplify the writing of the hulls, the limits of admissible defects in 2D modelling will be studied in the plane   A, z .

1,1/1,2 D can be written::

    2 2 1,1 1,1 1,1 A,1,1/1,2 A,1,1/1,2 1,1/1,2 2 1,1 1,1 2 1,1 B,1,1/1,2 B,1,1/1,2 . . 2 2 2 . . 2 2 2 t t t t t t                                                 ε y ε y ε y ε y D (7)
According to the property of the small displacement fields [START_REF] Clément | A Study of Optimal-Criteria Identification Based on the Small-Displacement Screw Model[END_REF] (see figure 4) it is possible to deduce:

B,1,1/1,2 A,1,1/1,2 1,1/1,2    ε ε BA ρ with 1,1 0 0 L       BA (8)
Hence with (8) 1,1/1,2 D is written at point A:

1,1 1,1 A,1,1/1,2 1,1/1,2 1,1 1,1 A,1,1/1,2 1,1 1,1/1,2 22 . 2 2 y y z tt t t L                         D   (9)
Figure 5 shows a graph of the deviation hull 1,1/1,2 D defined by (9). The translation deviation at point A following y is projected on axis Ay  and the rotation deviation following z is projected according to z  . The interior of the hull (boundary included) represents the geometric deviation in conformity with the coaxiality specification (see figure 4) [START_REF] Teissandier | Operations on polytopes: application to tolerance analysis In proc[END_REF].

1,1 1,1 t L  1,1 1,1 t L  1,1 2 t  1,1 2 t  Ay  z  1,1 1,1 t L  1,1 1,1 t L  1,1 2 t  1,1 2 t  Ay  z  Figure 5 : Deviation hull 1,1/1,2 D at point A (coaxiality of part 1).
By the same reasoning, the deviation hull 2,2/ 2,1 D is determined at point A, which characterises the coaxiality specification for part 2, which is shown following axes Ay  and z  in figure 6.

2,1 2,1 t L  2,1 2,1 t L  2,1 2,1 1,1 2,1 2 L t L L          2,1 2,1 1,1 2,1 2 L t L L          Ay  z  2,1 2,1 t L  2,1 2,1 t L  2,1 2,1 1,1 2,1 2 L t L L          2,1 2,1 1,1 2,1 2 L t L L          Ay  z  Figure 6 : Deviation hull of coaxiality of part 2 at point A.
Moreover, using a similar method to that used for the deviation hulls, from figure 7, illustrating the clearance of the cylindrical pair joint between surfaces 1,2 and 2,2, the clearance hull 1,2/ 2,2 D is defined. The clearance hull 1,2/ 2,2 D expressed at point A is represented in figure 8, and given in equation (10).

    max 2,1 2,1 1,2 1,2 J D d D d     (10) x y z 1 1,2 1,3 1,1 2 2,2 2,3 2,1 A x y z x y z 1 1,2 1,3 1,1 2 2,2 2,3 2,1 A 1 1,2 1,3 1,1 2 2,2 2,3 2,1 A x y z A 1,2 2,2 J max J max L 2,2 L 1,1 x y z A 1,2 2,2 J max J max L 2,2 L 1,1
Nominal model Model by substituted surfaces

x y z 1 1,2 1,3 1,1 2 2,2 2,3 2,1 A x y z x y z 1 1,2 1,3 1,1 2 2,2 2,3 2,1 A 1 1,2 1,3 1,1 2 2,2 2,3 2,1 A x y z A 1,2 2,2 J max J max L 2,2 L 1,1 x y z A 1,2 2,2 J max J max L 2,2 L 1,1
Nominal model Model by substituted surfaces

J L  max 2,1  J L 2,1 max 1,1 2,1 2 L J L L          2,1 max 1,1 2,1 2 L J L L          Ay  z  max 2,1 J L  max 2,1  J L 2,1 max 1,1 2,1 2 L J L L          2,1 max 1,1 2,1 2 L J L L          Ay  z  Figure 8 : Clearance hull 1,2/ 2,2 D of joint 1,2/2,2 at point A.
Equation ( 4) is a relation between small displacement torsors. A similar equation ( 11) can be written for the hulls that correspond to point A. This gives:

1,1/2,1 1,1/1,2 1,2/2,2 2,2/2,1 =++ D D D D (11)
The sum of two hulls is a Minkowski sum [F1], [START_REF] Mujezinovi | A new mathematical model for geometric tolerances as applied to round faces[END_REF], [START_REF] Roy | Representation and interpretation of geometric tolerances for polyhedral objects[END_REF] and [Z1]. A representation of hull 1,1/2,1 D expressed at point A is shown in figure 9.

In particular, from this is deduced expression (12) formalising respect for the FC in the worst case in rigid behaviour.

1,2 2,1 2,1 max max 1,1 2,1 1,2 2,1 2,1 max min 1,1 2,1 22 22 t L t J dL L t L t J dL L                           (12) 1,2 2,1 2,1 max 1,1 2,1 22 t L t J L L             1,2 2,1 2,1 max 1,1 2,1 22 t L t J L L             1,2 2,1 max 1,2 2,1 2,1 tt J L L L    1,2 2,1 max 1,2 2,1 2,1 tt J L L L    Ay  z  1,2 2,1 2,1 max 1,1 2,1 22 t L t J L L             1,2 2,1 2,1 max 1,1 2,1 22 t L t J L L             1,2 2,1 max 1,2 2,1 2,1 tt J L L L    1,2 2,1 max 1,2 2,1 2,1 tt J L L L    Ay  z  Figure 9 : Minkowski sum illustrating equation (9).

2.2-Thermomechanical behaviour: dimension-chain integrating thermomechanical strains

2.2.1-Physical hypotheses

Let us consider figure 10 showing a system with rigid behaviour and thermomechanical behaviour.

The rigid behaviour of this system was studied in §2.1 where respect for the FC is characterised by (10). Henceforth, it must be ensured that the FC is respected under thermomechanical behaviour. The hypothesis that part 2 remains undeformable is stated: only part 1 is subject to thermomechanical strain.

An other hypothesis is stated: the topological structure of the joint graph in rigid behaviour and all thermomechanical behaviours remains the same. This means that there is no supplementary joint and no suppression of a joint between the two behaviours. The main consequence is that the joints that are influent on the FC remain the same for both behaviours. The influent cycle identified in figure 3 remains the same. On the other hand, the different parameters that characterise a joint (minimal clearance, maximal clearance and nature of contact) may change. As a result, the thermomechanical strains on the joints and the influent parts will modify the different tolerances which ensure that the FC is respected. By tolerances are more precisely the dimensions of the tolerance zones of coaxiality and the tolerance intervals of the diameters: see figure 11. A

2 2,2 2,3 2,1 1 1,2 1,3 1,1 x y z A thermomechanical behaviour x y z A 1 1,2 1,3 1,1 2 2,2 2,3 2,1 rigid behaviour 2 2,2 2,3 2,1 1 1,2 1,3 1,1 x y z A 2 2,2 2,3 2,1 1 1,2 1,3 1,1
D 1, 2 d 1, 2 th t 1,1th A L 1 , 1 t 2,2th B D 2,2 d 2,2th B L 2, 1 A D 1, 2 d 1, 2 th t 1,1th A L 1 , 1 t 2,2th B D 2,2 d 2,2th B L 2, 1 t 2,2th B D 2,2 d 2,2th B t 2,2th B D 2,2 d 2,2th B L 2, 1

2.2.2-Integration of thermomechanical strains into the joints

A joint can be defined as a set of parameters. Several studies have been carried out on this subject: [START_REF] Defazio | A prototype of feature-based design for assembly[END_REF], [W1], [START_REF] Shen | Analysis neutral data structure for GD&T[END_REF] and [START_REF] Dantan | Tolerance synthesis: quantifier notion and virtual boundary[END_REF]. In this article the definition proposed in [START_REF] Dufaure | A tolerancing framework to support geometric specifications traceability[END_REF], which is a direct application of that proposed in [START_REF] Dantan | Tolerance synthesis: quantifier notion and virtual boundary[END_REF] is used. Thus a joint can be defined by the following parameters:

-type: planar pair, cylindrical pair, ball and cylinder pair, etc.

-situation elements: plane, line, point -nature: fixed, sliding or floating -clearance: minimal and maximal clearance: J min , J max .

-forced contact direction: optional

The type of joint conforms to one of the types defined in [I2].

The situation elements position the joint in space, and hence define the axes of the degrees of freedom [START_REF] Clément | TTRS Declarative information model[END_REF], [START_REF] Dufaure | A tolerancing framework to support geometric specifications traceability[END_REF].

The nature of the contact determines the behaviour of the joint [DM1]:

-Fixed: no displacement is possible (the degrees of freedom are suppressed).

-Sliding: displacements corresponding to the different degrees of freedom are possible, but the surfaces remain in permanent contact one with another. Clearance of the joint is null.

-Floating: displacements corresponding to the different degrees of freedom are possible, and other displacements are limited by clearance.

In the case of a turning pair, cylindrical pair, ball and cylinder pair, spherical pair, or prismatic pair joint clearance is defined as being the difference between the diameter of the hole and the diameter of the shaft. Clearance may be positive, null or negative. In the case of a ball and plane pair, cylinder and plane pair or planar pair, clearance is the distance between the two surfaces potentially in contact. Clearance can then be positive or null: negative clearance has no meaning in physical terms.

Features :

-Name : 1,2/2,2 -Type : Cylindrical Pair (CP)

-Situation element : line (P, x)

-Nature : floating contact -Clearance : J max >J min  0 x y z 1 1,2 1,3 1,1 2 2,2 2,3 2,1 P J L 1, 2 max 2,2 J L  max 2,2 J L  max 2 J  max 2 J  Py  z 

Clearance space : Features :

-Name : 1,2/2,2 -Type : Cylindrical Pair (CP) -Situation element : line (P, x) A joint is defined in rigid behaviour by a set of parameters that act as reference parameters. Each thermomechanical behaviour will be characterised by a change in the reference parameters for the joint. Only those changes where the type of joint remains invariant are considered.

-Nature : floating contact -Clearance : J max >J min  0 x y z x y z 1 1,2 1,3 1,1 2 2,2 2,3 2,1 P J L 1, 2 max 2,2 J L  max 2,2 J L  max 2 J  max 2 J  Py  z  Py  z  Clearance space :
Figure 12 shows rigid behaviour for the cylindrical pair joint used in figure 3. The situation element is the straight line (P, x).

The hypothesis that minimal clearance is positive or null is stated. The cylindrical surface 1,2 (shaft) is free to occupy the clearance space in the cylindrical surface 2,2 (hole): the nature of the joint is therefore: floating. The corresponding clearance hull (in the shape of a lozenge at point P, centre of the joint) is also illustrated in figure 12.

In the case of a cylindrical pair joint, all calculations of relative displacements between the two surfaces are carried out from a position of reference where the two surfaces are coaxial.

Let us consider a thermal flux, characteristic of thermomechanical behaviour 1, which expands the diameter of surface 1,2: see figure 13. This expansion is such that maximal clearance is negative or null: surface 1,2 is therefore clamped to surface 2,2. The nature of the joint is therefore: fixed. Surface 1,2 is constrained to be coaxial to surface 2,2: the clearance hull is therefore reduced to a point centred on the origin.

Let us consider a thermal flux, characteristic of thermomechanical behaviour 2, which fixes surface 1,2 by strain in a particular position on surface 2,2: see figure 14. This case defines a forced contact. The minimal clearance of the joint remains strictly positive but surface 1,2 occupies a fixed position in relation to surface 2,2: the two cylindrical surfaces are in contact according to their generatrices so that the greatest distance separating the two surfaces is equal to J following y-axis. y-axis characterises the direction of forced contact. In the case of a forced contact, a supplementary parameter is added to the description of the joint specifying the direction of the forced contact (y-axis in this case). The corresponding clearance hull illustrated in figure 14 is therefore a point that corresponds to null rotation and to a translation of an amount J/2 following y-axis.

The thermomechanical behaviour corresponding to figure 10 corresponds to thermomechanical behaviour 1 in figure 13.

An example of the use of thermomechanical behaviour 2 will be shown in §3 in an industrial application.

2.2.3-Integration of thermomechanical strains in the parts

The surfaces that have undergone thermomechanical strain using substituted surfaces are modelled [START_REF] Pierre | Analyse des tolérances géométriques dans un contexte multi-expertises, application à une turbine de moteur d'hélicoptère[END_REF]. A thermomechanical simulation tool based on a technique for calculating by finite elements discretises the nominal model of a part into points and determines the displacement of each point. It is thus straightforward to describe a surface strained by a set of points. The substituted surface of a thermomechanically strained surface is a surface obtained by approximation using the least squares technique for this set of points.

Thus the relative position of the thermomechanically strained substituted surface and its nominal model using a small displacements torsor is determined.

From equation (1) and figure 2, equation (13) can be written as follows:

1,1/1,2 1,1/1,1th 1,1th /1,1n 1,1n /1,2n 1,2n /1,2th 1,2th /1,2 d d d d dd                             (13)
Surfaces 1,1th and 2,2th are the substituted surfaces for the thermomechanically strained surfaces. Torsors

1,1th /1,1n d   and 1,2th /1,2n
d   characterise the deviations of the relative position of surfaces 1,1th and 2,2th respectively in relation to the nominal surfaces 1,1n and 1,2n.

Given that the nominal model of a part is by definition without defect and without strain, (14) can be written as follows:

1,1/1,2 1,1/1,1th 1,1th /1,1n 1,2n /1,2th 1,2th /1,2 d d d dd                         (14)
(14) can be expressed as follows:

1,1/1,2 1,1/1,2 1,1/1,2 ma th d d d               (15)
where:

1,1/1,2 1,1/1,1th 1,2th /1,2 ma d d d              (16) 1,1/1,2 1,1th /1,1n 1,2n /1,2th th d d d              (17)
Equation ( 16) characterises the geometric deviations of surface 1,1 in relation to surface 1,2 that derive from defects in manufacture. Equation (17) characterises the geometric deviations of surface 1,1 in relation to surface 1,2 that derive from thermomechanical strains. When considering equation (2) the following hypothesis is stated:

1,1/1,2 1,1/1,2 ma d d defined in rigid behaviour          (18)
Let us suppose that:

1,1,/1,2th 1,1/1,2 th,A A,1,1,/1,2th 1,1/1,2th ,1,1/1,2th 1,1,/1,2th 1,1/1,2th 1,1,/1,2th ,1,1/1,2th 1,1/1,2th ,1,1/1,2th d with :
and

x A x y A y z A z                              ρ ε ρ  (19)
According to the deviation hull defined at point A in §2.1.3 for rigid behaviour, the deviation hull 1,1/1,2 D for thermodynamic behaviour described in figure 10 by the Minkowski sum (20) is defined as follows:

1,1/1,2 1,1/1,2th 1,1/1,2ma

1,1/1,2ma 1,1/1,2 = + with : = defined in rigid behaviour

D D D D D (20) 1,1/1,2th
D is a hull reduced to a point corresponding to the thermomechanical strain defined in (17): see figure 15.

Deviation induced by thermomechanical strains

Ay  z  1,1 ,1,1/1,2th 2 A y t    1,1 ,1,1/1,2th 2 A y t    ,1,1/1,2th A y  1,1/1,2thz  1,1 1,1/1,2th 1,1 z t L    1,1 1,1/1,2th 1,1 z t L   

Deviation induced by thermomechanical strains

Ay  z  1,1 ,1,1/1,2th 2 A y t    1,1 ,1,1/1,2th 2 A y t    ,1,1/1,2th A y  1,1/1,2thz  1,1 1,1/1,2th 1,1 z t L    1,1 1,1/1,2th 1,1 z t L   
Figure 15 : Deviation hull integrating thermomechanical strains.

2.2.4-Synthesis: formalisation of relations between functional condition, geometric specifications and contact specifications

Equation (11) characterises the dependency between the functional condition, the geometric specifications and the contact specifications by a Minkowski sum for hulls under rigid behaviour. This equation remains valid for thermomechanical behaviour, given the conservation hypothesis for the topological structure of the contact graph shown in §2.2.1. In particular, it is possible to deduce the expression (21) formalising respect for the worst case FC with thermomechanical behaviour in figure 10.

Ay  z  1,2 2,1 2,1 0 1,2 2,1 2 2 Ax t L t L L                    1,2 2,1 0 1,2 2,1 x t t L L    1,2 2,1 0 1,2 2,1 x t t L L    1,2 2,1 2,1 0 1,2 2,1 2 2 Ax t L t L L                    ,1,1/1,2th A y  1,1/1,2thz  Ay  z  1,2 2,1 2,1 0 1,2 2,1 2 2 Ax t L t L L                    1,2 2,1 0 1,2 2,1 x t t L L    1,2 2,1 0 1,2 2,1 x t t L L    1,2 2,1 2,1 0 1,2 2,1 2 2 Ax t L t L L                    ,1,1/1,2th A y  1,1/1,2thz 
1,2 2,1 2,1 max ,1,1/1,2th 1,2 2,1 1,2 2,1 2,1 min ,1,1/1,2th 1,2 2,1 22 22 Ay Ay t L t dL L t L t dL L                           (21)

3-Application to a high-pressure turbine

3.1-Description of the turbine

The turbine consists of two sub-units: a rotor 12 and a stator (made up of parts numbered 1, 2, 3, 4, and 7) in a turning pair joint along x by means of ball bearings, numbered 13 and a cylindrical roller bearing, numbered 14: see figure 17. Part 12 has a revolution shape where the largest diameter corresponds to the diameter of the heads of the blades. Clearance at the top of blade between part 12 and part 1 of the stator is characterised by the difference in diameter between part 1 and part 12: see figure 17.

In order to maximise turbine power, this clearance should be minimised while ensuring that in all turbine operational phases part 12 does not touch part 1. When the turbine is operational, a flux of hot gases (about 1000 °C) is created in the combustion chamber (not shown in figure 17) and the consequence of this is to transform parts 1, 4 and 12 into thermomechanically deformable parts: see figure 17. In this application, the study is limited to a sub-unit of the turbine made up of parts 1 and 2: see figure 18.

The functional condition guaranteeing control of clearance C at the top of the blades between the rotor and the stator (see figure 17), corresponds to an FC for the alignment of the axes of surfaces 1,8 and 2,8: see figure 18.

FC is defined by equation ( 22) where e represents the distance between the axes of surfaces 1,8 and 2,8 following a direction orthogonal to the x-axis at point A:

min max ee e  ( 22
)
In this example, when the turbine is operational, only thermomechanical strains in part 1 will be taken into account. Part 2 is considered infinitely rigid and geometrically perfect: see figure 18.

The specifications required to satisfy the FC with two different behaviours have to be studied. In a first phase the turbine is modelled under rigid behaviour where only manufacturing defects and clearances in the different contacts will be taken into account. Next, the thermomechanical behaviour is considered where strains due to the flux of hot gases from the combustion chamber are integrated, at a point where the turbine is in stationary operating regime.

y z=v 1 v 2 v 3 v 4 v 5 1,5  2,5 2,7 x  2 N 1 N 2 N 3 N 4 N 5 A 1,1 2,1 1,7 1,8 1,4  2,4 1,3  2,3 1,2  2,2 y z x A 2,8 C  D y z=v 1 v 2 v 3 v 4 v 5 1,5  2,5 2,7 x  2 N 1 N 2 N 3 N 4 N 5 A 1,1 2,1 1,7 1,8 1,4  2,4 1,3  2,3 1,2  2,

3.2-Rigid behaviour

Modelling 3D dimension chains of the turbine

Part 1 is placed in position by a support plane contact with part 2, via surfaces 1,7 and 2,7. Contacts between surfaces 1,1/2,1; 1,2/2,2; 1,3/2,3; 1,4/2,4; 1,5/2,5 modelled by ball and cylinder pair joints, complete the positioning between these two parts. Five clamping screws ensure that parts 1 and 2 are held in position. The graph of the joints corresponding to the system described in figure 18 is shown in figure 19 with the parameters of the different joints according to the model shown in §2.2.2. -Nature : floatting contact Note the deviations in position between surfaces 1,8 and 2,8 by a small displacement torsor 1,8 / 2,8 d   . From the joint graph, the six relations that characterise the relative position of surfaces 1,8 and 2,8 are deduced:

-Clearance : 0  J min  J max FC 2 2,1 2,8 2,2 2,7 2,5 2,4 2,3 1 1,3 1,4 1,2 1,5 1,7 1,8 1,1 2 2,1 2,8 2,2 2,7 2,5 2,4 2,3 2 2,1 2,8 2,2 2,7 2,5 2,4 2,3 1 1,3 1,4 1,2 1,5 1,7 1,8 1,1 1 1,3 1,4 1,2 1,5 1,7 1,8 1,1

Features of PP

1,8 / 2,8 1,8 /1,7 1,7 / 2,7 2,7 / 2,8 1,8 / 2,8 1,8 /1,i 1,i / 2,i 2,i / 2,8 d d d d d d d d with i [1;5]                                             (23)
Since part 2 is considered as infinitely rigid and without any geometric defect, torsors

2,7 / 2,7 d   2,i / 2,8 d i [1;7]    are null
torsors. Equations (23) then become:

1,8 / 2,8 1,8 /1,7 1,7 / 2,7 1,8 / 2,8 1,8 /1,i 1,i / 2,i d d d d d d with i [1;5]                                  (24)
In part 1 it is necessary to control the displacement of surface 1,8 simultaneously in relation to surfaces 1,1 to 1,7. To do this, a system of AB references according to ISO standards is defined [I4], [I1] and [I3]. The plane surface 1,7 designated A is taken as primary reference and the group of five cylindrical surfaces (1,1 to 1,5) designated B is taken as a secondary reference: see figure 20. AB designates the system of references that make up primary reference A and secondary reference B.

As for part 1, the CD reference system is defined on part 2. The primary reference designated as C corresponds to surface 2,7. The secondary reference designated as D corresponds to the group of 5 surfaces 2,1 to 2,5. The clearance hull characterising the displacement limits of surface 1,8 in relation to surface 2,8 is expressed by the following Minkowski sum:

y z 1,5 x N 1 N 2 N 3 N 4 N 5 1,1 1,7 1,8 1,4 1,3 1,2 y z x A t 1 B A 5D 1 d 1 t 1, 8 B D 1,8 d1,8 A  D L 1, 8 L 1, 7 B C y z 1,5 x N 1 N 2 N 3 N 4 N 5 1,1 1,7 1,8 1,4 1,3 1,2 y z x y z x A t 1 B t 1 B A 5D 1 d 1 t 1, 8 B D 1,8 d1,8 A  D  D L 1, 8 L 1, 7 B C
i iy i i ix ii i iz i i ix D J D i                                        D   (28)
i i i i i i i ix i ix i i i iy i i i iy i i iz i iz D D                                             i i ε ε N C ρ ε ρ N C          (29)

Modelling the location of surface 1,8

By the same method as that presented in §2.1.3, the deviation hull 1,8/ AB D is defined at point A. This hull characterising the location specification (see figure 20) is shown in figure 26. This hull remains invariant in relation to any translation orthogonal to the straight line  

1,8 1,8  t L 1,8 1,8  t L 1,8 2  t 1,8 2  t Az  y  1,8 1,8  t L 1,8 1,8  t L 1,8 2  t 1,8 2  t Az  y  Figure
A, x given that the tolerance zone for the location specification is a cylinder [I1]. It is axisymmetrical.

Relations between functional condition, geometric specifications and contact specifications

The Minkowski sum defined by equation ( 25) is represented by figure 27. In particular, expression (30) formalising respect for the worst case FC in rigid behaviour can be deduced.

3.3-Thermomechanical behaviour

The influence of the thermomechanical strains caused by a flux of gases from the combustion chamber will be modelled: see figure 17.

Integration of thermomechanical strains into the joints

As was described in §2.2.2, the nature of the contacts and the clearances can change under thermomechanical behaviour. Similarly, a contact may become a forced contact.

A preliminary thermomechanical study was carried out to determine changes in contacts. Thermal marginal conditions were provided by expert engineers working in the energy domain. They consisted of two fluxes: one convection flux due to hot gas from the combustion chamber, situated on the internal face of part 1 ; and convective cooling flux on the external face of part 1 : see figure 28. As mechanical marginal conditions only contact conditions of the planar pair 1,7/2,7 are taken into account. This will enable first to analyse the behaviour of surfaces 1,1 to 1,5. This first thermomechanical study shows that the radial displacement of the slugs is much greater than the maximal clearance present in each BCP type joint. A consequence is that the nature of the contact of the five ball and cylinder pair joints changes from floating to fixed. Moreover, the strain on part 1 adds a forced contact to each BCP type joint according to direction i v : see figure 29.

y z 1,5 x N 1 N 2 N 3 N 4 N 5 2,3 1,4 1,3 1,2 V 3 V 4 V 5 V 1 V 2 1,1 2,2 2,1 2,4 2,5 2,7 1,7 y z 1,5 x N 1 N 2 N 3 N 4 N 5 2,3 1,4 1,3 1,2 V 3 V 4 V 5 V 1 V 2 1,1 2,2 2,

Integrating thermomechanical strains into the parts

A thermomechanical analysis of part 1 was carried out with the thermal marginal conditions defined in figure 28 and with mechanical marginal conditions that conformed to the parameters of the different joints under thermomechanical behaviour. Next, the geometry of surface 1,8 which was deformed by thermomechanical stresses has been studied. After analysing the deformed cylindrical surface the results for figures 31 and 32 using Samcef software have been obtained. The deviations that derive from the thermomechanical stress on part 1 (see figure 31) give result (31) :



Figure 31 shows that the thermomechanical strain on surface 1,8 is not axisymmetrical.

Figure 33 shows hull In particular, from this expression (32) formalising respect for the worst case FC with thermomechanical behaviour is deduced.

4-Conclusions and future studies

This article shows how deviation hulls and clearance hulls can take thermomechanical strains into account in 3D dimension chains. This paper demonstrated an application using a high pressure turbine from a helicopter engine. For specific thermomechanical behaviour, thermomechanical strains are taken into account from a reference definition of the variability of geometrical defects corresponding to rigid behaviour. Future work is centred on studying the distribution of geometrical defects arising from manufacturing dispersions and thermomechanical strains, for all behaviour modes of a turboshaft engine. As well as rigid behaviour, several different thermomechanical behaviours should be considered to represent the true operational cycle of a helicopter engine. One of the main objectives of this study will be the possibility of proposing qualification criteria for the technical solutions envisaged in the design cycle of a high pressure turbine for helicopter engines assisted by geometric specifications.

Figure 1 :

 1 Figure 1: Geometric modelling of a real surface.

Figure 2 :

 2 Figure 2 : Relative position of two substitution surfaces of a part.

Figure 4 :

 4 Figure 4 : Definition of a geometric specification with substituted surfaces according to ISO 1101.

Figure 7 :

 7 Figure 7 : Definition of a joint with clearance by substituted surfaces

Figure 10 :

 10 Figure 10 : System behaviour.

Figure 11 :

 11 Figure 11 : Specifications for thermomechanical behaviour.

Figure 12 :

 12 Figure 12 : Cylindrical pair joint in rigid behaviour.

-Figure 13 :--Figure 14 :

 1314 Figure 13 : Cylindrical pair joint in thermomechanical behaviour 1.

Figure 16 :

 16 Figure 16 : Minkowski sum illustrating equation (11) with thermomechanical behaviour in figure 10. Hull 1,1/ 2,1 D at point A corresponding to the thermomechanical behaviour represented in figure 10 will be determined, taking into account the deviation hull 1,1/1,2 D and the clearance hull 1,2/ 2,2 D integrating the thermomechanical strains illustrated in figures 15 and 13 respectively and the deviation hull 2,2/ 2,1 D illustrated in figure 6, since part 2 is considered as being infinitely rigid in the hypotheses formulated in §2.2.1. The hull shown in figure 16 is obtained.

Figure 17 :

 17 Figure 17 : Simplified geometric model of a high-pressure turbine.

Figure 18 :

 18 Figure 18 : Sub-unit of high-pressure turbine.

  Ball and Cylinder Pair (BCP) -Situation element : point N i , line (N i , x)

Figure 19 :

 19 Figure 19 : Graph of joints.

Figure 20 :

 20 Figure 20 : Drawing showing definition of part 1.

Figure 21 :Figure 23 :Figure 24 :Figure 25 :

 21232425 Figure 21 : Clearance hull 1,i / 2,i D

26 :

 26 Figure 26 : Deviation hull of location of surface 1,8.

Figure 28 :

 28 Figure 28 : Thermal marginal conditions.

Figure 30 :

 30 Figure 30 : Clearance hull of complete joint in two dimensions.

Figure 31 :Figure 32 :

 3132 Figure 31 : Thermomechanical strain in part 1.

Figure 34 :

 34 Figure 34 : Minkowski sum illustrating relation (25) with thermomechanical behaviour in figure 28.

  It is also possible to determine the translation deviation 1,8() b  from a point on the cylindrical surface 1,8 in relation to its nominal position. For rigid behaviour, equation (33) characterizes the translation deviation 1,8() b  and the limits of ,1,8 / 2,8 Az  are given by the hull shown in figure 27. For thermomechanical behaviour, equation (34) characterizes the translation deviation 1,8 () b  and the limits of ,1,8 / 2,8 Az  are given by the hull in figure 34.

  represents the deviation hull characterised by the location of surface 1,8 in relation to the reference system AB.

	1,8/ 2,8 D	1,8/ AB  D	D	AB/CD	(25)
	D	1,8/ AB
	D	AB/ CD	represents the clearance hull defined by the support plane joint (1,7/2,7) and the five ball and cylinder pair joints
	(1,i/2,i). Hull AB/ CD D		is formalised by the intersection of the clearance hulls of the joint between surfaces 1,7/2,7 and the five
	joints 1,i/2,i:
	AB/ CD DD 1,7 / 2,7 
						1,1/ 2,1 D		1,2 / 2,2 D		1,3/ 2,3 D		1,4 / 2,4 D		1,5 / 2,5 D		(26)
			3.2.2 Expression of the clearance hull between parts 1 and 2
	For each slug the clearance is maximal: it corresponds to the difference between the maximum diameter of the hole and the
	minimum diameter of the shaft (cf. §2.1.1). Then the hull 1,i / 2,i D	at point i N according to (27) is defined: see figure 21. This
	clearance is open in accordance with the four projections	, , et y z     x	i N x	.	corresponding to the four degrees of freedom of a
	ball and cylinder pair joint.
	D	1, / 2, i i		    		ε	Ni,1, / 2, i i	  2 .. Ni,1, / 2, i i  y ε z		2		max 2  J  	2	    	(27)
	with 1	i 	5
	Equation (28) deduces from (27) hull 1,i / 2,i D	expressed at point C by considering (29).
																2
									C,1, / 2,		2	sin .	1, / 2,
	1, / 2,													2	2
									C,1, / 2,		2	cos .	1, / 2,	max 2
	with 1					5

  Hull defining displacement limits of surface 1,8 in relation to surface 2,8 at point A.

					 	Az Az
					1,8 22 max  t J 1,8 max 22  t J
				 	1,8 1,8 t L 1,8 1,8 t L	 	1,8 1,8 t L 1,8 1,8 t L	 	y y
					1,8 22 max  t J 1,8 max 22  t J
				Figure 27 :
	max e	 	1,8 22 max t J 	   	(30)
	min e	 	1,8 22 max t J 	  

  1,8/ 2,8th D at point A in a process similar to that used in §2.2.3. Hull strain in surface 1,8 in relation to the nominal surface. 3.3.3 Synthesis: formalisation of relations between functional condition, geometric specifications and contact specifications Using a similar method to that described in §2.2.4, based on relation (25), the Minkowski sum of hull 1,8/ 2,8 D characterising the displacement limits for surface 1,8 in relation to surface 2,8 in relation to thermomechanical behaviour for the high pressure turbine can be calculated: see figure 34. The strain on surface 1,8/2,8 is not axisymmetrical.

		 	Az Az
			 	y y
	  Figure 33 : 1,8 t 1,8  t 1,8 2  t  1,8 t 1,8  t 1,8 2  t Az Az
	1,8 L 1,8 L		1,8 L 1,8 L	
				y
	 	1,8 2 t 1,8 2 t
				22

Acknowledgements: We would like to thank Jean Luc Breining and Rainer Thomas from Turbomeca for their advice and technical assistance with this project.


Elisa Rebolini

Julien Toulouse

Andrew M Teale

Trygve Helgaker

Andreas Savin

Excited states from range-separated density-functional perturbation theory

We explore the possibility of calculating electronic excited states by using perturbation theory along a range-separated adiabatic connection. Starting from the energies of a partially interacting Hamiltonian, a first-order correction is defined with two variants of perturbation theory: a straightforward perturbation theory, and an extension of the Görling-Levy one that has the advantage of keeping the ground-state density constant at each order in the perturbation. Only the first, simpler, variant is tested here on the helium and beryllium atoms and on the dihydrogene molecule. The first-order correction within this perturbation theory improves significantly the total ground-and excited-state energies of the different systems. However, the excitation energies are mostly deteriorated with respect to the zeroth-order ones, which may be explained by the fact that the ionization energy is no longer correct for all interaction strengths. The second variant of the perturbation theory should improve these results but has not been tested yet along the range-separated adiabatic connection.

I. INTRODUCTION

In density-functional theory (DFT) of quantum electronic systems, the most widely used approach for calculating excitation energies is nowadays linear-response time-dependent density-functional theory (TDDFT) (see, e.g., Refs. 1, 2). However, in spite of many successes, when applied with the usual adiabatic semilocal approximations, linear-response TDDFT has serious limitations for describing systems with static (or strong) correlation [3], double or multiple excitations [4], and Rydberg and charge-transfer excitations [5,6]. Besides, the Hohenberg-Kohn theorem [7] states that the timeindependent ground-state density contains all the information about the system so that time dependence is in principle not required to describe excited states.

Several time-independent DFT approaches for calculating excitation energies exist and are still being developed. A first strategy consists in simultaneously optimizing an ensemble of states. Such an ensemble DFT was pioneered by Theophilou [8] and by Gross, Oliveira and Kohn [9] and is still a subject of research [10][11][12][13], but it is hampered by the absence of appropriate approximate ensemble functionals. A second strategy consists in selfconsistently optimizing a single excited state. This approach was started by Gunnarsson and Lundqvist [14], who extended ground-state DFT to the lowest-energy state in each symmetry class, and developed into the pragmatic (multiplet-sum) ∆SCF method [15,16] (still in use today [17]) and related methods [18][19][20]. Great efforts have been made by Nagy, Görling, Levy, Ayers * Electronic address: erebolini@kjemi.uio.no and others to formulate a rigorous self-consistent DFT of an arbitrary individual excited state [21][22][23][24][25][26][27][28][29][30][31][32][33] but a major difficulty is the need to develop approximate functionals for a specific excited state (see Ref. 34 for a proposal of such excited-state functionals). A third strategy, first proposed by Grimme, consists in using configurationinteraction (CI) schemes in which modified Hamiltonian matrix elements include information from DFT [35][36][37][38].

Finally, a fourth possible approach, proposed by Görling [39], is to calculate the excitation energies from Görling-Levy (GL) perturbation theory [40,41] along the adiabatic connection using the non-interacting Kohn-Sham (KS) Hamiltonian as the zeroth-order Hamiltonian. In this approach, the zeroth-order approximation to the exact excitation energies is provided by KS orbital energy differences (which, for accurate KS potentials, is known to be already a fairly good approximation [42]). It can be improved upon by perturbation theory at a given order in the coupling constant of the adiabatic connection. Filippi, Umrigar, and Gonze [43] showed that the GL first-order corrections provide a factor of two improvement to the KS zeroth-order excitation energies for the He, Li + , and Be atoms when using accurate KS potentials. For (nearly) degenerate states, Zhang and Burke [44] proposed to use degenerate second-order GL perturbation theory, showing that it works well on a simple one-dimensional model. This approach is conceptually simple as it uses the standard adiabatic connection along which the ground-state density is kept constant (in contrast to approaches introducing generalized adiabatic connections keeping an excited-state density constant [21,22,24,29]). In spite of promising early results, this approach has not been pursued further, perhaps because it can be considered an approximation to TDDFT [45].

In this work, we explore further this density-functional perturbation-theory approach for calculating excitation energies, introducing one key modification in comparison to the earlier work of Refs. 39, 43: As a zerothorder Hamiltonian, instead of using the non-interacting KS Hamiltonian, we use a partially interacting Hamiltonian incorporating the long-range part of the Coulomb electron-electron interaction, corresponding to an intermediate point along a range-separated adiabatic connection [START_REF] Savin | Recent Dev. Appl. Mod. Density Funct. Theory[END_REF][START_REF] Yang | [END_REF][48][49][50][51]. The partially interacting zeroth-order Hamiltonian is of course closer to the exact Hamiltonian than is the non-interacting KS Hamiltonian, thereby putting less demand on the perturbation theory. In fact, the zeroth-order Hamiltonian can already incorporate some static correlation.

The downside is that a many-body method such as CI theory is required to generate the eigenstates and eigenvalues of the zeroth-order Hamiltonian. However, if the partial electron-electron interaction is only a relatively weak long-range interaction, we would expect a faster convergence of the eigenstates and eigenvalues with respect to the one-and many-electron CI expansion than for the full Coulomb interaction [50,[START_REF] Franck | basis convergence of range-separated density-functional theory[END_REF], so that a small CI or multi-configuration self-consistent field (MCSCF) calculation would be sufficiently accurate.

When using a semi-local density-functional approximation for the effective potential of the range-separated adiabatic connection, the presence of an explicit longrange electron-electron interaction in the zeroth-order Hamiltonian also has the advantage of preventing the collapse of the high-lying Rydberg excitation energies. In contrast to adiabatic TDDFT, double or multiple excitations can be described with this density-functional perturbation-theory approach, although this possibility was not explored in Refs. 39, 43. Finally, approximate excited-state wave functions are obtained, which is useful for interpretative analysis and for the calculation of properties.

We envisage using this density-functional perturbation theory to calculate excited states after a rangeseparated ground-state calculation combining a longrange CI [START_REF] Leininger | [END_REF]54] or long-range MCSCF [55,56] treatment with a short-range density functional. This would be a simpler alternative to linear-response range-separated MCSCF theory [57,58] for calculations of excitation energies. In this work, we study in detail the two variants of range-separated density-functional perturbation theory and test the first, simpler variant on the He and Be atoms and the H 2 molecule, using accurate calculations along a range-separated adiabatic connection without introducing density-functional approximations.

Both variants of the range-separated perturbation theory are presented in Section II. Except for the finite basis approximation, no other approximation is introduced and the computational details can be found in Section III. Finally, the results obtained for the He and Be atoms, and for the H 2 molecule are discussed in Section IV. Section V contains our conclusions.

II. RANGE-SEPARATED DENSITY-FUNCTIONAL PERTURBATION THEORY

A. Range-separated ground-state density-functional theory

In range-separated DFT (see, e.g., Ref. 50), the exact ground-state energy of an N -electron system is obtained by the following minimization over normalized multi-determinantal wave functions Ψ

E 0 = min Ψ Ψ| T + Vne + Ŵ lr,µ ee |Ψ + Ēsr,µ Hxc [n Ψ ] , (1)
where we have introduced the kinetic-energy operator T , the nuclear.attraction operator Vne = v ne (r)n(r)dr written in terms of the density operator n(r), a longrange (lr) electron-electron interaction

Ŵ lr,µ ee = 1 2 w lr,µ ee (r 12 )n 2 (r 1 , r 2 )dr 1 dr 2 , (2)
written in terms of the error-function interaction w lr,µ ee (r) = erf(µr)/r and the pair-density operator n2 (r 1 , r 2 ), and finally the corresponding complement short-range (sr) Hartree-exchange-correlation density functional Ēsr,µ

Hxc [n Ψ ] evaluated at the density of Ψ: n Ψ (r) = Ψ|n(r)|Ψ . The parameter µ in the error function controls the separation range, with 1/µ acting as a smooth cut-off radius.

The Euler-Lagrange equation for the minimization of Eq. ( 1) leads to the (self-consistent) eigenvalue equation

Ĥlr,µ |Ψ µ 0 = E µ 0 |Ψ µ 0 , (3)
where Ψ µ 0 and E µ 0 are taken as the ground-state wave function and associated energy of the partially interacting Hamiltonian (with an explicit long-range electronelectron interaction)

Ĥlr,µ = T + Vne + Ŵ lr,µ ee + V sr,µ Hxc , (4)
which contains the short-range Hartree-exchangecorrelation potential operator,

V sr,µ Hxc = vsr,µ Hxc [n 0 ](r)n(r)dr, (5)
where vsr,µ Hxc [n](r) = δ Ēsr,µ Hxc [n]/δn(r), evaluated at the ground-state density of the physical system n 0 (r) = Ψ µ 0 |n(r)|Ψ µ 0 for all µ. For µ = 0, the Hamiltonian Ĥlr,µ of Eq. ( 4) reduces to the standard non-interacting KS Hamiltonian, Ĥlr,µ=0 = ĤKS , whereas, for µ → ∞, it reduces to the physical Hamiltonian Ĥlr,µ→∞ = Ĥ. Therefore, when varying the parameter µ between these two limits, the Hamiltonian Ĥlr,µ defines a range-separated adiabatic connection, linking the non-interacting KS system to the physical system with the ground-state density kept constant (assuming that the exact short-range Hartreeexchange-correlation potential vsr,µ Hxc (r) is used).

B. Excited states from perturbation theory

Excitation energies in range-separated DFT can be obtained by linear-response theory starting from the (adiabatic) time-dependent generalization of Eq. ( 1) [59].

Here, the excited states and their associated energies are instead obtained from time-independent many-body perturbation theory. In standard KS theory, the singledeterminant eigenstates and associated energies of the non-interacting KS Hamiltonian,

ĤKS |Φ KS k = E KS k |Φ KS k , (6)
where ĤKS = T + Vne + VHxc , give a first approximation to the eigenstates and associated energies of the physical Hamiltonian. To calculate excitation energies, two variants of perturbation theory using the KS Hamiltonian as zeroth-order Hamiltonian have been proposed [39,43]. We here extend these two variants of perturbation theory to range-separated DFT. As a first approximation, it is natural to use the excited-state wave functions and energies of the long-range interacting Hamiltonian

Ĥlr,µ |Ψ µ k = E µ k |Ψ µ k , (7)
where Ĥlr,µ = T + Vne + V sr,µ Hxc + Ŵ lr,µ ee is the same Hamiltonian that appears in Eq. ( 4) with the shortrange Hartree-exchange-correlation potential V sr,µ Hxc = vsr,µ Hxc [n 0 ]n(r)dr evaluated at the ground-state density n 0 . These excited-state wave functions and energies can then be improved upon by defining perturbation theories in which the Hamiltonian Ĥlr,µ is used as the zerothorder Hamiltonian.

First variant of perturbation theory

The simplest way of defining such a perturbation theory is to introduce the following Hamiltonian dependent on the coupling constant λ

Ĥµ,λ = Ĥlr,µ + λ Ŵ sr,µ , (8)
where the short-range perturbation operator is

Ŵ sr,µ = Ŵ sr,µ ee -V sr,µ Hxc , (9)
with the short-range electron-electron interaction Ŵ sr,µ ee = (1/2) w sr,µ ee (r 12 )n 2 (r 1 , r 2 )dr 1 dr 2 defined with the complementary error-function interaction w sr,µ ee (r) = erfc(µr)/r. When varying λ, Eq. ( 8) sets up an adiabatic connection linking the long-range interacting Hamiltonian at λ = 0, Ĥµ,λ=0 = Ĥlr,µ , to the physical Hamiltonian at λ = 1, Ĥµ,λ=1 = Ĥ, for all µ. However, the ground-state density is not kept constant along this adiabatic connection.

The exact eigenstates and associated eigenvalues of the physical Hamiltonian can be obtained by standard Rayleigh-Schrödinger perturbation theory-that is, by Taylor expanding the eigenstates and eigenvalues of the Hamiltonian Ĥµ,λ in λ and setting λ = 1:

|Ψ k = |Ψ µ k + ∞ n=1 |Ψ µ,(n) k , (10a)
E k = E µ k + ∞ n=1 E µ,(n) k , (10b)
where

Ψ µ k ≡ Ψ µ,(0) k and E µ k ≡ E µ,(0) k
act as zerothorder eigenstates and energies. Using orthonormalized zeroth-order eigenstates Ψ µ k |Ψ µ l = δ kl and assuming non-degenerate zeroth-order eigenstates, the first-order energy correction for the state k becomes

E µ,(1) k = Ψ µ k | Ŵ sr,µ |Ψ µ k . ( 11
)
As usual, the zeroth-plus-first-order energy is simply the expectation value of the physical Hamiltonian over the zeroth-order eigenstate:

E µ,(0+1) k = E µ k + E µ,(1) k = Ψ µ k | Ĥ|Ψ µ k . ( 12
)
This expression is a multi-determinantal extension of the exact-exchange KS energy expression for the state k, proposed and studied for the ground state in Refs. 60-62.

The second-order energy correction is given by

E µ, (2) k
= - l =k | Ψ µ l | Ŵ sr,µ |Ψ µ k | 2 E µ l -E µ k , (13)
whereas the first-order wave-function correction is given by (using intermediate normalization so that

Ψ µ k |Ψ µ,(n) k = 0 for all n ≥ 1) |Ψ µ, (1) k
= - l =k Ψ µ l | Ŵ sr,µ |Ψ µ k E µ l -E µ k |Ψ µ l . (14)
For µ = 0, this perturbation theory reduces to the first variant of the KS perturbation theory studied by Filippi et al., see Eq. ( 5) of Ref. 43.

We now give the behaviors of the zeroth+first-order energies with respect µ near the KS system (µ = 0) and near the physical system (µ → ∞), which are useful to understand the numerical results in Section IV. The total energies up to the first order of the perturbation theory are given by the expectation value of the full Hamiltonian over the zeroth-order wave functions in Eq (11). Using the Taylor expansion of the wave function

Ψ µ k = Φ KS k + µ 3 Ψ (3)
k + O(µ 5 ) around the KS wave function [51], it implies that the zeroth+first-order energies are thus given by

E µ,(0+1) k = Φ KS k | Ĥ|Φ KS k + 2µ 3 Φ KS k | Ĥ|Ψ (3) k + O(µ 5 ), (15)
where Ψ

k is the contribution entering at the third power of µ in the zeroth-order wave function.

¿From the asymptotic expansion of the wave function

Ψ µ k = Ψ k + µ -2 Ψ (-2) k + O(µ -3
), which is valid almost everywhere when µ → ∞ (the electron-electron coalescence needs to be treated carefully) [51], the first correction to the zeroth+first-order energies enter at the fourth power of µ E µ,(0+1) k

= E k + 1 µ 4 E (0+1,-4) k + O 1 µ 6 , (16)
where

E (0+1,-4) k
is the contribution entering at the fourth power of 1/µ.

Second variant of perturbation theory

A second possibility is to define a perturbation theory based on a slightly more complicated adiabatic connection, in which the ground-state density is kept constant between the long-range interacting Hamiltonian and the physical Hamiltonian, see Appendix A. The Hamiltonian of Eq. ( 8) is then replaced by

Ĥµ,λ = Ĥlr,µ + λ Ŵ sr,µ -V sr,µ,λ c,md , (17)
where Ŵ sr,µ is now defined as

Ŵ sr,µ = Ŵ sr,µ ee -V sr,µ Hx,md , (18)
in terms of a short-range "multi-determinantal (md) Hartree-exchange" potential operator

V sr,µ Hx,md = δE sr,µ Hx,md [n 0 ] δn(r) n(r) dr, (19)
and a short-range "multi-determinantal correlation" potential operator

V sr,µ,λ c,md = δE sr,µ,λ c,md [n 0 ] δn(r) n(r) dr, (20)
that depends non-linearly on λ so that the ground-state density n 0 is kept constant for all µ and λ. The density functionals E sr,µ Hx,md [n] and E sr,µ,λ c,md [n] are defined in Appendix A.

One can show that, for non-degenerate ground-state wave functions Ψ µ 0 , the expansion of V sr,µ,λ c,md in λ for λ → 0 starts at second order:

V sr,µ,λ c,md = λ 2 V sr,µ,(2) c,md + • • • . (21)
Hence, the Hamiltonian of Eq. ( 17) properly reduces to the long-range Hamiltonian at λ = 0, Ĥµ,λ=0 = Ĥlr,µ , whereas, at λ = 1, it correctly reduces to the physical Hamiltonian, Ĥµ,λ=1 = Ĥ. This is so because the short-range Hartree-exchange-correlation potential in the Hamiltonian Ĥlr,µ can be decomposed as V sr,µ Hxc = V sr,µ Hx,md + V sr,µ c,md ,

where V sr,µ c,md = V sr,µ,λ=1 c,md is canceled by the perturbation terms for λ = 1. Equation ( 22) corresponds to an alternative decomposition of the short-range Hartree-exchangecorrelation energy into "Hartree-exchange" and "correlation" contributions based on the multi-determinantal wave function Ψ µ 0 instead of the single-determinant KS wave function Φ KS 0 [60][61][62], which is more natural in range-separated DFT. This decomposition is especially relevant here since it separates the perturbation into a "Hartree-exchange" contribution that is linear in λ and a "correlation" contribution containing all the higher-order terms in λ.

The first-order energy correction is still given by Eq. ( 11) but with the perturbation operator of Eq. ( 18), yielding the following energy up to first order:

E µ,(0+1) k = E µ k + E µ,(1) k = Ψ µ k | Ĥ + V sr,µ c,md |Ψ µ k . ( 23
)
The second-order energy correction of Eq. ( 13) becomes

E µ,(2) k = - l =k | Ψ µ l | Ŵ sr,µ |Ψ µ k | 2 E µ l -E µ k -Ψ µ k | V sr,µ, (2)
c,md

|Ψ µ k , (24)
whereas the expression of the first-order wave function correction is still given by Eq. ( 14) but with the perturbation operator of Eq. ( 18). For µ = 0, this density-fixed perturbation theory reduces to the second variant of the KS perturbation theory proposed by Görling [39] and studied by Filippi et al. [Eq. ( 6) of Ref. 43], which is simply the application of GL perturbation theory [40,41] to excited states. In Ref. 43, it was found that the first-order energy corrections in density-fixed KS perturbation theory provided on average a factor of two improvement on the KS zerothorder excitation energies for the He, Li + , and Be atoms when using accurate KS potentials. By contrast, the first-order energy corrections in the first variant of KS perturbation theory, without a fixed density, deteriorated on average the KS excitation energies.

The good results obtained with the second variant of KS perturbation theory may be understood from that fact that, in GL perturbation theory, the ionization potential remains exact to all orders in λ. In fact, this nice feature of GL theory holds also with range separation, so that the second variant of rangeseparated perturbation theory should in principle be preferred. However, it requires the separation of the short-range Hartree-exchange-correlation potential into the "multi-determinantal Hartree-exchange" and "multideterminantal correlation" contributions (according to Eq. ( 22)), which we have not done for accurate potentials or calculations along the double adiabatic connection with a partial interaction defined by Ŵ lr,µ ee + λ Ŵ sr,µ ee (cf. Appendix A). We therefore consider only the first variant of range-separated perturbation theory here but note that the second variant can be straightforwardly (in hartree) of the helium atom as a function of µ (in bohr -1 ). The zeroth-order energy E µ 0 is recalled for the ground state in plain line and the total energies of the physical system E k are plotted as horizontal dotted lines.

applied with density-functional approximations-using, for example, the local-density approximation that has been constructed for the "multi-determinantal correlation" functional [60,63].

III. COMPUTATIONAL DETAILS

Calculations were performed for the He and Be atoms and the H 2 molecule with a development version of the DALTON program [START_REF] Angeli | DALTON, a molecular electronic structure program[END_REF], see Refs. 65-67. Following the same settings as in Ref. 51, a full CI (FCI) calculation was first carried out to get the exact ground-state density within the basis set considered. Next, a Lieb optimization of the short-range potential v sr,µ (r) was performed to reproduce the FCI density with the longrange electron-electron interaction w lr,µ ee (r 12 ). Then, an FCI calculation was done with the partially-interacting Hamiltonian constructed from w lr,µ ee (r 12 ) and v sr,µ (r) to obtain the zeroth-order energies and wave functions according to Eq. (7). Finally, the zeroth+first order energies were calculated according to Eq. ( 12). The basis sets used were: uncontracted t-aug-cc-pV5Z for He, uncontracted d-aug-cc-pVDZ for Be, and uncontracted d-aug-cc-pVTZ for H 2 .

IV. RESULTS AND DISCUSSION

A. Helium atom

The ground-and excited-state total energies to first order of helium along the range-separated adiabatic connection are shown in Figure 1. In the KS limit, when µ = 0, the total energies are significantly improved with respect to the zeroth-order ones. In fact, as shown for the ground-state energy, the zeroth-order total energies were off by approximately 1.2 hartree with respect to the energies of the physical system. When the first-order correction is added, the error becomes smaller than 0.06 hartree for all states. Moreover, the singlet and triplet excitedstate energies are no longer degenerate. With increasing the range-separation parameter µ, a faster convergence towards the total energies of the physical system is also observed for all states.

The description of the total energies is therefore much improved with the addition of the first-order correction. The linear term in µ present in the zeroth-order total energies [51] is no longer there for the zeroth+first order total energies, which instead depend on the the third power of µ for small µ (cf. Eq. ( 15)). At large µ, the error relative to the physical energies enters as 1/µ 4 rather than as 1/µ 2 in the zeroth-order case, explaining the observed faster convergence of the first-order energies.

The excitation energies of the helium atom correct to zeroth and first orders are plotted in Figure 2. As previously noted, at µ = 0, the degeneracy of the zeroth-order singlet and triplet excitation energies is lifted by the firstorder correction, However, the excitation energies correct to first order overestimate the physical excitation energies by 0.1-0.2 hartree such that the error is actually larger than at zeroth order. For the 1 1 S → 1 3 P excitation energy, the correction is even going in the wrong direction and the singlet-triplet splitting is too large by about a factor 1.5.

When the very long-range part of the Coulombic interaction is switched on with positive µ close to 0, the initial overestimation is corrected. In fact, for small µ, all the excitation energies decrease in the third power of µ which is in agreement with Eq. (15). When µ ≃ 0.5-1, this correction becomes too large and the excitation energies of the partially interacting system become lower than their fully interacting limits. As µ increases further so that more interaction is included, the excitation ener- Excitation energies in hartree 1 1 S → 1 3 P 0th-order 1 1 S → 1 3 P 0th+1st order 1 1 S → 1 1 P 0th-order 1 1 S → 1 1 P 0th+1st order (in hartree) of the helium atom as a function of µ (in bohr -1 ). The excitation energies of the physical system ∆E k = E k -E0 are plotted as horizontal dotted lines.

gies converge toward their fully interacting values from below. The zeroth-order excitation energies, which do not for small µ, converge monotonically toward their physical limit and are on average more accurate than the zeroth+first order excitation energies. In short, the first-order correction does not improve excitation energies, although total energies are improved.

The inability of the first-order correction to improve excitation energies should be connected to the fact that, since the ground-state density is not kept constant at each order in the perturbation, the ionization potential is no longer constant to first order along the adiabatic connection. This behavior results in an unbalanced treatment of the ground and excited states. Moreover, high- Excitation energies in hartree 1 1 S → 1 3 P 0th order 1 1 S → 1 3 P 0th+1st order 1 1 S → 1 1 P 0th order 1 1 S → 1 1 P 0th+1st order energy Rydberg excitation energies should be even more sensitive to this effect, as observed for transitions to the P state. The second variant of perturbation theory should correct this behavior by keeping the density constant at each order, as shown in the KS case [41,43].

B. Beryllium atom

When the first-order perturbation correction is applied to the ground-state and valence-excited states of beryllium, total energies are again improved (although not illustrated here). In Figure 3, we have plotted the zerothand first-order valence excitation energies of beryllium against the range-separation parameter µ.

Since valence excitation energies should be less sensitive to a poor description of the ionization energy than Rydberg excitation energies, the first-order correction should work better for the beryllium valence excitations than for the helium Rydberg excitations. However, although the singlet excitation energy of beryllium is improved at µ = 0 at first order, the corresponding triplet excitation energy is not improved. In fact, whereas the triplet excitation energy is overestimated at zeroth order, it is underestimated by about the same amount at first order.

As the interaction is switched on, a bump is observed for small µ for the singlet excitation energy but not the triplet excitation energy, which converges monotonically to its physical limit. The convergence of the excitation energies with µ is improved by the first-order excitation energies, especially in the singlet case. Excitation energies in hartree

1 1 Σ + g → 1 3 Σ + u 0th order 1 1 Σ + g → 1 3 Σ + u 0th+1st order 1 1 Σ + g → 1 1 Σ + u 0th order 1 1 Σ + g → 1 1 Σ + u 0th+1st order -0.1 0 0.1 0.2 0.3 0 1 2 3 4 5
µ in bohr -1 Excitation energies in hartree (dashed line) excitation energies of the hydrogen molecule (in hartree) as a function of µ in bohr -1 at the equilibrium distance (left) and three times the equilibrium distance (right). The excitation energies of the physical system ∆E k are plotted as horizontal dotted lines.

1 1 Σ + g → 1 3 Σ + u 1 1 Σ + g → 1 3 Σ + u 1 1 Σ + g → 1 1 Σ + u 1 1 Σ + g → 1 1 Σ + u 1 1 Σ + g → 2 3 Σ + g (σ + u ) 2 1 1 Σ + g → 2 3 Σ + g (σ + u ) 2

C. Hydrogen molecule

In Figure 4, we have plotted the excitation energies of H 2 as a function of µ at the equilibrium distance R eq and at 3R eq . At the equilibrium geometry, the first-order correction works well. At µ = 0, the correction is in the right direction (singlet and triplet excitation energies being raised and lowered, respectively); for nearly all µ > 0, the error is smaller than for the zeroth-order excitation energies.

Unfortunately, when the bond is stretched, this is no longer the case. At the stretched geometry, the first excitation energy 1 1 Σ + g → 1 3 Σ + u becomes negative for small values of µ and the error with respect to the physical excitation energy is higher than in the zeroth-order case. Moreover, the ordering of the two singlet excitation energies is incorrect at small µ and they present a strong oscillation when the interaction is switched on. In this case, therefore, the zeroth-order excitation energies are better approximations to the physical excitation energies.

V. CONCLUSION

We have considered two variants of perturbation theory along a range-separated adiabatic connection. The first and simpler variant, based on the usual Rayleigh-Schrödinger perturbation theory, was tested on the helium and beryllium atoms and on the hydrogen molecule at equilibrium and stretched geometries. Although total energies are improved to first order in the perturbation, excitation energies are not improved since the theory does not keep the density constant along the adiabatic connection at each order of perturbation. It would be interesting to examine the evolution of the ionization potential to understand better the effect of this variant of the perturbation theory on our systems of interest.

The second variant of the perturbation theory, based on Görling-Levy theory, should improve the results significantly by keeping the ground-state density constant at each order in the perturbation [41], as already observed on the KS system [43]. However, this more complicated theory has not yet been implemented for µ > 0.

An interesting alternative to perturbation theory is provided by extrapolation, which make use of the behavior of the energies with respect to µ near the physical system to estimate the exact energies from the energy of the partially interacting system at a given µ and its first-or higher-order derivatives with respect to µ [68,69]. Work using this approach will be presented elsewhere.

1 S 1 S → 2 1

 111 → 2 3 S 0th-order 1 1 S → 2 3 S 0th+1st order 1 1 S → 2 1 S 0th-order 1

Figure 2 :

 2 Figure 2: Zeroth-order (plain line) excitation energies ∆E µ k = E µ k -E µ 0 and zeroth+first-order (dashed line) excitation energies ∆E µ,(0+1) k = E µ,(0+1) k -E µ,(0+1) 0

Figure 3 :

 3 Figure 3: Valence excitation energies of the beryllium atom (in hartree) at zeroth order ∆E µ k (plain line) and zeroth+first order ∆E µ,(0+1) k(dashed line), as a function of µ (in bohr -1 ). The excitation energies of the physical system ∆E k are plotted as horizontal dotted lines.

Figure 4 :

 4 Figure 4: Zeroth-order ∆E µ k (plain line) and zeroth+first-order ∆E µ,(0+1) k

Acknowledgements

This work was supported by the Norwegian Research Council through the CoE Centre for Theoretical and Computational Chemistry (CTCC) Grant No. 179568/V30 and and through the European Research Council under the European Union Seventh Framework Program through the Advanced Grant ABACUS, ERC Grant Agreement No. 267683.

Appendix A: Double adiabatic connection with a constant density

We here present a double adiabatic connection, depending on two parameters, that keeps the groundstate density constant. It is the basis for the perturbation theory presented in Section II B 2. A different density-fixed double adiabatic connection was considered in Refs. 70, 71.

The Levy-Lieb universal density functional for the Coulomb electron-electron interaction Ŵee is given by [72][73][74]

We here generalize it to the interaction Ŵ lr,µ ee + λ Ŵ sr,µ ee , where Ŵ lr,µ ee and Ŵ sr,µ ee are long-range and short-range electron-electron interactions, respectively, that depend on both a range-separation parameter µ and on a linear parameter λ:

The total universal density functional F [n] is then decomposed into F µ,λ [n] and a (µ, λ)-dependent shortrange Hartree-exchange-correlation density functional Ēsr,µ,λ Hxc [n],

giving the following expression for the exact ground-state energy of the electronic system

where the minimization is over normalized multideterminantal wave functions. The Euler-Lagrange equation corresponding to this minimization is is the short-range Hartree-exchange-correlation potential operator, evaluated at the ground-state density of the physical system at µ and λ, n 0 (r) = Ψ µ,λ 0 |n(r)|Ψ µ,λ 0 . The Hamiltonian Ĥµ,λ thus sets up a double adiabatic connection with a constant ground-state density.

The range-separated ground-state DFT formalism of Section II A is recovered in the limit λ = 0. To set up a perturbation theory in λ about 0, we rewrite Ĥµ,λ of Eq. (A6) as the sum of the noninteracting Hamiltonian Ĥlr,µ = Ĥµ,λ=0 and a perturbation operator. For this purpose, the Hartree-correlation-exchange functional is written as Ēsr,µ,λ Hxc [n] = Ēsr,µ,λ=0 is the short-range Hartree-exchange-correlation potential operator associated with E sr,µ,λ Hxc [n]. The dependence on λ of E sr,µ,λ Hxc [n] can be made more explicit. It is easy to show that

which leads to the following decomposition

where

is a multi-determinantal (md) generalization of the usual short-range Hartree-exchange functional [60][61][62]. Using the variational property of the wave function Ψ µ,λ [n], and for non-degenerate wave functions Ψ µ,λ=0 [n], the expansion of E sr,µ,λ c,md [n] in λ about 0 starts at second order:

as in standard GL perturbation theory [40,41]. The Hamiltonian of Eq. (A9) can now be rewritten as Ĥµ,λ = Ĥlr,µ + λ Ŵ sr,µ -V sr,µ,λ c,md ,

where the perturbation operator Ŵ sr,µ = Ŵ sr,µ ee -V sr,µ Hx,md and V sr,µ Hx,md = δE sr,µ Hx,md [n 0 ] δn(r) n(r) dr (A16) has been introduced to collect all the linear terms in λ, the remaining perturbation operator V sr,µ,λ c,md = δE sr,µ,λ c,md [n 0 ] δn(r) n(r) dr (A17) containing all higher-order terms in λ.