B Burgardt

P Cartraud

Continuum modeling of beamlike lattice trusses using averaging methods

Keywords: Truss structures, Continuum modeling, Beam-like structures, Homogenization

Introduction

Lattice structures have been used for many years in civil engineering, and they are candidates for future applications in space, due to their low cost, light weight and high stiness. A discrete ®nite element model of a large repetitive structure can become computationally expensive, and inappropriate to study the overall structural response and the vibration modes. This is why simpli®ed models have to be developed.

Several methods have been discussed in the literature, and they are reviewed in the survey paper [START_REF] Noor | Continuum modeling for repetitive lattice structures[END_REF]. Among them, the substitute continuum approach has several advantages, and has proved to be suitable for the analysis of large repetitive lattices [START_REF] Noor | Continuum modeling for repetitive lattice structures[END_REF]. In this ®eld, the approaches used can be grouped into two main classes.

The ®rst one includes methods in which the continuum behaviour is assumed a priori (e.g. a Timoshenko-beam for a beam-like lattice). Then, the problem amounts to ®nding the equivalent continuum properties. The equivalence can be de®ned on the energies [1±5], or on the constitutive relations, by transforming the force±displacement relationships in the repetitive lattice into constitutive equations in the continuum model [START_REF] Sun | Global±local approach to solving vibration of large truss structures[END_REF].

In the second class of methods, the governing laws of the continuum are obtained from the discrete nodal equations. Then, either an asymptotic expansion and a double-scale writing [START_REF] Verna | Mode lisation des structures discreÁ tes par homoge ne isation: cas des treillis[END_REF], or a ®nite dierence operator relating the displacements of two consecutive lattice sections [START_REF] Renton | The beam-like behavior of space trusses[END_REF] are used, or both [START_REF] Tollenaere | Mode les bidimensionnels de tisse s. Homoge ne isation des treillis en vibrations libres[END_REF].

In this paper, the general framework of the energy equivalence between the lattice and the continuum model is used to derive the equivalent properties. Usually, the main steps of such an approach consist of [1±4]:

1. Isolating a typical (repeating) cell. 2. Expanding the nodal displacements in a Taylor series. 3. Identifying the strain parameters for the continuum. 4. Transforming the energies expressed in the nodal quantities into continuum displacements and strains. 5. Computing the eective continuum properties.

The method presented herein follows this procedure, but with changes in steps 2, 3 and 5. A new displacement ®eld expression is introduced based on piecewise linear functions in order to describe the exact displacement ®eld in the trusses more accurately. Besides, in connection with the homogenization theory of periodic materials [10], the continuum strain and stress parameters are de®ned by their average values over the cell. The characterization of the continuum stresses represents the most signi®cant improvement to the methods based on the energy equivalence proposed in the literature [1±5]. Then, a rigorous partitioning of the strain parameters is possible, so as to derive the equivalent continuum properties, and take the lattice periodicity into account. Moreover, from the matrix formulation obtained, it will be shown that the static condensation which is widely used [2±5] to obtain the eective continuum properties, is an approximation.

In the present paper, the method is applied to plane pin-jointed beam-like lattices in static analysis, but it can also be used for three-dimensional lattices in dynamic analysis.

The expressions for the displacement ®eld and the de®nition of the continuum strain and stress parameters are presented in Section 2. From the matrix formulation in Section 3, dierent methods to determine the constitutive relations for the Timoshenkobeam model are proposed in Section 4, and comparisons between the methods are made in Section 5. Finally, numerical examples are given in Section 6 and prove the reliability of our approach.

Preliminary developments

Expressions of displacements

The starting point for the continuum modeling of lattice trusses is usually the introduction of a displacement ®eld, which is herein based on the representation of the deformed shape of the bay boundary. In order to illustrate our approach, let us ®rst consider a singlebay lattice (only one bay on the cell length). The case of the double-bay lattice will be treated afterwards.

Single-bay lattice

A study is made of the planar rectangle of length l and height h shown in Fig. 1, with the origin of the coordinate system at the cell center.

All the truss members are pin-jointed, and the cell can either be a single-or a double-laced lattice truss (only the vertical and the horizontal trusses have been represented). Then, a ®nite element model introduces 8 degrees of freedom.

It can be seen that the following displacement ®eld expressions ((x,y ) are the coordinates of a truss point):

ux,y a 0 a 1 x a 2 y a 3 xy vx,y b 0 b 1 x b 2 y b 3 xy 1
gives an exact representation of the deformed shape of the bay boundary. Therefore, the number of parameters introduced is equal to the total number of degrees of freedom of the cell. Moreover, in each truss, the displacement ®eld is linear, and if we consider the two vertical (respectively horizontal) trusses, the axial displacement v (respectively u ) varies linearly in y (respectively in x ). As a result, Eq. ( 1) gives an exact description of the displacement ®eld for these trusses parallel to the coordinate axes. However, for the diagonal trusses, the axial displacement is found from Eq. ( 1) to be an approximation in the form of a polynomial of second degree in x', the local abscissa of the truss.

If the continuum which completely covers the cell is considered, the displacement ®eld in it is also given by Eq. ( 1), and is extended to the continuum domain, with (x,y )

$ S c =[Àl/2,l/2] Â [Àh/2,h/2].
After dierentiation, one can easily obtain (the subscript ,x indicates a partial derivative with respect to x ):

ux,y u 0 xE x y 1 2 g xy À r 0 xyE x,y vx,y v 0 x 1 2 g xy r 0 yE y xyE y,x 2
where the rigid body displacement (u 0 ,v 0 ), the rotation r 0 =(v ,x Àu ,y )/2 and the strains are evaluated at the cell origin.

Double-bay lattice

Let us now add a second bay (of the same length for the sake of simplicity) to the previous single-bay lattice, as shown in Fig. 2, where O i is the origin and (x i ,y ) the coordinates of an ith bay point.

The length of the cell is again noted, l, with each bay of length l/2. The displacement ®eld of this double-bay lattice is simply written from the results of the previous section, considering each bay separately, hence Eq. ( 2), with eight parameters in each bay. The number of parameters is easily reduced to 12 (the number of degrees of freedom of the cell) and the continuity of the displacement ®eld at the junction of the two bays is taken into account. The following equations are then obtained:

u 1 0 l 4 E 1 x u 2 0 À l 4 E 2 x v 1 0 l 4 1 2 g 1 xy r 1 0 v 2 0 À l 4 1 2 g 2 xy r 2 0 1 2 g 1 xy À r 1 0 l 4 E 1 x,y 1 2 g 2 xy À r 2 0 À l 4 E 2
x,y

E 1 y l 4 E 1 y,x E 2 y À l 4 E 2 y,x 3
where the superscript i stands for the bay i.

The previous displacement ®eld expression has the main advantage of exactly representing the deformed shape of the bay boundary and the displacement in the trusses parallel to the coordinate axes, as is also the case for the single-bay lattice. This displacement is approximate only for the other trusses.

Comparison with other methods

In the literature, the problem of ®nding the displacement ®eld is generally treated by means of Taylor series expansion in the axial direction of the equivalent beam. This expansion can be made directly on the displacement ®eld [3±4], or on the strain ®eld [1±2]. In both cases, the method consists of representing the displacements by polynomials, the order of which is dependent on the cell con®guration.

For a single-bay lattice, such an approach leads to a linear displacement ®eld, equivalent to that given in Eq. (1) (see Ref. [START_REF] Noor | Anisotropic continuum models for beamlike lattice trusses[END_REF], where the rigid body rotation is expressed with a parameter dierent from r 0 ).

In the case of a double-bay lattice, the Taylor series expansion method leads to quadratic terms in x in the displacement ®eld expression [START_REF] Dow | Continuum models of space station structures[END_REF]. It turns out that the deformed shape of the bay boundary is approximated, while it is exactly described by the piecewise linear functions proposed in this paper.

Therefore, it can be concluded that, for a multiplebay lattice, our approach gives a more accurate approximation of the displacement ®eld.

Macroscopic strains

The next step to building a substitutive continuum model is to de®ne the strains associated with the overall behaviour of the cell. As in the homogenization methods of periodic materials [10], we choose to work with the average values of the local strains over the continuum cell.

Single-bay lattice

Since the origin of the coordinate system is at the cell center, the average value of each component of strains or strains gradient in Eq. ( 2) is equal to its value at the cell origin. Then, if h Á i stands for the averaging operator de®ned over the continuum cell, we have:

hai 1 jS c j S c a dS c 4
where v S c v is the surface of the continuum cell, i.e. v S c v =lh in our case (a unit thickness is assumed). Eq. ( 2) is now written as follows:

ux,y u 0 xhE x i y 1 2 hg xy i À r 0 xyhE x,y i vx,y v 0 x 1 2 hg xy i r 0 yhE y i xyhE y,x i 5

Double-bay lattice

The average value over the cell is calculated from the intermediate results obtained for each bay. Let h Á i i denote the averaging operator over the ith bay; then, for two bays of the same length l/2 we have:

hai hai 1 hai 2 =2 6
De®ning a as:

a hai À hai 1 hai 2 À hai 7
Eq. ( 3) becomes:

l 4 hE x i v 0 l 4 1 2 hg xy i hr 0 i 1 2 g xy À r 0 l 4 hE x,y i E y l 4 hE y,x i 8
In each bay, Eq. ( 7) can be used to express the displacement ®eld as a function of hai and a , where a is one of the following items: {u 0 , v 0 , r 0 , E x , E y , g xy , E x,y , E y,x }. The number of parameters can be reduced from 16 to 12 by eliminating u 0 , v 0 , 1 2 g xy À r 0 and E y using Eq. ( 8).

Comparison with other methods

For the methods based on a Taylor series expansion, the strain and strain gradient parameters for the continuum are de®ned by their value at the cell center. In the case of the single-bay lattice, our approach is equivalent to the preceding ones, both for the displacement ®eld (see Subsection 2.1.3) and for the continuum strain parameters.

For the multiple-bay lattice, the comparison is more dicult because of the dierence between the two displacement ®elds. However, it can be noticed that the macroscopic strains used in our approach are not equal to their values at the cell center. Since we are interested in the overall behaviour of the continuum, the choice of averaging the macroscopic strain gradients seems to be better, as in the homogenization theory [10].

Macroscopic stresses

The eective stiness of the lattice relates the macroscopic strains (de®ned in the previous section) to the macroscopic stresses. Two means can be used to calculate these macroscopic stresses.

Beam analogy

Let us consider a typical cell of the (single-or multiple-bay) lattice, as an element of a repetitive lattice subjected to external loads. The cell is in equilibrium with concentrated loads applied to the nodes at x= 2 l/2. The loads are assumed to be zero at the other nodes of the boundary y= 2 l/2, as shown in Fig. 3.

For such a beam-like lattice, the continuum model is a one-dimensional beam. Then, if the cell is considered as a beam of length l, the extensional force N and the transverse shear force T are constant along the beam, and the bending moment M is linear. Therefore, one directly obtains the average value of the integrated stresses over the beam length:

hNi X 3 X 4 ÀX 1 X 2 1 l X i x i X i hTi 1 l X i x i Y i hMi 1 2 h 2 X 2 À X 1 h 2 X 3 À X 4 À 1 l X i x i y i X i 9
where the summation is made from 1 to 4.

It should be noticed that, after assuming the continuum to behave like a beam, the macroscopic integrated stresses are averaged over the beam length.

Energy equivalence

Another means of obtaining the macroscopic stresses is to equate the strain energy of the repeated discrete lattice and that of the continuum element. The energy equivalence is classically used in the substitute continuum approach, and also in the homogenization methods. It will enable us to characterize the macroscopic stresses as the dual quantities of the macroscopic strains.

Single-bay lattice

The ®nite element analysis of the lattice cell shown in Fig. 3 gives a system of the form:

