Jana Brunner
email: brunner@zas.gwz-berlin.de

Phil Hoole

Pascal Perrier

ARTICULATORY OPTIMISATION IN PERTURBED VOWEL ARTICULATION

Keywords: optimisation, motor control, perturbation

A two-week perturbation EMA-experiment was carried out with palatal prostheses. Articulatory effort for five speakers was assessed by means of peak acceleration and jerk during the tongue tip gestures from /t/ towards /i, e, o, y, u/. After a period of no change speakers showed an increase in these values. Towards the end of the experiment the values decreased. The results are interpreted as three phases of carrying out changes in the internal model. At first, the complete production system is shifted in relation to the palatal change, afterwards speakers explore different production mechanisms which involves more articulatory effort. This second phase can be seen as a training phase where several articulatory strategies are explored. In the third phase speakers start to select an optimal movement strategy to produce the sounds so that the values decrease.

INTRODUCTION

Several findings from the area of speech motor control and neural processing have led to the assumption that speech motor control involves an internal model in the speaker's brain. In its most general form this model can be seen as some kind of an internal image of the speech production system which serves to predict the acoustic output for a certain motor command input.

Different proposals have been made with regard to the exact structure of the internal model. Most of the proposals include two main mappings, one between motor commands and articulatory configurations and another one between articulatory configurations and sounds. Both these mappings are set up during speech acquisition when the speech production system provides outputs for a given input, namely an articulatory configuration for a given set of motor commands and an acoustic output for a given articulatory configuration.

In this context, speech production planning then would inverse the process by finding a set of mo-tor commands for a given articulatory configuration and an articulatory configuration for a given acoustic output. Since, however, both these mappings are one-to-many this process involves finding a solution to an ill-posed problem.

There are several proposals as to how speakers solve this problem. Most of them can be seen as a kind of optimisation in the sense that the "best" mapping is preferred over the others. They differ in the way in which "best" is defined. One approach (Jordan e.g. [START_REF] Jordan | Indeterminate motor skill learning problems[END_REF]) argues that the search for the best mapping involves an optimisation of the articulatory movement via minimisation of the articulatory effort.

In the study presented here we look for evidence for the existence of such an optimisation process. Under perturbation, the speech production system is changed; as a consequence parts of the internal model should also change. Just as during the speech acquisition process, the system is confronted with new mappings of motor commands to articulatory configurations, possibly also with new mappings between articulatory configurations and acoustic outputs. Since the internal model lacks information about the perturbed condition, one can therefore hypothesise that speakers at first will just shift the motor commands in response to the perturbation without changing other movement characteristics, as for example the smoothness of the movement. Since the acoustic output will not in all cases be perfect, one can furthermore hypothesise that in a second phase they will try out other strategies. This will involve increased articulatory effort. In a final phase, they will try to optimise the movement. This should result in a drop of articulatory effort.

Support for this hypothesis of on adaptation process involving several stages comes from previous experimental results for example for the development of accuracy in pointing tasks for children. [START_REF] Favilla | Reaching Movements in Children: accuracy and reaction time development[END_REF] found that this development is not linear but there are stages where accuracy decreases. These stages are comparable to our second stage where new strategies are tried out, which are not necessarily successful in all cases.

METHODS

Palatal prostheses were made for five German native speakers (more speakers are currently being analyzed). For three of them (represented by solid lines in Fig. 1) these prostheses moved the alveolar ridge posteriorily, for the other two (represented by dotted lines) the palatal vault was filled out so that the palate became flatter. Subjects wore the palates for 14 days and were recorded regularly over this period via EMA (Carstens AG 500 for the speaker represented by a black solid line, Carstens AG 100 for the others) in different conditions:

• Day 1:

-session 1: without artificial palate -session 2: with artificial palate, with auditory feedback masking due to white noise over headphones (session missing for one speaker) 1 -session 3: with artificial palate with full auditory feedback • Day 8:

-session 4: with artificial palate • Day 15:

-session 5: with artificial palate 2 The target sounds /e, i, o, y, u/ were embedded in the nonsense words /'te:ta/, /'ti:ta/, /'to:ta/, /'ty:ta/ and /'tu:ta/ which were produced in the German carrier phrase Ich sah ... an. (I looked at ...). The sentences were repeated 20 times per session in randomised order. Three sensor coils were placed on the tongue, one at about 1cm behind the tongue tip, one opposite the border between the hard and the soft palate and a third one in between the two. One further sensor was placed below the lower incisors in order to track jaw movements, one at the upper lip and one at the lower lip. Two coils at the upper incisors and the bridge of the nose served as reference sensors to compensate for head movements. 3 A parallel acoustic recording was carried out on a DAT recorder. The vocalic gestures of the tongue tip (downward movement from the consonant to the vowel) was labeled on the tangential velocity signal using a 20% threshold criterion. The tongue tip sensor was taken because the tongue tip is not involved in fulfilling the task and optimisation should therefore be more evident. Articulatory effort was assessed by means of peak tangential acceleration and tangential jerk ( [START_REF] Nelson | Physical Principles for Economies of Skilled Movements[END_REF]). These measures were chosen because they describe kinematic optimisation and not for example an optimisation in terms of a reduction of the Euclidean distance between two sounds following each other.

Articulatory effort at three points in time was of interest:

• During the unperturbed session. Here one can expect maximally optimised movements since the speaker uses articulatory strategies acquired a long time ago. • During the perturbation at the point of maximal effort. Maximal effort signals that the speaker tries out new strategies at this stage. Since the effort decreases afterwards this can be seen as the onset of the optimisation. • During the last perturbed session. Here again one can expect optimised movements since the speakers can be expected to have chosen the most efficient strategy among all the strategies they tried out in phase 2. Possibly optimisation takes longer than our experiment so the optimisation observed here might only be partial. Repeated measures ANOVAs have been calculated for data split by speaker (SPSS 13, Windows XP). Since the standard error varies for the different sessions, a Greenhouse-Geisser correction of the degrees of freedom has been carried out. Bonferroni posthoc tests for the difference between the unperturbed value and the maximal value, and between the maximal value and the last perturbed value for each of the two parameters were calculated.

RESULTS

The following figure shows the results for the measurements of articulatory effort during the tongue tip gesture towards /o/. The upper panel presents the results for peak acceleration, the lower one the ones for jerk. Different colours and line styles stand for different speakers. The abscissa presents the number of the session as given in the list above. The ordinate presents the value of the parameter. Plots give means and standard error. If the difference between the unperturbed value and the highest value is significant the part of the line from session 1 to the highest value is plotted in bold. If the difference between the highest value and the last session is significant the part of the line from the highest value to the last value is plotted in bold.

For all speakers a significant increase of articulatory effort can be found (difference between session 1 and the session with the maximal values is significant). A significant decrease can be found in all but one case (jerk for the speaker represented by the black dotted line).

The temporal location of the peak points is speaker dependent (e.g. session 3 for the speaker represented by the light grey solid line, but session 4 for two other speakers). Phase 1, where we expected constant values, can, if at all, only be found for the speaker represented by black dotted lines. For the others it is so short (one or two repetitions) that it cannot be seen in the plot. Phase 2 (with maximal values) can be found for all speakers. Phase 3, the optimisation process (starting at the peak points), can be found for four speakers only, for the fifth the difference between the point with maximal effort (session 4) and the final session is not significant (signaled by the thin line). For /u/ the results are similar. A significant increase in peak acceleration has been found for all speakers and a significant increase in jerk for four of the five speakers. A significant decrease has been found for all but one speaker for both parameters. Phase 1 can again be found for only for one speaker. For /y/ the results are less clear. Again phase 1 could be found for only one speaker. A significant increase of both, peak acceleration and jerk, has been found for four speakers, a significant decrease for three speakers if measured as peak acceleration and four speakers if measured as jerk. For /i/ a significant increase has been found for only one (peak acceleration) and two (jerk) speakers, a decrease has been found for four speakers if measured as peak acceleration and for three speakers if measured as jerk.

For /e/ a significant increase has been found for one speaker only. A decrease has been found for three speakers if measured as peak acceleration and two if measured as jerk.

CONCLUSION

To summarise the results, evidence for the existence of phase 1 where the speaker shifts the motor commands but articulatory effort stays constant is difficult to find. This is because this phase is very short, maybe one or two repetitions. For one speaker, however, a rather long phase 1 could be found. Evidence for our hypothesised phase 2, where new strategies are explored and the articulatory effort increases is rather clear. The same is true for phase 3, where the effort decreases because the movements become more and more optimised.

