Daniel Hirst
email: daniel.hirst@lpl-aix.fr

Saandia Ali
email: saandia.ali@lpl-aix.fr

Optimizing the automatic functional annotation of English intonation

One of the fundamental aims of prosodic analysis is to provide a reliable means of extracting functional information (what prosody contributes to meaning) directly from prosodic form. It has been argued that an explicit model of the mapping from prosodic function to prosodic form could provide an objective way of approaching this task. In this presentation we look specifically at some of the problems of optimizing this mapping in order to extract the functional information automatically from the formal representation, hence ultimately directly from the acoustic data.

Introduction

Phonological models of intonation and prosody in general are notoriously theory dependent. Such models obviously need to address the problem of validation. What criteria can be used to decide that one particular representation of is more adequate than another.

There has, in the past, been little or no consensus on either the way to represent the forms of prosodic patterns, or on how to annotate prosodic functions, the way that speech prosody contributes to the interpretation of an utterance. In many systems of representation, formal and functional aspects are in fact conflated. This is the case for example with the most widely used annotation system, ToBI [START_REF] Silverman | ToBI : a Standard for Labelling English Prosody[END_REF], originally designed as a standard for the representation of American English intonation, but later adapted to a number of other languages [START_REF]Prosodic Typology and Transcription: A Unified Approach[END_REF]. It has been argued [START_REF] Hirst | Form and function in the representation of speech prosody[END_REF] that maintaining a clear separation between prosodic form and prosodic function would provide a more effective means of validating or evaluating specific proposals.

In an earlier study [START_REF] Ali | Analysis by synthesis of English intonation patterns: generalising from form to function[END_REF], we suggested that a solution to the problem of evaluating models could be found by a process of analysis by synthesis -generating prosodic forms by an explicit mapping from a representation of prosodic functions and comparing the output with observable data.

We showed specifically how five successively more complex representations of English intonation patterns could be converted by rule into surface phonological representations of intonation using the INTSINT alphabet [5][3]. These can in turn be converted into phonetic representations by means of the Momel algorithm and the ouput can then be compared directly to the original recordings of 40 continuous 5 sentence passages from the Eurom1 corpus.

The five successive models we used to illustrate the process were : a. none. The intonation of the whole passage was modelled as a rising then falling pattern represented with INTSINT as [M-T…B] (where "-" indicates a fixed duration between the two tones and where [ and ] correspond to the beginning and end of the passage).

b. IU (Intonation Unit). Taking the boundaries of intonation units as functional annotation, each IU was modelled with the sequence [M-T…B] (where here [ and ] correspond to the beginning and end of the intonation unit). f. emphasis in this model a distinction was made between non-emphatic patterns, coded as in model (e) and emphatic nuclei and heads, which were assigned more complex patterns, described in detail in [START_REF] Ali | Analysis by synthesis of English intonation patterns: generalising from form to function[END_REF]. This final pattern corresponds to the IF (intonation functions) annotation system originally proposed in [START_REF] Hirst | Intonative Features. A Syntactic Approach to English Intonation[END_REF] and more recently in [START_REF] Hirst | A Praat plugin for Momel and INTSINT with improved algorithms for modelling and coding intonation[END_REF]. Figure 1 shows an example of one of the Eurom1 passages with the functional information corresponding to model (f). In fact all the other models were derived from this functional representation by adapting the mapping rules appropriately. By comparing the output of the successive models with the original recordings, it was shown that progressively increasing the complexity of the functional information leads to a corresponding increase in the correlation between the predicted values and the observed values, as can be seen in Table 1. We noted that while this result is encouraging, none of the mean correlations obtained were particularly good (although several individual values for specific sentences were much higher, often over 0.9). It should be noted, though, that in this preliminary study no attempt was made to optimise the parameters of the model. The mapping rules were formulated by hand, as was the specific functional annotation.

Nevertheless, the study demonstrated that it is possible, using a technique of this type, to quantify the contribution of a particular functional annotation to the quality of the model.

In the rest of this paper we present some further research on this problem and in particular we make some preliminary steps towards deriving a functional representation such as that of Figure 1 directly from the acoustic signal rather than annotating it by hand as in our earlier presentation.

Optimizing the representation of prosodic form

The first stage of the optimization of the mapping between function and form was carried out using model (b) described above ie "IU" (Intonation Unit). The sentences recorded by four different speakers from the Eurom1 corpus were analyzed. The only functional information taken into account at this stage was the pause, hand-labelled and represented by "#". These pauses coincide with the five sentences read by each speaker within one sound file.

Optimizing the coding of INTSINT tones

It was decided to represent prosodic form between two pauses, using just four INTSINT tones: [t1-t2….t3-t4] where t1 and t2 were aligned respectively with the left and right boundaries and t3 and t4 were aligned at a fixed duration (200ms) after or before the initial or final boundary. By means of a Praat script, and for each extract, all of the combinations of four INTSINT tones were tested except that the first tone of a unit was necessarily an absolute tone, either T, M or B, since a relative tone presupposes that there is a preceding target. For the other tones there were 8 possible tones at each point and we also included the possibility of no tones at all (_) at each point except the last. This gave a total of 1944 (=3*9*9*8) possible sequences of tones.

Results of the experiment

The F0 curve was calculated for each possible combination of tones and the correlation with the hand-corrected Momel curve was calculated. The combination of tones which obtained the best correlation on all of the extracts was then selected as shown in table 2.

Optimizing the alignment of the tones

Once an optimal sequence of tones was defined for each speaker: Fe, FF (MULL), Fa (TTLB), Fg (T_LB), we tried to optimize the alignment of these tones. Only the alignment tones t2 and t3 were optimized at this stage leaving t1 and t4 at a fixed offset. Thus for speaker Fe the tones were fixed as follows:

M fixed at an offset of 10ms from the left boundary U from 100 to 800 ms after the left boundary with steps of 50ms (15 iterations) L from 100 to 800 ms before the right boundary with steps of 50ms (15 iterations) L fixed at an offset of 10 ms from the right boundary

The pair of points obtaining the best correlation with the original Momel curve was selected using the same method as in the optimization of tones.

Optimizing the representation of prosodic function

The ultimate aim of this experiment is to extract the functional representation automatically from the formal representation. Hence, using the optimization of the sequence of tones and their alignment, our next attempt is to automatically detect intermediate boundaries of intonation units between two pauses.

Automatic segmentation of speech using melody only: from form to function

The minimal functional information used for this experiment was the word boundaries, on the fairly uncontroversial assumption that an Intonation Unit boundary must coincide with a word boundary.

A Praat script was used to test every possible location for an Intonation Unit boundary. In each case, the optimal sequence of tones, as determined from the first experiment, was applied and the result of the correlation saved for comparison with other possible boundary locations.

The ultimate aim was to obtain a better correlation with the introduction of a new boundary and so the simple criterion for a new intonation unit boundary was that the correlation should be higher with two groups than with just one.

To begin with, the correlation was calculated for the whole passage, then each possible location was tested and the boundary corresponding to the highest correlation was retained provided that this correlation was higher than that obtained with no boundary. The process was then iterated until no further boundaries improved the correlation.

Two different series of detection of boundaries were run, the first using fixed parameters of tones and alignment as mentioned above and the second one allowing a variation of the alignment of the tones between the values giving the best correlations. Thus for Speaker Fe, the first experiment was run with the sequence of tones [M-U…L-L] and a fixed alignment of t2 and t3 (where t2=0.2 and t3=0.7) and in the second one t2 varied from 0.2 to 0.3 with three iterations and t3 from 0.65 to 0.8 with three iterations. It was expected that a model of detection allowing more flexibility in terms of target alignment would improve the detection of boundaries.

Evaluation of the detection of boundaries

In order to evaluate, the automatic detection of boundaries, the results were compared with the annotation by the authors of this paper.(table 3). Boundaries were added by hand and qualified as terminal or non terminal. The results of the new experiments of boundary detections were classified in keeping with the type of experiment and the type of boundary (terminal/ non terminal).

Table3. Results of the detection of boundaries. IU= hand-labeled Intonation

Recurrent errors

As can be seen in table 3, the second experiment of detection of boundaries (DB2) found 20 more boundaries than the first one. Most terminal boundaries were detected in both cases (46 out of 61/ 45 out of 61) while most non-terminal boundaries were missed (48 out of 328/ 70 out of 301). This suggests that a further optimisation of the model should be developed taking into account the distinction between terminal and nonterminal boundaries. Indeed, all of the optimised sequences of tones used for our speakers so far contain a final falling tone usually associated with terminal boundaries.

Furthermore, boundaries were frequently identified just one or two words after the one annotated by a linguist. (noted NextW/O in table 3). These were generally grammatical words like pronouns or articles, which are unaccented and sometimes also reduced. It seems that no major change of the shape of F0 can be noticed around these points so that a better correlation is obtained by adding them to the preceding intonation unit.

Examples for file Feo1074 : In these examples a boundary was detected just before what was interpreted as the nucleus. The higher value of F0 at the beginning of the nucleus was interpreted as the beginning of a new IU. This suggests that a more complex model taking into account accents might provide a better fit.

I

Results

In this section, we summarize the results obtained after the different steps of optimisation. After the automatic detection of intermediate boundaries, the correlations between the original Momel curves and the output of our model were evaluated anew.

Conclusions

In this paper, we make a first attempt at optimizing the parameters of a previously described model of the mapping between prosodic function and prosodic form. Initially, the only functional information taken into account was the pause, which, in read speech, can be fairly safely be assumed to correspond to a prosodic boundary [START_REF] Cho | Empirical evidence for prosodic phrasing: pauses as linguistic annotation in Korean read speech[END_REF]. The prosodic form was optimized as a sequence of four INTSINT tones the alignment of which was in turn optimized for four speakers. Our next endeavor was to use these optimal forms to try to automatically detect intermediate intonation unit boundaries. The results of this experiment were evaluated by comparison with the annotation of boundaries by hand. These were generally quite promising. A closer look at the errors and missing boundaries suggested the need to develop a more complex model which would take into account the type of boundary (terminal /non terminal) and prominences thus getting closer to the more complex models of synthesis described in section one.

  c. terminal Intonation boundaries were marked as either terminal or non-terminal and modelled respectively as [M-T…B-B] and [M-T…B-H]. d. accent In addition to the tones assigned in model (c), in this model the T tone is aligned with the first accent and each subsequent accent is assigned a D tone. e. nucleus In this model, a distinction was made between prenuclear accents coded as in model (d). and post-nuclear accents which are assigned a B tone.