K d fu d g fF d g 10
where [K d ] is the stiness matrix, {u d } the vector of nodal displacements and {F d } the vector of equivalent nodal forces (X i ,Y i ) (the subscript `d' refers to the discrete lattice). Eq. ( 5) can be written at each node i = 1,4 as follows: where {u c } T ={u 0 ,v 0 ,r 0 ,{E } T ,{E 1 } T ,{E 2 } T } (the subscript `c' refers to the continuum) with fE g T fhE x i,hg xy i,hÀE x,y ig

E 1 hE y i; E 2 hE y,x i 12
and where the matrix [G i ] is de®ned as:

G i 1 0 Ày i x i y i =2 Àx i y i 0 0 0 1 x i 0 x i =2 0 y i x i y i 13
This matrix diers from the matrix [T i ] given in Ref. [START_REF] Dow | Continuum models of space station structures[END_REF] where quadratic terms in x i and y i are added because E x,y and E y,x appear in the coecients of the polynomials representing the displacement ®eld respectively in v ,xx and u ,yy . Since the origin of the coordinate system is at the cell center, these terms are constant, and it can be derived that transformation of Eqs. ( 11) and ( 13) is equivalent to that proposed in Ref. [START_REF] Dow | Continuum models of space station structures[END_REF].

Then, with [G ] T =[[G 1 ] T ,[G 2 ] T ,[G 3 ] T ,[G 4 ] T ],
we can write: fu d g Gfu c g 14 and substituting Eq. ( 14) in Eq. (10) yields:

G T K d Gfu c g G T fF d g 15
Denoting {F c }=[G ] T {F d }, this vector is found to be:

fF c g T ( X i X i , X i Y i , X i Ày i X i x i Y i , X i x i X i , X i y i X i x i Y i =2 À X i x i y i X i , X i y i Y i , X i x i y i Y i ) 16
Next, the strain energy on the discrete cell expressed in terms of the continuum parameters is given by:

U d 1 2 fu c g T fF c g 17
On the other hand, the strain energy on the continuum (a unit thickness is assumed) is expressed as:

U c 1 2 S c s:E dS c 18
where S c is the continuum cell surface. Using Eq. ( 5), U c can be re-written as:

U c 1 2 fu c g T fF s g 19
where the vector {F s } is found to be:

fF s g T S c f0,0,0,

20

In order to calculate the components of {F s }, given that the continuum is a beam of axis x, we write:

S c s xx dS c l=2 Àl=2 h=2 Àh=2 s xx dy ! dx l=2 Àl=2
Nx dx lhNi 21

In the same way, it can easily be seen that the ®fth component of {F s } is equal to lhTi. In order to show that the sixth component of {F s } is equal to lhMi, let us prove that:

S c xs xy dS c 0 22
Since the center of the cell is at x = 0, it is sucient to show that f h/2 Àh/2 s xy dy is independent of x. To this end, the continuum is considered as a bidimensional domain in plane stress state, subjected to external loads on its boundary x= 2 l/2 and with free edges at y=2l/2. Then, using the equilibrium equation, we have:

@ @ x h=2 Àh=2 s xy dy ! À h=2 Àh=2
s yy,y dy 23 and this last integral is zero because of the boundary conditions.

From the energy equivalence U d =U c , and using Eqs. ( 17) and ( 19), one can derive:

fF c g fF s g 24
As a result, it can be checked that the ®rst three components of {F c } are zero, which corresponds to the cell equilibrium. The next three relationships show that the energy equivalence enables us to de®ned the average integrated stresses as the dual quantities of the average strains. The last two components of {F s } will be ident-i®ed later.

It follows that {F s } can be put in the form:

fF s g T lf0,0,0,fSg

X i Ày i X i x i Y i 0
we have:

X i y i X i x i Y i =2 X i x i Y i 2.3.2.

Double-bay lattice

The modi®cations to the previous development for the case of a double-bay lattice are presented now. From the results of Section 2.2.2, the vector {u c } can be noted: Then, expressing Eq. ( 5) in the global coordinate system (the origin of which is at the cell center), the matrix [G ] de®ned in Eq. ( 14) can be found as follows:

fu c g T ffu à g T ,fE g T ,fE à g T ,fE 1 g T ,
. for the ®rst-bay nodes (the superscript indicates the bay number):

u 1 i v 1 i G 1 i fu c g 29 with .
for the second-bay nodes, similar expressions hold, with the four terms including À(x i +l/4) changed into +(x i Àl/4).

We obtain

[G ] T =[[G 1 1 ] T , [G 1 2 ] T , [G 2 3 ] T , [G 2 4 ] T , [G 1 5 ] T , [G 1 6 ] T ],
where the nodes 1±4 are shown in Fig. 3 and the nodes 5 and 6 are located at the junction of the two bays. Since there is no external load on these two interior nodes, the calculation of {F c }=[G ] T {F d } yields:

fF c g T ( X i X i , X i Y i X i Ày i X i x i Y i , X i x i X i , 1 2
X i y i X i x i Y i , À X i x i y i X i , l 4
X i X i , l 4
X i Y i , X i y i Y i , l 4
X i y i X i , l 4
X i y i Y i , X i x i y i Y i )
where the summation is made from 1 to 4.

The energy equivalence gives the eight relations in Eqs. ( 16), (20), ( 24), (26) and four new relationships which are not reported here. The vector {F s } can then be de®ned as: fF s g T lf0,0,0,fSg T ,fF Ã g T ,fF 1 g T ,fF 2 g T g 32 with fSg T fhNi,hTi,hMig; fF Ã g T f0,0g

lfF 1 g T X i y i Y i , l 4
y i X i , l 4
y i Y i lF 2 X i x i y i Y i 33

Comments on the continuum stresses

In Section 2.3, it is shown that the energy equivalence makes it possible to calculate the macroscopic integrated stresses. It must be noticed that in the literature, the developments are usually focused on the strain and kinetic energies [1±6] and no special attention is paid to the integrated stresses. It will be seen in Section 4 that the knowledge of these eorts allows a better understanding of the assumptions used for the computation of the equivalent continuum parameters.

G 1 i 1 0 Ày i x i y i =2 Àx i y i Àx i l=4 0 0 Àx i l=4y i 0 0 0 1 x i 0 x i =2 0 0
Àx i l=4 y i 0 Àx i l=4y i x i y i 30

Matrix formulation

After suppressing the rows and columns corresponding to the rigid body terms, Eqs. (15), ( 24) and (25) lead to:

2 4 K 1 11 K 1 12 K 1 13 K 1 21 K 1 22 K 1 23 K 1 31 K 1 32 K 1 33 3 5 8 < : fE g E 1 E 2 9 = ; l 8 < : fSg F 1 F 2 9 =
; 34 for a single-bay lattice, and in the same way, the relation:

2 6 6 4
K 2 11 K 2 12 K 2 13 K 2 14 K 2 21 K 2 22 K 2 23 K 2 24 K 2 31 K 2 32 K 2 33 K 2 34 K 2 41 K 2 42 K 2 43 K 2 44 3 7 7 5 8 > > < > > : fE g fE Ã g fE 1 g fE 2 g 9 > > = > > ; l 8 > > < > > : fSg f0g fF 1 g fF 2 g 9 > > = > > ; 35
is obtained for a double-bay lattice.

Computation of beam parameters

The problem is to ®nd [D hom ] such that: fSg D hom fE g 36

Three dierent methods are proposed herein, consisting of eliminating the strain and stress terms in Eqs.

(35) and (36) that are not involved in the Timoshenkobeam theory.

The periodic solution

Due to the periodicity of the repetitive lattice, the strains and internal eorts are locally periodic in the lattice. This also holds for the strains and stresses in the equivalent continuum which is periodic too. Then:

fX 1 ,X 2 ,Y 1 ,Y 2 g ÀfX 3 ,X 4 ,Y 3 ,Y 4 g 37
and E y will take the same value on the opposite sides of the cell. Thus, it turns out that {E 2 }={F 1 }=0. As a result:

D hom Per 1 l K 11 À K 12 K À1 22 K 21 38
where K ij =K 1 ij , i,j = 1,2 for the single-bay lattice, and

K 11 K 2 11 ; K 12 K 2 12 K 2 13 ; K 21 K 2 21 K 2 31 ; K 22 K 2 22 K 2 23 K 2 32 K 2 33 39
for the double-bay lattice. This method of computation of the beam parameters is equivalent to that proposed in Ref. [START_REF] Noor | Continuum models for beam-and platelike lattice structures[END_REF]. In the above mentioned paper, the strain components corresponding to {E 2 } are set to zero so as to ensure the compatibility between the repeated elements, which is of course equivalent to the strain periodicity. On the other hand, {F 1 } is set to zero according to the Timoshenko-beam theory, because this force is associated with the strain component E y .

Therefore, we have presented a rigorous way of obtaining the eective coecients of the beam, and the lattice periodicity is taken into account. However, it is interesting to show that approximate methods can be derived from the matrix formulation of Section 3, and compare them with other methods discussed in the literature.

A force method

Another way of obtaining the overall beam behaviour is to assume that {F 1 }={F 2 }={0}, which gives:

D hom F 1 l K Ã 11 À K Ã 12 K Ã 22 À1 K Ã 21 40
where K Ã ij , i,j = 1,2 is given by a formula similar to Eq. (39) for the single-bay lattice, and for the double bay lattice:

K Ã 11 K 2 11 ; K Ã 12 K 2 12 K 2 13 K 2 14 ; K Ã 21 2 4 K 2 21 K 2 31 K 2 41 3 5 ; K Ã 22 2 4 K 2 22 K 2 23 K 2 24 K 2 32 K 2 33 K 2 34 K 2 42 K 2 43 K 2 44 3 5 41
Such a method amounts to eliminating, through static condensation, all the strain parameters which do not appear in the macroscopic strains {E } of the continuum. This approach has already been applied in Refs. [3±4]. Another means of determining the matrix D hom F is to submit the lattice (with the help of the matrix [G ] de®ned in Subsection 2.3.2.1 and 2.3.2.2) to three dierent loading cases corresponding to a unit load, respectively, in hNi, hTi or hMi with {F 1 }={F 2 }={0}. Then, from the inverse of the matrix [G ] and the vector of the nodal displacements, {E } is calculated and yields the inverse of D hom F , column by column. For the single bay lattice, the three loading cases are shown in Fig. 4.

These loading cases can be found in Ref. [START_REF] Sun | Global±local approach to solving vibration of large truss structures[END_REF]. One of the merits of the formulation used herein is to show that, because of the regularity of [G ], there is only one loading case for a given {S } under the assumption that {F 1 }={F 2 }={0}. For a double-bay lattice and for the loading case correspond to a unit hTi, there are external loads exerted on the nodes at the junction of the two bays.

Lastly, it must also be mentioned that such a method based on the loading cases is cumbersome compared to Eq. ( 40). However, it is very easy to use with a standard ®nite element program and it gives a physical interpretation to the hypothesis {F 1 }={F 2 }={0}, which amounts to assuming the load distribution on the lattice boundary for a given {S }.

A displacement method

The dual standpoint of the previous force method is to assume that {E 1 }={E 2 }={0}, which yields:

D hom D 1 l K 1 11 42
for the single-bay lattice and

D hom D 1 l K 2 11 À K 2 12 K 2 22 À1 K 2 21 43
for the double-bay lattice. Again,the matrix D hom D can be computed after sol-ving a local problem on the lattice cell with displacement boundary conditions corresponding to a unit hE x i, hg xy i or hÀE x,y i with {E 1 }={E 2 }={0}. The vector {S } is then calculated from the reactions (X i ,Y i ) at each node (using Eqs. ( 16) and ( 24)). For a single-bay lattice, the boundary conditions for each loading case are easily derived from the matrix [G ]; they are shown in Fig. 5. For a single-bay lattice, with {E 1 }={E 2 }={0}, we take six parameters to represent the displacement ®eld. Moreover, the left and right sections of the beam remain plane, with E y =0. Then, the kinematical assumptions are equivalent to the approach presented in Ref. [START_REF] Lee | Dynamic continuum modeling of beamlike space structures using ®nite element matrices[END_REF]. However, the method proposed in the paper mentioned above is restricted to a one-bay lattice, and the equivalent beam parameters are not obtained directly, but through the comparison of ®nite element stiness matrices of the lattice and the extended Timoshenko beam models.

For the double-bay lattice, if {EÃ}=0 was assumed, the displacement boundary conditions shown in Fig. 5 would be determined, but this assumption is useless, because {FÃ}=0 holds. We have therefore:

fE Ã g ÀK 2 22 À1 K 2 21 fE g 44
from which boundary conditions for a given {E } can be obtained. It can be checked that a zig-zag pattern can occur for each loading case. The last remarks in Section 4.2 are also valid for the displacement method.

Comments on the previous methods

First, one should bear in mind that in the case of the multiple-bay lattice, the initial system of equations given in Eqs. (34) or ( 35) is dierent from that proposed in Refs. [1±4]. With our approach, not only are the displacements written dierently, but also the strains and integrated stresses which are averaged over the cell (these dierences vanish for a single-bay lattice, see Section 2.2.3).

As mentioned previously, the periodic approach for the computation of the eective beam properties is equivalent to the method initially proposed in Ref. [START_REF] Noor | Continuum models for beam-and platelike lattice structures[END_REF] by Noor and his coworkers. The static condensation was used later by Noor et al. [START_REF] Noor | Anisotropic continuum models for beamlike lattice trusses[END_REF], and Dow [START_REF] Dow | Continuum models of space station structures[END_REF], and is equivalent to the procedure of Sun et al. [START_REF] Sun | Global±local approach to solving vibration of large truss structures[END_REF], and to the force methods presented in this paper. Lastly, there are similarities between the approach of Lee [START_REF] Lee | Dynamic continuum modeling of beamlike space structures using ®nite element matrices[END_REF] and the displacement method given herein, but the latter is more general.

Moreover, it is worth noticing that the characterization of the right-hand member of the system of equations (i.e. the macroscopic integrated stresses) enables us: . to give an estimate of the internal truss forces in terms of the resultant forces acting on the equivalent beam section (from Eqs. (34) or ( 35) and the matrix [G ]); . to obtain a unifying approach for deriving dierent methods for computing the eective beam properties; . to gain insight into the assumptions of each method.

Thus, it can be concluded that the periodic approach is the most accurate way of obtaining the eective stiness of the continuum, because it is the only method which really takes the periodicity of the repetitive lattice into account (in this way, it can be said that {E 2 }={F 1 }=0 are not assumptions). Consequently, the force and displacement methods are approximate, and their drawback is to lead to equivalent beam properties which depend on the cell choice.

Besides, the approaches used in the present work can be compared with homogenization methods for periodic materials. Indeed, the force method, the periodic approach and the displacement method are equivalent to the localization problems found in the homogenization theory, with respectively uniform stress, periodic, or uniform strain boundary conditions [10]. Furthermore, as can be seen in Section 5 hereafter, the strain energies corresponding to each method are classi®ed in the same order as in homogenization.

Comparison between the three methods

In Section 4, three dierent methods for the computation of beam parameters have been proposed. The aim of this section is to examine the possibility of comparing the results given by each method.

Let us ®rst prove the following lemma: Lemma 5.1. let K ij , i,j = 1,2 be submatrices of a positive de®nite symmetric stiness matrix K. Then the following inequality holds:

VfE g,fE g T K 11 À K 12 K À1 22 K 21 fE gRfE g T K 11 fE g 45
Proof. the result is a direct consequence of the symmetry and the de®nite positiveness of the matrix K.

One can write: For the second inequality, the lemma is used after the static condensation of {E 1 } (and {EÃ} for the doublebay lattice). A system of equations is then obtained with {E } and {E 2 } and the submatrix of indices 11 is found to be lD hom Per . As shown in Section 4.4, the periodic approach is the most accurate method for calculating the equivalent continuum properties; it follows from Eq. (47) that the force method underestimates the stinesses, while the displacement method overestimates it.

fE g T K 12 K À1 22 K

Numerical examples

The methods presented previously are now illus- Let us ®rst consider a single-laced single-bay lattice shown in Fig. 6; its dimensions and material properties are given in Table 1. This lattice is analyzed with the three methods presented in Section 4, using a cell with two vertical trusses of area S/2. The results are shown in Table 2.

It can be checked that the force method results coincide with those given in Ref. [START_REF] Sun | Global±local approach to solving vibration of large truss structures[END_REF] (the sign of C 2 is dierent because of the rotation convention). We can notice that the force method and the periodic approach give the same results for this example. Examining the stiness matrix of Eq. (34) shows that K 1 13 =K 1 31 =0 and that the submatrix K 1 33 is diagonal, hence this equality. Besides, the results are in agreement with the inequalities of Eq. ( 47).

If a cell with only one vertical truss of area S is used, the periodic and displacement methods provide the same results as those of Table 2. On the other hand, the force method gives very dierent results from those obtained with the cell with two vertical trusses of area S/2, especially for the second row of the stiness matrix, which is found to be zero, due to the fact that the loading case corresponding to the unit shear force leads to a vertical load applied at a node which is only connected with a horizontal bar.

The latter result emphasizes the importance of the approximations made in the force method to derive the equivalent beam properties. However, from the results of Table 2, it is not possible to judge which is better between the periodic and the displacement approaches. Therefore, a direct analysis of the actual lattice structure with 11 bays is made, with loadings corresponding to a unit extensional force, a transverse shear force and a bending moment, through the boundary conditions shown in Fig. 4. Then for each loading case, the macroscopic strains and then a compliance matrix are derived from the nodal displacements of the central bay of the lattice; and the stiness matrix giving the eective beam properties is ®nally computed. The results show that the non-zero stiness coecients obtained with the periodic approach are within 0.02% of those of the discrete models, which con®rms that this method is the most rigorous means of deriving the equivalent beam properties. Moreover, for a cantilever repeated element loaded with an axial load, a shear load, and a pure moment, the continuum model based on the periodic solution is found to be the most accurate.

6.2. Double-bay lattice

A test problem

The ®rst step consists of showing the reliability and the accuracy of the displacement ®eld proposed herein by studying a double-bay lattice made of two identical one-bay lattices described in Section 6.1.

Then, with the expression of displacements given in Section 2.2.2, the equivalent beam parameters are calculated. It is possible to verify that for the three methods, the results exactly match those of Table 2.

As a second step, the example is taken to illustrate a method in which a Taylor series expansion is used to write the displacement ®eld [START_REF] Dow | Continuum models of space station structures[END_REF] given by:

where all the rigid body displacements, strains and strain gradients are taken at the cell center. As mentioned in Section 4.4, one can easily check that, for a single-bay lattice, the method based on the Taylor series expansion is equivalent to that proposed here. Then, this method also gives the results of Table 2.

For a double-bay lattice, from Eq. ( 49), we can de-®ne a transformation matrix between the discrete nodal displacements and the vector of the rigid body displacements, strains and strain gradients at the cell center. Then, the strain components are partitioned as in Eq. ( 27) with: fu à g T fu 0 ,v 0 ,r 0 g; fE à g T fE x,x ,g xy,x À E x,y g fE g T fE x ,g xy , À E x,y g fE 1 g T fE y ,E x,xy ,E y,xx g; fE 2 g T fE y,x g 50

It is therefore possible to apply the methods proposed in Section 4 to compute the equivalent beam parameters.

The results obtained for the particular double-bay lattice considered herein are shown in Table 3. As can be seen, except for the displacement method, these results deviate from the results of Table 2, with a dierence which reaches 10% on GS and C 1 .

This discrepancy is not very large, but illustrates that the Taylor series expansion method gives results which depend on the way this single-bay lattice is treated (with a method developed for a one-or a doublebay lattice). The error results from the expression used in Eq. (49) for the displacements, where terms in x 2 and x 2 y appear.

The previous drawback (dependence of the results on the cell choice) does not appear with the displacement ®eld proposed herein, based on piecewise linear functions.

A single-laced double-bay lattice

A single-laced double-bay lattice as shown in Fig. 7 is now taken as an example. Following the procedure presented in Section 2.2.2, the equivalent beam properties shown in Table 4 are obtained. It is noticed that the only dierence between the three methods is the equivalent shear stiness, which is classi®ed as in Eq. (47).

We also present the results found with the method based on the displacement given in Eq. ( 49). The results are shown in Table 4, except for the equivalent shear stiness which is: 5.459 Â 10 5 N, 5.701 Â 10 5 N, and 6.745 Â 10 5 N, for the force, periodic and displacement methods, respectively.

In order to evaluate the accuracy of the stiness coecients obtained, the procedure described in Section 6.1 to identify the eective beam parameters was carried out. The results are given in the last column of Table 4, and show that the periodic method gives the most accurate solution.

Conclusions

A simple approach to the representation of repetitive beam-like lattice trusses by an equivalent continuum has been presented. The procedure is based on the energy equivalence between lattice and continuum models, and attempts are made to improve the accuracy of the displacements approximation in methods using polynomials. A piecewise linear displacement ®eld is introduced.

Moreover, some features of the homogenization theory for periodic materials are used, e.g. the averaged strains and integrated stresses are taken as continuum parameters, and the lattice periodicity is used to determine the equivalent beam properties. The general formulation presented in this paper extends methods which have already been given in the literature to calculate the continuum properties. The assumptions made in these methods are pointed out, especially for the widely used method based on static condensation.

The numerical results given illustrate the reliability of the present method. The method has only been applied for the planar static analysis of beam-like lattice trusses, but it also can be used for three-dimensional beam-or plate-like lattice trusses in dynamic analysis. A similar method can also be developed for the study of rigid-joint lattices, with a displacement ®eld written with third-order polynomials.

thesis (in French), Institut National Polytechnique de Grenoble 1994.

[10] Suquet P, Elements of homogenization for inelastic solid mechanics. In: Sanchez-Palencia E, Zaoui A, editors. Homogenization techniques for composite media. Lecture Notes in Physics, 272. Berlin: Springer, 1985. p. 194±278.

Fig. 1 .

 1 Fig. 1. The single-bay lattice.

Fig. 2 .

 2 Fig. 2. The double-bay lattice.

Fig. 3 .

 3 Fig. 3. The cell isolated from the repetitive lattice.

Fig. 4 .

 4 Fig. 4. The force method loading cases (single-bay lattice).

Fig. 5 .

 5 Fig. 5. Boundary conditions of the displacement method (single-bay lattice).

Fig. 6 .

 6 Fig. 6. The single-laced single-bay lattice.

Fig. 7 .

 7 Fig. 7. A single-laced double-bay lattice.

  s xx ,s xy , À ys xx xs xy ,s yy ,xs yy ys xy g dS c

  T ,fF 1 g T ,fF 2 g T g

					25
	with
	fSg T fhNi,hTi,hMig
	lF 1	s yy dS c	X y i Y i
	S c		i
	lF 2	xs yy ys xy dS c	X x i y i Y i	26
	S c			i
	It can be seen that the results are in agreement with
	those of Subsection 2.3.2.1 where the energy equival-
	ence was not used: taking into account that

  fE 2 g T g 27

	with
	fu	à g T fhu 0 i,hv 0 i,hr 0 ig; fE	à g t	E x ,	1 2	g xy	r 0
	fE g

T fhE x i,hg xy i,hÀE x,y ig fE 1 g T fhE y i, E x,y , E y,x g; fE 2 g T fhE y,x ig 28

  21 fE g K 21 fE g T K À1 22 K 21 fE g 46 which proves the lemma, because the second term of the equality is clearly positive.

	Using the lemma, one can show that:
	VfE g,fE g T D hom Per fE gRfE g T D hom D	fE g
	VfE g,fE g T D hom F	fE gRfE g T D hom Per fE g	47

Table 1

 1 Material properties and geometric dimensions of the single-laced single-bay lattice

	Lattice members	E (N/m 2 )	L (m)	S (m 2 )	Designation
	Long bar	71.7 Â 10 9	7.5	8 Â 10 À5	1
				18 Â 10 À5
	Diagonal bar	71.7 Â 10 9	9.0	4 Â 10 À5	Ð
	Batten	71.7 Â 10 9	5.0	6 Â 10 À5	0

---

Table 2

 2

	Equivalent beam properties
	Methods	Force	Periodic	Displacement
	ES (10 7 N)	2.013	2.013	2.029
	GS (10 5 N)	6.592	6.592	7.343
	EI (10 7 Nm 2 )	11.651	11.651	11.651
	C 1 (10 5 N) C 2 (10 7 Nm) C 3 (10 7 Nm)	9.889 1.793 0.000	9.889 1.793 0.000	11.014 1.793 0.000

Table 4

 4 Equivalent beam properties

Methods


F Dupoirieux

A Vincent

N Bertier

A Banh

Numerical simulation of a premixed CH4-air burner for comparison of RANS and LES methodologies

Introduction

The manufacturers of gas turbine are facing tough challenges linked to the reduction of the specific fuel consumption, pollutant emission and noise, and to the improvement of the reliability. The reduction of the specific fuel consumption is tackled by an increase of the Operating Pressure Ratio (OPR) which leads to higher pressure and temperature in the combustion chamber. The counterpart of this temperature increase is that pollutant formation such as NO x is enhanced and that wall and turbine blades are submitted to a higher thermal stress. To cope with this problem, gas turbine manufacturers have to design new combustor concepts. To do so, they rely on the results obtained on realistic experiments but also more and more commonly on results yielded by numerical tools. These tools are based on the resolution of the Navier-Stokes (NS) equations for reactive flows which can be of Reynolds Averaged type (or RANS for Reynolds Averaged Navier-Stokes), or unsteady type such as Large Eddy Simulation (LES). Whatever the resolution type, physical models have to be introduced to close the NS equations. Moreover the spatial and time discretization necessary to solve the NS equations introduces numerical errors, the control of which can be difficult. For those reasons, the numerical tools must be validated by comparison with results of experiments reproducing the real physics encountered in gas turbine combustors. Some of these experiments aim at reproducing only the turbulent mixing between hot and fresh gas [START_REF] Borghi | Turbulent combustion in a premixed flow[END_REF] while others have in addition the purpose to mimic the swirl induced by the injection devices either with gaseous fuel [START_REF] Meier | Detailed characterization of the dynamics of thermoacoustic pulsations in a lean premixed swirl flame[END_REF] [3] or with liquid fuel [START_REF] Dupoirieux | Methodology for the numerical prediction of soot formation in turbulent reactive flows and application to aircraft engine combustors[END_REF]. Some experiments reproduce the process of flame holding by recirculation of hot gases behind an obstacle or a sudden pipe expansion [START_REF] Pitz | Combustion in a turbulent mixing layer formed at a rearward-facing step[END_REF] [START_REF] Sainte-Rose | A DES method applied to a backward facing step reactive flow[END_REF]. Despite its apparent simplicity this type of configuration is a severe validation case for the numerical tools because the hot gas recirculation and the flame developing from this recirculation to the opposite wall are not confined in a closed room, so that the correct simulation of their extension depends on the quality of the physical models and the numerical resolution. For that reason the configuration described hereafter, where a premixed flame is stabilized by recirculation of hot gases behind a rearward facing step, has been considered to assess different resolution methodologies (RANS, LES). Section 2 of this paper gives an overview of the experimental set-up and associated optical diagnostics. In section 3, the computation results concerning the simulation of the rearward facing step burner are compared to experimental results for different numerical methodologies, and key points deduced from these comparisons are summarized in conclusion section 4.

Experimental set-up and diagnostics

The numerical simulations presented in this paper are related to the experiments achieved at ONERA's LAERTE facility in the case of a premixed methane-air flame stabilized by a backward facing step [START_REF] Magre | Further studies by CARS of premixed turbulent combustion in a high velocity flow[END_REF].

The combustor is fed with preheated air. Methane and heated air mass flow are measured by use of choked nozzles. The preheated air passes through the inlet (premixing) duct before entering the combustor. The inlet duct includes several rectangular-section pipes (100 mm × 65 mm cross section). The methane injector is made of five vertical bars pierced with a total of 55 holes. It is placed about 1.2 m upstream from the step. Downstream of the injector is a honeycomb grid (5 mm × 5 mm × 100 mm channels), which straightens the flow, reduces the turbulence level and improves the homogeneity of the airflow. The thickness of the boundary layer at the step section is about 13 mm. Velocity profile corresponds to a fully-developed turbulent boundary layer and is well-described by a 1/7 power law. The step is 35 mm high, and the 1.4 m long chamber has a square section of 100 mm × 100 mm.

The combustor walls are water cooled allowing long duration tests and large optical accesses (100 mm × 260 mm quartz windows) are used to achieve optical diagnostics. A flush-mounted spark igniter is inserted through the lower combustor wall.

An exit nozzle equipped with a throttling plug is connected downstream of the combustor to control the chamber pressure (moving this plug also modifies the acoustic impedance).

TEMPERATURE MEASUREMENTS BY CARS (A3C CONFIGURATION)

The temperature of the flow has been measured at different positions of the combustor by use of Coherent Anti-stokes Raman Scattering technics. These measurements have been obtained only for the reacting case. Because of the low dynamics of the detector, most of the acquisitions have been done twice, once for "cold" values and once for "hot" values.

Most of the histograms presents a bimodal distribution, one centered on T = 550 K (fresh gases), the other around 2000-2200 K (burnt gases). As a consequence, the average temperature has a low probability to be effectively measured. The lack of intermediate values can be explained by the very thin thickness of the flame front.

The conditions of temperature measurement by CARS spectroscopy are the following:

• time-resolution of 12 ns for a single-shot acquisition (i.e. laser pulse duration),

• spatial resolution of 3 mm × 40 µm (Boxcars arrangement of the laser beams),

• measurement accuracy of 3-5 % for the temperature.

Average test conditions were:

• ṁ air = 0.194 ± 0.012 kg.s -1

• P = 0.103 ± 0.003 MPa • T air inlet = 532 ± 19 K,

• methane fuel* • E.R. = 0.833 ± 0.022.

VELOCITY MEASUREMENTS BY LDV (A3C CONFIGURATION)

The velocity has been measured in the combustor by use of Laser Doppler Velocimetry (LDV) technics. For each measurement location, a selection of 650 to 5000 single-shot measurements has been used to build the current database (improper measurements have been removed). Sub-micron particles used for seeding were made of zirconium dioxide.

For the non-reacting case, test conditions were:

• ṁ air = 0.196 ± 0.009 kg.s -1

• P = 0.100 ± 0.002 MPa • T air inlet = 529 ± 28 K.

For the reacting case, test conditions were:

• ṁ air = 0.193 ± 0.007 kg.s -1

• P = 0.103 ± 0.002 MPa • T air inlet = 521 ± 18 K • methane fuel • E.R. = 0.855 ± 0.036.

Numerical simulation of the rearward facing step burner

Conditions of the calculations

The rearward facing step burner described in section 2 is simulated with the help of the solver CHARME of CEDRE software [START_REF] Refloch | CEDRE Software[END_REF]. The resolution of the 3D Navier-Stokes equations is achieved over a computational domain including a small part of the duct upstream of the step and a much longer part downstream of it. Since the turbulent boundary layer upstream of the step is fully developed and has been experimentally characterized, the influence of the upstream turbulence level on the downstream mixing is introduced through the inlet boundary conditions in terms of velocity profile. Therefore the calculation of the evolution of the upstream boundary layer along the lower wall is not required, so that the length of duct upstream of the step considered in the calculation is reduced to 200 mm. On the other hand, the length of duct taken into account downstream of the step is much longer, because it must encompass all the flame which is attached to the step corner and slowly propagates to the upper wall. This length is equal to 1300 mm.

For the RANS simulation the transverse direction is not considered and the calculation grid is 2D. For the LES the calculation grid is of course 3D but, to offset the great length of the calculation domain and to keep the computation time reasonable, only a part of the duct is considered in the transverse horizontal direction: The width of the calculation domain is only 70 mm while the real one is 100 mm, and periodic boundary conditions are applied on the lateral walls instead of usual friction conditions, in order to avoid the mesh refinement necessary to deal with the lateral boundary layers. Knowing that the measurements have been achieved in the lengthwise symmetry plane, this simplification does not alter the comparison between numerical and experimental results as long as the lateral boundary layers are not too thick.

Even if the CHARME solver may deal with any type of mesh, it has been decided to generate a calculation grid with rectangular (in 2D) or parallelepiped (in 3D) meshes to optimize the resolution accuracy. The RANS grid includes 501000 rectangular cells. The LES grid is made up of 3 domains: -a non refined domain in the most upstream part of the step including 30000 cells, -a refined domain slightly upstream and downstream of the step with 2500000 cells, -a third domain located far downstream of the step, with a progressive loss of refinement, including 236000 cells. The grid is refined near the walls and a law of the wall is used both in 2D and 3D to correctly reproduce the evolution of the wall boundary layers. Figure 2 highlights the features of the LES calculation grid. The spatial discretization relies on upwind uncentered convective fluxes with a MUSCL type limiter of order 2. The time integration is semi-implicit of order 2 as well. The boundary conditions are: -Inlet: Bulk velocity of 50 m/s for a CH 4 -air mixture at equivalence ratio 0.855, synthetic turbulence rate of 5 %, temperature of 520 K; -Outlet: Pressure of 101300 Pa combined with non reflexive conditions for acoustic waves;

-Lower wall: adiabatic conditions before the step, temperature of 1000 K downstream of the step. -Upper wall: adiabatic conditions A two-equations turbulence model is used in RANS simulations while the Smagorinsky subgrid scale model is used for the LES to account for the non resolved turbulent structures. Concerning combustion, an EBU (Eddy Break-Up) model with a limitation of the reaction rate to account for the chemical equilibrium temperature [START_REF] Dupoirieux | The Models of Turbulent Combustion in the CHARME Solver of CEDRE[END_REF] is used in RANS simulations. For the LES, the combustion model is based on the Thickened Flame Model [START_REF] Colin | A Thickened Flame Model for Large Eddy Simulations of Turbulent Premixed Combustion[END_REF] including a two-steps mechanism [START_REF] Franzelli | Large Eddy Simulation of combustion instabilities in a lean partially premixed swirl flame[END_REF] aiming at reproducing the correct laminar propagation speed of premixed CH 4 -air flames.

The time step used for the LES is 10 -6 s, which makes 5 to 10 times the CFL (Courant -Friedrichs -Lewy number) for a characteristic cell size of 10 -4 m. This does not alter the accuracy of the simulation since, in our configuration, the time evolution is driven by the Kelvin-Helmoltz instabilities and the flame pulsation induced by the large turbulent structures, but not by the acoustics. 2×10 5 time steps are required to obtain a statistically converged flowfield, which corresponds to a CPU time of 20 hours on 128 Nehalem processor cores of the ONERA super-computer. Until x=60 mm, the LES profiles match very well with the measurements. From x=80 mm to x=250 mm, the agreement between numerical and experimental profiles remains acceptable but, from x=340 mm to x=910 mm, some discrepancies must be noted. Two reasons can explain the deterioration of the numerical results with the distance from the step. Firstly the calculation grid is less and less refined with increasing x, which worsens the accuracy of the calculation; secondly, the lateral boundary layers, which are not considered in the calculation, tend to thicken, which can accelerate the flow in the lengthwise symmetry plane of the burner, where the velocities are measured. This last point could explain the higher measured velocities obtained for x=340 mm and x=460 mm. We note that the RANS results are much less satisfactory than the LES ones: the RANS simulation greatly overestimates the reverse velocities in the recirculation zone taking place behind the step and, in x=210 mm and x=250 mm, underestimates the velocity in the part of the flow located above the step. More downstream RANS and LES deteriorate in a similar way. x=910 mm, the matching of which being pure coincidence, LES results are much more satisfactory than RANS results.

Non reactive flow

Reactive flow

Figure 5 gives the average temperature field obtained by LES and an example of instantaneous field yielded by LES. The angle of the average turbulent flame seems quite small despite the existence of large turbulent structures evidenced by the instantaneous field. However the corresponding propagation speed of the turbulent flame is of the order of 5 m/s, i. e. 10 to 15 times the propagation speed of the laminar flame. We also note that the temperature inside the recirculation zone immediately downstream of the step is lower than the flame temperature. This lower temperature is induced by the thermal boundary condition (1000 K on the wall) which has a great impact on recirculating gases associated with a high residence time. Let us notice that we obtain a good agreement between LES and experiment near the lower wall, where we noted some discrepancies concerning the velocity fluctuations. This is not surprising: the temperature in the lower part of the combustor (except in the very thin thermal boundary layer close to the wall) must be equal to the chemical equilibrium temperature, whatever the value of the velocity, and this temperature is recovered both in the calculations and the experiment. Figure 9 gives the profiles of temperature fluctuations for the experiment and the LES (this information is not available in the RANS calculations). We observe a very satisfactory agreement until x=340 mm. Further downstream the fluctuations given by the LES naturally vanish because the lower part of the burner is filled with burnt gases for which no temperature fluctuation can be produced by the calculation. The fluctuations observed in the lower part of the last experimental profiles can be explained either by intermittent local extinctions or by huge turbulent fluctuations introducing intermittently some fresh gases in the lower part of the combustor. These phenomena can not be captured without a sophisticated combustion model and/or a very refined calculation grid.

Conclusion

The calculation of the rearward facing step burner is very challenging because the location of the flame is not imposed by the geometry but by the physical models, which induce the turbulent flame propagation, and the accuracy of the numerical resolution. We compare in this paper RANS and LES simulations of this burner. As expected more accurate results are obtained with LES as long as the calculation grid remains refined. More sophisticated chemical mechanisms could help improve the prediction of the turbulent flame propagation to the upper wall. We also note that an effort must be paid to better describe the behavior of the boundary layer around the reattachment point, in order to correctly predict the velocity fluctuations in this region.

In conjunction with a grid refinement in the downstream part of the burner, the calculation of the lateral boundary layers is another way of improving the quality of the simulation.

Figure 1

 1 Figure 1 presents a schematic of the facility.

Figure 1 :

 1 Figure 1: Schematic of the experimental set-up

Figure 2 :

 2 Figure 2: Calculation grid for the LES of the backward facing step

Figure 3 :

 3 Figure 3: Non reactive case: profiles of average axial velocity; x: distance from the step; red dots: measurements; blue lines: LES; purple lines: RANS

Figure 4 :

 4 Figure 4: Non reactive case: profiles of axial velocity fluctuation; x: distance from the step; red dots: measurements; blue lines: LES; purple lines: RANS Figure 4 gives the profiles of axial velocity fluctuations. Again, the agreement between LES and experiment is very satisfactory until x=60 mm. It remains acceptable until x=120 mm but worsens more downstream, probably because of the less refined spatial discretization. The RANS values are deduced from the turbulent momentum fluxes given by the 2 equations turbulence model. For all profiles except

Figure 5 :Figure 6 :

 56 Figure 5: Average temperature field obtained by LES (top); instantaneous temperature field given by LES (bottom)

Figure 7 :

 7 Figure 7: Reactive case: profiles of axial velocity fluctuation; x: distance from the step; red dots: measurements; blue lines: LES; purple lines: RANS Figure 7 gives the profiles of axial velocity fluctuations. One notes some discrepancies between LES and experiment, especially near the lower wall. It seems that the law of the wall used in the calculation is not damping enough the large turbulent structures interacting with the lower wall, especially around the reattachment point (x=120 mm) where the LES gives a large peak of fluctuation not observed in the experiment. The RANS calculation behaves better, probably because it has not to deal with this interaction, since the effect of the large structures is included in the turbulent viscosity.

Figure 8 :

 8 Figure 8: Profiles of average temperature; x: distance from the step; red dots: measurements; blue lines: LES; purple lines: RANS

Figure 9 :

 9 Figure 9: Profiles of temperature fluctuations; x: distance from the step; red dots: measurements; blue lines: LES; purple lines: RANS

Corresponding author : axel.vincent@onera.fr


C Olivier

G Poignand
email: gaelle.poignand@univ-lemans.fr

G Penelet

P Lotton

SIMPLIFIED MODELING OF A THERMO-ACOUSTO-ELECTRIC ENGINE FORCED BY AN EXTERNAL SOUND SOURCE

Keywords: Thermoacoustics, Nonlinear dynamics, Active Control

Introduction

Thermoacoustic wave-generators are usually designed and optimized using software tools [START_REF] Ward | Design Environement for low Amplitude ThermoAcoustic Engines[END_REF][START_REF] Swift | Thermoacoustics: A unifying perspective for some engines and refrigerators[END_REF] based on the linear thermoacoustic theory. These tools enable to predict the operating point of a thermoacoustic engine from the balance between the thermoacoustic amplification process and the numerous nonlinear effects saturating wave amplitude growth, the latter effects being very difficult to describe properly. A non exhaustive list of these effects includes (apart from the thermoacoustic heat pumping accompanying wave amplification) different kinds of streaming, nonlinear acoustic propagation or aerodynamical and thermal entrance effects. These nonlinear effects are responsible for both dissipation of acoustic power and perturbations of the temperature field in the thermoacoustic core, that work together to limit the overall performances of thermal-to-acoustic conversion. Therefore, the development of adequate simulation tools is still needed to describe the high level of complexity of the processes involved above the onset of thermoacoustic instability. Direct numerical simulation seems to be the only way to reproduce quantitatively the effects mentioned above, but it is still limited by large computation times inherent to the complicated physics and multiple time and space scales involved in the description of thermoacoustic engines. Simplified analytical models can help in getting a deeper physical insight about the processes involved, but they are also based on substantial approximations. In this study, we use a simplified modeling of thermoacoustic engines to help understanding recent experimental observations dealing with the external forcing of thermoacoustic oscillations [START_REF] Desjouy | Active control of thermoacoustic amplification in an annular thermoacoustic engine[END_REF][START_REF] Olivier | Active control of thermoacoustic amplification in a thermo-acousto-electric engine[END_REF][START_REF] Poignand | Â ńÂȃHysteretic behavior induced by an electroacoustic feedback loop in a thermo-acousto-electric generatorÂȃÂ ż[END_REF].

Aims and scope of the study

Recent proof of concept studies [START_REF] Desjouy | Active control of thermoacoustic amplification in an annular thermoacoustic engine[END_REF][START_REF] Olivier | Active control of thermoacoustic amplification in a thermo-acousto-electric engine[END_REF] performed at LAUM have already shown that the external forcing of thermoacoustic oscillations could be an interesting way to increase the efficiency of a thermoacoustic prime-mover. In particular, we have shown that the adjunction of an adequately tuned electroacoustic feedback loop to a thermo-acousto-electric transducer can lead to the increase of its global efficiency [START_REF] Olivier | Active control of thermoacoustic amplification in a thermo-acousto-electric engine[END_REF]. Using a reference pressure signal from the engine as an input, an auxiliary acoustical source is used to control self-sustained thermoacoustic oscillations in the engine. By tuning the properties of this feedback (in terms of amplification gain and phase-shift relative to the reference signal), the efficiency of the transducer can either be reduced down to the death of the auto-oscillations, or increased, thus generating a greater output power for the same heat input. More recently, we have shown on the same device that the electroacoustic feedback loop is also responsible for the occurrence of hysteretic behavior [START_REF] Poignand | Â ńÂȃHysteretic behavior induced by an electroacoustic feedback loop in a thermo-acousto-electric generatorÂȃÂ ż[END_REF] of both the temperature difference across the regenerator and the efficiency versus the heat input. However, at the present time the multiplicity of the involved nonlinear effects and their complexity does not allow to understand why and how the mechanisms at stake in this process of "active control" enable to increase the global efficiency of the thermoacoustic engine. It is there-fore worth considering to make theoretical studies on the topic. This task is however arduous because a pre-requisite is to understand deeply and to describe properly the self-sustained oscillator itself before describing its nonlinear coupling with external sound sources. However, we plan as a first step to find a qualitative description of the dynamics of the system which would be well-suited for the analysis of its bifurcation diagram and its response to external forcing. Such kind of approach is very common in nonlinear dynamics in which complicated systems are sometimes summarized by a few ordinary differential equations.

In this study, a simplified model of a thermo-acousto-electric generator is proposed. This reduced-order model is based on a lumped element approach for the description of the acoustic propagation in the engine, coupled with a resolution of unsteady heat transfer through the thermoacoustic core based on a finite difference scheme. This description allows the study of both transient and steady-state operation of the engine. Applying an external forcing to that model also allows to get a better overview of the effects of the feedback loop on the engine. Notably, by gradually increasing its degree of complexity, it enables to study various parameters, such as the coupling of the auxiliary electrodynamic source with the engine, and the influence of the feedback on both the acoustic field and on the distribution of the temperature field within the thermoacoustic core. This simplified model of a forced thermoacoustic transducer will help to get a better understanding of the phenomena responsible for the complex behaviors induced by the electroacoustic feedback loop on the actual prototype: death of oscillations, efficiency improvement accompanied by a temperature gradient reduction, as well as a hysteretic behavior induced by the electroacoustic feedback loop during the processes of gradual heat input increase followed by a gradual decrease.


Philippe Bernardoff
email: philippe.bernardoff@univ-pau.fr

Laplace copulas of multifactor gamma distributions are new generalized Farlie-Gumbel-Morgenstern copulas

Keywords: Archimedean copula, cumulative distribution function, copula, exponential families, infinitely divisible distribution, generalized Farlie-Gumbel-Morgenstern copulas, generalized hypergeometric function, generalized Lauricella functions, Horn function, Kendall's tau, Laplace copula, Laplace transform, multi-factor gamma distribution, multivariate gamma distribution, Spearmann's rho

This paper provides bifactor gamma distribution, trivariate gamma distribution and two copula families on [0, 1] n obtained from the Laplace transforms of the multivariate gamma distribution and the multi-factor gamma distribution given by [P (θ)] -λ and [P (θ)] -λ n i=1 (1 + piθi) -(λ i -λ) respectively, where P is an affine polynomial with respect to the n variables θ1, . . . , θn.

These copulas are new generalized Farlie-Gumbel-Morgenstern copulas and allow in particular to obtain multivariate gamma distributions for which the cumulative distribution functions and the probability distribution functions are known.

Introduction

This paper is motivated by the opportunity to produce explicit cumulative distribution functions (cdf) and probability distribution functions (pdf) of multivariate gamma distributions and multi-factor gamma distributions. In this way, we present two new classes of multivariate copulas generalizing the Farlie-Gumbel-Morgenstern Copulas. From Sklar [START_REF] Sklar | Fonctions de répartition à n dimensions et leurs marges[END_REF] who states that the cdf F of a random vector X = (X 1 , . . . , X n ) with continuous marginal cdfs can be uniquely written in the form

F (x 1 , . . . , x n ) = C [F 1 (x 1 ) , . . . , F n (x n )] , x = (x 1 , . . . , x n ) ∈ R n , ( 1
)
where n is the dimension of the random vector X, where C : [0, 1] n → [0, 1] is a copula and where F 1 , . . . , F n are the marginal cdfs of X.

If f is the pdf of the random vector X and f 1 , . . . , f n the marginal pdfs of X, then

c (u 1 , . . . , u n ) = ∂ n ∂u 1 . . . ∂u n C (u 1 , . . . , u n ) ,
and we have the following equality

f (x 1 , . . . , x n ) = c [F 1 (x 1 ) , . . . , F n (x n )] f 1 (x 1 ) . . . f n (x n ) . (2)
From Equality (1) it is possible to express the cdf F of X. This expression cannot simply explicitly give the copula C for the multivariate gamma distribution associated with (P, λ) and the multi-factor gamma distribution associated with (P, Λ) . But with Joe [START_REF] Joe | Multivariate Models and Dependance Concepts[END_REF], we can give copulas deduced of the Laplace transform of the multivariate gamma distribution associated with (P, λ) and the Laplace transform of the multi-factor gamma distribution associated with (P, Λ) . By applying the formulas [START_REF] Balakrishnan | Continuous Multivariate Distributions[END_REF] and [START_REF] Bar-Lev | The diagonal multivariate natural exponential families and their classification[END_REF] we can give a new multivariate gamma distribution associated with (P, λ) and a new multi-factor gamma distribution associated with (P, Λ) for which we have an explicit formula for its cdf and its pdf.

The paper is organized as follows. Section 2 gives definitions of multivariate gamma distributions and multi-factor gamma distributions for which Laplace transform is given, and considers the bidimentional and tridimentional case. Section 3 defines the Laplace copula. Section 4 states the two main results.

Section 5 applies the main cases to bidimensional and tridimensional cases. In particular, the Joe's family BB10 is generalized. For ease of the fluent exposition of the paper, proofs are collected in the Appendix.

2 Multivariate gamma distributions and multi-factor gamma distributions

In the literature, the multivariate gamma distributions on R n have several non-equivalent definitions.

Many authors require only that the marginal distributions are ordinary gamma distributions (Balakrishnan et al., 1997). In the present paper, we use the extension of the classical one-dimensional definition to R n obtained as follows: we consider an affine polynomial P (θ) in θ = (θ 1 , . . . , θ n ) where 'affine' means that, for j = 1, . . . , n, ∂ 2 P/∂θ 2 j = 0. We also assume that P (0) = 1. For instance, for n = 2, we have

P (θ 1 , θ 2 ) = 1 + p 1 θ 1 + p 2 θ 2 + p {1,2} θ 1 θ 2 .
We denote by P n = P ([n]) the family of all subsets of [n] and P * n the family of non-empty subsets of [n] = {1, . . . , n} . For simplicity, if n is fixed and if there is no ambiguity, we denote these families by P and P * , respectively.

We denote by N the set of non-negative integers.

If z = (z 1 , . . . , z n ) ∈ R n and α = (α 1 , . . . , α n ) ∈ N n , then α! = α 1 ! . . . α n !, |α| = α 1 + . . . + α n , a α = a α1,...,αn and z α = n i=1 z αi i = z α1 1 . . . z αn n .
For T in P n , we simplify the above notation by writing z T = t∈T z t instead of z 1T where

1 T = (α 1 , . . . , α n ) with α i = 1 if i ∈ T and α i = 0 if i / ∈ T.
We also write z -T for t∈T 1/z t . For a mapping a : P → R, we shall use the notation a : P → R, T → a T . In this notation, an affine polynomial with constant term equal to 1 is P (θ) = T ∈P p T θ T , with p ∅ = 1. For simplicity, if T = {t 1 , . . . , t k } , we denote a {t1,...,t k } = a t1...t k . The indicator function of a set S is denoted by ½ S , that is, ½ S (x) = 1 for x ∈ S and 0 for x / ∈ S.

We fix λ > 0. If a random vector X = (X 1 , . . . , X n ) on R n with probability distribution (pd) µ X is such that its Laplace transform is

E {exp [-(θ 1 X 1 + • • • + θ n X n )]} = [P (θ)] -λ , (3)
where E denotes the expectation, for a set of θ with non-empty interior, then we denote µ X = γ (P,λ) , and γ (P,λ) will be called the multivariate gamma distribution associated with (P, λ). These multivariate gamma distributions occur naturally in the classification of natural exponential families in R n (Bar-Lev

et al., 1994).
The marginal distributions of the multivariate gamma distribution associated with (P, λ) are ordinary gamma distributions of parameters (p i , λ) for i = 1, . . . , n, with Laplace transform (1 + p i θ i ) -λ , and pd

γ (pi,λ) (d x) = x λ-1 p -λ i /Γ (λ) exp (-x/p i ) ½ (0,∞) (x)d x.
We extend the first definition to the multi-factor gamma distribution associated with (P, Λ) where Λ = (λ, λ 1 , . . . , λ n ) and λ i λ > 0 for all i = 1, . . . , n by its Laplace transform

E {exp [-(θ 1 X 1 + • • • + θ n X n )]} = [P (θ)] -λ n i=1 (1 + p i θ i ) -(λi-λ) . (4)
The marginal distributions of the multi-factor gamma distribution associated with (P, Λ) are ordinary gamma distributions of parameters (p i , λ i ) for i = 1, . . . , n.

We state first a proposition, whose proof is obvious.

Proposition 1 A random vector X with distribution γ P,Λ can be obtained in the following way:

Let Y be a vector with distribution γ P,λ . Let Z = (Z 1 , . . . , Z n ) be a random vector constituted of independent margins for which its pds are γ (pi,λi) , and such that Z and Y are independent. Then the vector X = Y + Z has Laplace transform (4), and consequently is a multi-factor gamma distribution associated to (P, Λ) .

For the bidimensional case, Dussauchoy and Berland, [START_REF] Dussauchoy | Lois gamma à deux dimensions[END_REF]) consider the random vector X = (X 1 , X 2 ) with Laplace transform