Our results therefore support the hypothesis of the solution of the inverse problem via an optimisation of articulatory effort. As a response to a perturbation speakers compensate in three phases: At first, there is a shift in the system (phase 1). Different articulatory positions are produced but the articulatory effort stays the same. Afterwards, new strategies are explored and the effort increases. Finally, the speakers select the most optimal movement strategies so that the effort decreases.

There is a high dependency on the vowel, the results are clearer for rounded than for unrounded vowels and clearer for back vowels than for front vowels. Furthermore, there is a dependency on the speaker. The speakers represented by the black solid and the dark grey solid line show clearer results than the others. The temporal location of phase 2 (peak point) also depends on the speaker more than on other factors. A possible dependency on the kind of prosthesis can be suspected: The speakers represented by the solid lines (with alveolar prostheses) present clearer results. However, this has to be verified by further experiments with more speakers being carried out at the moment.

For the differences between the vowels several explanations can be found. A reason for the better results for /u/ and /o/ as compared to the front vowels could be that the speakers use the palate as an upper limit ( [START_REF] Stone | How the tongue takes advantage of the palate during speech[END_REF]) during the high front vowels so that the possibilities to produce other strategies which possibly involve more articulatory effort are rather limited. This would mean that most speakers do not leave stage 1 during the production of /i/ and /e/.

Another reason for the clearer results for the back vowels is the fact that the acoustics of the two back vowels is not disturbed by the prosthesis since the formants originate in resonances of the back tube ( [START_REF] Apostol | A model of acoustic interspeaker variability based on the concept of formant-cavity affiliation[END_REF]). This means that for the two back vowels the second phase can start immediately after perturbation onset.

A reason for the better results for rounded as compared to unrounded vowels could be that the rounded sounds allow for more variation in tongue position (and therefore more strategies to explore) because acoustic consequences of the variation can be compensated for by different degrees of lip rounding. The investigation of lip protrusion for the vowels is not finished yet, the investigation of /S/ for the same experiment, however, has shown that under pertur-bation speakers indeed use different degrees of lip rounding for this sound ( [START_REF] Brunner | Temporal development of compensation strategies for perturbed palate shapes in German /S/production[END_REF]). One can assume that they use similar strategies for rounded vowels resulting in a greater number of articulatory strategies.

The results could also be interpreted according to a gestural approach. If one assumes that the target words involve coordinative structures the increases in acceleration and jerk could indicate initial difficulty in coordinating C and V gestures which later diminishes with more and more practice. This then results in less articulatory effort.

One could therefore argue that every model, not only the Jordan approach, predicts optimisation. However, it is important to note that the optimisation process here optimises kinematic parameters and not only Euclidean distance towards the target as for example in the approach by Guenther ([4]). This kind of optimisation process is well illustrated by Jordan's model.

It has been remarked that the inverse problem from articulatory configurations to acoustics has been overstated. Evidence has been provided that in normal speech speakers do not use different vocal tract shapes in order to produce the same acoustic output (e.g. [START_REF] Hogden | Accurate Recovery of Articulator Positions From Acoustics: New Conclusions Based on Human Data[END_REF] and [START_REF] Yehia | A study on the speech acousticto-articulatory mapping using morphological constraints[END_REF] who managed each to recover vocal tract shapes measured on a subject from the acoustic output). However, rather than claiming that the inverse problem has been overstated, our study rather opts for another explanation of this finding. Speakers do not use different articulatory strategies because, once they have found the most efficient strategy, they stay with it. 1 The influence of feedback masking is not subject of this paper and will therefore not be discussed. During the session speakers could view a sound level display and were instructed to keep their level about equal. Articulatory effort, if it is increased for this session, does therefore not increase as a result of a higher sound level. 2 A postperturbed session was recorded as well. Since, however, the topic of this paper is optimisation, aftereffects will not be discussed. 3 [9] report a slight and inconsistent perturbation of speech due to the sensors in EMA recordings resulting in tongue retraction and tongue lowering for some speakers. Since, however, the influence is small and we are furthermore comparing EMA-recordings with each other and not EMA recordings with non-EMA recordings, the effect found between sessions cannot be attributed to the EMA sensors but only to the prosthesis.

Figure 1 :

 1 Figure 1: Optimisation during the gesture towards /o/. Solid lines show results for speakers with alveolar prostheses, dashed lines results for speakers with central prostheses. Different colours and line styles present different speakers. Error bars show standard error.


Bruno Gaume
email: gaume@univ-tlse2.fr

Karine Duvignau
email: duvignau@univ-tlse2.fr

Laurent Prévot
email: prevot@univ-tlse2.fr

Yann Desalle
email: desalle@univ-tlse2.fr

Toward a cognitive organization for electronic dictionaries, the case for semantic proxemy

We compare a psycholinguistic approach of mental lexicon organization with a computational approach of implicit lexical organization as found in dictionaries. In this work, we associate dictionaries with 'small world' graphs. This multidisciplinary approach aims at showing that implicit structure of dictionaries, mathematically identified, fits the way young children categorize. These dictionary graphs might therefore be considered as 'cognitive artifacts'. This shows the importance of semantic proximity both in cognitive and computational organization of verbs lexicon.

Introduction

According to [START_REF] Dik | Functional grammar[END_REF] a linguistic theory should be compatible with psycholinguistic research on language acquisition, treatment, production, interpretation and memorization of linguistic expressions. We agree with this view and postulate that elaborating electronic dictionaries on the ground of a linguistic theory, satisfying Dik's principle, will confer them good ergonomics that will increase their usability. Our approach is to some extent comparable to WordNet initiative [START_REF] Fellbaum | WordNet: An Electronic Lexical Database[END_REF], in the sense that we are trying to characterize speakers' mental lexicon.