Figure 1 .

 1 Figure 1. Passage fao30072 with the final functional annotation used in the analysis.

  Units. DB1= detection of boundaries with fixed alignment of t2-t3. DB2= variation of t2-t3's alignment. Nterm= Non terminalboundary/ Term= terminal boundary.

Figure 2 .

 2 Figure 2. Different models of detection of boundaries for the sentence : Unfortunately, she's going to leave us soon. , tier 3= IU, tier8=DB2, Tier 9= DB1

Table 1

 1

		mean and standard deviation of correlation
	coefficients for each Intonation Unit between the
	output of the model and the output from the hand-
	corrected Momel targets
	model	none	IU	terminal accent emphasis
	mean	0.1238 0.5714 0.5215 0.5902 0.6573
	sd	0.2929 0.2533 0.3386 0.2754 0.2238

Table 2 .

 2 Best sequences of tones and target alignment after optimization

	Speaker Form	Left align t1	Left align t2	Right align t3	Right align t4	Mean Correlation
	FA	T T L B	0.01	0.2	0.7	0.01	0.666
	FE	M U L L 0.01	0.2	0.7	0.01	0.683
	FF	M U L L 0.01	0.2	0.4	0.01	0.650
	FG	T _ L B	0.01	0.1	0.4	0.01	0.642

Table 4 .

 4

	mean and standard deviation of correlation
	coefficients for each Intonation Unit between the
	output of the model and the output from the hand-
	corrected Momel targets
	DB1	DB2	No bound
	correl NB Mean SD NB Mean SD NB Mean SD
	>= 0.75 146 0.833 0.06 175 0.840 0.06 141 0.834 0.06
	>= 0.50 163 0.638 0.07 146 0.656 0.06 147 0.627 0.07
	< 0.50 44 0.197 0.40 31 0.176 0.42 65 0.251 0.34
	Total 353 0.664 0.249 352 0.705 0.232 353 0.640 0.26
	In both experiments (DB1 and DB2), the mean correlation for
	all speakers increased, rising from 0.57 to 0.66 and 0.705.
	(see table 1 and

table 4 )

 4 . It must be added, though, that these results are not entirely comparable, since the first table sums up results taken from just one speaker as opposed to four speakers in the second. Furthermore, model (b) was based on the Intonation unit labelled by hand as opposed to an automatic detection in the second case.


Dominique Leguillon
email: dol@ccr.jussieu.fr

Strength or toughness? A criterion for crack onset at a notch

Both energy and stress criteria are necessary conditions for fracture but neither one nor the other are sufficient. Experiments by Parvizi et al. on transverse cracking in cross-ply laminates corroborate this assumption. Thanks to the singularity at the tip of the notch, the incremental form of the energy criterion gives a lower bound of admissible crack lengths. On the contrary, the stress criterion leads to an upper bound. The consistency between these two conditions provides a general form of a criterion for crack nucleation. It enjoys the desirable property of coinciding with the usual Griffith criterion to study the crack growth and with the stress criterion for the uniform traction along a straight edge. Comparisons with experiments carried out on homogeneous notched materials and on bimaterial structures show a good agreement.

Introduction

The problem of a crack onset at a notch is still an open problem if no assumption is made on a hypothetic pre-existing flaw [START_REF] Grenestedt | Crack initiation from homogeneous and bimaterial corners[END_REF][START_REF] Leblond | Crack propagation from a pre-existing flaw at a notch root -I: Introduction and general from of the stress intensity factors at the initial crack tip -II: Detailed form of the stress intensity factors at the initial crack tip and conclusions[END_REF]. When employing classical tools it leads to a paradox as will be seen below. Let us first recall briefly the energy and stress criteria that are commonly used in brittle fracture mechanics.

The energy criterion

We consider the initial state of a loaded structure to be elasto-static. The equilibrium state is characterized by a potential energy W p and a zero kinetic energy W k = 0. Next, we consider the same structure after the onset of a new crack or the growth of a preexisting one. The start point of the energy criterion is an unquestionable balance between these two states:

δW p + δW k + G c δS = 0.
(1)

Here, δW p and δW k are the changes respectively in potential and kinetic energy. The newly created crack surface is denoted δS and G c is the fracture energy per unit surface, the so-called toughness. Since the initial state is static δW k 0 and a necessary condition for fracture derives from (1):

- δW p δS G c . (2)
This incremental form of the energy criterion is the foundation of Finite Fracture Mechanics (FFM). It requires the knowledge of the crack increment surface δS. If the crack grows continuously, the above condition must hold for any small surface change δS, then considering the limit δS → 0 leads to the differential (Griffith) form of (2):

- ∂W p ∂S = G G c , (3)
1
where G is the energy release rate. Nevertheless, there are some contra-indications to the use of the differential form (3) as explained at the end of Section 1.2.

The stress criterion

The stress criterion is based on the data of a critical tension σ c (or shear τ c ), the so-called strength, that a material can bear before it breaks. The fracture of a surface occurs if:

σ σ c (or τ τ c ), (4)
where σ and τ are the tension and shear components of the stress tensor acting on the surface. Such a criterion sounds like a necessary and sufficient condition for failure. A counter example will evidence that in fact it is only a necessary one. Applied to a crack onset at a notch, these two criteria lead to the announced paradox. The Griffith criterion (3) is unable to predict such a mechanism. The energy release rate G vanishes and thus can never reach the critical toughness G c . On the other hand, the notch tip is singular, the stress field tends to infinity when approaching this point and then tension and shear stress components are always above the material strength. The stress criterion would thus conclude systematically to a crack onset whatever the applied load. The two criteria are contradictory. Moreover, neither one nor the other conclusion agree with the experiments. It can be observed that such a notch is a privileged site for crack nucleation and fracture does occur at this point but not for any small applied load.

Failure of a bar

A first example can be obtained from the simple problem of the failure of a homogeneous isotropic bar submitted to an increasing applied strain ε a . At a given load level the bar breaks in two parts and the potential energy vanishes. Then, the change in potential energy to take into account in ( 2) is simply the potential energy of the bar before failure:

-δW p = W p = 1 2 SLσ a ε a = 1 2 SL σ 2 a E , (5)
where S is the cross-section area of the bar, L its length and E is the Young's modulus. The energy criterion (2) gives:

1 2 SL σ 2 a E G c S ⇒ σ a 2EG c L . (6)
This inequality holds for a sufficiently large L leading to the contradictory conclusion that the bar breaks whatever the applied load provided it is sufficiently long! It is clear that it is the stress criterion that prevails in that case, omitting it and fulfilling only the energy condition induces an erroneous statement.

The Parvizi et al. experiments

Although the situation departs from the original topic of this paper, the experiments by [START_REF] Parvizi | Constrained cracking in glass fibre-reinforced epoxy cross-ply laminates[END_REF] proposed an answer to the above paradox that will be a key point for the present work. They carried out traction tests on a cross-ply glass fiber reinforced composite made of two outer layers oriented at 0 • (fibers are parallel to the traction direction) and a 90 • inner one. As first damage, they observed that transverse cracks appear in the inner ply and plotted the applied strain vs. the inner ply thickness at the first occurrence of a transverse crack (Fig. 1).

As long as the inner ply keeps sufficiently thick, the first transverse crack appears for a constant applied strain. It is the stress criterion that predicts correctly this mechanism (and it will be checked that the energy criterion automatically holds true). Below a critical thickness e 0 , the applied strain must be increased in order to trigger inner ply cracking. Obviously, the stress criterion no longer prevails.

Using the shear-lag model, they showed that the change in potential energy between the initial uncracked state and a new one containing just one transverse crack is written:

-δW p = Aσ 2 a e 2 d, (7)
where σ a is the applied stress, e is the inner ply thickness, d the constant width of the specimen (plane strain elasticity) and A a scaling coefficient. The energy balance (2) gives:

Aσ 2 a e 2 d G c ed ⇒ σ a G c Ae . (8)
Condition (8) imposes that if e<e 0 , the applied stress triggering a transverse crack increases like 1/ √ e, a condition that matches with the experiments (Fig. 1). The energy criterion now governs the process while the stress criterion is trivially fulfilled. The two conditions are complementary.

Conclusion

The main conclusion to draw from these examples is that, when fracture occurs the two criteria are fulfilled simultaneously, even if one often hides the other. Both are necessary conditions and together they seem to form a sufficient one. Based on this ascertainment, a criterion for crack onset at a notch is derived in the next section using the singular stress field around the notch tip. It ensures that the two criteria hold true. The answer to the question suggested in the title is then "strength and toughness". Giving both the toughness G c and the strength σ c brings us to define a characteristic length for crack onset. The failure is assumed to be a sudden and quasi-spontaneous mechanism as proposed by [START_REF] Aveston | Theory of multiple fracture of fibrous composites[END_REF], [START_REF] Parvizi | Constrained cracking in glass fibre-reinforced epoxy cross-ply laminates[END_REF], [START_REF] Wang | Initation and growth of transvers cracks and edge delamination in composite laminates, part 1. An energy method[END_REF] and [START_REF] Nairn | The strain energy release rate of composite microcracking: a variational approach[END_REF]. Other recent works are based on neighbouring arguments: [START_REF] Hashin | Finite thermoelastic fracture criterion with application to laminate cracking analysis[END_REF], [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF], Leguillon (1999Leguillon ( -a, 1999-b)-b), Leguillon andLacroix (2000-a, 2000-b).

Crack onset at a notch

The singular stress field

In plane linear elasticity, the displacement field in the vicinity of the tip of a notch with angle ω (see Fig. 2) in a homogeneous material is given by:

U (x 1 ,x 2 ) = U (0, 0) + kr λ u(θ) +•••. ( 9
)
The origin is at the tip of the notch, x 1 ,x 2 and r, θ are respectively the Cartesian and polar coordinates. The singularity exponent is such that 1/2 λ 1. It is 1/2 for a crack (ω = 0) and it is 1 (no longer singular) for a straight edge (ω = π ). The function u(θ) is an angular shape function and k is the generalized intensity factor. The exponent λ and the function u(θ) are solutions to an eigenvalue problem [START_REF] Leguillon | Computation of Singular Solutions in Elliptic Problems and Elasticity[END_REF]. The first term in ( 9) is constant and present for consistency, it is the irrelevant rigid translation. Moreover, it is assumed that there is a single real governing singularity term, it is the case for a symmetric loading for instance. In the particular case of the crack (9) is written:

U (x 1 ,x 2 ) = U (0, 0) + k I √ ru I (θ) +•••, (10)
where k = k I is the usual mode I stress intensity factor and √ ru I (θ) the opening mode I. For the straight edge the expansion is written: where k = T is the uniform traction parallel to the edge. As a consequence of ( 9), the stress field, which behaves like r λ-1 , tends to infinity as r decreases to 0. The stress criterion is thus trivially fulfilled at the tip of such a notch.

U (x 1 ,x 2 ) = U (0, 0) + Trt(θ) +•••, (11)

The crack initiation length

The potential energy change at a crack onset in the direction θ 0 is written [START_REF] Leguillon | Calcul du taux de restitution de l'énergie au voisinage d'une singularité[END_REF]:

-δW p = k 2 K(ω, θ 0 )ℓ 2λ d +•••, ( 12
)
where ℓ is the small newly created crack length and d the width of the specimen (plane elasticity). The coefficient K(ω, θ 0 ) is a scaling term depending on the local geometry (ω) and on the direction of fracture (θ 0 ) (see [START_REF] Leguillon | Computation of Singular Solutions in Elliptic Problems and Elasticity[END_REF] and the Appendix). The condition (2) becomes:

-δW p G c ℓd ⇒ k 2 K(ω, θ 0 )ℓ 2λ-1 G c . (13)
The intensity factor k is proportional to the applied load:

k = κσ a , (14)
and ( 13) is a lower bound for the increment lengths ℓ:

ℓ 2λ-1 G c K(ω, θ 0 )κ 2 σ 2 a (2λ -1 > 0). ( 15
)
Since the applied load at onset cannot be infinitely large, the increment length ℓ cannot be infinitely small. At onset, there is a jump from 0 to ℓ which is an illustration of the FFM. Of course, this lower bound must be compatible with the asymptotic framework, it has to be small with respect to a characteristic length of the structure. However, it is essential to recall that the energy balance is only a necessary condition.

In the two previous examples, the failure is assumed to occur spontaneously, which seems reasonable since the critical traction acts uniformly all over the fracture surface. This condition will be extended to the present case to provide an upper bound for the crack extension length. The singular tension σ θ (i.e. the component σ θθ of the stress tensor) at a distance ℓ from the tip in the direction θ 0 reads:

σ θ (ℓ, θ 0 ) = kℓ λ-1 s θ (θ 0 ) +•••. ( 16
)
It is a decreasing function of ℓ. If the condition (4) holds at any point between 0 and ℓ, it becomes an upper bound for ℓ:

σ θ (ℓ, θ 0 ) σ c ⇒ ℓ 1-λ κσ a s θ (θ 0 ) σ c (1 -λ>0). ( 17
)
Once again, this bound must be small in order to be sure that higher order terms in (16) are negligible. For a small applied load σ a , (15) leads to a high lower bound while (17) defines an incompatible low upper bound, thus for a monotonically increasing load the solution is achieved when equality holds in both ( 15) and ( 17). The increment length derives from these two equalities:

ℓ 0 = G c s θ (θ 0 ) 2 K(ω, θ 0 )σ 2 c . ( 18
)
The structure embedding the micro-crack with length ℓ 0 is in equilibrium in the sense that the elastic solution is characterized by the absence of kinetic energy. However, it is highly unstable from the point of view of the growing crack. The energy release rate at the tip of the newly created crack is an increasing function of its length and moreover is still above the critical toughness G c :

G(ℓ 0 ) = 2λG c (λ > 1/2), (19)
where G(ℓ 0 ) is computed using ( 12) and considering a small increment δℓ to ℓ 0 :

G(ℓ 0 ) =-lim δℓ→0 W p (ℓ 0 + δℓ) -W p (ℓ 0 ) δℓ . (20)
Here, W p (ℓ 0 ) denotes the potential energy of a structure embedding a crack with length ℓ 0 . In a first step the crack length jumps from 0 to ℓ 0 and then grows continuously. There are two particular cases. If λ = 1 (the straight edge, ω = π ) the stress criterion does not provide any upper bound:

ℓ G c s θ (θ 0 ) 2 K(π, θ 0 )σ 2 c . ( 21
)
If λ = 1/2 (the crack tip, ω = 0), it is the energy criterion that does not impose any lower bound:

ℓ G c s θ (θ 0 ) 2 K(0,θ 0 )σ 2 c . ( 22
)
The crack increment length can be taken as small as needed and the differential approach of Griffith is permitted. Indeed, another conclusion would definitely condemn the approach.

The crack onset criterion

It is now commonly admitted that the intensity factor k is the relevant parameter to define a crack onset criterion at a notch [START_REF] Seweryn | Brittle fracture criterion for structures with sharp notches[END_REF][START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF][START_REF] Qian | An experimental investigation of failure initiation in bonded joints[END_REF][START_REF] Reedy | Comparison of butt tensile strength data with interface corner stress intensity factor prediction[END_REF], 1997, 1999). It takes the Irwin-like form:

k k c , (23)
where k c is the critical value of the intensity factor. Experiments are required for each notch angle to determine the critical value k c = k c (ω). Herein, it will be derived in terms of material toughness G c and strength σ c (or τ c ). Replacing for (18) in ( 13) or (17) leads to a condition for a crack onset in the direction θ 0 :

k G c K(ω, θ 0 ) 1-λ σ c s θ (θ 0 ) 2λ-1 . ( 24
)
The direction of fracture θ c can be determined by the minimum value of the right-hand side of (24). In a homogeneous material, if the fracture properties are isotropic, i.e. independent of the direction of fracture, G c and σ c are constant and θ c is characterized by:

K(ω, θ c ) 1-λ s θ (θ c ) 2λ-1 K(ω, θ 0 ) 1-λ s θ (θ 0 ) 2λ-1 , ∀θ 0 , 0 <θ 0 < 2π -ω. (25)
This condition coincides with the G-max branching criterion for a crack (λ = 1/2) [START_REF] Leguillon | Fracture in heterogeneous materials, weak and strong singularities[END_REF]. Once θ c is known, the eigenmode r λ u(θ) can be normalized in such a way that:

s θ (θ c ) = 1 (note that usually it is s θ (θ c ) = 1 √ 2π for a crack), (26)
and the criterion takes the simplified form:

k G c K(ω, θ c ) 1-λ σ 2λ-1 c . ( 27
)
It enjoys the following very nice property, it coincides with the Griffith criterion for a crack (λ = 1/2) and with the strength criterion for a straight edge (λ = 1) (see ( 10) and ( 11)).

For a notch in a homogeneous isotropic material under symmetric loading the fracture direction is known θ c = πω/2, thus we can define a critical value k c of k as a function of the material properties σ c and G c and of the notch angle ω (through λ and K):

k c = G c K(ω) 1-λ σ 2λ-1 c , ( 28
)
where K(ω) holds for K(ω, πω/2).

The outer wedge

It can be considered as a particular case entering in the same framework with ω>π (Fig. 3).

The first exponent to consider is λ>1 (it remains one term corresponding to λ = 1 but it is the irrelevant rigid rotation). As a consequence the stress field in the vicinity of the wedge vanishes, it is now an increasing function of the distance to the corner tip:

σ θ (r, θ 0 ) = kr λ-1 s θ (θ 0 ) → 0a s r → 0 (λ -1 > 0). ( 29
)
Thus the stress criterion can never be fulfilled and then, as expected, a crack cannot initiate from an outer wedge.

Comparison with Seweryn's formula

Seweryn (1994) proposed a similar formulation that is written with the present notations (and especially with ( 26)):

k (S) c = 2 2(1-λ) λ G c K(0) 1-λ σ 2λ-1 c = 2 2(1-λ) λ K(ω) K(0) 1-λ k c . (30)
It has been derived from the singular stress field, assuming an average stress condition (31) (similar arguments are employed in the average stress criterion of Whitney and Nuismer (Isupov and Mikhailov, 1998;[START_REF] Whitney | Stress fracture criteria for laminated composites containing stress concentrations[END_REF]) without any reference to the energy criterion (except for ω = 0, see below the calibration of d 0 ):

1 d 0 d 0 0 σ dr σ c . ( 31
)
The distance d 0 is calibrated for a crack (ω = 0) and is expressed in terms of k Ic (which is equivalent to G c through the Irwin formula):

d 0 = 4k 2 Ic σ 2 c . ( 32
)
Thus the length d 0 keeps constant for any notch angle ω, moreover it is four times larger than the distance ℓ 0 computed using ( 22) with ω = 0. c . Unit for k c is MPa m 1-λ as a consequence of (26).

ωk (D)

c from [START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF] k c from (28)k The average stress criterion ( 31) is somewhat artificial, it looks like a non local law and is much more questionable than (17). All along ℓ 0 (18) the tension is above the limit strength of the material, while it is not along d 0 (32), the necessary condition (4) is not fulfilled. Moreover, the increment length ℓ 0 changes with the notch angle, it is typically a structural effect, while Seweryn tried to identify d 0 to a material property, making many references to Novozhilov's work.

Nevertheless the two formulations give numerical results close to each other, the deviation between ( 28) and ( 30) does not exceed 10% in the example of the next section (see Tables 1 and2). [START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF] performed three-point flexure experiments on a specimen of PMMA for different notch angle ω and depth (measured by the ratio a/h, see Fig. 4). They also used a finite element simulation to identify the intensity factor k (9) of the singularity. They showed, by varying the notch depth, that it is the relevant parameter to define a crack onset criterion at the notch under the Irwin-like form: Table 1 compares the experimental critical values of Dunn et al. and the estimates derived from ( 28) and (30) using the following material data for PMMA: E = 2.3GPa,ν = 0.36, σ c = 12 MPa and G c = 394 J m -2 . Table 2 exhibits the singularity exponent λ and the scaling coefficient K(ω) and also compares the crack increment lengths ℓ 0 (18) (the present analysis) and d 0 (32) (Seweryn approach). The exponent λ and the coefficient K(ω) involved in (28) are numerically determined (using procedures described in [START_REF] Leguillon | Computation of Singular Solutions in Elliptic Problems and Elasticity[END_REF][START_REF] Leguillon | Fracture in heterogeneous materials, weak and strong singularities[END_REF], see also the Appendix) even for ω = 0 where both are analytically known. This explains the very slight discrepancy observed in line 1 of Table 1. The same eigenmode normalization ( 26) is used in both approaches, thus the values on line 1 (the crack) differs from the classical one by √ 2π (see ( 26)). Although there is not a perfect coincidence, a good agreement is found (Fig. 5) and the energy balance seems to play a predominant role even for notch openings quite different from the crack.

Experimental results

Initiation at a notch in a homogeneous body

Failure of a butt joint

The failure of a butt joint is another example that shows the efficiency of the proposed criterion. Although there is no notch, the situation is very similar at the singular points. [START_REF] Reedy | Comparison of butt tensile strength data with interface corner stress intensity factor prediction[END_REF], 1997, 1999) made many experiments on epoxy (E = 3.5GPa,ν = 0.35) joints of different thicknesses in aluminum (E = 69 GPa, ν = 0.33) and steel (E = 207 GPa, ν = 0.30) specimens. Their main goal was to bring into evidence the influence of the joint thickness e on the intensity factor of the singularity located at the ends of the interface between the substrates and the joint (see Fig. 6).