L X (θ 1 , θ 2 ) = (1 + p 1 θ 1 ) -λ1 (1 + p 2 θ 2 ) -λ2 1 - r 12 p 1 p 2 θ 1 (1 + p 1 θ 1 ) (1 + p 2 θ 2 ) -λ = (1 + p 1 θ 1 ) -(λ1-λ) (1 + p 2 θ 2 ) -(λ2-λ) (1 + p 1 θ 1 + p 2 θ 2 + p 12 θ 1 ) -λ , (5)
where r 12 = 1 -p 12 / (p 1 p 2 ) > 0 and p 1 , p 2 , p 12 > 0.

Unfortunately the Laplace transform of these pds is simple, but its pdfs and cdfs are unknown, except for the case n = 2 for the multivariate gamma distribution associated with (P, λ). Let F p m be the generalized hypergeometric function [START_REF] Slater | Generalized Hypergeometric Functions[END_REF] defined by

F p m (α 1 , . . . , α p ; β 1 , . . . , β m ; z) = ∞ k=0 (α 1 ) k • • • (α p ) k (β 1 ) k • • • (β m ) k z k k! , (6)
where (a) k = Γ (a + k) /Γ (a) for a > 0 and k ∈ N, or more generally by (a) 0 = 1, (a) n+1 = (a + n)(a), ∀n ∈ N, ∀a ∈ R, is the Pochhammer's symbol. For simplification, we denote F 0 m by F m . Bernardoff [START_REF] Bernardoff | Which multivariate gamma distributions are infinitely divisible?[END_REF] gives the Proposition 2 Let P (θ 1 , θ 2 ) = 1 + p 1 θ 1 + p 2 θ 2 + p 12 θ 1 θ 2 be an affine polynomial where p 1 , p 2 > 0 and p 1,2 > 0. Let µ = γ (P,λ) be the gamma distribution associated to (P, λ) . The measure µ exists if and only if c = (p 1 p 2 -p 12 ) /p 2 12 = p 1 p 2 /p 2 12 r 12 > 0. Then we have