In this paper, we focus on verb lexical organization through the examination of verbal pivot metaphorical utterances (VPMU). Such utterances involve an understudied structural aspect of the lexicon: interdomain co-hyponymy (Duvignau, c 2008. Licensed under the Creative Commons Attribution-Noncommercial-Share Alike 3.0 Unported license (http://creativecommons.org/licenses/by-nc-sa/3.0/). Some rights reserved.

2002; [START_REF] Duvignau | Between words and world: Verbal 'metaphor' as semantic or pragmatic approximation[END_REF]. In this context, we take semantic proximity as a central principle for cognitive ergonomics influencing dynamic lexical acquisition and adult lexical organization. VPMU generally consists in substituting elements from different semantic domains. They are usually considered as deviants while they might constitute a linguistic illustration of the categorial flexibility advocated in [START_REF] Piaget | La formation du symbole chez lŠenfant[END_REF][START_REF] Ny | La sémantique psychologique[END_REF][START_REF] Hofstadter | Fluid concepts and creative analogies[END_REF]. They might therefore reveal an early lexical structuring mode that may form a ground for improving electronic dictionaries.

This paper presents a mathematical method able to discover the areas in which this structuring mode appears in dictionaries. Our approach is to take advantage of the mathematical structure of the network generated by verb definitions. This structure has been mentioned in [START_REF] Watts | Collective dynamics of small-world networks[END_REF], studied for WordNet by [START_REF] Sigman | Global organization of the wordnet lexicon[END_REF], refined in [START_REF] Gaume | Forms of meaning, meaning of forms[END_REF] and exploited in the current proposal.

The paper is organized as follows. The next section brings evidence of categorization by semantic proximity from early lexicon acquisition experiments. Section 3 presents the computational model, hereafter 'proxemy'. Section 4 details our work on lexical graphs while section 5 compares the results of experimental studies with those of the computational model.

Toward a categorization by semantic proximity: evidences from early lexicon acquisition

In order to show the importance of semantic approximation, we have chosen to support our claim with productions observed at the crucial period of lexical construction (between 2 and 4 years-of-age) and to compare these with adult speakers that have a stabilized lexicon.

Inter-domains vs. intra-domain semantic approximations

Studies in this field are almost exclusively limited to nominal utterances. [START_REF] Duvignau | Semantic approximations intraconceptȃvs. interconcepts in early verbal lexicon: flexibility against error[END_REF] established the existence of the production of metaphor-like utterances with a verbal pivot in 2-4 years-old children and proposed to consider them, at this stage of language development, as semantic approximations and not as mistakes or true metaphors. Duvignau distinguished two kinds of semantic approximations: Inter-domains proximity and intra domain proximity between verbs [START_REF] Duvignau | La métaphore berceau et enfant de la langue[END_REF].

-Inter-domains proximity / co-hyponymy between verbs : a 'linguistic approximation'

(1) Elle déshabille l'orange (She undresses the orange) [Age: 3 years] [movie: a lady peels an orange]

In this category of approximation, the verb used by the speaker constitutes a reference to a semantic domain different from the one of element it is combined to ('undress' / 'orange'). For this reason, the approximate character of the verb is understandable independently of the context of the utterance: detecting the approximation occurs at the linguistic level. We call this type of production 'semantic approximation'. They might constitute a metaphor or an 'analogic surextention'. When someone has a conventional verb in the mental lexicon ('to peel') and use a non conventional but relevant verb like 'to undress the orange' for the action [to peel the orange his verbal semantic approximation constitutes a metaphor. On the contrary when someone does not have a conventional verb in the mental lexicon but manages to use a non conventional but relevant verb in saying 'to undress' for this action, his verbal semantic approximation constitutes a 'surextension' but not an error because of the lexical relation that links these verbs. In fact, according to [START_REF] Duvignau | Between words and world: Verbal 'metaphor' as semantic or pragmatic approximation[END_REF] 'to undress' and 'to peel' are related by an inter-domains synonymic relation.

-Intra-domain proximity / co-hyponymy between verbs: a 'pragmatic approximation' In this category, illustrated by (2) the approximate character of the verb comes only from a non-correspondence between the verb used and the reality it designates. This happens with utterances in which the use of the verbal form does not create any semantic tension within the utterance but designates a way of carrying out an activity that does not correspond precisely to the action undertaken.

(2)

Elle coupe l'orange (She cuts the orange)[age: 3 years][movie: a lady peels an orange]

We propose an experimental study of the production of verbal semantic approximations like (2) or (1) by way of a naming task of 17 action-movies with young children (from 2 to 4 years old). We compare their performances with adult's ones.

Experimental Design

In order to elicit the production of semantic approximations we proposed to all our participants an action-video naming task. The population sample consisted of:

• 54 non-disturbed children (2-4 years old), monolingual in French

• 77 non-disturbed young adults (18-40 years old), monolingual in French

The action movies sequences are coming from the Approx protocol [START_REF] Duvignau | Semantic approximations intraconceptȃvs. interconcepts in early verbal lexicon: flexibility against error[END_REF]. The material consists in 17 action-movies sequences described in table 1.

The 17 action movies are presented in random order to each participant. Instructions were given at the time the action in the movie was completed and its results were visible (e.g when the glass is broken). At that moment a question was asked to the participant:'What did the woman do? (just now)'

Results

Each of the children produced between 2 and 5 approximations: 'Elle casse la tomate' -She breaks a tomato ' [action = to squash], 'Elle épluche le bois' -She peels the wood' [action = to strip the bark off a log]. Globally, children produced semantic approximations for 34 % of the naming tasks, which were distributed as follows: 24 % intra-domain semantic approximations, 10 % inter-domains semantic approximations. They produced them significantly more frequent than adults : 5 % with 4% intra-domain semantic approximations and 1 % inter-domains semantic approximations. The student Test shows the difference between children and adults in terms of production of semantic approximation is very significant: here p < 0, 01 while p < 0, 05 is enough.

These results signal the importance of semantic approximations and of semantic proximity between verbs in the cognitive organization of verbs lexicon.

In the rest of the paper we present a computational model of semantic proximity and then compare this model with the experimental data obtained from the children.

3 Proxemy: a computational approach A theory of language useful for computational work must account for language statistical regularities. Zipf law [START_REF] Zipf | Human behavior and the principle of least effort[END_REF] satisfy this observation but provides little insight on lexical structural organization. More recent graph theory studies [START_REF] Ferrer-I-Cancho | The small world of human language[END_REF][START_REF] Sigman | Global organization of the wordnet lexicon[END_REF], capitalizing on results in other scientific domains, provided interesting contributions to the establishment of such a theory of language. All structures discovered in this field research satisfy the "hierarchical small word" (HSW) definition (see section 3.1). Our approach takes place in this general framework. Our specificities are:

• a new linguistic and psycholinguistic insight that guides us and help us on our results validation;

• the kind of objects studied (dictionaries);

• our analysis of graph structure resulting in a computational model of semantic proximity among vertices (here vertices are French verbs).

The study by [START_REF] Resnik | Measuring verb similarity[END_REF] signaled that although existing models for verb similarity performed reasonably well against human judgments, none managed to handle certain types of metaphorical pairs such as to undress / to peel off that are nonetheless declared to be rather similar by speakers. We aim to develop a model addressing this issue.

Small World Networks

Networks corresponding to structures found in real world (henceforth real world networks) are sparse: in a graph with n nodes, the maximum number of possible edges is O(n 2 ) while the number of edges in real networks is generally inferior to O(nlog(n)). Watts and Strogatz [START_REF] Watts | Collective dynamics of small-world networks[END_REF] proposed two indicators to characterize a large sparse graph G:

• L : the characteristic path length, i.e the mean of the shortest path between two nodes of G

• C : the clustering coefficient, C ∈ [0, 1], it measures the graph tendency to host zones very dense in edges. (The more clustered the graph is, the more the graph's C approaches 1, whereas in random graphs C is very close to 0).

In applying these criteria to different types of graphs, Watts and Strogatz found that:

• real world networks have a tendency to have a small L: generally there is at least one short path between any two nodes ;

• real world networks have a tendency to have a large C: this reflects a relative tendency for two neighbors on the same node to be interconnected;

• random graphs have a small L: If someone builds a graph randomly with a density of edges comparable to real world networks, it will obtain graphs with a small L;

• random graphs have a small C: They are not composed of aggregates. In a random graph there is no reason why neighbors of a same node are more likely to be connected than any two other nodes, hence their poor tendency to form aggregates.

Watts and Strogatz proposed to call the graphs having these two characteristics (a small L and a large C) small worlds (SW). They recognized these SW in all the real world networks they observed, and therefore postulated for being a SW was an universal property of real world networks. A complete presentation of Small Words can be found, for example, in [START_REF] Newman | The structure and function of complex networks[END_REF].

More recent research has shown that most SW also have a hierarchical structure (hereafter hierarchical small worlds, HSW ). The distribution of the vertices incidence degrees follows a power law. The probability P (k) that a given node has k neighbors decreases as a power law, P (k) ≈ k -λ , where λ is a constant characteristic of the graph [START_REF] Barabási | Emergence of scaling in random networks[END_REF], while random graphs conforms to a Poisson Law.

In the next section, we present 'proxemy', a semantic proximity measure based on a distance we define. A interesting particularity of this distance is to calculate the distance between two vertices on the ground of the complete graph, and not only on their direct neighbors.

The mathematical model

PROX (PROXemy) is a stochastic method designed for studying "Hierarchical Small Worlds".1 This method takes graph as input and transform them in a Markov chain whose states are graph vertices. Metaphorically, energy particles wander randomly from vertex to vertex through the edges of the graph. It is their trajectory dynamics that give us the structural properties of the graph.

PROX takes a graph in input and output a similarity measure between the vertices of the graph. Our problem is therefore the opposite than the one of Pathfinder networks (PFNETs see [START_REF] Schvaneveldt | Graph theoric foundations of pathfinder networks[END_REF]). PFNETs take a full proximity matrix in input and output a sparse graph. Their goal is to minimize the number of edges required in the sparse graph to be able to approximate the full distance matrix corresponding to the initial full proximity matrix.

PROX build a similarity measure between the vertices. The hypothesis is that areas having a high density in edges (hereafter, these areas will be called aggregates) correspond to closely related verb meanings (in a graph of verbs).

Given a graph with n vertices, G = (V, E), we will note [G] the matrix n × n such that ∀r, s ∈ V , [G] r,s = 0 if {r, s} ∈ E and 1 otherwise. [G] is called the adjacency matrix of G.

Given G = (V, E) a reflexive graph with n vertices.

[ Ĝ] is a n × n matrix defined by ∀r, s ∈ V , [ Ĝ] r,s = [G]r,s x∈V {[G]r,x} . [ Ĝ] is the Markovian matrix of G.
[ Ĝ] is the n × n matrix is a transition matrix of the homogeneous Markov chain whose states are the vertices of the graph such that the probability of going from one vertex r ∈ V at an instant t onto another s ∈ V at the instant t + 1 is equal to:

• 0 if {r, s} ∈ E (s is not neighbor of r) • 1/D if {r, s} ∈ E and r has D neighbors (s is a neighbor of r) 2
Given G = (V, E) a reflexive graph with n vertices and [ Ĝ] its Markovian matrix, ∀r, s ∈ V, ∀t ∈ N * , P ROX(G, t, r, s) = [ Ĝt ] r,s P ROX(G, t, r, s) is therefore the probability for a particle departing from r at the instant zero to be on s at the instant t.

Therefore when, P ROX(G, t, r, s) > P ROX(G, t, r, u), the particle has more probability to be, at instant t on s than on u and it is graph structure that determine these probabilities.

For the rest of this paper we will set the value of t to 4 since L is less than 4 in the kind of graph we are concerned with. Therefore, we take into account the global graph simply by calculating P ROX(G; 4; r; s). Now we have defined our model we will present lexical graphs on which we apply it.

Lexical graphs

Several types of lexical graphs can be built according to the type of the semantic relation used for defining the graph's edges. The two principal types of relations used are:

• Syntagmatic relationships, like co-occurrence relationships: they define edges between nodes corresponding to words found near to each other in a corpus.

• Paradigmatic relationships, like synonymy: they define, on the ground of lexical databases such as WordNet [START_REF] Fellbaum | WordNet: An Electronic Lexical Database[END_REF], edges between nodes of words being in a synonymy relationship in such resource.

Moreover, we are interested into less specific relations, called semantic proximity relations or semantic relatedness, and which covers both paradigmatic and syntagmatic dimensions.

Dictionary graphs

Meaning in dictionary definition is at least partially brought by the relations they create between the words constituting the entries. Our approach consists in exploiting the small word properties of the graphs corresponding to dictionaries. More precisely, we are taking advantage of our hypothesis that aggregates correspond to areas of closely related senses. We illustrate our approach on two kinds of dictionary, two traditional dictionaries, Le Grand Robert3 and TLFi4 , and an synonym dictionary (Dicosyn) made of compilation of synonym relations extracted from seven other dictionaries (Bailly, Benac, Du Chazaud, Guizot, Lafaye, Larousse et Robert). 5We create a graph from a dictionary in the following way. The entries constituted the vertices. Edges between two vertices A and B were added if and only if B appears in A's lemmatized definition 6 as illustrated in Figure 4.1

We proceed in this way for each entry and obtained a graph of the dictionary. By extracting the subgraph composed only of verbs, the 'neighborhood' we get for the verb 'écorcer' is illustrated by Figure 4.1. Then we render the graph symmetric and reflexive. These modifications on the graphs are allowed thanks to its paradigmatic nature. Graphs created in this way are typical small before being corrected at CRISCO laboratory (http://elsap1.unicaen.fr/dicosyn.html).

6 Lemmatization has been realized with TreeTagger (http://www.ims.unistuttgart.de/projekte/corplex/TreeTagger/). [START_REF] Duvignau | La métaphore berceau et enfant de la langue[END_REF] has shown that co-hyponymy verb lexical organization according fits with a power law distribution of incidence degrees. In our opinion, (i) the hierarchical organization of dictionaries is a consequence of the special role of the hypernymy relation together with the polysemy of some specific vertices; (ii) the strong C reflects the role of interdomain co-hyponyny [START_REF] Duvignau | La métaphore berceau et enfant de la langue[END_REF][START_REF] Duvignau | Linguistic, psycholinguistic and computational approaches to the lexicon: Contributions to early verb-learning[END_REF]. For example, in French language, 'casser (to break)' appears in many definitions: 'émietter (to crumble)', 'fragmenter (to fragment)', 'détériorer (to damage)', 'révoquer (to dismiss)', 'abroger (to abrogate)'. This results in a very high incidence for the vertex 'casser (to break)'. Moreover, many triangles exist ( {casser, émietter, fragmenter}, {casser, révoquer, abroger}...,) and they help to create aggregates. These areas that are bringing co-hyponyms closer in the resulting graph.

Disambiguization for creation dictionary graphs

Word Sense Disambiguation is a general issue for natural language processing that we need to address when we build our graphs. We need to disambiguate the verbs we found in the definition facing a similar problem as [START_REF] Harabagiu | Wordnet 2 -a morphologically and semantically enhanced resource[END_REF]. For example, in French dictionary Le Grand Robert, there are two distinct entries for the verb 'causer': to cause (3) and to chat (4).

(3) CAUSER-1: être la cause de. (to be the cause of)

(4) CAUSER-2: S'entretenir familièrement avec qqn. to chat with

Of course, the word 'causer' may appear in other definitions like 'bavarder' (to chat) . Although a French speaker knows that the 'causer' in ( 5) refers to the definition (4) our system for building the graph cannot disambiguate. The solution we propose is to (i) first create a fictive vertex which is not a dictionary entry and then (ii) adds two edges {CAUSER, CAUSER-1} and {CAUSER, CAUSER-2 }. When 'causer' is found in another definition like (5), we add the edge { BAVARDER, CAUSER } as illustrated in Figure ( 5). More generally, let suppose we found a word with k entries in a definition, we will then have S 1 , . . . , S k vertices corresponding to the entries a fictive vertex S. In case there is an edge {A, S} it is replaced by {A, S i } where S i is such that [ Ĝ4 ] A,S i = M AX 0<i≤k { Ĝ4 ] A,S i }. Then we remove all fictive vertices from the graph to get a disambiguated graph.

We can then apply PROX a last time on the disambiguated graph in order to get the closest word of a word according to our proxemy measure. For example, the PROX-closest words of écorcer (to bark -a tree-), calculated with t = 6 are: 1 ECORCER (to bark), 2 DÉPOUILLER (strip), 3 PELER (peel), 4 TONDRE (mow, shear), 5 ÔTER (remove), 6 ÉPLUCHER (peel, pare), 7 RASER (shave), 8 DÉMUNIR (divest), 9 DÉ-CORTIQUER (decorticate), 10 ÉGORGER (slit the throat of), 11 ÉCORCHER (skin), 12 ÉCALER (husk), 13 VOLER (steal), 14 TAILLER (prune), 15 RÂPER (grate), 16 PLUMER (pluck), 17 GRAT-TER (scrape), 18 ENLEVER (remove), 19 DÉ-SOSSER (bone), 20 DÉPOSSÉDER (dispossess), 21 COUPER (cut), 22 BRETAUDER (shear sloppily), 23 INCISER (incise), 24 GEMMER (tap), 25 DÉMASCLER (remove first layer of cork)7

Proxemy and Experimental studies

Prox is a robust method: changing randomly a few edges does not change significantly the results. The repartition of aggregates is not strongly affected by a random redistribution of some edges. However the relevance of our proxemy approach of lexical networks is tied to the linguistic representativity of the networks we use. Therefore, we tested the PROX model of four different dictionary graphs and we compared them to the psycholinguistic experimental results presented in section 2. The graph we compared were: 1. Graph.TLFI.Verb, a graph built as explained in 4.1 from TLFi 8 dictionary, 2. Graph.Robert.Verb, a graph built as explained in 4.1 from Le Grand Robert dictionary, 3. Graph.DicoSyn.Verb, in which there is a edge between two verbs if there are given as syn-onyms by one of the synonym dictionary composing DicoSyn 4. Graph.DicoSyn_20 built from Graph.DicoSyn but in which 20% of the edges are randomly removed and re-added.

For each of these graphs we looked at two variables to be related with the psycho-linguistics experiments: the answers incidence and the proximity of answers to a 'reference verb' Answers incidence We compare in the graph the average incidence degree between adult (ID adult ) and children answers (ID children ). The proximity of answers to a 'reference verb' Three linguist judges determined together for each movie which was the most appropriate verb to describe the action performed in the movie (hereafter R i is the reference verb for the movie M i ). For a given movie M i , an answer may therefore be ranked according to its proxemy according to R i .

For a lexical graph G = (V, E) composed of n words, and for a reference verb R i ∈ V , one can define rank Ri for ranking all the vertices of V in decreasing order resulting from a PROX iteration P ROX(G, t, Ri, •) on V (see section 3.2).

Table 3: Proximity between answers and reference verb

Our first hypothesis was that ID adult < ID children . According to the hypothesis children would learn first words corresponding to high incidence vertices. Then they would use them for talking about an large lexical area (e.g 'casser' (to break) is used by children while adults use a more precise verb like 'déchirer' (to tear) which has a lower incidence in dictionary graphs).

Our second hypothesis was that the mean of the rank of the children answers according to the reference verb is higher that the adult ones. When a child is attempting to communicate an event (e.g déchirer un livre, to tear a book) for which he does not have an already constituted verbal category, he would do an analogy with a past event (e.g to break a glass) and use this verb for describing the current event (e.g casser un livre, to break a book). The adult could use a number of more accurate verbs but their proxemic rank, with regard to the reference verb, is generally lower than the children ones.

The table 2 shows the results concerning answers incidence. Although some variability is observed across the graphs, our first hypothesis is validated for the 4 graphs. On the three first graphs the average incidence of answers is roughly twice as the adults one.

The table 3 illustrates the results concerning proxemic rank of answers according to the reference verb. Again, in spite of some variability across the graphs our second hypothesis is validated as well. Moreover, having in mind that the graph has about 10 000 vertices, we observe that although less close that adults answers, the children answers remain relatively close to the reference verb according to our proxemic measure.

Conclusion

Our psycholinguistic approach allows us to establish that semantic proximity between verbs play a fundamental role during the period of early lexical acquisition. We signaled the existence in the organization of the lexicon of a relation of cohyponymy between verbs. Based on these first observations we consider that productions based on semantic proximity are particularly interesting: they manifest the existence, at the surface level of discourse, of a lexical relation of inter-domain 'semantic proximity' between verbs not yet considered in linguistics.

Moreover we have seen that semantic approximations for verbs appear to fit the proximity values calculated by PROX. On the ground of these first results, we postulate that constructing electronic dictionaries on the ground of linguistic theory of lexical semantic organization that fits with early lexicon acquisition as well with adult lexical organization will provide them interesting ergonomics properties. This should increase their usability and might be taken into account for normalizing electronic dictionaries.

For example, we are developing a 'proxemic electronic dictionary' from TLFi. Such dictionaries enable to find an uncommon but precise verb like 'to bark' by using (i) a common verb like 'to undress' which is related to 'to bark' by semantic proximity and (ii) a word (e.g 'tree') bringing a relevant semantic domain. Moreover, in the definition of 'to bark' one can find: 'tree', 'grain' 'fruit' which are close from each other according to PROX ran on nouns. Finally, when we look for verbs that are close from both 'to undress' and 'tree', PROX provides the verbs: 'to cut, to ring, to peel, to notch, to bark, to incise,...' which constitute relevant verbs. Such a dictionary can be useful for didactic studies where it can complements approaches like and NLP for word sense desambiguization [START_REF] Gaume | Désambiguisation par proximité structurelle[END_REF] or de-metaphorization.

Figure 1 :

 1 Figure 1: Sub-graph near 'écorcer (to bark -a tree-)' from Le Grand Robert

Figure 2 :

 2 Figure 2: Sub-graph of the verbs near 'écorcer' from Robert

  Figure 3: Disambiguation: 'Causer', fictive vertice

Table 1 :

 1 Approx 17 action movies

Table 2 :

 2 Results for 'Answers incidence'

In this paper we will use the term 'proxemy' to refer to the obtained by PROX algorithm. It corresponds to some kind of semantic proximity.

In the context of this presentation of the model we do not consider weighted graphs. However when building the graphs we do consider information, such as the position of the word in the definition, for giving weight to the edges.s

A significant amount of work has been done to encode 'Le Grand Robert in a graph.

We would like to thank ATILF for making the TLFi resource available to us.

5 Dicosyn has been first realized at ATILF (Analyse et Traitement Informatique de la Langue Française),

Proposing a translation for such fine grained and sometimes polysemous words is impossible since proposing the translation include a certain form of disambiguisation as it is suggested by the work of[START_REF] Gale | A method for disambiguating word senses in a large corpus[END_REF].

8 http://atilf.atilf.fr/tlf.htm


Zahra Afrassiabian

Willy Leclerc

Mohammed Guessasma

Khashayar Saleh
email: khashayar.saleh@utc.fr

Numerical simulation of mechanical resistance of wet and dry powder cakes

Keywords: DEM simulation, yield stress, caking, agglomeration, capillary condensation, solid bridges, beam model

This paper presents a multi-scale modelling approach to describe the mechanical behaviour of caked powders as an aid to properly design and analyse the results of laboratory caking tests to assess the propensity of powders to cake. This approach is based on the integration of particle-particle interactions due to liquid or solid bridges over the total volume of a cake using Discrete Element Method (DEM). The overall objective of the study was to simulate the behaviour of agglomerated (caked) samples subjected to a mechanical, compressive or tractive, stress. First, simulations were performed on an assembly of particles subjected to capillary forces. The capillary forces being amenable to a theoretical description, the tensile strength of wet agglomerates was calculated theoretically and used as a reference to validate the numerical results. Simulations were based on the combination of conventional DEM method and the phenomenological model of capillary condensation. The stress-strain curves were then simulated and compared to the theoretical values predicted by Rumpf's model for wet agglomerates. After validation, numerical simulations were used to study the mechanical behaviour of wet or dry agglomerates subjected to compression. In order to take into account the cohesive effect of solid bridges at the level of elementary contacts in discrete modelling, an original Euler-Bernoulli type beam model between particles in contact was established. The simulated stress-strain results was in very good agreement with experimental data.

Introduction

Powdered products are a significant part of finished or intermediate products encountered in the various industrial sectors [START_REF] Charpentier | Managing complex systems: some trends for the future of chemical and process engineering[END_REF]. However, despite the numerous benefits of powders, their handling and storage often cause some severe problems (i.e., explosion, clogging, dust generation, etc.). This could lead to a destruction of products and even their corresponding manufacturing units or to deterioration of the expected functionalities of products. Among these issues, caking of powders is undoubtedly one of the most detrimental phenomena in terms of their end-use properties. Generally speaking, caking corresponds to spontaneous and undesired formation of a coherent mass from individual grains [START_REF] Afrassiabian | An overview of the role of capillary condensation in wet caking of powders[END_REF]. Caking phenomenon occurs as a result of formation of material bonds at the contact points between particles. These links, whose fundamental mechanisms of appearance are still insufficiently understood, are reinforced by pressure, the migration of matter and variations in humidity and temperature. Note also that this unwanted agglomeration of particles of a powder is often very slow, irreversible and difficult to predict. In this context, it is essential to have predictive tools (caking tests or numerical models) to characterize, as reliably as possible, the propensity of the powders to cake in a given environment.

A multitude of methods has been reported in the literature to characterize the propensity of powders to cake [START_REF] Zafar | A review of bulk powder caking[END_REF][START_REF] Carpin | Caking of lactose: A critical review[END_REF]. The test methods reported in the literature range from very basic tests (drop test) to much more sophisticated tests such as shear tests.

The common point of all these methods is that they include two stages [START_REF] Afrassiabian | A study on the caking behaviour of binary mixtures of lactose due to solid-state crystallisation of the amorphous phase[END_REF][START_REF] Afrassiabian | Multiscale investigation of caking phenomenon of polymorphic powders: from physico-chemical aspects to industrial applications[END_REF]:

-The first step consists in placing the powder under conditions such that bonds will form between the particles. This can be done by applying a stress (formation of van der Waals bonds or solid bonds), by humidifying the powder (liquid bridges) or by subjecting the powder to a moisture and/or temperature cycle (formation of solid bridges) -Once the powder has been caked, the sample can then be characterized mechanically by methods such as shear, compression, traction or indentation tests.

One of the main challenges in the characterization of caking remains the determination of the mechanical strength of the cake. Zafar et al. [START_REF] Zafar | A review of bulk powder caking[END_REF] and Carpin et al. [START_REF] Carpin | Caking of lactose: A critical review[END_REF] reviewed the different tests reported in the literature. These authors made a classification of tests based on the method used to measure the mechanical resistance of cakes. They distinguish mechanical test methods including shear cell, uniaxial compression, tensile, ICI, penetration and powder rheometer testings as well as alternative tests (sensory, sieving, sticky point measurement, blow tester, microscopic observations, etc.). Because the caking is a slow process, the preparation of cakes is generally accelerated by prior wetting of powders or by using an air current instead of stagnant air. Moreover, in a recent work, [START_REF] Samain | Caractérisation multi-échelle de l'efflorescence et du mottage du saccharose[END_REF] Samain provided an extensive review of the caking tests by crossing the methods to prepare the caked samples and those used to characterize the cakes. These works clearly showed that, the existing tests suffer from a poor repeatability, a part of which is often inherent in the system, but also because of inadequate testing methods or non-optimized cake geometry. On the other hand, there is still no reliable theoretical model to describe the dynamic behavior of a powder bed during mechanical tests under compression, traction or shear. The most used model in this field is that of Rumpf but this model corresponds to the tensile strength of a static stack of grains when in practice, it is practically impossible to perform the tensile tests (but compression or shearing) and in addition, tests are dynamic rather than static. There is therefore a need to establish new models in order to better describe the behavior of cakes during different solicitations and to be able to simulate this behavior. This allows the experiments to be well designed and the data to be better interpreted for a more reliable assessment of caking behavior. This is crucial for a better understanding of the phenomena and decoupling them to extract the intrinsic properties of the cakes. However, modelling of dispersed media in order to link microscopic phenomena occurring at the particle scale to their overall macroscopic behaviour remains still a major challenge in powder science and technology. In the case of wet agglomeration or caking, for example, given the complexity of phenomena occurring simultaneously at particle's scale (i.e., adsorption, absorption, capillary condensation, crystallization, diffusion, desorption, liquid or solid bridging, etc.) on the one hand, and the heterogeneity of the system (shape, size, composition, etc.), on the other hand, a detailed description of the media, even on a small scale, remains almost inaccessible. To this must be added the problem of the excessive number of particles present in real systems, which makes a complete description of the system almost impossible until now.

In the present work, a method based on DEM simulation of granular media based on use of multiple modelling approaches across length scales is developed to describe the mechanical resistance of wet or dry agglomerated (caked) samples. use of. The purpose of simulations was to provide a better understanding and interpretation of the experimental results taking into account the phenomena related to the heterogeneity of the medium as well as to have a reliable numerical tool for the prediction of the behaviour of cakes in the conditions other than those experimentally tested.

From a numerical cost point of view, both CPU time and memory depend strongly on the density of the simulated granular packing. Regarding the numerical simulations carried out in this study, the maximum number of particles that composed the agglomerated powder is in the order of 10,000 which does not lead a large computational time. Typically, one day in length for one simulation.

Firstly, the simulations were carried out on the tensile strength of wet cakes due to the presence of liquid bridges. The results for tensile strength of wet cakes calculated by a phenomenological model were used as a reference case to validate the consistency of simulation results. After validation, DEM simulations were used to describe the yield stress of cakes resulting from different attractive forces (liquid bridge, solid bridge, etc.). The validation of the model was carried out based on the experimental data of the authors published elsewhere [START_REF] Afrassiabian | A study on the caking behaviour of binary mixtures of lactose due to solid-state crystallisation of the amorphous phase[END_REF][START_REF] Afrassiabian | Multiscale investigation of caking phenomenon of polymorphic powders: from physico-chemical aspects to industrial applications[END_REF].

Brief reminder of previous works 1.Experimental data

Experimental data used in this study regarding the procedure of cake preparation and the post-caking tests in order to validate the model have already been published elsewhere [START_REF] Afrassiabian | A study on the caking behaviour of binary mixtures of lactose due to solid-state crystallisation of the amorphous phase[END_REF][START_REF] Afrassiabian | Multiscale investigation of caking phenomenon of polymorphic powders: from physico-chemical aspects to industrial applications[END_REF] and will not be detailed here. Very briefly, the caking experiments were conducted in a multi-cell caking device allowing to prepare cylindrical compacts (40mm ID, 26 mm height) under accelerated (humid air crossing the samples) and controlled conditions (T, RH, P, air flow rate). The mechanical resistance of cakes was then measured using different tests: uniaxial compression, diametral compression (Brazilian test) as well as shearing (translational) test.

Modelling

In a previous work, Afrassiabian et al. [START_REF] Afrassiabian | An overview of the role of capillary condensation in wet caking of powders[END_REF] established a phenomenological model of capillary condensation in granular media under humid conditions. This approach is based on the integration of water condensation and liquid bridge formation within binary contacts in regular or random packing of particles. The model was based on the combination of Kelvin law of capillary condensation and Laplace law for capillary forces. Briefly, the model predicts first the water activity, corresponding to a liquid bridge formed between two spherical particles with radii (Figure 1). Knowing the filling angle, , and the interparticle gap, ℎ, the geometry (principal radii of curvature, and and consequently the volume, ) of the liquid bridge can be obtained by the following set of equations:

= -(1 - ) + ℎ ( + ℎ) (1) = - [1 - ( + ℎ) ] (2)
= + - 1 2 -2 2 1 - 2 3 3 (2 -3 -cos 3 ) (3)
All symbols are illustrated in Figure 1.

Knowing and , the equivalent capillary (Kelvin) radius and the corresponding water activity can be calculated as follow:

1 % & = 1 2 ' 1 + 1 ( (4) ln
= -2 + • - % & • ℛ • / (5)
whereis the vapor-liquid surface tension, ℛ the universal gas constant, / the temperature, and + the molar volume of liquid.

Furthermore, according to Laplace theory for capillary forces, the total attractive force between the two particles including both the capillary and the surface tension components is

given by [START_REF] Rumpf | The strength of granules and agglomerates[END_REF]:

0 &12 = 4 -sin( ) 5sin( + 6) + 2 sin ' 1 + 1 (7 (6)
Finally, the tensile strength of a random packing of monosize particles can be estimated by Rumpf's model:

8 9 = 1 -: ; 2 • 0 (7)
where < is the mean void fraction of the bed, = is the average coordination number of the packing and 0 is the average force between two single particles. In the case of liquid bridges 0 = 0 &12 .

DEM simulation of mechanical strength of agglomerates 2.1 Basic concepts and equations

Discrete numerical modelling was first used to study the effect of capillary condensation in granular media at a Representative Element of Volume (REV). This amounts to solving the classical equations of dynamics for a particle assembly. This method has already been used to model the behaviour of granular materials in the presence of capillary forces [START_REF] Delenne | Compressive strength of an unsaturated granular material during cementation[END_REF][START_REF] Kim | DEM simulation of collapse behaviours of unsaturated granular materials under general stress states[END_REF]. The advantage of a discrete approach lies particularly in its ability to integrate physical quantities by being closer to the microstructure of the granular medium (at contact scale). Expressed in the global coordinate system, the set of equations originally developed by Cundall and Strack [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Golshan | Modelling Methods for Gravity Flow of Granular Solids in Silos[END_REF][START_REF] Norourzi | Coupled CFD-DEM Modeling: Formulation, Implementation and Application to Multiphase Flows[END_REF] is the following:

> ? @ ? A B C DE = 0 C FGH + I 0 J→C J L C 6 E = M C FGH + I M J→C J (8)
where DE and 6 E are linear and angular acceleration vectors, respectively, B C and L C are the mass and the moment of inertia matrix of the particle , 0 J→C and M J→C the force and moment vectors due to the action of the particle N on the particle , and 0 C FGH and M C FGH the external force and moment vectors.

Mechanical contact interactions are calculated with an explicit contact law [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. The particles were assimilated to rigid spheres. The contact forces in the absence of capillary effect were presented with a linear contact model depending on elastic force displacement law, Coulomb friction and viscous damping coefficient. For a pair of particles in contact ( ; N) the normal and tangential forces, 0 P,H , expressed in the coordinate system related to the contact plane were provided by the Kelvin-Voigt model:

> ? @ ? A 0 P,H = -= P,H R P,H -P,H S P,H T Uℎ ∶ 0 H = -B (|0 H |, X|0 P |) × Z (S H ) = H = P ∈ 5 2 3 , 17
P,H = 2\= P,H B * (9)
where = P,H and P,H are the normal and tangential stiffness of contact and damping coefficients, respectively, and μ is the dry friction coefficient of Coulomb. R P,H and S P,H are the relative displacements and velocities, respectively, according to the normal and tangential directions. The repulsive force, i.e., normal force, is obtained according to Hertz theory where the normal stiffness = P is function of the elastic moduli and the effective radius of both particles in contact [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. The normal stiffness

= P ^=P = 2 _ ``\ F``RP a is proportional to R P / inducing a nonlinear normal force ^0P ∝ R P c/ a, where _ F`` ' d eff = gν h i d h + gν j i d j
( is the effective modulus of elasticity function of Young's moduli, _ C and _ J and the Poisson's ratios,

ν C and ν C , of the two contacting particles. F`` is the effective radius ' F``= k h k j k h lk j (, with C
and J the radii of particles in contact and N, respectively. The tangential force is calculated by means of the tangential stiffness = H ^ = H = 8 n F``\ F``RP a, where n F`` is the effective shear modulus.

The mechanical system requires an explicit time integration based on the velocity-Verlet scheme:

o D(U + ∆U) = D(U) + ∆U Dq (U) + ∆U 2 DE (U) Dq (U + ∆U) = Dq (U) + ∆U 2 (DE (U) + DE (U + ∆U) (10)
DEM simulations were carried out using the MULTICOR-3D and MULTICOR3D++ discrete element codes [START_REF] Martin | Simulation of sintering using a Non Smooth Discrete Element Method. Application to the study of rearrangement[END_REF][START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF]. The method allows describing both the case of liquid capillary and solid bridges.

Discrete simulation of liquid bridges

In the presence of capillary bridges, the normal component of the contact force 0 P is the resultant of repulsive and attractive interactions:

0 P = -= P R P -P S P -0 &12 (11)
The attractive interaction, 0 &12 , due to the capillary effect is obtained according to Rumph's model [START_REF] Rumpf | The strength of granules and agglomerates[END_REF] taking into account the filling angle, , and the principal radii of curvature and . More specifically, at each pair of particles with interparticle gap 2ℎ, the filling angle is calculated iteratively using the Newton-Raphson algorithm with less than 50 iterations and an accuracy of 10 -6 rad. Indeed, since the number of capillary bridges may change during the simulation, this leads to run the contact detection algorithm at each time step. This of course means that the filling angle for each contact has also to be updated appropriately.

Nevertheless, the numerical cost of the capillary bridge updating can be considered as minor in comparison to the contact detection cost. By running the Newton-Raphson algorithm, the filling angle is obtained with less than 50 iterations and an accuracy of 10 gs % . Moreover, as well-known, the contact detection process is by far the costliest step in discrete element method. Once the filling angle is calculated, the principal radii of curvature are determined and the attractive force is then calculated. It should be pointed out that ℎ and R P are the same and both express the difference between the centre to centre distance separating two particles and 2 (Figure 1).

In this study, it is assumed that the tangential effect of the capillary bridge is negligible and that, in addition, the wettability is perfect (contact angle 6 = 0). The viscous effect of liquid was also neglected in the modelling because the main liquid involved in caking is water, which is not very viscous. The capillary force at each contact point was then calculated using Eq. (6).

Discrete simulation of solid bridges (Beam model)

This section deals with the mechanical behaviour of powders containing solid bridges.

In order to take into account the cohesive effect of solid bridges at the level of elementary contacts in discrete modelling, we have introduced an Euler-Bernoulli type beam model between particles in contact. This approach has been developed by LTI/UPJV and have been described in details elsewhere [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF][START_REF] Haddad | A DEM-FEM coupling based approach simulating thermomechanical behaviour of frictional bodies with interface layer[END_REF][START_REF] Leclerc | On the suitability of a Discrete Element Method to simulate cracks initiation and propagation in heterogeneous media[END_REF].

In this model, the solid bridges are assimilated to beams linking the particles together.

The local (or microscopic) parameters that characterize the beam model are _ t , n t , u t and v t , respectively the Young's modulus, the shear (or Coulomb) modulus, the cross-section of the beam and the length of the beam (Figure 2). In a matrix form, the generalized forces of cohesion induced by the beam are calculated via a matrix of stiffness (Eq. 12) [START_REF] Gagneux | Analytic Calculation of Capillary Bridge Properties Deduced as an Inverse Problem from Experimental Data[END_REF]. At the end of each beam element ( ; N) are affected the degrees of freedom in translation, D P C,J , D H C,J , D C,J , and in rotation 6 P C,J , 6 H C,J , 6 C,J , expressed in the local base of the beam (D wx P , D wx H , D wx `). The subscripts , U and y are related to the normal, tangential and bending effects, respectively.

The cohesion forces between two particles, and N, are represented by the following system: we saw in the previous section, such bridges cause an attractive force between particles, leading to the formation of a cake with intrinsic physicochemical and mechanical properties.

                          - - - -                                   - - - =                     → → → → → → j f i f j t i t j n i n j f i f j t i t j n i n 2 2 2 2 i j f i j t i j n i j f i j t i j n u u u u u u 6 3 0 0 0 0 2 0 0 0 6 3 0 2 0 0 0 0 0 0 0 0 0 0 0 2 2 0 0 0 2 2 0 0 0 0 0 0 0 0 0 0 0 0 M M M µ µ µ µ L K L K L K L K L K L K S L K L K K L K L K K K t t
By considering a Representative Elementary Volume (REV) of 10.000 particles, the DEM was then performed by means of MULTICOR3D software taking into account the properties of the cake (degree of saturation) in order to establish relationships between the microscopic parameters and the macroscopic behavior (tensile strength) of cakes. The software uses the explicit (also called smoothed) formulation [START_REF] Lubachevsky | Geometric properties of random disk packings[END_REF] which enables to fulfil the assumptions of the Random Close Packing (RCP).

DEM simulations were based on the combination of conventional DEM method and the phenomenological model of capillary condensation, presented in the previous section. At each step of time, capillary forces at each binary contact were calculated by assuming that the system is in thermodynamic equilibrium (aw=RH). The volume of liquid bridges and the extent of capillary forces within a known arrangement of particles at a desired water activity was calculated by the combination of the Kelvin law of capillary condensation and the Laplace law for capillary forces, as explained earlier:

-Initially, the degree of saturation in water was set at a given value, the corresponding volume of liquid was distributed equitably over the potential contacts (0 < 2ℎ * < 0.3), with ℎ * = " k . -At each time step, the geometric parameters of each capillary bridge ( , , ) were obtained by an iterative calculation.

-Finally, the redistribution of the liquid volume for the following time step was calculated based on the evolution of the contacts.

Figure 3 presents a snapshot from MULTICOR3D software, showing four binary contacts for different dimensionless separating distances (2ℎ * = 0.02; 0.08; 0.20; 0.30). At each contact, the capillary bridge shape is drawn with respect to their respective geometric parameters, namely the principal radii, , and the filling angle, .

Densification of regular packings under the effect of capillary bridges

Initially, the implementation of DEM allowed to model the effect of capillary condensation on regular packing of particles (CS, CC and CFC), with a particle radius of R=1.5

μm. The objective was to simulate the densification of stacks initially generated with noncontact particles (0 < 2ℎ * < 0.3). After having fixed the degree of saturation, the creation of capillary bridges between particles mechanically leads to a densification of the stack. At the end of this step, the theoretical values of compactness are verified respectively for the different studied cases, namely 52% for the CS, 68% for the CC and 74% for the CFC packing (Figure 4 provided by the discrete model.

Mechanical behaviour of powder compacts in the presence of capillary bridges

In a second time, the DEM was used to characterize the behaviour of an assembly of monodisperse particles in order to study the influence of the particle size and the degree of saturation on randomly generated packing of particles. The characteristics of the powder compacts and the tensile test parameters are summarized in table 1. The granular packing was generated by the efficient Lubachevsky-Stillinger Algorithm (LSA) [START_REF] Lubachevsky | Geometric properties of random disk packings[END_REF]. Considering a constant box volume, (10 g… ) c B c , six packings are generated with a number of particles varying in the interval [2 × 10 c ; 2 × 10 … ], a particle radius ranging between 2 μm and 4 μm, a constant compactness (≈50%) and a density of 1500 kg.m 3 . The degree of saturation, z, varies within the interval [0.64%; 12.8%].

The numerical characterization of the behaviour was carried out through a simple tensile test with an imposed speed of 10 -4 m.s -1 and a time step of 10 -7 s. The six stacks were submitted to simple tension in order to highlight a possible influence of their behaviour as a function of particle size. It should be reminded here that the initial volume of the domain is constant regardless of the number of particles in the stack. At the beginning, when the capillary links at the contacts are created, the packing is in a state of agitation (Figure 5-a).

Therefore, before applying the tensile stress, a stabilization phase is necessary in order to reach the equilibrium state. Figure (5-a) provides information about the densification process of the stack with a kinetic energy close to zero when the equilibrium state is reached.

Figure (5-b) and (c) present the capillary bridges and displacement field, respectively, at the end of the densification process. Once equilibrium state is achieved, the packing is subjected to a quasi-static tensile stress. During the test, the tensile stress, σ (Pa), is evaluated in the middle section of the packing to avoid edge effects, and the strain, ε (%), is estimated with respect to the particle radius.

The curves plotted in Figure 6 show the variation of the time-weighted tensile stress as a function of the relative elongation of the packing with a water saturation rate of 10.4%. The first observation concerns the two phases that characterize the first moments of the tensile test. Indeed, the traction curve has a very steep first phase, which reflects the reaction of the stack to the load. This sharp increase is followed very quickly by a second phase of low stiffness.

The second observation is related to the relative sensitivity of the maximum stress to particle size, even though the tensile curve is less and less erratic as the number of particles increases. Finally, the stack containing 20.000 particles is the one whose tensile curve shows a decrease in the slowest stress. This behaviour is partly explained by a high density of capillary links in the right section of the stack. This results in a wider flow area beyond the maximum stress.

Figure 6-b shows the evolution of a stack of 8000 particles at different times of the numerical tensile test. The wall placed at the middle of the stack allows computing with a good accuracy the tensile stress during the simulation in order to avoid a possible disturbance near the upper and lower boundaries. We note that there is no interaction at all between the particles and the middle wall as shown in figure 6-b. To avoid the boundary effects, a fictitious wall is considered, the role of which is only to localize those particles at the middle of the sample that are subjected to the tensile force.

The second study focused on the influence of the water saturation rate on the tensile strength. The number of particles was 10.000, which is a good compromise between accuracy and calculation time. The reference saturation rate was S = 0.64%. The test protocol is identical to that described above. Figure 7 shows a clear influence of the saturation rate on the tensile strength of the stack. The results show the evolution of the coordination number Z (mean number of capillary bridges by particle). This result was obviously expected, since a greater quantity of water increases both the intensity and the number of bridges. Indeed, this results in a coordination number, ;, which increases as a function of the saturation rate (Figure 7-b).

This evolution is in agreement with the results of tensile tests. Moreover, during the upward phase of the test, the parameter Z decreases until the tensile force reaches its maximum. In fact, during this phase, the particles may lose capillary bridges because of the increasing of the gap between some particles and its neighbours (2ℎ * > 0.3). Conversely, in the second phase of the test, when the sample is broken, the particles are no longer subjected to the tensile force, thereby facilitating the reactivation of new capillary bridges, thus increasing slightly the coordination number Z.

3.4 Mechanical behaviour of powder compacts in the presence of solid bridges to link the macroscopic properties of a powder compact to the local parameters of the strong interparticle links. Thus, it is possible from numerically obtained abacuses to identify the microscopic properties of the beam element in order to reproduce the behavior of a caked powder whose macroscopic properties have been experimentally characterized [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF][START_REF] Leclerc | On the suitability of a Discrete Element Method to simulate cracks initiation and propagation in heterogeneous media[END_REF]. At the end of the calibration process, local properties are thus assigned to the solid bonds of a caked powder, in order to reproduce numerically the desired macroscopic behavior.

Mechanical tests of indirect tensile (Brazilian test) and shear tests were carried out using the MULTICOR3D++ discrete element simulation tool [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF] on compacts of anhydrous lactose powder. The experimental characterization tests are detailed elsewhere [START_REF] Afrassiabian | Multiscale investigation of caking phenomenon of polymorphic powders: from physico-chemical aspects to industrial applications[END_REF].

The damage occurred during the mechanical tests in the specimen is taken into account by means of an appropriate failure criterion depending on the stress state in the particle. The cohesive links of a given particle fail when the material strength is reached. Two criteria are considered for the simulated mechanical tests. For Brazilian and shear tests, we have considered criteria based on hydrostatic stress and Mohr Coulomb model, respectively. In the case of discrete material, the virial stress does not characterize the measure for mechanical force between material points and cannot be considered as a measure for mechanical stress.

Hence, we determine the stress state at the particle scale by using an equivalent Cauchy stress tensor 8 ‡ 2 C . This tensor was proposed by Zhou [START_REF] Zhou | A new look at the atomic level virial stress: on continuum-molecular system equivalence[END_REF] for molecular dynamic simulations, and was used in DEM lately by Jebahi et al. [START_REF] Jebahi | Simulation of Vickers indentation of silica glass[END_REF]. The stress tensor is given by the following equation:

8 ‡ 2 C = 1 2 C I 1 2 ^ CJ ⊗ y CJ + y CJ ⊗ CJ a J ∈ ˆh (13)
with C the volume of the particle , ; C the coordination number of the particle , CJ the position vector connecting the center of the particles and N and y CJ the cohesive force vector acting on the particles and N.

The mechanical tests were conducted on samples of cylindrical shape, of diameter D=40 mm, thickness ‰ = 16 BB and mean void fraction = 0.64 . The simulated compacts consist of 10 4 particles, with a mean diameter of dm=675 μm. The indirect tensile and shear tests are controlled at a speed of amplitude of 1 mm.min -1 . The sample sizes and the test conditions were then identical to those imposed during the experimental campaign. The characteristics of the caked powder and the Brazilian and shear test parameters are summarized in table 2.

Brazilian test

Before performing numerically the brazilian test, we have adjusted _ t , for a given, % t ^%t ≈ 0.72a to calibrate the macroscopic mechanical behavior of the caked powder, in terms of macroscopic Young's modulus, _ • and Poisson ratio, ν • . The solid bridges at particle scale fail when the positive hydrostatic stress exceeds the uniaxial tensile strength of the material [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF]. However, the value of the tensile strength of the particle should be first adjusted so that the simulated macroscopic strength is close to that obtained experimentally. This first trial demonstrated the ability of the simulation tools developed to reproduce the mechanical behaviour of a cohesive multicontact system. The quality of the results obtained is, however, conditioned by the calibration step discussed previously, hence the importance of a rigorous identification of the mechanical and geometric parameters of the solid bridges. Taking into account the breaking of solid bridges, using a rupture criterion based on the hydrostatic stress measured at the particle scale [START_REF] Leclerc | On the suitability of a Discrete Element Method to simulate cracks initiation and propagation in heterogeneous media[END_REF], it was possible to simulate and follow the state of mechanical degradation of the sample at the particle scale.

Shear test

The second type of simulated behavior was that obtained during a shear test on a cylindrical compact as explained in [START_REF] Samain | Characterization of caking for crystalline materials: comparison and statistical analysis of three mechanical tests[END_REF]. The diameter of the sample is identical to that considered in the case of indirect traction. The particle size of the caked powder also remains unchanged, namely an average diameter of the order of 675 microns. The shear test was simulated with a speed of 1 mm.mn -1 during the mechanical test. In the case of shear test, the solid bridges created between particles fail when the Mohr Coulomb criterion is not fulfilled.

Here we just recall that the Mohr Coulomb criterion is expressed by means of a combination of shear and normal stresses, depending on the minor and major principal stresses. In order to obtain a good agreement between the numerical curve and the experimental one, during the upward phase of the shear test, the shear strength at the particle scale should first be adjusted so that the simulated macroscopic strength is close to that obtained experimentally.

Moreover, when the failure mechanism starts in the shear plane, the damaged particles are removed as well as the cohesive links between damaged particles and their neighbours, as proposed in Jebahi et al. [START_REF] Jebahi | Simulation of Vickers indentation of silica glass[END_REF] and Leclerc et al. [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF][START_REF] Leclerc | On the suitability of a Discrete Element Method to simulate cracks initiation and propagation in heterogeneous media[END_REF]. It should be noted that the failure criterion (Mohr-Coulomb criterion) depends on the stress field in the particle determined by means of the tensor formula suggested by Zhou [START_REF] Zhou | A new look at the atomic level virial stress: on continuum-molecular system equivalence[END_REF].

The results are shown in Figure 9. It appears that the shear behaviour obtained numerically partially reproduces that observed experimentally (Figure 9

,.3 Results and discussion 3 . 1

 31 and M J→C are the components of the vector of generalized forces of cohesion expressed in the local base of the beam and = P , = H and z P of the matrix are the local stiffnesses, with: { P = d | } | | L • is the quadratic moment. Note that the radius • of the cross section u • of the beam depends on a dimensionless parameter, % • , such that: • = % • k h lk j . The behaviour of the beam is therefore a function of the local Young's modulus _ • and the dimensionless parameter, % • ∈ ]0,1]. Discrete simulations of capillary effects at a REV scale First, we focus on the role of capillary condensation and subsequent liquid bridge formation within a granular medium exposed to fluctuations of ambient relative humidity. As

  -a). The results of Figure(4-b) show a homogeneous densification of the stack. The kinetics of densification of all the particles was similar to the behaviour of two isolated particles. These test cases are, of course, not representative of the actual packings, but they allow to be placed in an ideal case without rearrangement and thus to validate the results

3. 4 . 1

 41 Calibration of beam properties and stress state at particle scale The identification of the local properties of the beam, namely its normal and transverse rigidities (_ • ; n • ) are obtained by means of a calibration technique, which makes it possible

Figure 8 -

 8 a compares the numerical and experimental curves obtained in a Brazilian test. The mechanical response of the compact powder follows fairly well the experimental curve. In addition, the ultimate load reached by the sample is substantially close to that obtained during the experimental test. Also, Figure 8-b shows clearly a concentration of the tensile stress in the central zone of the sample. This type of mechanical test makes it possible to characterize the tensile strength of brittle materials. Finally, Figure 8-c qualitatively shows the presence of particles (particles coloured blue) whose interparticle bonds are broken.

2 CFigure 2 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 9 .

 22345679 Figure1. Schematic representation of a pendular liquid bridge between two equal-sized particles

Table 1 :

 1 Characteristics of the powder compacts for the tensile tests

	Sample dimensions

Table 2 :

 2 Characteristics of the caked powder for the Brazilian and shear tests

	Sample diameter

Acknowledgements

Join support by Région Haut-de-France (CR Picardie) and EU funding (FEDER) for MOTTAMORPH project are gratefully acknowledged.

Nomenclature

Degrees of freedom in rotation of the beam