Using a stretching of the domain by 1/e allows, at the leading order, rewriting of the problem in an unbounded domain undergoing a uniform vertical traction Trt(θ) (see ( 11) and the Appendix) at infinity:

U e (x 1 ,x 2 ) = U e (ey 1 ,ey 2 ) = Cte + Te ρt(θ) + V (y 1 ,y 2 ) +•••, (34)
where y 1 ,y 2 are the stretched variables, y i = x i /e and ρ = r/e (Fig. 6). In this domain the joint thickness is now unity, V (y 1 ,y 2 ) is the solution to a well-posed problem and is singular at the two ends of the interfaces separating the substrates and the joint (see ( 9)):

V (y 1 ,y 2 )) = V (0, 0) + kρ λ u(θ) +•••. ( 35
)
Note that the well-posedness of the problem is not trivial since ρu (θ) is not admissible within the layer and through the interfaces between the layer and the substrates [START_REF] Leguillon | Mode III near and far fields for a crack lying in or along a joint[END_REF]. Table 3 Comparison with [START_REF] Reedy | Comparison of butt tensile strength data with interface corner stress intensity factor prediction[END_REF], 1997, 1999) The fictitious intensity factor k in ( 35) is independent of the applied loads, it depends only on the elastic properties of the substrates and the joint. Using (34) and the stretching allows us to write the actual intensity factor k of the singularity (expressed in the unstretched, i.e. physical, coordinates):

k = Te 1-λ k, (36)
as already proposed by [START_REF] Reedy | Comparison of butt tensile strength data with interface corner stress intensity factor prediction[END_REF] (with another method to introduce the coefficient k). The thickness e is a parameter of the experiments, they checked many different values from 0.25 mm to 2 mm. They also measured the ultimate strength of the specimen σ ult , it corresponds to the critical intensity factor T = T c = σ ult we are interested in. Thus, it is sufficient for our purpose to compute once for all the scaling coefficient k to get the critical intensity factors for different joint thicknesses and load levels:

k (RG) c = σ ult e 1-λ k. ( 37
)
Gathering informations from the three articles leads to the results of Table 3. The toughness G c and strength σ c involved in (28) are taken to be the bulk properties of the epoxy: G c = 45 J m -2 , σ c = 45 MPa, although it is observed that failure occurs at the interface between the joint and one of the substrates where the stress concentration takes place. Experimental results are scattered and the predicted values lie well within the range of measures. The effect of thermal residual stresses can also be taken into account, predictions are in good agreement with results by Qian and Akisanya (1998) on butt and scarf joints.

Conclusion

As in other models (CZM, damage), the fracture onset criterion requires two parameters, the toughness (in general a characteristic energy) and the strength (in general a peak stress). These two data are equivalent to one of them plus a characteristic length. This length derives herein from a consistency argument between the classical energy and stress criteria. It is not a material parameter but a structural one, it depends for instance on the notch opening.

Fig. 1 .

 1 Fig. 1. The Parvizi et al. (1978) results. Theoretical (1) denotes the energy criterion, theoretical (2) the stress criterion, they are compared to the experimental results.

Fig. 2 .

 2 Fig. 2. The cracked notch in a homogeneous material.

Fig. 3 .

 3 Fig. 3. The unbreakable outer wedge.

  of the intensity factor beyond which fracture occurs.

Fig. 4 .

 4 Fig. 4. The three-point bend test on the notched PMMA specimen.

Fig. 5 .

 5 Fig. 5. Comparison of the theoretical critical intensity factor values with Dunn et al. (1997) and Yosibash (2001) experimental measures.

Table 1

 1 Comparison of theoretical critical intensity factor values k c and k

	(S) c (Seweryn, 1994) with Dunn et al.

Table 2

 2 Singular exponent, scaling coefficient and fracture lengths (from (18) and (32)). Unit for K is somewhat meaningless and not specified here.

	ωλK ( ω ) (×10 2 )	ℓ 0 from (18) (µm)	d 0 from (32) (µm)
	0 • 30 • 45 • 60 • 90 • 120 • 150 • 165 • 180 •	0.50 .248 0.502 0.243 0.506 0.242 0.513 0.237 0.545 0.212 0.616 0.176 0.752 0.128 0.857 0.098 10 .079	10.34 1 .3 10.54 1 .3 10.64 1 .3 10.84 1 .3 12.14 1 .3 14.64 1 .3 20.04 1 .3 26.14 1 .3 32.44 1 .3

  results on the critical intensity factor value for butt joint failure. Here k c is expressed in MPa mm 1-λ .

		λk	(RG) c	from (37)	k c from (28)
	Steel	0.70	11.4 k c 15.0	14.3
	Aluminum	0.73	14.0 k c 19.6	16.2

The main assumption made to get these results lies in the truncated expansion (9) of the elastic solution in terms of singular modes in the vicinity of the notch tip. The omitted part is supposed to be negligible. In the case of the butt joint there is a single singular exponent, the others are larger than 1 and then the assumption holds true. It is also true for a notch in a homogeneous body under symmetric loading. The first mode (retained) is symmetric while the second one (omitted) is anti-symmetric and then not relevant (the corresponding intensity factor vanishes). Moreover, in that last case, the results extend obviously to a pure anti-symmetric loading considering the shear strength.For a notch in a bimaterial wedge, the above arguments based on symmetries fail. Two singular terms must in general be accounted for, they are characterized by a pair of generalized intensity factors k 1 and k 2 (or even a complex one). The consistency between the two fracture conditions can no longer be explicitly solved. Nevertheless, a new expression of the criterion can be derived, it involves a mixity parameter evaluating the mixture of singular modes[START_REF] Leguillon | Finite fracture mechanics -Application to the onset of a crack at a bimaterial corner[END_REF]. Results are compared to some experiments worked out byMohammed and Liechti (Hutchinson, 2000;[START_REF] Mohammed | Cohesive zone modeling of crack nucleation at bimaterial corners[END_REF] on a notched bimaterial.In this respect, the only case that is easily accessible concerns a multiple singular eigenvalue. The definition of the mode mixity is no longer troublesome. It is for instance trivially the case for a crack in a homogeneous body (multiplicity of λ = 1/2 equals 2 in plane elasticity) but it applies also to a crack impinging orthogonally on an interface[START_REF] He | Crack deflection at an interface between dissimilar elastic materials[END_REF] Leguillon and Lacroix, 2000-a, 2000-b) where the double root λ is no longer 1/2. It will be the topic of a forthcoming paper dealing with new conditions for crack deflection by an interface derived from the above criterion.

Appendix

This appendix is dedicated to a brief recall of the matched asymptotics procedure that leads to some of the expressions used within the paper.

In plane linear elasticity, we consider a domain in which a corner is slightly perturbed by a flaw (a short crack or a small cavity for instance, Fig. 7(a)). The dimensionless diameter of this flaw is denoted ε. The solution U ε (x 1 ,x 2 ) to an elasticity problem in this domain can be expressed as the unperturbed solution U 0 (x 1 ,x 2 ) defined in 0 ( 0 is the limit of ε as ε → 0, Fig. 7(b)) plus a small correction:

where f 1 (ε) → 0asε → 0. Such an expansion (38) a so-called 'outer', is valid in the whole domain 0 (or ε ) except near the corner point where the geometry is perturbed. The solution U 0 (x 1 ,x 2 ) is singular at the corner point and can be expanded as:

See ( 9) and Fig. 2 for the meaning of the different terms. For simplicity, it is assumed that the leading singular term is a single real one.

In order to have a description of the near fields, the domain ε is stretched (×1/ε)andasε → 0 it leads to the unbounded 'inner' domain in spanned by the stretched variables y 1 = x 1 /ε and y 2 = x 2 /ε. The perturbation diameter is now 1. The solution can be expanded in this domain as:

where

It is the 'inner' expansion. Conditions at infinity are missing to define well-posed problems for the unknown functions V 0 (y 1 ,y 2 ) and V 1 (y 1 ,y 2 ). They derive from the matching conditions based on the existence of an intermediate area where both expansions ( 38) and ( 40) hold. In other words, the behaviour of the outer terms (in ( 38)) when approaching the singular point must match with the behaviour of the inner terms (in ( 40)) at infinity. The asymptotic behaviour of U 0 (x 1 ,x 2 ) being defined by ( 39), it leads to:

In addition the behaviour of V 1 (y 1 ,y 2 ) at infinity is prescribed: The function V 1 (y 1 ,y 2 ) is independent of the applied loads, it depends only on the local geometry and on the shape of the singular mode defining the condition at infinity (42). Using a superposition principle:

allows us to conclude that the second term of the inner expansion ( 40) is well-defined.

Let us now consider the change in potential energy δW p between the solutions in the unperturbed U 0 (x 1 ,x 2 ) and perturbed U ε (x 1 ,x 2 ) situations for unchanged boundary conditions (unchanged applied loads). Betti's theorem gives:

where Ŵ is any contour surrounding the corner and n its normal pointing toward the origin. The integral can be taken either in 0 (or ε )orin in . Selecting the inner domain and substituting the asymptotics ( 40)-( 43) in ( 44) leads finally to:

with (here d = 1, see ( 12)):

Clearly K depends on the shape of the perturbation through V 1 but not on its actual size since it has been stretched to 1. Moreover it is independent of the applied loads that appear in the inner expansion ( 40)-( 41) only through the coefficient k.The function V 1 can be computed once for all by finite elements and K can be derived using the contour-independent integral (46).