γ (P,λ) (dx 1 , dx 2 ) = p -λ 12 Γ (λ) 2 e -p 2 p 12 x1- p 1 p 12 x2 (x 1 x 2 ) λ-1 F 1 (λ; cx 1 x 2 ) ½ (0,∞) 2 (x) d (x) . (7)
For the case Λ = (λ, λ, λ 2 ), the multi-factor gamma distribution associated with (P, Λ) is named by Chatelain et al. [START_REF] Chatelain | Change detection in multisensor SAR images using bivariate gamma distributions[END_REF] the multisensor gamma distribution associated with (P, λ, λ 2 )

and they have proved that its pd is given by the equality

γ (P,Λ) (dx 1 , dx 2 ) = p -λ 12 p -(λ2-λ) 2 1 Γ (λ) Γ (λ 2 ) x λ-1 1 x λ2-1 2 e - p 2 p 12 x1- p 1 p 12 x2 Φ 3 (λ 2 -λ; λ 2 ; c p 12 p 2 x 2 ; cx 1 x 2 ) × ½ (0,∞) 2 (x 1 , x 2 ) dx 1 dx 2 , (8)
where

Φ 3 (a; b; x, y) = m,n 0 (a) m (b) m+n x m m! y n n! (9)
is the Horn function.

For the bidimensional general case, we have the following Theorem. Let F I be the function defined by

F I (a, b, c, z) = ∞ m1,m2,m3=0 (a) m1 (b) m2 (c) m3 (a + c) m1+m3 (b + c) m2+m3 z m1 1 m 1 ! z m2 2 m 2 ! z m3 3 m 3 ! ; ( 10
)
it is a particular generalized Lauricella function defined, by example, in [START_REF] Panda | Some integrals associated with the generalized Lauricella functions[END_REF].

Theorem 3 The pd of γ (P,(λ,λ1,λ2)) is given by the equality

γ (P,(λ,λ1,λ2)) (dx 1 , dx 2 ) = p -λ 12 p -(λ1-λ) 1 p -(λ2-λ) 2 Γ (λ 1 ) Γ (λ 2 ) x λ1-1 1 x λ2-1 2 e - p 2 p 12 x1+ p 1 p 12 x2 × F I λ 1 -λ, λ 2 -λ, λ, p 12 p 1 x 1 , p 12 p 2 x 2 , cx 1 x 2 ½ (0,∞) 2 (x 1 , x 2 ) dx 1 dx 2 , (11)
If we get λ 1 = λ in the equality [START_REF] Letac | Laplace transforms which are negative powers of quadratic polynomials[END_REF], we obtain Chatelain and Tourneret's result [START_REF] Joe | Multivariate Models and Dependance Concepts[END_REF] because