Comment on the article "Probability of radial anisotropy in the deep mantle" by [START_REF] Visser | Probability of radial anisotropy in the deep mantle[END_REF] Based on the inversion of a large surface waves dataset, including both fundamental and higher modes, [START_REF] Visser | Probability of radial anisotropy in the deep mantle[END_REF] constructed new maps of the probability of surface wave polarization (or radial) anisotropy. These maps bring new constraints on the variation of radial anisotropy with depth as a function of the geodynamic environment. However, the interpretation of the seismic anisotropy distribution in terms of mantle deformation is oversimplified, since it is essentially based on data of the evolution of olivine crystal preferred orientations and seismic anisotropy at high temperature, but low pressure conditions, neglecting recent findings on the deformation of mantle minerals under high-pressure in-situ mantle conditions.

Three points, in particular, deserve a throughout discussion, since they highlight the limitations of using basic concepts that have been established from the analysis of shallow mantle rocks and low pressure deformation experiments, as for instance: (i) lack of anisotropy implies deformation by diffusion creep or (ii) fast polarization (or propagation for P or SV waves) directions mark systematically the flow direction, to interpret seismic anisotropy data below depths of 200-250 km in the mantle.

Absence of anisotropy does not necessarily imply diffusion creep

Estimation of seismic anisotropy using polycrystal plasticity models based on recent experimental data on the deformation of olivine, wadsleyite, and perovskite highlights that development of a crystal preferred orientation does not necessarily result in significant seismic anisotropy, even if the single crystal has strongly anisotropic elastic constants [START_REF] Tommasi | Strain-induced seismic anisotropy of wadsleyite polycrystals and flow patterns in the mantle transition zone[END_REF][START_REF] Mainprice | Pressure sensitivity of olivine slip systems and seismic anisotropy of the Earth's upper mantle[END_REF][START_REF] Mainprice | Predicted slip systems and crystal preferred orientation of polycrystalline silicate Mg-Perovskite at mantle pressures: Implications for the seismic anisotropy of the lower mantle[END_REF].

Strong anisotropy at the rock scale is produced when the crystal preferred orientations have strong single maxima concentrations of the crystallographic directions corresponding to the highest and lowest propagation velocities or S-wave splitting of the single crystal. This is the case, for instance, for olivine under shallow mantle conditions, where dislocation creep with dominant activation of the (010)[100] slip system leads to a strong orientation of the [100] axis, which is the fastest P-wave propagation and the lowest S-wave splitting direction, and of the [010] axis, which is the slowest P-wave propagation direction [START_REF] Tommasi | Viscoplastic self-consistent and equilibrium-based modeling of olivine lattice preferred orientations. Implications for upper mantle seismic anisotropy[END_REF]. In contrast, dominant activation of {hk0}[001] systems in the deep upper mantle results in a point concentration of the [001] axis that has intermediate values for both P-wave propagation velocities and S-wave splitting and in girdle distributions of the [100] and [010] axes in a plane normal to the flow direction and hence in averaging between high and low values [START_REF] Mainprice | Pressure sensitivity of olivine slip systems and seismic anisotropy of the Earth's upper mantle[END_REF]. Similarly, for wadsleyite, seismic velocity and anisotropy distributions of a polycrystalline aggregate do not correlate in a simple manner with the single-crystal properties [START_REF] Tommasi | Strain-induced seismic anisotropy of wadsleyite polycrystals and flow patterns in the mantle transition zone[END_REF]. In the single crystal, P-wave propagation is fastest parallel to the [010] axis and slowest parallel to [001]. The maximum S-wave polarization anisotropy is observed for propagation parallel to b110N, the fast S-wave being polarized parallel to the [010] axis. In wadsleyite-rich polycrystals, the slowest P-wave propagation is still parallel to the [001] maximum, but the fastest propagation direction and the polarization of the fast S-wave are not parallel to the [010] maximum, but at low angle to the [100] maximum, i.e., at low angle to the shear direction. The dispersion of [010] axes in the polycrystal results in a weak seismic anisotropy for wadsleyite-rich rocks (N 2%), even if the single-crystal is clearly anisotropic (∼ 13%, [START_REF] Mainprice | The seismic anisotropy of the Earth's mantle: from single crystal to polycrystal[END_REF] and the rock has a well-developed wadsleyite CPO. Finally, in the lower mantle, coupling of ab initio and polycrystal plasticity models suggests that even if deformation by dislocation creep produces clear perovskite CPO, the resulting seismic anisotropy is weak [START_REF] Mainprice | Predicted slip systems and crystal preferred orientation of polycrystalline silicate Mg-Perovskite at mantle pressures: Implications for the seismic anisotropy of the lower mantle[END_REF]. The maximum polarization anisotropy is N 2% at the low temperatures and pressures that prevail just below DOI of original article: 10.1016/j.epsl.2008.03.041. ⁎ Corresponding author.

E-mail address: deia@dstu.univ-montp2.fr (A. Tommasi).

1

the transition zone and it decreases with increasing temperature and pressure.

2. Variations in sign of radial anisotropy (from V SH N V SV to V SV N V SH ) with increasing depth cannot be univocally interpreted as a change in flow direction

This interpretation is based on the premise, derived from experiments at high temperature and moderate pressure and extensive data on naturally deformed shallow upper mantle rocks, that olivine deforms by dominant slip on {0kl}[100] slip systems. Under these conditions, the deformation-induced CPO results in the fastest S-waves being polarized parallel to the flow direction. It also leads to faster propagation of P and Rayleigh waves parallel to the flow direction. However, if the relation between finite strain axes and olivine CPO changes, this relation is no longer valid. Recent experimental data show that at the high pressures that prevail in the deep upper mantle [START_REF] Couvy | Shear deformation experiments of forsterite at 11 GPa -1400 °C in the multianvil apparatus[END_REF][START_REF] Raterron | Pressure induced slip-system transition in forsterite: single-crystal rheological properties at mantle pressure and temperature[END_REF]o ra t high water contents [START_REF] Jung | Water-induced fabric transitions on olivine[END_REF] olivine deforms by dominant glide on {hk0}[001] systems. There is still debate on the relative importance of these two parameters, but they are not independent, since water solubility in olivine increases with pressure [START_REF] Bolfan-Casanova | Water partitioning between nominally anhydrous minerals in the MgO-SiO2-H2O system up to 24 GPa: implications for the distribution of water in the Earth's mantle[END_REF]. Moreover, both pressure and water have the same effect, resulting in the same olivine CPO evolution and hence in fast polarization directions normal to the flow direction. Polycrystal plasticity models based on these data show indeed that the change in intensity and sign (V SH b V SV ) of anisotropy with increasing depth in the upper mantle observed in previous studies [START_REF] Montagner | How to reconcile body-wave and normal-mode reference Earth models[END_REF][START_REF] Panning | A three-dimensional radially anisotropic model of shear velocity in the whole mantle[END_REF]Zhou et al., 2006) and confirmed by [START_REF] Visser | Probability of radial anisotropy in the deep mantle[END_REF] is perfectly explained as resulting from the change in olivine deformation with depth from dominant [100] slip to dominant [001] slip [START_REF] Mainprice | Pressure sensitivity of olivine slip systems and seismic anisotropy of the Earth's upper mantle[END_REF].

Anisotropy in the transition zone cannot be interpreted based on olivine deformation data