F I (0, λ 2 -λ, λ, z 1 , z 2 , z 3 ) = ∞ m2,m3=0 (b) m2 (b + c) m2+m3 z m2 2 m 2 ! z m3 3 m 3 ! = Φ 3 (b; b + c; z 2 , z 3 ) .
Bernardoff [START_REF] Bernardoff | Which multivariate gamma distributions are infinitely divisible?[END_REF] gives the following Proposition:

Proposition 4 Let µ be a multivariate gamma distribution on R n associated to (P, λ) . Assume that µ is not concentrated on a linear subspace of R n of the form {x ∈ R n ; x k = 0} for some k in [n] = {1, . . . , n}.

Then:

(i) For all i ∈ [n] , p i = 0. (ii) If p 1 , . . . , p k < 0 and p k+1 , . . . , p n > 0, then Supp (µ) ⊂ (-∞, 0] k × [0, ∞) n-k . (iii) If p 1 , . . . , p n > 0 then p [n] 0.
Bernardoff [START_REF] Bernardoff | Which multivariate gamma distributions are infinitely divisible?[END_REF] gives a necessary and sufficient condition for infinite divisibility of the multivariate gamma distribution associated with (P, λ), in the sense that the Laplace transform of γ (P,λ)

power t for all positive t is still the Laplace transform of a positive measure, by the following theorem:

Theorem 5 Let µ = γ P,λ be a gamma distribution associated with (P, λ) , where λ > 0 and P (θ) =

T ∈Pn p T θ T is such that p i > 0, for all i ∈ [n] , and p [n] > 0. Let P (θ) = T ∈Pn p T θ T be the affine polynomial such that p T = -p T /p [n] for all T ∈ P n , where

T = [n] T . Let b S = b S ( P ) = |S| k=1 (k -1)! T ∈Π k S T ∈T p T ,
where |S| denotes the cardinality of the set S, and Π k S denotes the set of all partitions of S into k non-empty subsets of S. Then the measure µ is infinitely divisible if and only if

p i < 0 for all i ∈ [n] , (12)
and

b S 0 for all S ∈ P * n such that |S| 2. ( 13
)
Corollary 6 By the properties of infinite divisible distributions we conclude that the necessary and sufficient conditions for infinite divisibility of a multivariate gamma distribution associated to (P, λ) of theorem [START_REF] Chatelain | Change detection in multisensor SAR images using bivariate gamma distributions[END_REF], are also necessary and sufficient conditions for infinite divisibility of multivariate multifactor gamma distribution associated to (P, Λ) .

To illustrate the difficulty to calculate the multivariate gamma distribution associated with (P, λ) we give, for the tridimensional case, the following theorem. Let F II be the function defined by

F II (λ 1 , λ 2 , z 1 , z 2 , z 3 , z 4 ) = ∞ m1,...,m4=0 1 (λ 1 ) m1+m2+m3 (λ 2 ) 2m1+m2+m4 z m1 1 m 1 ! z m2 2 m 2 ! z m3 3 m 3 ! z m4 4 m 4 ! ; (14)
it is still a particular generalized Lauricella function.

Theorem 7 In the case n = 3,

p i > 0 for i ∈ [3] , p ij > 0 for (i, j) ∈ [3] 2 , b ij = -b k p123 + p jk p ik p 2 [Γ (λ)] 3 exp( p 1 x 1 + p 2 x 2 + p 3 x 3 ) (x 1 x 2 x 3 ) λ-1 × F II (λ, λ, b 13 x 1 x 3 b 23 x 2 x 3 , b 123 x 1 x 2 x 3 , b 12 x 1 x 2 , b 13 x 1 x 3 + b 23 x 2 x 3 )½ (0,∞) 3 (x) dx. ( 15
)
Remark 8 The case p 123 = 0 is solved by Letac and Weso lowski (Letac and Weso lowski, 2008).

γ (λ,P ) (dx) = p -λ 123 [Γ (λ)] 3 exp( p 1 x 1 + p 2 x 2 + p 3 x 3 ) (x 1 x 2 x 3 ) λ-1 F 2 λ, λ; b 123 x 1 x 2 x 3 ½ (0,∞) 3 (x) dx,
and if we put λ = 1 in this last equality, we obtain the Kibble and Moran distribution given in (Balakrishnan et al., 2000).

Laplace copula

We recall the following theorem [START_REF] Marshall | Families of Multivariate Distributions[END_REF]

(i = 1, . . . , n). Let K be a copula. If F i (x) = exp -φ -1 i [H i (x)] (i = 1, . . . , n) , then H (x 1 , . . . , x n ) = • • • K [F 1 (x 1 )] θ1 , . . . , [F n (x n )] θn dG (θ 1 , . . . , θ n ) (16)
is an n-variate cdf with marginals H 1 , . . . , H n .

Corollary 11 Under the conditions of the Theorem 10, with

K (x 1 , . . . , x n ) = n i=1 x i , (0 x i 1; i = 1, . . . , n) , H (x 1 , . . . , x n ) = φ φ -1 1 [H 1 (x 1 )] , . . . , φ -1 n [H n (x n )] ( 17
)
defines an n-variate cdf with marginals H 1 , . . . , H n .

Therefore copula for H is associated to G and is given by the formula

C (u 1 , . . . , u n ) = φ φ -1 1 (u 1 ) , . . . φ -1 n (u n ) . (18)
Remark 12 If we know φ and φ i , for i = 1, . . . , n, we can find the copula C, that is the case for the multivariate gamma distributions and for the multivariate multi-factor gamma distributions. For the case H i = G i , we can write :

H (x 1 , . . . , x n ) = φ φ -1 1 [G 1 (x 1 )] , . . . φ -1 n [G n (x n )] , (19)
and this last formula defines n-variate cdf gamma H with margin distributions G 1 , . . . , G n .

Then the copula associated to H is

C (u 1 , . . . , u n ) = φ φ -1 1 (u 1 ) , . . . , φ -1 n (u n ) .
Now, from the Formula 18, we can give the following Definition

Definition 13 Let X = (X 1 , X 2 , . . . , X n ) be a random vector in ]0, +∞[ n with Laplace transform ϕ X defined by L X (θ) = L X (θ 1 , . . . , θ n ) = ϕ X (θ 1 , . . . , θ n ) = E {exp [-(θ 1 X 1 + • • • + θ n X n )]}
, and let ϕ Xi be the Laplace transform of the random variable X i , defined by

L Xi (θ i ) = ϕ Xi (θ i ) = E [exp (-θ i X i )] .
C L X (u 1 , . . . , u n ) = ϕ X ϕ -1 X1 (u 1 ) , . . . , ϕ -1 Xn (u n )
is a copula. We call this copula C L X the Laplace copula associated to the random vector X. If X has pd

µ X (dx) , then we denote still C L X (u 1 , . . . , u n ) = C L µ X (u 1 , . . . , u n ) .
Proposition 14 Let F 1 , . . . , F n , be n univariate cdfs on ]0, +∞[ . The relation

F (x 1 , . . . , x n ) = ϕ X ϕ -1 X1 [F 1 (x 1 )] , . . . , ϕ -1 Xn [F n (x n )] = C L X [F 1 (x 1 ) , . . . , F n (x n )]
defines a cdf F with marginal cdfs F 1 , . . . , F n . Then the copula associated to F is C L X .

If we choose, for i = 1, 2, . . . , n, F i = F Xi , the cdf of X i , then F is the cdf of a random vector X for which the marginal cdfs are the cdfs of X i , and we have F = F X where F X calls the cdf of X. In the case of multivariate gamma distributions and multi-factor gamma distributions, we get other multivariate gamma distributions and multi-factor gamma distributions with given copula.

Laplace copula can be seen as a generalization of Archimedean copula [START_REF] Joe | Multivariate Models and Dependance Concepts[END_REF][START_REF] Joe | Dependence Modeling with Copulas[END_REF].

Main results

Now, we are in the capacity to give the two main results of this paper Theorem 15 Let P an affine polynomial in the n variables θ i , i = 1, . . . , n, with P (0) = 1, then the Laplace copula of the multivariate gamma distribution γ (P,λ) associated to (P, λ) is

C Lγ (P,λ) (v) = v [n] [1 + T ⊂[n],|T |>1 (-1) |T | P (- 1 p 1 T ) t∈T (1 -v 1 λ t )] -λ , (20)
where v = (v 1 , . . . , v n ) , |T | is the cardinality of T and the vector 1 p 1 T is defined by

1 p 1 T i = 1 pi if i ∈ T , 1 p 1 T i = 0 if i / ∈ T , for i ∈ {1, 2, . . . , n}.
This family is given, by example, for the simpler case P (θ) = n i=1 (1+p i θ i )]-βp [n] θ [n] , with 0 β < 1, corresponding to the family gived by [START_REF] Fang | A family of bivariate distributions with non-elliptical contours[END_REF], namely

C (v) = v [n] [1 -β n i=1 (1 -v 1 λ i )] -λ .
Theorem 16 Let P an affine polynomial in the n variables θ i , i = 1, . . . , n, with P (0) = 1, let Λ = (λ, λ 1 , . . . , λ n ) with λ i > λ > 0, for i ∈ {1, 2, . . . , n} , then the Laplace copula of the multivariate multifactor gamma distribution γ (P,Λ) associated to (P, Λ) is

C Lγ (P,Λ) (v) = v [n] [1 + T ⊂[n],|T |>1 (-1) |T | P (- 1 p 1 T ) t∈T (1 -v 1 λ t t )] -λ . ( 21
)
This copula family is a new generalization of Farlie-Gumbel-Morgenstern Copulas.

Bekrizadeh et al., 2012 propose a similar formula only for λ being a negative integer.

We note that, if the conditions of Theorem 5 are checked, formulas (20) and (21) define copulas.

For these copulas, if we inject univariate gamma distributions γ (pi,λ) in the first case or γ (pi,λi) in the second case, then we obtain other multivariate gamma distributions and other multivariate multifactor gamma distributions in the sense that their marginal distributions are respectively γ (p1,λ) and γ (p1,λi) . Their cdf and pdf are given by ( 1) and ( 2) respectively, and have a link with multivariate gamma distribution associated with (P, λ) and multi-factor gamma distribution associated with (P, Λ) respectively.

The bidimensional and tridimensional cases

For the bidimensional case, Theorem [START_REF] Slater | Generalized Hypergeometric Functions[END_REF] and Theorem (16) give the following corollary Corollary 17 For the bivariate gamma distribution such that L γ (P,λ) (θ) = (P (θ))

-λ , we have

C Lγ (P,λ) (v 1 , v 2 ) = v 1 v 2 [1 -r 12 (1 -v 1 λ 1 )(1 -v 1 λ 2 )] -λ ,
where, for X = (X 1 , X 2 ) with pd γ (P,λ) , r 12 = -P (-p -1 1 , -p -1 2 ) = 1 -p 12 / (p 1 p 2 ) is the linear correlation coefficient of the random variables X 1 , X 2 ; it checks 0 r 12 1. It is the copula of the BB10 family p. 154 in the Joe's book [START_REF] Joe | Multivariate Models and Dependance Concepts[END_REF]). This result is due to the following equalities :

E (X i ) = λp i, V ar (X i ) = λp 2 i , Cov(X 1 , X 2 ) = λ (-p 12 + p 1 p 2 )
obtained by derivating L γ (P,λ) (θ) at 0. For the bivariate multi-factor gamma distribution such that L γ (P,Λ) (θ) = (P (θ)) -λ 2 i=1 (1 + p i θ i ) -(λi-λ) , we have a more general family

C Lγ (P,Λ) (v 1 , v 2 ) = v 1 v 2 [1 -r 12 (1 -v 1 λ 1 1 )(1 -v 1 λ 2 2 )] -λ , ( 22
)
where, for

X = (X 1 , X 2 ) = (Y 1 , Y 2 ) + (Z 1 , Z 2 ) = Y + Z
as defined in Proposition 1 with pd γ (P,Λ) , and

0 r 12 = -P (-p -1 1 , -p -1 2 ) = 1 -p 12 / (p 1 p 2 ) 1, is the linear correlation coefficient of the random variables Y 1 , Y 2 in the bivariate gamma distribution such that L Y (θ) = [P (θ)] -λ .
This bivariate family is not given in the Joe's books [START_REF] Joe | Multivariate Models and Dependance Concepts[END_REF][START_REF] Joe | Dependence Modeling with Copulas[END_REF]. Bekrizadeh et al., 2012 give a similar formula for λ being a negative integer.

We recall the formulas [START_REF] Joe | Dependence Modeling with Copulas[END_REF] for the computation of τ, the Kendall's tau, and ρ S , the Spearmann's rho

τ = 1 -4 [0,1] 2 ∂C ∂u (u, v) ∂C ∂v (u, v) dudv, ( 23
)
ρ S = 12 [0,1] 2 C (u, v) dudv -3. ( 24
)
Then, we can give the following result Proposition 18 The Kendall's tau and the Spearman's rho of the Laplace copula of the bivariate multifactor gamma distribution are given by the following formulas

τ = 1 -F 3 2 (2λ, 1, 1; 2λ 1 + 1, 2λ 2 + 1; r 12 ) + 4λ (2λ 1 + 1) (2λ 2 + 1) r 12 F 3 2 (2λ + 1, 1, 2; 2λ 1 + 2, 2λ 2 + 2; r 12 ) - λ 2 (2λ 1 + 1) (2λ 2 + 1) (λ 1 + 1) (λ 2 + 1) r 2 12 F 3 2 (2λ + 2, 2, 2; 2λ 1 + 3, 2λ 2 + 3; r 12 ) , (25)
and

ρ S = 3 F 3 2 (1, 1, λ; 2λ 1 + 1, 2λ 2 + 1; r 12 ) -1 (26) = 3λ (2λ 1 + 1) (2λ 2 + 1) r 12 F 3 2 (1, 2, λ + 1; 2λ 1 + 2, 2λ 2 + 2; r 12 ). ( 27
)
For the tridimensional case, Theorem [START_REF] Slater | Generalized Hypergeometric Functions[END_REF] and Theorem (16) give the following corollary

Corollary 19 For the trivariate gamma distribution γ (P,λ) such that L γ (P,λ) (θ) = (P (θ))

-λ , we have

C Lγ (P,λ) (v 1 , v 2 , v 3 ) = v 1 v 2 v 3 [1 -r 12 (1 -v 1 λ 1 )(1 -v 1 λ 2 ) -r 13 (1 -v 1 λ 1 )(1 -v 1 λ 3 ) -r 23 (1 -v 1 λ 2 )(1 -v 1 λ 3 ) + 2r 123 (1 -v 1 λ 1 )(1 -v 1 λ 2 )(1 -v 1 λ 3 )] -λ ,
where, for X = (X 1 , X 2 , X 3 ) with pd γ (P,λ) , we denote by r ij = -P (-

1 p 1 {i,j} ) = 1 -p ij /(p i p j ), 1 i = j
3, the linear correlation coefficient of the random variables X i , X j , they check 0 r ij 1, and we denote by r 123 the number defined by

r 123 = E{ 3 i=1 [X i -E (X i )]}/ 3 i=1 (E{[X i -E (X i )] 3 }) 1/3 (28) = -1 2 P (-p -1 1 , -p -1 2 , -p -1 3 ) (to compare with r 12 = E{ 2 i=1 [X i -E (X i )]}/ 2 i=1 (E{[X i -E (X i )] 2 }) 1/2 ).
This last result is due to the following equalities :

E (X i ) = λp i, E{[X i -E (X i )] 3 } = 2λp 3 i , E{ 3 i=1 [X i -E (X i )]} = -λp 1 p 2 p 3 P (-p -1 1 , -p -1 2 , -p - 1
3 ) obtained by derivating L γ (P,λ) (θ) at 0.

For the trivariate multi-factor gamma distribution such that L γ (P,Λ) (θ) = [P (θ)]

-λ 3 i=1 (1 + p i θ i ) -(λi-λ) (λ i λ, i = 1, 2, 3 ), we have C Lγ (P,Λ) (v 1 , v 2 , v 3 ) = v 1 v 2 v 3 [1 -r 12 (1 -v 1 λ 1 1 )(1 -v 1 λ 2 2 ) -r 13 (1 -v 1 λ 1 1 )(1 -v 1 λ 3 3 ) -r 23 (1 -v 1 λ 2
2 )(1 -v

1 λ 3
3 ) + 2r 123 (1 -v

1 λ 1 1 )(1 -v 1 λ 2
2 )(1 -v

1 λ 3
3 )] -λ , where, for 6 Appendix 6.1 Proof of Theorem 3

X = (X 1 , X 2 , X 3 ) = (Y 1 , Y 2 , Y 3 ) + (Z 1 , Z 2 , Z 3 ) = Y + Z
Let X = (X 1 , X 2 ) and Y = (Y 1 , Y 2 ) be independent random variables. The pdf of Z = (X 1 + Y 1 , X 2 + Y 2 )
where X i has pd γ (p i , λ i -λ) , i = 1, 2 and Y has pd γ (P, λ), is obtained by convolution. By changing

variables 1 -v i /z i = u i , i = 1, 2, and with the notation c = b 12 ( P ) = -p -1 12 + p 1 p 2 p -2 12 , d i = cp 12 p -1 i , i = 1, 2 we obtain the pdf f Z of Z by f Z (z 1 , z 2 ) = 0 if (z 1 , z 2 ) / ∈ (0, ∞) 2 and for (z 1 , z 2 ) ∈ (0, ∞) 2 f Z (z 1 , z 2 ) = z1 0 z2 0 1 [Γ (λ)] 2 p λ 12 e - p 2 p 12 v1+ p 1 p 12 v2 (v 1 v 2 ) λ-1 F 1 (λ; cv 1 v 2 ) × e - z 1 -v 1 p 1 Γ (λ 1 -λ) z 1 -v 1 p 1 λ1-λ e - z 2 -v 2 p 2 Γ (λ 2 -λ) z 2 -v 2 p 2 λ2-λ dv 1 z 1 -v 1 dv 2 z 2 -v 2 = z λ1-1 1 z λ2-1 2 e - p 2 p 12 z1+ p 1 p 12 z2 (Γ (λ)) Γ (λ 1 -λ) Γ (λ 2 -λ) p λ 12 p λ1-λ 1 p λ2-λ 2 k 0 (cz 1 z 2 ) k Γ (λ + k) k! × 1 0 e d1z1u1 u λ1-λ-1 1 (1 -u 1 ) λ+k-1 du 1 1 0 e d2z2u2 u λ2-λ-1 2 (1 -u 2 ) λ+k-1 du 2 . (29)
As we have, with the notation B (α, β) =

1 0 u α-1 (1 -u) β-1 du = Γ (α) Γ (β) /Γ (α + β), α, β > 0 for the Euler's Beta function, 1 0 e δu u α-1 (1 -u) β-1 du = B (α, β) ∞ n=0 (α) n (α + β) n δ n n! , we obtain f Z (z) = z λ1-1 1 z λ2-1 2 e - p 2 p 12 z1+ p 1 p 12 z2 Γ (λ) Γ (λ 1 -λ) Γ (λ 2 -λ) p λ 12 p λ1-λ 1 p λ2-λ 2 k 0 (cz 1 z 2 ) k Γ (λ + k) k! Γ (λ 1 -λ) Γ (λ + k) Γ (λ 1 + k) × ∞ n=0 (λ 1 -λ) n (λ 1 + k) n (d 1 z 1 ) n n! Γ (λ 2 -λ) Γ (λ + k) Γ (λ 2 + k) ∞ m=0 (λ 2 -λ) m (λ 2 + k) m (d 2 z 2 ) m m! = z λ1-1 1 z λ2-1 2 e - p 2 p 12 z1+ p 1 p 12 z2 Γ (λ 1 ) Γ (λ 2 ) p λ 12 p λ1-λ 1 p λ2-λ 2 ∞ k=0 ∞ n=0 ∞ m=0 (λ 1 -λ) n (λ 2 -λ) m (λ) k (λ 1 ) k+n (λ 2 ) k+m (d 1 z 1 ) n n! (d 2 z 2 ) m m! (cz 1 z 2 ) k k! .
Hence, we have proved the formula (11).