It is essential to consider the phases that are actually present at this depth range: wadsleyite or ringwoodite and pyrope-rich garnet. New high pressure, high temperature experimental data show that garnet [START_REF] Li | Plastic flow of pyrope at mantle pressure and temperature[END_REF], wadsleyite (Thurel and Cordier, 2003;Thurel et al., 2003a;[START_REF] Wenk | In situ observation of texture development in olivine, ringwoodite, magnesiowus-tite and silicate perovskite at high pressure[END_REF], and ringwoodite [START_REF] Wenk | In situ observation of texture development in olivine, ringwoodite, magnesiowus-tite and silicate perovskite at high pressure[END_REF] deform by dislocation creep under high-pressure conditions. TEM analysis of experimentally deformed samples [START_REF] Karato | Plastic deformation of silicate spinel under the transition zone conditions of the Earth's mantle[END_REF]Thurel et al., 2003b) and ab initio modelling of dislocation plasticity [START_REF] Carrez | Slip systems and plastic shear anisotropy in Mg2SiO4 ringwoodite: insights from numerical modelling[END_REF] allow constraining the active slip systems and their relative strength. Based on these data, one may use polycrystal plasticity models to simulate the evolution of crystal preferred orientations in garnet [START_REF] Mainprice | Crystal preferred orientations of garnet: comparison between numerical simulations and electron back-scattered diffraction (EBSD) measurements in naturally deformed eclogites[END_REF], wadsleyite [START_REF] Tommasi | Strain-induced seismic anisotropy of wadsleyite polycrystals and flow patterns in the mantle transition zone[END_REF], and ringwoodite polycrystals [START_REF] Carrez | Slip systems and plastic shear anisotropy in Mg2SiO4 ringwoodite: insights from numerical modelling[END_REF] during plastic flow and estimate the resulting anisotropy. For the upper part of the transition zone, which is essentially composed of wadsleyite and garnet, these models predict a weak seismic anisotropy for a polycrystal of pyrolitic composition (60% wadsleyite, 40% garnet) at transition zone conditions: ∼2% for P and ∼ 1% for S-waves for a shear strain of 1 with fastest propagation of P-waves at low angle to the shear direction and the maximum shear wave splitting in the shear plane, perpendicular to the shear direction [START_REF] Tommasi | Strain-induced seismic anisotropy of wadsleyite polycrystals and flow patterns in the mantle transition zone[END_REF]. Swave anisotropy is characterized by faster propagation of waves polarized at low angle to the shear direction. Horizontal shearing results therefore in higher velocities for horizontally-propagating Pwaves (PH) and horizontally-polarized S-waves (SH), as well as in weak azimuthal variation of SV and SH velocities. On the other hand, vertical flow leads to higher velocities for vertically-propagating Pwaves (PV) and vertically-polarized S-waves (SV) and to a weak azimuthal variation of SV velocity, but to a roughly constant SH velocity. In the light of these models, the higher probability of V-SV N V SH in the transition zone below convergent or divergent plate boundaries observed by [START_REF] Visser | Probability of radial anisotropy in the deep mantle[END_REF], [START_REF] Beghein | Probability density functions for radial anisotropy: implications for the upper 1200 km of the mantle[END_REF], and [START_REF] Panning | A three-dimensional radially anisotropic model of shear velocity in the whole mantle[END_REF] may imply dominant vertical flow beneath plate boundaries, whereas the global azimuthal anisotropy patterns inferred from Love-wave overtone data by [START_REF] Trampert | Global azimuthal anisotropy in the transition zone[END_REF] are better explained by horizontal shearing, which may predominate in most of the transition zone, away from plate boundaries. Similar models using lower transition zone, ringwoodite-rich compositions show that anisotropy below 520 km depth cannot be explained by strain-induced CPO of ringwoodite [START_REF] Carrez | Slip systems and plastic shear anisotropy in Mg2SiO4 ringwoodite: insights from numerical modelling[END_REF]. Thus, if, as suggested by [START_REF] Visser | Probability of radial anisotropy in the deep mantle[END_REF] probability maps, seismic anisotropy is higher in the lower transition zone, another mechanism, such as alternating layers with highly contrasting seismic properties has be invoked explain the observations. Such a layering may, for instance, be associated with stagnant hydrated slabs in the transition zone as proposed by [START_REF] Ritsema | Global transition zone tomography[END_REF] and [START_REF] Mainprice | Predicted elastic properties of the hydrous D phase at mantle pressures: Implications for the anisotropy of subducted slabs near 670-km discontinuity and in the lower mantle[END_REF].

  EPSL 270:241-250 Andréa Tommasi ⁎, David Mainprice Geosciences Montpellier, CNRS & Université Montpellier 2, CC.060, Pl. Eugène Bataillon, Montpellier, F-34095, France


Nicolas Lupoglazoff

Recent progress in numerical simulations for jet noise computation using LES on fully unstructured meshes

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

Recent progress in numerical simulations for jet noise computation using LES on fully unstructured meshes

Nicolas Lupoglazoff [START_REF] Lupoglazoff | Comprehensive 3D unsteady simulations of subsonic and supersonic hot jet flow-fields. Part 1: aerodynamic analysis[END_REF] and François Vuillot2 Onera -The French Aerospace Lab, F-92322 Châtillon, France Numerical simulations of a hot turbulent jet are performed on full unstructured grids with the aim to compute the noise generated in the far field. The jet flow is computed with the Onera CEDRE code using the LES (Large Eddy Simulation) approach and provide the unsteady inputs to the Onera KIM acoustic solver, based on the integral surface FW-H (Ffowcs Williams and Hawkings) approach which is used to reconstruct the acoustic far field. Several unstructured grids of increasing resolution are considered and the results obtained for both the computed jet flow and far field acoustics are compared to measurements carried out in the Onera CEPRA19 anechoic wind tunnel and to previous computations with the same codes but based on structured grids. Results show that unstructured grids allow a better grid tailoring work with an evident ease of construction, compared to the structured grid approach. Very good results are obtained and compare well with available measurements with an unstructured grid of roughly 130 millions cells. These results open the way to the use of unstructured grids in more complex geometries such as those encountered in real life configurations dealing with noise reduction devices or installation effects.

Nomenclature

I. Introduction

ET noise contributes largely to the nuisance of aircraft in take-off situations. In order to have the means to control and reduce jet noise, it is desirable to develop reliable and easy to use simulation methods for noise prediction. In that regard, the use of unstructured meshes is seen as an important step to improve on the difficulty to construct structured meshes for real life nozzle geometries and aircraft environment.

This article describes the aeroacoustic computations carried out during 2013 and 2014 on the hot jet issued from Onera Φ80 nozzle at Mach 0.7 with fully unstructured, predominantly tetrahedral meshes. These computations follow many works on the jet noise, e.g. see references [START_REF] Lupoglazoff | Comprehensive 3D unsteady simulations of subsonic and supersonic hot jet flow-fields. Part 1: aerodynamic analysis[END_REF][START_REF] Biancherin | Comprehensive 3D unsteady simulations of subsonic and supersonic hot jet flow-fields. Part 2: acoustic analysis[END_REF][START_REF] Lupoglazoff | Application of the CEDRE unstructured flow solver to jet noise computations[END_REF][START_REF] Muller | Experimental and Numerical Investigation of the Near Field Pressure of a High Subsonic Hot Jet[END_REF][START_REF] Ahuja | Tone-excited jet: theory and experiments[END_REF][START_REF] Bridges | Establishing consensus turbulence statistics for hot subsonic jets[END_REF][START_REF] Piot | Investigation of the PSE approach for subsonic and supersonic hot jets. Detailed comparisons with LES and Linearized Euler Equations results[END_REF][START_REF] Bogey | Large-eddy simulation of the flow and acoustic fields of a Reynolds number 10 5 subsonic jet with tripped exit boundary layers[END_REF][START_REF] Huet | Influence of boundary layers resolution on heated, subsonic, high Reynolds number jet flow and noise[END_REF][START_REF] Lorteau | Experimental analysis of the near field pressure of a single stream subsonic jet[END_REF][START_REF] Prieur | Aeroacoustic integral methods and efficient numerical implementation[END_REF][START_REF] Rahier | Investigation of integral surface formulations for acoustic post-processing of unsteady aerodynamic jet simulations[END_REF][START_REF] Labbé | A CFD/CAA coupling method applied to jet noise prediction[END_REF][START_REF] Lorteau | Recent progress in LES computation for aeroacoustics of turbulent hot jet. Comparison to experiments and near field analysis[END_REF][START_REF] Bogey | Influence of initial turbulence level on the flow and sound fields of a subsonic jet at a diameter-based Reynolds number of 10 5[END_REF] . More specifically, the present works aim at evaluating the recent versions of the Onera CEDRE code on predominantly tetrahedral meshes on a test case which has been previously computed with the same code on hexahedral and hybrid meshes [START_REF] Huet | Influence of boundary layers resolution on heated, subsonic, high Reynolds number jet flow and noise[END_REF][START_REF] Lorteau | Recent progress in LES computation for aeroacoustics of turbulent hot jet. Comparison to experiments and near field analysis[END_REF] . The present results will contribute to revisiting the first results obtained on tetrahedral meshes with earlier versions of the CEDRE code, which concluded on the need to insert structured patches into the unstructured mesh [START_REF] Lupoglazoff | Application of the CEDRE unstructured flow solver to jet noise computations[END_REF] . Avoiding the need of structured patches will greatly ease the mesh generation work in complex configurations. The main advantage is that tetrahedral meshes can be easily generated with most grid generation softwares (present grids have been generated with CENTAUR 16 ) in any complex geometries.

It has been demonstrated [START_REF] Haider | Stability of the cell centered finite-volume MUSCL method on unstructured grids[END_REF] that tetrahedral meshes need some form of gradient limitations to insure stability. Earlier versions of gradient limiters used in CEDRE turned out to be too dissipative to properly capture flow instabilities on tetrahedral grids. The overall methodology was then adapted to include hexahedral zones in the mesh, at least in the vicinity of the jet mixing layers, and it provided satisfactory results [START_REF] Lupoglazoff | Application of the CEDRE unstructured flow solver to jet noise computations[END_REF] . However, for the hybrid grids, the inclusion of the structured zone is disadvantageous in terms of number of stitches (from hexahedra zone to tetrahedra zone), which produce unnecessary small cells, and requires greater sophistication in the construction of the mesh. The other major drawback is that the structured zones cannot produce isotropic grids without strong penalties on the surrounding regions, in terms of number of cells, while in the early jet development region isotropy is desirable. On the contrary this is easily achievable with tetrahedral meshes.

Predicting jet noise was first tested and validated on hexahedral structured meshes, for single and double stream nozzles. The methodology that we are using is the hybrid approach from previous computations at Onera 1-4, 7, 9, 14 which associates a LES unsteady flow solution, using the Smagorinsky subgrid scale model, and an acoustic code for far field reconstruction. The LES solution is computed with the Onera CEDRE code and is stored on particular surfaces surrounding the jet. The stored solution on the surfaces is then fed to the Onera KIM acoustic code, developed by Rahier et al. [START_REF] Prieur | Aeroacoustic integral methods and efficient numerical implementation[END_REF][START_REF] Rahier | Investigation of integral surface formulations for acoustic post-processing of unsteady aerodynamic jet simulations[END_REF] that relies on several surface integration methods, including the Ffowcs-Williams & Hawkings analogy (FW-H) that is used in the present work.

II. Simulation parameters A. Jet conditions

In this paper we are considering the Onera PHI80 nozzle, an academic single stream nozzle used for several years to validate numerical simulations [START_REF] Lupoglazoff | Comprehensive 3D unsteady simulations of subsonic and supersonic hot jet flow-fields. Part 1: aerodynamic analysis[END_REF][START_REF] Biancherin | Comprehensive 3D unsteady simulations of subsonic and supersonic hot jet flow-fields. Part 2: acoustic analysis[END_REF][START_REF] Muller | Experimental and Numerical Investigation of the Near Field Pressure of a High Subsonic Hot Jet[END_REF][START_REF] Piot | Investigation of the PSE approach for subsonic and supersonic hot jets. Detailed comparisons with LES and Linearized Euler Equations results[END_REF][START_REF] Huet | Influence of boundary layers resolution on heated, subsonic, high Reynolds number jet flow and noise[END_REF][START_REF] Lorteau | Experimental analysis of the near field pressure of a single stream subsonic jet[END_REF][START_REF] Labbé | A CFD/CAA coupling method applied to jet noise prediction[END_REF][START_REF] Lorteau | Recent progress in LES computation for aeroacoustics of turbulent hot jet. Comparison to experiments and near field analysis[END_REF] . It has a 80mm exit diameter and has been tested in the Onera CEPRA19 anechoic wind tunnel for hot or isothermal conditions. Figure 1 illustrates the nozzle geometry.

The hot conditions are the one used in the present work and are summarized in table 1.

The ambient pressure is set to 101 325 Pa and the nozzle stagnation temperature is 900 K Finally, the free field velocity is set to Ua = 5 m/s. The resolution of the Navier-Stokes equations is made with a finite volume approach for the conservative variables on generalized unstructured meshes. The simulations were performed following two steps, the first step is performed with higher time steps and is intended to initialize the jet flow across the computational domain, starting from rest conditions. The second step is performed with a reduced time step and is intended to provide unsteady solutions for averaging mean flow, computing flow statistic and perform acoustic computations from solutions stored during this second step. The time steps for the second computational phase are chosen so that, depending on grid resolution, the CFL values is less than 1 in the jet near field and of the order of unity or slightly larger in the vicinity of the nozzle lips and shear layer regions. This second step represents a physical time in excess of 50 ms, that is to say about 250 D/Uj convective time units.

Inside the nozzle, uniform stagnation pressure and temperature profiles are imposed, thus the boundary layers develop freely. Outside the nozzle, static pressure pa = 101 325 Pa is imposed at the outflow conditions. On the lateral and upstream boundaries, static temperature Ta = 280 K and velocity Ua = 5 m/s conditions are imposed. The grid is stretched from a refined zone, in which the jet flow is calculated, to the simulation domain boundaries in order to damp acoustic waves before they reach the borders and thus avoid spurious reflections. The nozzle walls are assumed to be adiabatic.

The computation was performed at the experimental Reynolds number using the Smagorinsky subgrid scale model (Cs = 0.1). Due to the high temperature of the jet flow we decided to rely on thermodynamics properties which are allowed to vary with the temperature.

C. Grid definition

Great care was exercised in defining the unstructured grids used in this work. The idea is to get rid of the structured grid approach or the mixed grid approach which was used in the preceding years and which consisted of imbedding structured patches into an overall unstructured grid. Although the structured grid approaches (either fully structured or mixed structured/unstructured) did work well it implies complex grid construction work that is not well suited to complex configurations as those encountered in practice. Since the recent versions of the code propose advanced gradient limitation procedures it was decided to evaluate them.

Consequently, a consistent unstructured grid hierarchy was constructed. Grids were created with the Centaursoft grid generation tools 16 for a cylindrical computational domain of radius 80D and extending from -29D to 100D in the axial direction (with origin chosen at the nozzle exit plane). These grids are composed of tetrahedra, completed by prisms along nozzle walls. Prisms are intended to produce some degree of mesh anisotropy in the boundary layers, although these are not fully resolved at this stage. The meshes are referenced as 2GG for a coarse version, and 2GM and 3GM for two versions of the medium size meshes. The main properties of these grids are summarized in table 2 below. The 2GG is illustrated on Figure 2. For these meshes, the refined region is a truncated cone of length 25 D whose radius varies from 2 D (x / D = 0) to 4D (x / D = 25). The largest cell size at the edge of this refined zone leads to an estimated cut-off frequency for the finest 3GM grid in excess of 4000 Hz, based on a 20 point per wave length criterion. A further refinement involves the wake region of the nozzle lip, where the mixing layer develops, and extends over 7.5 D. An additional buffer zone, extending from 25 D to 35 D is used to progressively coarsen the mesh before reaching the computational domain boundaries to avoid spurious wave reflection.

Except for the prismatic zones along the nozzle walls, the tetrahedron grids are mostly isotropic. This is particularly important in regions close to the nozzle exit. Table 2 provides the main characteristics of the grids in terms of spatial resolution. In particular, the grid resolution, , at the nozzle exit can be easily set in the grid generation software and values of /D range from 1.25% for the 2GG grid to 0.375% for the 3GM grid. This results in an azimuthal resolution with cell count ranging from 251 to 838 around the nozzle lip. Recent works [START_REF] Bogey | Large-eddy simulation of the flow and acoustic fields of a Reynolds number 10 5 subsonic jet with tripped exit boundary layers[END_REF][START_REF] Lorteau | Recent progress in LES computation for aeroacoustics of turbulent hot jet. Comparison to experiments and near field analysis[END_REF] have emphasized the importance of providing sufficient azimuthal resolution to properly resolve the early development of the jet shear layer. In this region, the resolution of the 3GM unstructured grid is almost twice as much as the one of the finest structured grids that relied on 480 azimuthal planes.

D. Acoustic far field computations

The far field noise radiation is computed using the Ffowcs Williams & Hawkings [START_REF] Williams | Sound generation by turbulence and surfaces in arbitrary motion[END_REF] (FW-H) porous surface formulation available in the KIM code developed at Onera, see Rahier et al. [START_REF] Rahier | Investigation of integral surface formulations for acoustic post-processing of unsteady aerodynamic jet simulations[END_REF] for more details. This formulation allows to compute time pressure histories at any observer location by integration of the flow fields on a control surface surrounding the jet and containing all the noise sources. It is preferred to the Kirchhoff method that was shown by Rahier et al. to be more exposed to the generation of spurious noise, especially for hot jets. In the present work, several control surfaces have been used. In particular, as discussed hereafter, the initial surfaces turned out to be too short and were completed by further extensions in the downstream direction. No closure of the surfaces has been used at both extremities to avoid pressure signals contamination by spurious noise generated by the crossing of the surface by turbulent spots. Moreover Rahier et al. showed that the closure does not affect much the results for a sufficiently long surface provided that the turbulence is rather weak on the downstream closing disc. In the present simulation, the far field microphones are located at 75D from the nozzle exit for the directions between 20 deg. and 150 deg., as in the experiments.

III. Flow solution

First the solutions obtained for the three unstructured grids are analyzed in term of jet flow characteristics. In particular, the benefits of a finer grid resolution are clearly illustrated. It must be stressed that for the present computations and contrary to recent works by Lorteau et al. [START_REF] Lorteau | Recent progress in LES computation for aeroacoustics of turbulent hot jet. Comparison to experiments and near field analysis[END_REF] , no turbulent triggering is used. The importance of turbulence triggering is still an open question and the present work will give attention to this aspect of the simulations.

A. Jet development

The computed jet development is visible, for the two extreme grids (2GG and 3GM), on figures 3-5. Note that large turbulent scales are less present for the 3GM grid than for the 2GG grid where they result in a faster and more abrupt dispersion of the jet, absent from the 3GM solution. This is especially visible for the axial velocities, figures 3 and 5. As mentioned above, due to the high temperature of the jet flow we decided to rely on variable thermodynamic properties as illustrated on figure 4 which presents the ratio γ of specific heat capacities. On this figure, the effective CFL number is also reported. For the second phase of simulation, the time step is set to 1 s for the 2GG, 0.5 s for the 2GM, and 0.3 s for the 3GM. The CFL values turned out to be less than unity almost everywhere and in the order of unity in the vicinity of the nozzle lips.

The mean and rms axial velocity profiles are shown in figures 6, 7 and compared to experimental measurements of references [START_REF] Muller | Experimental and Numerical Investigation of the Near Field Pressure of a High Subsonic Hot Jet[END_REF][START_REF] Ahuja | Tone-excited jet: theory and experiments[END_REF] . The improvement as the grid is refined is clearly visible with 3GM achieving very good comparisons with the measurements. These results also compare well with previous structured computations on the same jet, with the same code [START_REF] Huet | Influence of boundary layers resolution on heated, subsonic, high Reynolds number jet flow and noise[END_REF][START_REF] Lorteau | Recent progress in LES computation for aeroacoustics of turbulent hot jet. Comparison to experiments and near field analysis[END_REF] . Moreover, the 3GM results are found very comparable to those recently obtained with CEDRE by Lorteau et al. on a refined structured mesh of 240 million cells [START_REF] Lorteau | Recent progress in LES computation for aeroacoustics of turbulent hot jet. Comparison to experiments and near field analysis[END_REF] , see figure 11. Both computations produce very good comparisons with available measurements while the 3GM computation presents a clear advantage in term of mesh generation work and in term of CPU cost. For 3GM each job lasted an average of 13 hours CPU on 480 cores on ONERA cluster. It took 12 jobs to pass the transitional first step and 24 jobs for the second step to compute mean field and acoustics. Compared to the reference computation of Lorteau et al. the CPU cost, for computing 1 ms of flow, is found to be reduced by a factor of 3.5, when a second order implicit time scheme is selected.

2GG 3GM Figure 5. Comparison 2GG / 3GM axial velocity iso surfaces

The present results show a consistent trend both in terms of comparison with previous structured grid results and in term of grid dependency for the present tetrahedral grid family. This permits to validate the use of tetrahedral meshes for jet flow simulations. It is interesting to note that the progressive mesh refinement results in a constant improvement of the results with a very good comparison obtained for the finest mesh, 3GM, compared to available measurements for the jet aerodynamics, with a proper potential core length estimation and turbulence development. Indeed, the potential core length is estimated at Lc = 6.1 D, for both experiment and 3GM computation (with Lc being taken at Ux/Uj = 0.90 on the jet axis). Examining in more details the fluctuation of axial velocity, the finest 3GM mesh leads to a less abrupt turbulence development, which is in good agreement with the measurements (Figure 6).

Figure 6. Axial profiles of mean and rms axial velocity.

Green: 2GG, orange: 2GM, red: 3GM -symbols: experiments, square Onera/Φ80, triangles Ahuja et al., [START_REF] Ahuja | Tone-excited jet: theory and experiments[END_REF] .

Figure 7. Radial profiles of mean axial velocity.

Green: 2GG, orange: 2GM, red: 3GM -symbol: experiments, Onera/Φ80.

B. Nozzle exit flow

As recent works illustrated the importance of the early flow development, just at the nozzle exit, some attention was given to the computed solution in the nozzle exit region. As mentioned above, no turbulence triggering was used in the present work while Lorteau et al. [START_REF] Lorteau | Recent progress in LES computation for aeroacoustics of turbulent hot jet. Comparison to experiments and near field analysis[END_REF] showed that turbulence triggering is an important aspect for jet flow solution.

Figures 8-10 illustrate the importance of grid resolution on the early flow development just at the nozzle exit. From these figures it is apparent that coarser grids exhibit a retarded transition to turbulence, with large scale vortices which are generated at the nozzle lip, and lead to marked vortex pairing further downstream. On the contrary the finest resolution 3GM grid produces finer structures which are less prone to vortex pairing. As vortex pairing is seen as an additional source of acoustic emission to the far field, the benefit of an increased resolution grid, such as the 3GM grid, is clearly visible. This benefit is clearly apparent on figures 8 and 9, for both the vorticity and pressure fields. When looking at the normalized velocity fluctuations along the nozzle lip line it is clear that finest 3GM grid produces much earlier turbulence development which peaks at 20% at x/D=0.5 and then levels at 16% further downstream (see figure 10) which is in agreement with observed behaviors. Coarser grids 2GG and 2GM exhibit a retarded peak which is at a higher level, close to 20%. Green: 2GG, orange: 2GM, red: 3GM.

Further comparisons were made with the recent computations of Lorteau et al. [START_REF] Lorteau | Recent progress in LES computation for aeroacoustics of turbulent hot jet. Comparison to experiments and near field analysis[END_REF] on a structured grid of 240 million cells with turbulence tripping thanks to a small step inserted inside the nozzle. Figure 11 shows comparisons for the axial profiles of mean and rms axial velocity on the jet axis. This figure shows that very good agreement is obtained between these two computations, performed on very different grids. This demonstrates a satisfactory grid independence of the jet flow solution. To go further in this analysis, figure 12 shows the PSD (Power Spectral Densities) of the radial fluctuating velocity recorded at different axial locations along the nozzle lip line. The first plot is obtained from a numerical transducer located inside the nozzle (x/D=-0.25), downstream of the small step used in the work of Lorteau et al. [START_REF] Lorteau | Recent progress in LES computation for aeroacoustics of turbulent hot jet. Comparison to experiments and near field analysis[END_REF] . On that plot the effect of the step, absent from the present computations, is clearly visible as it produces a marked turbulent behavior which is absent from the present 3GM computation. On the next plot at x/D=0.125 one can see that the present 3GM computation is still not fully turbulent with a marked peak at Strouhal number (St = f D/Uj) close to 4, while the triggered computation shows a much broader spectrum. However just downstream of the nozzle, from x/D = 0.25 to x/D = 1.0, the present computation shows broadband spectra in good agreement with the triggered computation, indicating that a correct turbulent transition has occurred, thanks to the fine resolution of the 3GM grid. At this point, the role of turbulence triggering is still an open question and further works are needed to decide whether a refined grid is sufficient to provide the conditions to avoid having to consider turbulence triggering in the nozzle boundary layers. In analyzing these results, it is believed that the ability of the code to properly resolve the turbulence at the nozzle outlet is a crucial point for the simulation of the jet and the resulting noise. One of the features that is thought to be important is the proper reproduction of the early transition to turbulence, right at the nozzle outlet. This particular point seems to have a strong dependence on the azimuthal resolution of the grid.

The proper reproduction of this turbulent behavior seems to be essential to capture proper jet development and noise emisssion [START_REF] Bogey | Influence of initial turbulence level on the flow and sound fields of a subsonic jet at a diameter-based Reynolds number of 10 5[END_REF][START_REF] Lorteau | Recent progress in LES computation for aeroacoustics of turbulent hot jet. Comparison to experiments and near field analysis[END_REF] . In that regard, the refined 3GM grid clearly produces better results, although homogeneous nozzle inlet conditions and no turbulent tripping are yet used at the present stage of these computations. In particular, earlier turbulent transition is visible for the 3GM grid which reduces the formation of the vortex pairing phenomenon. This phenomenon is suspected to be the source of spurious acoustic emissions in less resolved LES computations. In that regard, the high degree of azimuthal resolution permitted by the tetrahedral mesh 3GM (see table 1) is seen as a favorable point to avoid excessive vortex pairing by a better description of isotropic resolved turbulent scales. The importance of azimuthal resolution is discussed by Bogey et al. [START_REF] Bogey | Large-eddy simulation of the flow and acoustic fields of a Reynolds number 10 5 subsonic jet with tripped exit boundary layers[END_REF] .

Acoustic far field A. Integral surfaces

Concerning the far field prediction, several FW-H surfaces were considered, as illustrated on figures 13 and 14. As shown on figure 13, older surfaces were found to be too close to the jet flow and new surfaces were considered, figure 14. However these new surfaces, used for the 3GM computations, appeared to be too short and were later on extended downstream, as illustrated on figure 14. The extension permitted to include low frequency sources, as illustrated on figure 15, that could not be properly taken into account with the shorter ones. To save computational resources the new extended surfaces were evaluated with the coarse grids only and the low frequency correction was simply added to the finer grids, assuming that the coarse grid did properly resolve this low frequency part of the radiated noise.

B. Far field OASPL

Figure 16 illustrates the far field pressure level (Overall Sound Pressure Level, OASPL) and PSD, reconstructed with the KIM code from LES solution stored on the FW-H surfaces. Levels and spectral contents compare favorably to measurements. It is also apparent on this figure that mesh refinement can significantly improve the prediction of both the spectral contents and the far field levels, with 3GM results being mostly in a ±2 dB band from the measurements. This result is comparable to the far field noise computed by Lorteau et al. [START_REF] Lorteau | Recent progress in LES computation for aeroacoustics of turbulent hot jet. Comparison to experiments and near field analysis[END_REF] on their "Very Fine" mesh (structured mesh of 240 million cells). The largest departure is observed in the downstream direction (low angles) and seems to correspond to an underestimation of the low frequency contents. It is anticipated that this discrepancy comes from the initial new FW-H surfaces which miss some low frequency acoustic contributions, although a crude correction was incorporated from coarse grid computation.

We also see the effect of storage in figure 15, which affects only the low frequencies. This helped to correct the low frequencies of the 3GM, but not for 2GM, because the old storage surfaces are different (compare Figure 13 and14). It is worthwhile to note that the mesh refinement, as expected, extends the cut-off frequency of the computations, as visible on figure 16.

These results, obtained on fully unstructured predominantly tetrahedral meshes, are in line with previous mesh dependency studies obtained, with the same code, on structured meshes, for the same jet, see references [START_REF] Huet | Influence of boundary layers resolution on heated, subsonic, high Reynolds number jet flow and noise[END_REF][START_REF] Lorteau | Recent progress in LES computation for aeroacoustics of turbulent hot jet. Comparison to experiments and near field analysis[END_REF] . This provides a good confidence in the present computations on tetrahedral meshes.

IV. Conclusion

The present tetrahedral mesh calculations with coarse and medium grids, with the 2013/2014 versions of the CEDRE code, have provided consistent results, similar to those obtained previously on structured meshes with the same code. For the finest unstructured grid considered in this work, which was composed of 130 millions cells, very good agreement with measurements was presented for both the jet flow solution and the radiated far field noise (levels and spectral contents). This agreement was found comparable to the one obtained from reference structured grid computation on 240 million cells. This provides an essential point of validation to the present numerical approach based on fully unstructured grid for predicting jet noise.

The question of turbulence triggering appears as an open question and further work may be needed to decide if this is a requirement that can be compensated with the higher degree of mesh isotropy and azimuthal resolution permitted by unstructured grids.

The advantage of the full unstructured approach is that coarse (test) grids, which are often sufficient to the overall setting of the computations, can be easily modified to produce more refined grids, in order to provide reliable computations of the solution.

This validation work provides a new opportunity to more freely choose and adapt mesh size, which is promising for applications to more complex geometries where mesh construction is often a limitation. This is particularly crucial for the computations of complex nozzles, including noise reduction devices, such as chevrons, mixers, or of installation effects.

Figure 2 .

 2 Figure 2. Overview and zoom of the 2GG tetrahedral based mesh.

2GG 3GM Figure 3 .Figure 4 .

 34 Instantaneous fields of axial velocity (top) and vorticity modulus (bottom).2GG 3GM Instantaneous fields of ratio specific heat capacities (top) and CFL number (bottom).

Figure 8 .Figure 9 .

 89 Comparison of instantaneous vorticity field in the jet at the nozzle exit2GG2GM 3GM Effect of refinement on instantaneous pressure field at nozzle exit.

Figure 10 .

 10 Figure 10. Evolution of normalized axial velocity fluctuations along the nozzle lip line.Green: 2GG, orange: 2GM, red: 3GM.

Figure 11 .

 11 Figure 11. Comparison of the axial evolution on the jet axis of the mean and rms axial velocity for the present 3GM simulation (red line) and Lorteau et al. 14 Very Fine computation (blue line) -symbols: experiments, square Onera/Φ80, triangles Ahuja et al. 5 .

Figure 12 .

 12 Comparison of PSD of radial fluctuating velocity normalized by Uj at r/D = 0.5 for x/D ∈ ∈ ∈ ∈ {-0.25; 0.125; 0.25; 0.5; 1.0} for the present 3GM computation (red line) and Lorteau et al. 14 Very Fine computation (blue line)

Figure 13 .

 13 Figure 13. Instantaneous fields of pressure, density gradient and divergence of mass flow, with old storage FW-H surfaces for acoustic KIM code.

Figure 14 .

 14 Figure 14. Instantaneous fields of pressure, density gradient and divergence of mass flow, with new storage FW-H surfaces and their downstream extension for acoustic KIM code.

Figure 15 .

 15 Figure 15. Effect of extended FW-H surfaces on low frequencies (30°, 60°, 90°, 120°) , blue line original surfaces and red line extended surfaces, as computed on the 2GG grid.

Figure 16 .

 16 Figure 16. Comparison of average spectra at 30, 60, 90 and 120 degrees, then OASPL.Green: 2GG, orange: 2GM, red: 3GM -black: experiment, Onera/Φ80 at Cepra19.

Table 1 . Hot jet characteristics

 1

	Uj (m/s)	Mj	Tj/Ta	Ma	Ttot/Ta	Ta (K)	Re		ptot/pa
	410	0.7	2.96	1.2	3.2	280	4. 10	5	1.4

Figure 1. Nozzle geometry B. Numerical settings and boundary conditions Present

  LES computations have been carried out with recent versions of Onera in-house CEDRE code (2013 and 2014 versions) relying on second order MUSCL space and one step implicit time schemes. The improved gradient limitation methods available on recent versions were activated. CEDRE is a multi-physics solver which is developed for industrial and research applications in the fields of energetic and propulsion and has been consistently used for jet noise studies.

Table 2 . Grid properties

 2

						Max cell	Number of
	Grid	# Mnodes	# Mcells	# Mfaces	Cell size at nozzle exit	size at edge of refined	cells in azimuthal direction at
						zone	nozzle exit
	2GG	3.1	17.1	34.4	1.0 mm	8.8 mm	251
	2GM	11.1	55.5	113.5	0.6 mm	5.5 mm	419
	3GM	23.5	130.2	262.3	0.3 mm	4.5 mm	838

Research scientist, Fundamental and Applied Energetics department

Research scientist, Fundamental and Applied Energetics department, francois.vuillot@onera.fr

Acknowledgments

The authors want to thanks colleagues from Onera, Mathieu Lorteau, Maxime Huet and Franck Cléro for numerous exchanges and comparisons with their computations of the same jet conditions with the same code on different grids. The presented work has been supported from internal grants from the Fluid Mechanics and Energetics Scientific Branch.