Proof of Theorem 7

We start from the Laplace transform of the trivariate gamma distribution γ (λ,P ) associated to (λ, P ).

First, we write

L γ (λ,P ) (θ) = (1 + p 1 θ 1 + p 2 θ 2 + p 3 θ 3 + p 12 θ 1 θ 2 + p 13 θ 1 θ 3 + p 23 θ 2 θ 3 + p 123 θ 1 θ 2 θ 3 ) -λ = (1 + p 3 θ 3 ) -λ 1 + p 1 + p 13 θ 3 1 + p 3 θ 3 θ 1 + p 2 + p 23 θ 3 1 + p 3 θ 3 θ 2 + p 12 + p 123 θ 3 1 + p 3 θ 3 θ 1 θ 2 -λ . Let Q (θ 1 , θ 2 ) = 1+q 1 θ 1 +q 2 θ 2 +q 12 θ 1 θ 2 where q 1 = (p 1 + p 13 θ 3 ) / (1 + p 3 θ 3 ) , q 2 = (p 2 + p 23 θ 3 ) / (1 + p 3 θ 3 ) ,
and q 12 = (p 12 + p 123 θ 3 ) / (1 + p 3 θ 3 ). With these notations, we have

L γ (λ,P ) (θ) = (1 + p 3 θ 3 ) -λ (1 + q 1 θ 1 + q 2 θ 2 + q 12 θ 1 θ 2 ) -λ . Let p T = -p T /p 123 , b ij = p ij + p i p j ,
( Q) = -q 12 + q 1 q 2 q 2 12 = b 12 + b 123 (θ 3 -p 3 ) + b 13 b 23 (θ 3 -p 3 ) 2 .
Let γ (λ,Q) be the bivariate gamma distribution associated to (λ, Q), its Laplace transform is

(1 + q 1 θ 1 + q 2 θ 2 + q 12 θ 2 θ 3 ) -λ . We have γ (λ,Q) (dx) = q -λ 12 [Γ (λ)] 2 e - q 2 q 12 x1- q 1 q 12 x2 (x 1 x 2 ) λ-1 F 1 [λ, b 1,2 ( Q)x 1 x 2 ]½ (0,∞) 2 (x 1 , x 2 ) dx 1 dx 2 .
Second, we are looking for the unidimensional distribution for which its Laplace transform with respect to the variable θ 3 is equal to

(1 + p 3 θ 3 ) -λ q -λ 12 (Γ (λ)) 2 e -q 2 q 12 x1- q 1 q 12 x2 (x 1 x 2 ) λ-1 F 1 [λ, b 1,2 ( Q)x 1 x 2 ]½ (0,∞) 2 (x 1 , x 2 ) dx 1 dx 2 = p -λ exp( b 13 x 1 + b 23 x 2 θ 3 ) (32)
From (Hadlik,1986), we have the following equality

L [Γ(λ)] -1 t λ-1 F1(λ,at)½ (0,∞) (t) dt (s) = s -λ exp a s (33)
We utilize the equality (33) in (32) and we obtain the following unidimensional distribution

∞ k=0 1 (λ) k (x 1 x 2 ) k ℓ+m+n=k 1 ℓ!m!n! b ℓ 13 b ℓ 23 b m 123 b n 12 [Γ (λ + 2ℓ + m)] -1 x λ+2ℓ+m-1 3 × F 1 λ + 2ℓ + m, ( b 13 x 1 + b 23 x 2 )x 3 ½ (0,∞) (x 3 ) dx 3 = [Γ (λ)] -1 x λ-1 3 ½ (0,∞) (x 3 ) × (ℓ,m,n)∈N 3 1 (λ) (ℓ+m+n) (λ) (2ℓ+m) ℓ!m!n! ( b 13 x 1 x 3 b 23 x 2 x 3 ) ℓ ( b 123 x 1 x 2 x 3 ) m ( b 12 x 1 x 2 ) n × F 1 (λ + 2ℓ + m, b 13 x 1 x 3 + b 23 x 2 x 3 ) dx 3 = [Γ (λ)] -1 x λ-1 3 ½ (0,∞) (x 3 ) × (ℓ,m,n,p)∈N 4 1 (λ) (ℓ+m+n) (λ) (2ℓ+m+p) ℓ!m!n!p! × ( b 13 x 1 x 3 b 23 x 2 x 3 ) ℓ ( b 123 x 1 x 2 x 3 ) m ( b 12 x 1 x 2 ) n ( b 13 x 1 x 3 + b 23 x 2 x 3 ) p dx 3 = [Γ (λ)] -1 x λ-1 3 F II (λ, λ, b 13 x 1 x 3 b 23 x 2 x 3 , b 123 x 1 x 2 x 3 , b 12 x 1 x 2 , b 13 x 1 x 3 + b 23 x 2 x 3 )½ (0,∞) (x 3 ) dx 3 (34)
We deduce from equality (35)

L γ (P,λ) (θ) = [ T ⊂[n] (-1) |T | P (- 1 p 1 T ) u -1 [n] T u T ] -λ = u [n] -λ [ T ⊂[n] (-1) |T | P (- 1 p 1 T )(1 [n] - 1 u ) T ] -λ = v [n] [1 + T ⊂[n],|T |>1 (-1) |T | P (- 1 p 1 T ) t∈T (1 -v 1 λ t )] -λ = C Lγ (P,λ) (v) .

Proof of Theorem 16

We start from the Laplace transform of the multivariate multi-factor gamma distribution associated to (P, Λ) , Λ = (λ, λ 1 , . . . λ n ), where P is an affine polynomial with respect to θ = (θ 1 , . . . , θ n ) and

λ i λ, for i ∈ [n]
. We use the following notations

u i = (1 + p i θ i ) -λ and v i = (1 + p i θ i ) -λi so that (1 + p i θ i ) -(λi-λ) = v i /u i and (1 + p i θ i ) -1 = u 1/λ i = v 1/λi i
. According to equality 20, we have

L γ (P,Λ) (θ) = [P (θ)] -λ i∈[n] (1 + p i θ i ) -(λi-λ) = u [n] [1 + T ⊂[n],|T |>1 (-1) |T | P (- 1 p 1 T ) t∈T (1 -u 1 λ t )] -λ i∈[n] v i u i = v [n] [1 + T ⊂[n],|T |>1 (-1) |T | P (- 1 p 1 T ) t∈T (1 -u 1 λ t )] -λ = v [n] [1 + T ⊂[n],|T |>1 (-1) |T | P (- 1 p 1 T ) t∈T (1 -v 1 λ t t )] -λ = C Lγ (P,Λ) (v) .
6.5 Proof of Proposition 18

Kendall's tau

By injecting the given copula in formula (22) in formula (23), we obtain

1 -τ 4 = [0,1] 2 u 1 u 2 [1 -r 12 (1 -u 1 λ 1 1 )(1 -u 1 λ 2 2 )] -2λ-2 {1 -r 12 [1 -(1 - λ λ 1 )u 1 λ 1 1 ](1 -u 1 λ 2 2 )}{1 -r 12 [1 -(1 - λ λ 2 )u 1 λ 2 2 ](1 -u 1 λ 1
1 )} du 1 du 2 . Finally, by using equality (6), we have proved equality (25).

By changing the variables t

i = 1 -u 1/λi i , i = 1, 2,

Spearman's rho

By injecting the given copula in formula (22) in formula ((24), we obtain

̺ S = 12 1 0 1 0 u 1 u 2 [1 -r 12 (1 -u 1 λ 1 1 )(1 -u 1 λ 2
2 )] -λ du 1 du 2 -3.

(36) By changing the variables v i = u 1 λ i i , i = 1, 2, we obtain

1 0 1 0 u 1 u 2 [1 -r 12 (1 -u 1 λ 1 1 ) (1 -u 1 λ 2
2 )] -λ du 1 du 2 This gives equality (27).

= 1 0 1 0 λ 1 λ 2 v 2λ1-1 1 v 2λ2-1 2 [1 -r 12 (1 -v 1 ) (1 -v 2 )] -λ dv 1 dv 2 = λ 1 λ 2 ∞ k=0 (λ) k k! r k 12 1 0 v 2λ1-1 1 (1 -v 1 ) k dv 1 1 0 v 2λ2-1 2 (1 -v 2 ) k dv 2 = λ 1 λ 2 ∞ k=0 (

Acknowledgement

I thank Gérard Letac for many helpful conversations.

Remark 9

 9 If b 12 = b 13 = b 23 = 0, Theorem 7 gives

The function ϕ

 ϕ Xi is a one and onto, and is decreasing of [0, +∞[ onto ]0, 1]. Its inverse function is denoted by ϕ -1Xi . Then the function C L X defined by

, 1 i = j 3 ,

 3 as defined in Proposition 1 with pd γ (P,Λ) , we denote by r ij = -P (-1 p 1 {i,j} ) = 1 -pij pipj the linear correlation coefficient of the random variables Y i , Y j in the trivariate gamma distribution such that L Y (θ) = (P (θ)) -λ , they check 0 r ij 1, and r 123 is defined by equality (28).

( 1 )( 1 ) 2 ( 12 ∞ k=0 ( 2 )

 11212k=02 With the equalityB (α, β) = Γ (α) Γ (β) /Γ (α + β) for α, β > 0, we obtain 1 -τ = ∞ k=0 k (1) k (2λ) k (2λ 1 + 1) k (2λ 2 + 1) k (2) k (2λ + 1) k (2λ 1 + 2) k (2λ 2 + 2) 2λ 1 + 1) (2λ 2 + 1) (λ 1 + 1) (λ 2 + 1) r 2 k (2) k (2λ + 2) k (2λ 1 + 3) k (2λ 2 + 3) k

  and b 123 = p 123 + p 1 p 23 + p 2 p 13 + p 3 p 12 + 2 p 1 p 2 p 3 , then we have b 1,2

  (1 -r 12 t 1 t 2 ) -2λ (1 -t 1 ) 2λ1-1 (1 -t 2 ) 2λ2-1 -λ λ 1 r 12 (1 -r 12 t 1 t 2 ) -2λ-1 (1 -t 1 ) 2λ1 t 2 (1 -t 2 ) -r 12 t 1 t 2 ) -2λ-1 t 1 (1 -t 1 ) 2λ1-1 (1 -t 2 ) -r 12 t 1 t 2 ) -2λ-2 t 1 (1 -t 1 ) 2λ1 t 2 (1 -t 2 ) 2λ2 dt 1 dt 2 .By using equality (1 -r 12 t 1 t 2 )

								we obtain
	1 -τ 4λ 1 λ 2	=		[0,1] 2
								2λ2-1
		-	λ λ 2	r 12 (1 2λ2
		+	λ 2 λ 1 λ 2 12 (1 -2λ = r 2	∞ k=0 (2λ) k	(r k 12 t k 1 t k 2 ) k!	in the last result, we obtain,
	1 -τ 4λ 1 λ 2	=	∞ k=0	(2λ) k	r k 12 k!	B (k + 1, 2λ 1 ) B (k + 1, 2λ 2 )
		-	λ λ 1	r 12	∞ k=0	(2λ + 1) k	r k 12 k!	B (k + 1, 2λ 1 + 1) B (k + 2, 2λ 2 )
		-	λ λ 2	r 12	∞ k=0	(2λ + 1) k	r k 12 k!	B (k + 2, 2λ 1 ) B (k + 1, 2λ 2 + 1)
		+	λ 2 λ 1 λ 2	r 2 12	∞ k=0	(2λ + 2) k	r k 12 k!	B (k + 2, 2λ

1 + 1) B (k + 2, 2λ 2 + 1) .

  By injecting equality (37) in equality (36), we obtain equality (26). We remark that we haveF 3 2 (λ, 1, 1; 2λ 1 + 1, 2λ 2 + 1; r 12 ) -1 =

				λ) k k!	r k 12 B (2λ 1 , k + 1) B (2λ 2 , k + 1)
	=	1 4	∞ k=0	(λ) k (1) 2 k (2λ 1 + 1) k (2λ 2 + 1) k	r k 12 k!
	=	1 4	F 3 2 (1, 1, λ; 2λ 1 + 1, 2λ 2 + 1; r 12 )	(37)
					λ (2λ 1 + 1) (2λ 2 + 1)	r 12	∞ k=0	(k + 1)!	(λ + 1) k (2λ 1 + 2) k (2λ 2 + 2) k	r k 12 .

Proof of Theorem 15

We utilize the following notations, for i ∈

We start from the Laplace transform of the multivariate gamma distribution associated to (P, λ).

From the definition of γ (P,λ) , we obtain

Hence, we have

where u = (u 1 , . . . , u n ) , and

form a basis of the vector space of affine polynomials of degree less or equal to n. Because

we have

consequently

).

Since P is an affine polynomial, we have

in particular α T = 0, if |T | = 1 and α ∅ = 1. Finally we obtain