Maël Guilleme

Simon Malinowski

Romain Tavenard

Xavier Renard

Mael Guillemé

Energiency

Localized Random Shapelets

Keywords: time series, machine learning, shapelets

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

Introduction

Time series classification has recently gained an increasing attention from researchers, due in particular to its possible application in various domains such as economics, agriculture, and health for instance. There are two main families of methods for that task: some deal with raw time series and use or design dedicated (dis-)similarity measures, while others rely on feature extraction to embed time series in metric spaces in which standard machine learning tools can be considered. The work presented in this paper is part of this latter category. Feature-based methods include works relying on hand-crafted features [START_REF] Bailly | Dense Bag-of-Temporal-SIFT-Words for Time Series Classification[END_REF][START_REF] Schäfer | The BOSS is concerned with time series classification in the presence of noise[END_REF] as well as learning-based approaches, among which the shapelet model plays an important role. Shapelets have first been introduced in [START_REF] Ye | Time series shapelets: a new primitive for data mining[END_REF]. They correspond to subsequences that are able to discriminate classes. The concept of shapelet Fig. 1: Comparison between shapelets extracted by the Learning Time-Series Shapelets (LS) algorithm and our Localized Random Shapelets (LRS) approach. This Figure has been generated using tslearn implementation of LS [START_REF] Tavenard | tslearn: A machine learning toolkit dedicated to time-series data[END_REF]. transform has then been proposed in [START_REF] Lines | A shapelet transform for time series classification[END_REF]. It consists in transforming time series into a vector whose components represent the similarity between the time series and shapelets that have been selected beforehand (or learned, as in [START_REF] Grabocka | Learning timeseries shapelets[END_REF]). After this transformation, time series are embedded in a Euclidean space where classifiers like Multi-Layer Perceptron or Support Vector Machines can be used. Techniques based on the shapelet transform are amongst the most accurate ones for time series classification (an interesting survey and comparison of time series classification methods can be found in [START_REF] Bagnall | The great time series classification bake off: a review and experimental evaluation of recent algorithmic advances[END_REF]). However, they have three main drawbacks. First, the step of generating (or learning) shapelets that lead to accurate classification is computationally demanding (especially when dealing with large time series datasets). Some works have been proposed in order to fasten this step, by relying on time series approximation [START_REF] Rakthanmanon | Fast shapelets: A scalable algorithm for discovering time series shapelets[END_REF] or by drawing random shapelets [START_REF] Wistuba | Ultra-fast shapelets for time series classification[END_REF][START_REF] Renard | EAST Representation: Fast Discriminant Temporal Patterns Discovery From Time Series[END_REF]. Second, no information about the localization of the shapelets in the time series is available after the transformation. Classical shapelet transform only makes use of the similarity between a shapelet and a time series. Information about when the shapelet occurs in the time series might be of importance to discriminate classes. This can be related to previous works that have shown that localization of extracted features in time series improve classification accuracy [START_REF] Tavenard | Efficient Temporal Kernels between Feature Sets for Time Series Classification[END_REF]. To better understand why, let us consider, for example, a sign language sentence recognition task. In this use case, being able to recognize salient patterns (gesture atoms) is key, but localizing them in the sentence is also important in order to understand the meaning of the sentence. Third, most shapelet-based models suffer from a lack of interpretability in the sense that (i) it is difficult to understand the impact of each shapelet on the final classification decision and/or (ii) extracted shapelets might not always be related to original time series. Figure 1 compares shapelets learned by the Learning Shapelet (LS) algorithm [START_REF] Grabocka | Learning timeseries shapelets[END_REF] on the EarthQuakes dataset and those extracted (at random) by our algorithm. Time series are rescaled to lie in the [0, 1] range before learning the shapelets. However, we can see that LS-extracted shapelets have very little in common with the time series and do not even fit in the [0, 1] range. In other words, shapelets cannot be seen as realistic time series patterns. For methods based on the shapelet transform [START_REF] Lines | A shapelet transform for time series classification[END_REF], shapelets are extracted from original time series. But they then feed an ensemble of classifiers, which makes it difficult to understand relationship between shapelets and class belongings.

In this paper, we propose a novel shapelet model that tackles these drawbacks. First, in order to reduce the computing cost of selecting the most discriminant shapelets, our model selects shapelets randomly from training data. Hence, extracted shapelets are, by definition, closely related to the time series. Moreover, this random shapelet framework allows us to easily take the shapelet localization information into account. Second, we propose a dedicated feature selection method called Semi-Sparse Group Lasso (SSGL) capable of either ignoring a shapelet, keeping only its distance information or using both the distance and localization information. Third, we show that obtained shapelets are meaningful and the resulting model can be easily analyzed to get insights about which are the important features (for the classification task) in the dataset and for each of these features, what kind of information (presence only or presence and localization) contributed to the decision. Overall, we show that we are able to reach competitive performance w.r.t. Learning Shapelets [START_REF] Grabocka | Learning timeseries shapelets[END_REF] (even outperforming this baseline when larger datasets are considered) with a more interpretable model. The rest of this paper is organized as follows. Our Localized Random Shapelet model is detailed in Section 2 and its interpretability is discussed in Section 3. Section 4 evaluates the benefit of the proposed model on time series classification.

Localized Random Shapelet Model

Background on shapelets and shapelet transform

A shapelet S = s 1 , . . . , s l is a temporal sequence (that can be extracted from existing time series or not). Given a time series T = t 1 , . . . , t L , the distance between s and T is defined as :

d(T, S) = min 1≤j≤L-l+1 l i=1 (s i -t i+j-1 ) 2 . (1)
In other words, euclidean distances between s and every subsequence of T (of length l) are computed and only the best match (minimum distance) is kept.

Given a set S = {S 1 , . . . , S K } of K shapelets, the shapelet transform of T is defined as the vector v 1 , . . . , v K such that v k = d(T, S k ) for all k. The original way to select a set of shapelets for a classification task is to evaluate the discriminatory power of a shapelet candidate by using for instance the information gain [START_REF] Ye | Time series shapelets: a new primitive for data mining[END_REF][START_REF] Lines | A shapelet transform for time series classification[END_REF], and to keep the ones with the higher gain. A strategy based on learning the shapelets that minimize an objective function was proposed in [START_REF] Grabocka | Learning timeseries shapelets[END_REF]. Other works consider the use of random shapelets and then rely on classical feature selection algorithms together with the classification step. In [START_REF] Renard | EAST Representation: Fast Discriminant Temporal Patterns Discovery From Time Series[END_REF], it has been shown that with this idea, a few thousands subsequences are enough to reach state-of-the-art classification performance on a standard benchmark [START_REF] Bagnall | The uea & ucr time series classification repository[END_REF].

Localized Random Shapelet model

In this framework, no information about the localization of the shapelets in the time series is used, while it has been shown that this kind of information helps improving classification performance [START_REF] Tavenard | Efficient Temporal Kernels between Feature Sets for Time Series Classification[END_REF]. In this section, we explain how we can integrate such information in the shapelet transform framework and derive a feature selection algorithm that keeps localization information only when needed.

In our Localized Random Shapelet (LRS) model, each shapelet S is drawn uniformly at random from the set of all training time series snippets. Each shapelet leads to two features for each time series T . The first feature is the same as in the classical shapelet transform, i.e. the shapelet distance1 d(T, S) between T and s as defined in Equation ( 1). The second feature corresponds to the first time instant at which this distance is reached. It is computed as

l(T, S) = argmin 1≤j≤L-l+1 l i=1 (s i -t i+j-1 ) 2 .
(

) 2
Let T = T 1 , . . . , T N be a dataset of N time series. This set is transformed by the localized random shapelet model into a feature matrix X, such that:

X =    d(T 1 , S 1 ) l(T 1 , S 1 ) • • • d(T 1 , S K ) l(T 1 , S K ) . . . . . . . . . . . . d(T N , S 1 ) l(T N , S 1 ) • • • d(T N , S K ) l(T N , S K )    . (3)
Once this feature matrix computed, we feed it to a standard multi-layer perceptron classifier as shown in Figure 2. This model takes as input a 2×K-dimensional vector and outputs probabilities for the different classes at stake. The number of hidden layers can be adapted to application needs. We denote the whole set of parameters of this model by θ, and inside this set of parameters, we isolate the parameters of the first layer in the model and denote them β.

β (k) β (k)
Fig. 2: Overview of our localized random shapelet model. Blue circles indicate distance features while orange ones correspond to location features. For each shapelet, a group is formed whose weights are denoted β (k) (where k is the shapelet index). Note that the number of hidden layers may vary from one application to the other.

Structured feature selection

As explained above, the proposed shapelet model relies on random shapelets (taken from the original time series) in order to fasten the shapelet generation step. When dealing with such shapelets, a feature selection strategy should be applied before (or jointly with) the classification step, in order to simplify the resulting representation, which tends to improve overall accuracy [START_REF] Renard | Random-shapelet : an algorithm for fast shapelet discovery[END_REF]. It can be seen in Equation (3) that the extracted features are structured. Distance features tell how well a shapelet matches a time series, while localization features inform about the location of the match. Classical feature selection strategies are hence not adapted to this kind of feature matrix. Indeed, it does not seem reasonable to exclude a distance feature related to a shapelet while keeping the associated localization feature. However, removing a localization feature and keeping the associated distance feature might be meaningful in cases where the localization of a shapelet does not impact the class belonging of the times series. In the following, we design a feature selection strategy adapted to the kind of features extracted from the localized random shapelet model. This strategy, based on regularization, is described below.

Regularization strategy In the following, we assume that we have a prediction problem with a loss function L to be minimized. This loss function is computed over a dataset X of N observations associated with a target vector y. A standard approach to bias the learning process towards sparse solutions is to derive a regularized loss function, as done for Lasso regression [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]:

L Lasso (X, y, θ) = L(X, y, θ) + λ β 1 (4)
Simon et al. [START_REF] Simon | A sparse-group lasso[END_REF] introduced Sparse-Group Lasso (SGL), a more structured regularization scheme that could take feature group information into account. The resulting regularized loss function for a K-group problem is:

L SGL (X, y, θ) = L(X, y, θ) + αλ β 1 + (1 -α)λ K k=1 √ p k β (k) 2 (5)
where p k is the number of features in group k and β (k) is the sub-vector of β made of features from group k. Here, the β 1 term enforces per-feature sparsity while β (k) 2 pushes towards sparsity at the group level. The α parameter hence acts as a trade-off between these two regularization terms.

We consider a slightly different setting in which, inside a group, only part of the features can be dropped. In our case, each group corresponds to a shapelet. For each shapelet, two features are available: one for the distance and one for localization of the match. We aim at designing a strategy that can, for each shapelet (or group), either:

drop all information related to that shapelet (if the shapelet is useless for prediction), -keep only the distance information (if the location of the shapelet does not help for prediction), -keep both the distance and localization information.

In order to meet these constraints, we introduce the Semi-Sparse-Group-Lasso (SSGL) framework, which consists in minimizing the following loss function:

L SSGL (X, y, θ) = L(X, y, θ) + αλ M ind θ 1 + (1 -α)λ K k=1 √ p k β (k) 2 (6)
where M ind is an indicator diagonal matrix made of ones and zeros, the latter corresponding to variables that should be kept as soon as the group is not zeroedout (i.e. variables that will not be considered for 1 regularization). In the case of our localized shapelet model, M ind has a diagonal that alternates between zeros for dimensions corresponding to distances and ones for those related to the localization information. Optimization in practice In practice, minimization of such a regularized loss function can be tackled by several means. First, building on [START_REF] Simon | A sparse-group lasso[END_REF], one can derive a block-wise procedure that considers feature groups one at a time and starts by deciding on whether it should be zeroed-out or not. In the case of an ordinary least square regression setting, it gives the following sufficient condition for zeroing-out group k (i.e. dropping all the information from the corresponding shapelet):

β (1) 1 β (1) 2 β (2) 1 β (2) 2 β (3) 1 β (3) 2 10 -4
M ind S X (k) r -k N , αλ 2 ≤ √ p k (1 -α)λ, (7)
where X (k) is the submatrix of X in which only variables from group k are kept, r -k is the partial residual of y, substituting all group fits other than group k and S(•, •) is the coordinate-wise soft thresholding operator:

(S(z, αλ)) j = sign(z j )(|z j | -αλ) + . (8)
If this condition is satisfied, all β (k) coefficients are set to zero and the process goes on to the next group. Otherwise, optimization of the coefficients inside group k should be performed, either using subgradient equations in a coordinate-wise algorithm or by performing gradient descent steps inside group k.

In all the experiments presented in this paper, we rather use full gradient descent on the regularized loss, in order not to face known limitations of blockwise gradient descent algorithms such as slow convergence and low parallelism capabilities. This strategy was already successfully applied in [START_REF] Scardapane | Group sparse regularization for deep neural networks[END_REF].

Toy example: structured linear regression Figure 3 presents a comparison of the three regularization schemes described in this section. For this comparison, we used the following model to draw samples: where X and ε are drawn from centered normal distributions of standard deviation 1 and 0.01 respectively and β = β

y = X • β + ε, (9)
(1) 1 , β (1) 2 , β (2) 1 , β (2) 2 , β (3) 1 , β (3) 2
has only three non-zero components:

β (2) 1 , β (3) 1 and β (3) 2
(cf. Figure 3). We assume to have a problem similar to the shapelet setting, i.e. we have groups of two variables among which only β (k) 2 is concerned by 1 -norm regularization (for group k). This structural information is used for SGL and SSGL variants in the experiments, while Lasso is blind to such structure.

For this example, we use a simple model without any hidden layer. First, SSGL outperforms both SGL and Lasso in terms of mean squared error (MSE), showing the benefit of taking variable structure into account in the model. Second, from a more qualitative perspective, the structure of the ground truth coefficients is better preserved with SSGL (all three null coefficients have low value estimators and β (1) 2 is better estimated thanks to the intra-group specific regularization scheme).

Model interpretability

In this section, we illustrate the interpretability of our method through a simple use-case. We consider the TwoPatterns dataset from the UCR & UEA time series classification repository. TwoPatterns is a synthetic dataset in which each time series is made of two pattern occurrences surrounded by noise. There are two different patterns (named A and B in the following) in the dataset and the classification problem is hence made of four classes, corresponding to all possible permutations of patterns A and B, as shown in Fig. 4. For this illustration, we rely on a simple variant of our model in which we have no hidden layer and draw K = 2, 000 shapelets to build our feature space. We show that, once trained, our model can be easily analyzed by scrutinizing its weights. More precisely, in the following, we aim at answering the following questions:

-Which shapelets are the most important for classification? -For each of these shapelets, is the localization information important for the decision ?

Note that, in this simple setting where no hidden layer is used, this information comes in a straight-forward manner by analyzing weights, but the same analysis could be performed in the multi-layer case through gradient descent in the shapelet transform space.

First, we can extract most important shapelets (in terms of classification power) by ranking them with respect to the 2 -norm of their associated weight matrix β (k) . Top-4 shapelets for dataset TwoPatterns are presented in Fig. 5. They fully match expectations since they focus on discriminative parts of the time series and cover both patterns A and B, as well as successions of these.

Then, for a shapelet that is considered discriminant, we can wonder whether its localization and distance are both used by the model or not. To assess the importance of the distance feature for a shapelet, we compute the 2 -norm of the coefficients of β (k) that are associated with its distance feature. We call this value distance coefficient. Similarly, to assess the importance of the localization feature for a shapelet, we compute the 2 -norm of the coefficients of β (k) that are associated with its localization feature. We call this value localization coefficient.

The first shapelet of Figure 5 has the 4 th highest localization coefficient amongst the 2,000 shapelets, together with the 14 th highest distance coefficient. This means that for this shapelet both the distance and the localization are important. It is coherent as this shapelet corresponds to a B pattern. The distance feature enables to discriminate class 1 from the others, while the localization features enables to discriminate class 2 and class 3. This is confirmed by the histograms of Figure 6. The second shapelet of Figure 5 has the 2 nd highest distance coefficient amongst the 2,000 shapelets. However, its localization coefficient is not as important as its distance coefficient (217 th out of 2,000). This means that for this shapelet the distance feature is important but not the localization. This is coherent as this shapelet corresponds to the end of a pattern A followed by the beginning of a pattern B, which only occurs in class 3, hence distance feature is enough to discriminate class 3 from the others, as confirmed by Figure 6. Similar conclusions can be drawn for shapelets 3 and 4.

We have seen in this section that our model allows easy extraction of the shapelets that have more contributed to discriminate classes. One can also easily analyze whether the localization of these shapelets was used by the model. This kind of information is very important as it brings a richer interpretability in the decision process. Conversely, state-of-the-art shapelet-based algorithms tend to suffer from lack of interpretability. Shapelet Transform is a weighted ensemble of several standard classifiers, which makes it difficult to assess the importance of shapelets in the decison process. In the Learning Shapelet (LS) framework, obtained shapelets are not constrained to be similar to training time series (cf. Figure 1), which limits the interpretability of such framework.

Experiments

In this section, we present experimental results to assess the performance of the proposed Localized Random Shapelet (LRS) model and compare it to state-ofthe-art shapelet-based methods.

Experimental setup

The proposed LRS is compared to the following state-of-the-art baselines: Shapelet Transform (ST) [START_REF] Lines | A shapelet transform for time series classification[END_REF], the Learning Shapelet (LS) model [START_REF] Grabocka | Learning timeseries shapelets[END_REF] and the Fast Shapelet (FS) one [START_REF] Rakthanmanon | Fast shapelets: A scalable algorithm for discovering time series shapelets[END_REF]. Experiments are based on the datasets from the UCR & UEA time series classification repository [START_REF] Bagnall | The uea & ucr time series classification repository[END_REF], a classical benchmark in the time series classification field. It is composed of 85 datasets with diverse time series classification problems. A description of the datasets can be found in [START_REF] Bagnall | The great time series classification bake off: a review and experimental evaluation of recent algorithmic advances[END_REF]. Classification performance for baseline methods is retrieved from the website associated to this dataset.

For each dataset, we draw K = 2, 000 shapelets uniformly at random from the time series training set. For the length of the shapelets, we followed the procedure described in [START_REF] Grabocka | Learning timeseries shapelets[END_REF]. All time series are represented by a feature vector which size is twice the number of shapelets (cf. Eq. ( 3)). This feature vector embeds information about shapelet distance and localization.

We use a 3-hidden-layer model with 256, 128 and 64 units in the hidden layers and each internal layer is followed by a batch normalization layer as suggested in [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF]. We add some dropout in the second and third layer. We use ReLU activations [START_REF] Glorot | Deep sparse rectifier neural networks[END_REF] for all internal layers and softmax activation for the final layer. The regularization in LRS requires two more hyper-parameters: λ that controls the regularization strength, and α that controls the trade-off between lasso and group-lasso terms. We cross-validated the dropout from the set {0.0, 0.3, 0.5, 0.7}, λ from the set {10 -1 , 10 -2 , . . . 10 -7 } and α from the set {0.1, 0.3, 0.5, 0.7, 0.9}. Models are learned using RMSProp optimizer [START_REF] Tieleman | Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude[END_REF] with a learning rate of 0.001.

Runtime and memory cost

Experiments are run on a Laptop with a Fedora 24 system, 16 GB DDR4 RAM and dual core CPU (i7-6600U 2.6 GHz). Python code used for these experiments is publicly available 2 . As an indication of the low computing and memory cost of our method, it requires 66 seconds and 1.5 GB of memory to learn a model (for a given set of hyper-parameters) on the Electric Devices dataset that is the largest in the UCR & UEA archive.

Evaluation of the impact of the SSGL regularization

Figure 7 shows the error rates on 85 datasets for two different regularization strategies: Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] and SSGL. It can be seen that the SSGL regularization strategy (that allows selection for each shapelet of either distance feature, or distance and localization features, or none of these) slightly improves classification performance over Lasso (that makes the selection for each features independently).

Performance comparison against baselines

Critical diagrams show the average ranks of the classifiers in order, and cliques, represented by a solid bar. A clique represents a group of classifiers within which there is no significant pairwise difference. Figure 8 presents critical diagrams of the performance against the baselines when considering all datasets, and when considering only datasets with more than 300 training instances. Overall, in terms of classification performance, only Shapelet Transform outperforms our approach. As already stated, ST is a weighted ensemble of several standard classifiers, hence suffering from lack of interpretability about the decision process. The performance of LRS is slightly better than LS (but not significantly). When we only consider large datasets (the 42 datasets with more than 300 training instances, which are better suited to the characteristics of our method), performance of LRS is better than LS and FS and gets closer to that of ST. This seems to indicate that taking the localization information into account enables to improve classification performance provided that sufficient training data is available to correctly fit the model.

Conclusion

In this paper, we have proposed a novel shapelet model that has low computing cost, competitive performance and better interpretability compared to tradi-tional shapelet-based methods. This model uses shapelet localization in addition to the traditional shapelet transform representation in a random shapelet framework. As such a framework needs many shapelets as input, we have designed a dedicated hierarchical regularization framework (SSGL) to fit our application needs. Experiments show that our model produces a ranking of realistic shapelets, and is able to explain their importance both in terms of distance and localization.This enables an easy and rich interpretability of the classification process. Moreover our classification performance is competitive with state-ofthe-art shapelet-based classifiers.

Fig. 3 :

 3 Fig. 3: Coefficients learned using different regularization schemes for a linear regression problem. Ground-truth coefficients are reported in blue.

Fig. 4 :

 4 Fig. 4: An example of each class of TwoPatterns dataset.

Fig. 5 :

 5 Fig.5: Four most important shapelets (in red) extracted by our method from the TwoPatterns dataset.

Fig. 6 :

 6 Fig. 6: The distribution of the localization (left) and distance (right) for the most important (first row) and second most important (second row) shapelet in TwoPatterns dataset.

Fig. 7 :

 7 Fig.7: Error rates comparison on 85 UCR datasets between LRS with ssgl regularization against LRS with lasso regularization.

Fig. 8 :

 8 Fig. 8: Critical diagrams of the performance against the baselines.

Note that the term distance is an abuse of notation since d(T, S) is not a distance, mathematically speaking.

https://github.com/rtavenar/localized_random_shapelets


T Xiong

H Ilmi

Y Hoarau

P Choquet

C Goetz

A Fouras

S Dubsky

M Braza

S Sainlos-Brillac

F Plouraboué

D Lo Jacono

Flow and particles deposition in anatomically realistic airways

Keywords: particle deposition, rat lungs, rabbit lungs, numerical simulations

Introduction

The aerosol particles present in our environment are identified as increasingly risk factors for the human health. In particular, small particles that reach the surface of gas exchange in the alveolar region of lungs are considered most harmful to health. Nevertheless, aerosols can also be used for diagnostic or therapeutic purposes to treat lung diseases. Inhalers, for example, are commonly used to administer drugs. In cases of lung disease, administration of medication by aerosol has the advantage of reaching the region directly affected and requires lower doses (side effects are minimized). The morphological descriptions of the bronchial tubes of the man appeared in the literature since 1963, when Weibel [START_REF] Weibel | Morphometry of the human lung[END_REF] introduced a physical symmetric model. The human lung is a sequence of bifurcations in two parts. It has an average of 23 generations of bifurcations which represents about 17 million of bifurcations. The last 5 generations are covered with cellular structures through which the exchange takes place with the blood gas [START_REF] Weibel | The pathway for oxygen, structure and function in the mammalian respiratory system[END_REF]. The understanding of the flow structures and the particle transport/deposition across the human bronchial system remains a challenge to achieve because of the complexity of the geometry of human lungs. This work relies a strong collaboration between physicians, medical imaging researchers, fluid mechanics researcher and CFD researchers.

Methods

Our study focuses on the numerical modeling of the flow and particles deposition in airways using the commercial package CFD-ACE. The objective of this study is to propose a reliable and complete numerical model of transport of particles in the pulmonary flow by taking into account the full breathing cycle. For the moment the movement of the airways imposed by the breathing are not taken into account.

This study is divided into four parts :

 First the CFD-ACE package is validated regarding flow [START_REF] Zhao | Steady inspiratory flow in a model symmetric bifurcation[END_REF][START_REF] Comer | Flow structures and particle deposition patterns in double bifurcation airway models, part 1. air flow fields[END_REF] and particles deposition in a simple model [START_REF] Kim | Deposition characteristics of aerosol particles in sequentially bifurcating airway models[END_REF][START_REF] Comer | Flow structures and particle deposition patterns in double bifurcation airway models. part 2. aerosol transport and deposition[END_REF]]  Simple geometries derived from the Weibel model are studied  An anatomically realistic model of rat lung is simulated  A more complex geometry of rabbit lung is simulated

The particles trajectories and deposition are obtained by solving a lagrangian transport equation where only the drag force is retained. The deposition efficiency is computed as the ratio of the amount of particles sticked at the wall over the total number of released particles. The rat geometry is obtained after an adult rat killed by IV administration of barbiturics, was infused with a iodine contrast agent (10% v/v, Visipaque320, GE Healthcare, Chalfont St Giles, UK). Imaging of the thorax was performed ex vivo using an X-Ray μCT scanner (eXplore Vision 120, GE Healthcare, Waukesha, USA). The protocol used involved 360 views over 360°, with 1 frame averages, at 100kV, 50mA. After reconstruction using a Feldkamp algorithm of backprojection, the imaging volume was made of cubic voxels of 49μm x 49μm x 49μm. Visualisation was performed using MicroView (GE Healthcare, Waukesha, USA). After invertion of grey levels, a surface rendering (marching cubes algorithm) was applied after thresholding. The surface was then exported in stl format and prepared in gmsh before being meshed in CFD-GEOM.

The rabbit geometry was obtained using synchrotron-based dynamic computed tomography. Details of the experimental technique and segmentation protocol are provided elsewhere [START_REF] Dubsky | Synchrotron-based dynamic computed tomography of tissue motion for regional lung function measurement[END_REF]. Briefly, newborn rabbit pups are ventilated and imaged using propagation-based phase-contrast imaging. The lungs are an ideal sample for phasecontrast imaging , as they exhibit large phase boundaries at air-tissue interfaces. *Corresponding author. Email: hoarau@unistra.fr This technology can provide rapid imaging of the lungs without contrast agent [START_REF] Fouras | Altered lung motion is a sensitive indicator of regional lung disease[END_REF]. By synchronising the image acquisition with the ventilation and rotating the sample, four-dimensional (4D=3D+time) images of the bronchial tree can be achieved. Again, one stl geometry was prepared in gmsh before being meshed in CFD-GEOM. Most of this work was performed during H. Ilmi PhD thesis [START_REF] Ilmi | Modélisation numérique des écoulements pulmonaires[END_REF], the remaining part was studied by T. Xiong during her master 2 internship.

Results and Discussion

Concerning the Weibel's generic models, the solutions generated by the solver CFD-ACE were in very good agreement with the literature and we evidenced the importance of the Dean vorticies on both the flow and the particles deposition. The deposition is essentially of inertial type and the deposition efficiency increases with both the Stokes number and the Reynolds number. The deposition occurs mainly at bifurcations. We have simulated different geometries where the successive bifurcations were rotated by 90° and we observed a strong influence of the geometry : at the first bifurcation the Dean vorticies are properly splitted in each bifurcation and the flow rate in the inner tube is much higher but at the second bifurcation the Dean vorticies are broken and the flow is more chaotic, the particle deposition is also significantly higher.

Fig. 1 : deposition in rat lungs

The simulations of the rat and the rabbit lungs raised the major difficulty of this study : the appropriate grid generation. Very fine grids were necessary to obtain grid convergence on the particles deposition and such grids is very difficult to obtain in this kind of complex geometries. The rat geometry is quite linear and no recirculations were found. The particles deposition is again mainly inertial, increases with the Stokes number and occurs in the lower bifurcations (Fig. 1) because they are pushed toward the wall after each bifurcation. The rabbit geometry is much more complex and a large and strong vortex appear near the inlet because of the intubation which causes abnormal particle deposition at the inlet (Fig. 2). Dean vorticies were evidenced in few bifurcations. The deposition occurs at the bifurcations in the lower part of the airway and again it increases with the Stokes number which is characteristic of the inertial process.

Conclusions

After a validation phase, we have simulated the flow and the particle deposition in many Weibel's generic models and in anatomically realistic geometries of rat and rabbit coming from medical imaging. The particle deposition is mainly inertial and increases with the Stokes number and the Reynolds number. One of the major difficulty of this study is the generation of the appropriate mesh for the realistic geometries.

Fig. 2 :

 2 Fig. 2 : deposition in rabbit lungs


G Martinat
email: martinat@ccpo.odu.edu

M Braza

Y Hoarau

G Harran

Turbulence modelling of the flow past a pitching NACA0012 airfoil at 10 5 and 10 6 Reynolds numbers

This paper provides a study of the NACA0012 dynamic stall at Reynolds numbers 10 5 and 10 6 by means of two and three dimensional numerical simulation. The turbulence effect on the dynamic stall is studied by statistical turbulence modelling. The results are compared with experiments concerning each test case. Standard URANS turbulence modelling have shown a quite dissipative character that attenuates the instabilities and the vortex structures related to the dynamic stall. The URANS approach OES, Organised Eddy Simulation, has shown an improved behaviour at the high Reynolds number range. Emphasis is given to the physical analysis of the three dimensional dynamic stall structure, for which there exist few numerical results in the literature, as far as the Reynolds number range is concerned. This study has shown that the downstroke phases of the pitching motion are subjected to strong three-dimensional turbulence effects along the span, whereas the flow is practically two-dimensional during the upstroke motion.

Introduction

The prediction of dynamic stall phenomenon at high Reynolds number is a crucial need in aeronautics and more specifically in rotorcraft dynamics. In this context, the forced unsteadiness (organised motion) interacts non-linearly with the fine-scale random turbulence and produces a strong irreversibility effect that usually leads to hysteresis loops in the aerodynamics coefficients versus the angle of incidence. This hysteresis effect is the result of a phase lag between the bodys motion and the near-region velocity-vorticity gradients. Under these above effects, the stall angle is found higher than the normal static stall angle. The dynamic stall creates conditions of strong non-equilibrium turbulence, where the turbulence kinetic energy production is not equal to the dissipation. In this case, standard turbulence modelling approaches are often insufficient to capture correctly the hysteresis effect and to predict the dynamic stall loads at high Reynolds number. The applications of these flows occur in turbomachinery and in helicopter rotorblades as well as in wind turbine airfoils. It is important to have a good prediction of the dynamic stall to ensure efficiency for design. A comprehensive review of the dynamic stall is described in [START_REF] Mccroskey | The phenomenon of dynamic stall[END_REF][START_REF] Mccroskey | Unsteady airfoils[END_REF], [START_REF] Chandrasekhara | Flow visualisations studies of the Mach number effets on dynamic stall[END_REF] and [START_REF] Piziali | An experimental investigation of 2D end 3D oscillating wing aerodynamics for a range of angle of attack including stall[END_REF]. Experimental work from Laboratoire d'Aérodynamique and Biomécanique du Mouvement (LABM) has given a detailed physical analysis of this phenomenon [START_REF] Berton | Embedded l.v. methodology for boundary-layer measurements on oscillating models[END_REF]. More recently, [START_REF] Tsang | Dynamic stall behavior from unsteady force measurements[END_REF] provides an experimental study of dynamic stall using force measurement and [START_REF] Akbari | Simulation of dynamic stall for a NACA 0012 airfoil using a vortex method[END_REF] provides a numerical study using vortex method to model the flow past a pitching airfoil.

In the context of a pitching airfoil, with a sufficiently high incidence angle that depends on flow parameters such as pitching amplitude, pitching frequency, average incidence and profile geometry, a high energy vortex develops on the leading edge of the profile and is advected near wall towards the trailing edge. Consequently, this vortex grows and leads to a boundary layer thickness increase and finally to separation. When the leading edge vortex is shed, the flow is fully stalled inducing a abrupt loss in the lift coefficient, an increase of drag coefficient and a negative peak on the moment coefficient. With the decrease of incidence, the flow finally reattaches from the leading edge. From the point of view of numerical modelling, the work of the european program UNSI (Unsteady Viscous Methods in the Context of Fluid-Structure Interaction, [START_REF] Haase | Progress in computational flow-structure interaction, results of the project UNSI, supported by the european union 1998 -2000[END_REF]) reports a comparison of CFD approaches for flows around pitching airfoils at high Reynolds number based on [START_REF] Piziali | An experimental investigation of 2D end 3D oscillating wing aerodynamics for a range of angle of attack including stall[END_REF] and [START_REF] Berton | Embedded l.v. methodology for boundary-layer measurements on oscillating models[END_REF] experiments. It is reported that stantard URANS methods are unable to properly capture experimental hysteresis cycle. A slight improvement can be achieved by using non-linear eddy viscosity modelling (NLEVM), explicit algebraic Reynolds stress model [EARSM, [START_REF] Wallin | An explicit algebraic Reynolds model for incompressible and compressible turbulent flows[END_REF]] or by calibrating a one-equation model based on the eddy viscosity transport equation of [START_REF] Spalart | A one equation turbulence model for aerodynamics flows[END_REF] model, specifically for a pitching airfoil [SALSA model, Rung et al. (2002)]. Finally, [START_REF] Benyahia | Detailed evaluation of cfd predictions against lda measurement for flow n an aerofoil[END_REF] have shown an improved behaviour by using k -ω SST model.

The aim of the present study is to evaluate the turbulence modelling performance by comparing classical and advanced URANS approaches in two-and three-dimensional simulation. Furthermore, dynamic stall is a peculiar case, less dependant on three-dimensional effects than "static" bluff body turbulence during the upstroke phase and more dependant on three-dimensionality during the downstroke phase [START_REF] Berton | Experimental methods for subsonic flow measurements[END_REF]. Therefore, an emphasis is given is given on three-dimensional simulation of the dynamic stall at high Reynolds number. Hybrid (URANS-LES) turbulence modelling is employed in this case by using a DDES, Delayed Detached Eddy Simulation [DDES, [START_REF] Spalart | Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[END_REF][START_REF] Spalart | A new version of detached eddy simulation, resistant to ambiguous grid densities[END_REF]]. The results are compared ith two-dimensional statistical turbulence modelling.

Turbulence modelling : macrosimulation for unsteady flows

Organised Eddy Simulation (OES)

Advanced URANS modelling aims at capturing non-equilibrium turbulence effects especially occuring in the dynamic stall context. Among different macrosimulation approaches, the Organised Eddy Simulation (OES) [START_REF] Braza | Turbulence properties in the cylinder wake at high Reynolds number[END_REF]) was developed by the research group EMT2 (Ecoulements Monophasiques Transitionnels et Turbulents) at CNRS/IMFT (Centre National de la Recherche Scientifique, Institut de Mécanique des Fluides de Toulouse). This approach consists of distinguishing the structures to be resolved from the structures to be modelled on the basis of their physical nature, organised or chaotic and not on their size as in LES approach. According to OES, the turbulence spectrum is extended in whole the range from low to high wavenumbers and therefore it can be modelled by statistical turbulence modelling consideration. However, due to strong non-equilibrium turbulence effects in the present class of inhomogeneous turbulence flows, shape and slope modification occur in the inertial range of the turbulence energy spectrum, as shown in detailed experimental studies by [START_REF] Braza | Turbulence properties in the cylinder wake at high Reynolds number[END_REF] using time-resolved PIV. As a consequence, production is not equal to dissipation like in URANS equilibrium turbulence modelling. Therefore, the statistical turbulence scales used to model the continuous part of the energy spectrum have to be properly reconsidered to capture the non-equilibrium effects. In OES, the organised part of the motion corresponds to distinguished peaks in the spectrum that correspond to the resolved turbulence. In the time-domain, due to the pronounced periodic character of the pitching motion (forced unsteadiness), the phase-averaged Navier-Stokes equations are considered according to the decomposition of [START_REF] Cantwell | An experimental study of entrainment and transport in the turbulent near wake of a circular cylinder[END_REF]. The periodic nature of the flow past an oscillating airfoil allows definition of phase-averaged quantities. The flow is decomposed into a periodic component and a random fluctuation : U i = Ūi + Ũi . A detailed presentation of the OES approach is provided in a companion paper by Bourguet et al., in the present dedicated volume of the the IUTAM Symposium "Unsteady Separated Flows and their Control". The reader is kindly invited to refer to this companion paper for more details concerning this approach. OES has been applied succesfully in a number of European aeronautics research programs, among which, ETMA [Efficient Turbulence Models in Aeronautics, [START_REF] Dervieux | Computation and comparison of efficient turbulence models for aeronautics[END_REF]], FLOMANIA [Flow Physics Modelling an Integrated Approach, [START_REF] Bouhadji | Turbulence modelling of unsteady flows with a pronounced periodic character[END_REF]] and DESIDER [Detached Eddy Simulation for Industrial Aerodynamics].

In the present study, a first-order OES modelling is employed [START_REF] Bourdet | Prediction and physical analysis of unsteady flows around a pitching airfoil with the dynamic mesh approach[END_REF], derived from a second order moment closure DRSM [START_REF] Launder | Progress in the development of a Reynolds stress turbulence closure[END_REF]. In this context, a modified two-equation model was derived, where the turbulence length scales were modified by means of evaluation of the Cµ eddy diffusion coefficient through DRSM and by reconsidering the turbulence damping law towards the wall [START_REF] Hoarau | Prediction of turbulent unsteady aerodynamic flows with a pronounced periodic character[END_REF][START_REF] Jin | A two-equation turbulence model for unsteady separated flows around airfoils[END_REF]. In the context of high-order stress-strain modelling, a tensorial OES eddy-viscosity model has been derived to capture the non-equilibrium turbulence [START_REF] Bourguet | Anisotropic eddy-viscosity concept for strongly detached unsteady flows[END_REF] where the Cµ eddy diffusion coefficient varies according to a directional criterion of stress-strain misalignment. The first order OES modelling yields a two-equation low-Reynolds that derives from the original low-Reynolds two-equation model from [START_REF] Chien | Predictions of channel and boundary-layer flows with a low-Reynolds number turbulence model[END_REF]. the present model is proven less dissipative than the standard two-equation first-order models and is described by the following system of equations :

Dk Dt = ∂ ∂y (ν + ν t ) ∂k ∂x j + ν t ∂Ui ∂x j 2 -- 2νk y 2 , ( 1
)
D Dt = ∂ ∂xj [(ν + ν t σ ) ∂ ∂xj ] + C 1 k ν t ( ∂U i ∂xj )2 -k [C 2 f + 2νk e -C 4 U * y ν y 2 ], (2)
where,

f = 1 - 0.4 1.8 exp(-(k 2 /6ν ) 2 ), ( 3
)
ν t = C µ f µ k 2 (4)
with Cµ = 0.02, ( 5)

fµ = 1 -exp(0.0002y + 0.000065y +2 ). (6)
The advantages of this approach are the robustness at high Reynolds number wall bounded flows and the fact that the method is not intrinsically three-dimensional.

Detached Eddy Simulation (DES)

A reported by [START_REF] Travin | Detached-Eddy Simulations past a Circular Cylinder[END_REF], the Detached-Eddy Simulation is a three-dimensional simulation using statistical turbulence modelling equations, where the turbulence length scale is chosen between its RANS expression and a sub-grid length scale as in Large Eddy Simulation. This choice operates as a RANS model in the near-wall region and as a LES model in the regions of flow detachment. The advantage of the present hybrid (RANS-LES) approach comparing to other hybrid approaches is that it uses a single set of statistical turbulence modelling equations and does not require the definition of an interface between RANS and LES regions. The DES approach is intrinsically 3-D.

The DES length scale is chosen according to the following equation :

d = min(d, CDES∆), ( 7
)
where CDES is the DES constant calibrated by means of homogeneous, isotropic turbulence spectrum, ∆ is the largest dimension of the elementary control volume cell, ∆ = max(∆x, ∆y, ∆z) and d is the URANS lengthscale. This is the distance to the nearest wall in the case of Spalart and Allmaras modelling and a distance homogeneous to k 3/2 / in two equation modelling. For example, for a k -ω turbulence model, d = √ k/βω. This replaces the length scale expression in the dissipative source term of the transport equation of the turbulence kinetic energy k.

It has been demonstrated that the DES approach leads to URANS in the near wall region and to a subgrid LES modelling in the farther region, where the vortices detachment occurs. DES is less dissipative than URANS thanks to the choice of the lengthscale that allows increasing the dissipation term in the turbulence kinetic energy equation comparing to URANS.

Delayed detached eddy simulation (DDES)

To avoid a transition from URANS to LES in the boundary layer that could produce non physical artefacts, [START_REF] Spalart | A new version of detached eddy simulation, resistant to ambiguous grid densities[END_REF] introduced a modification of the r d parameter of the Spalart-Allmaras model, involved in near-wall damping, as follows :

r d = ν + ν t Sijκ 2 d 2 . ( 8
)
The definition of terms and constants can be found in appendix A. By adding the viscosity ν on the numerator as a small parameter (in high Reynolds number flows), it is ensured that r d remains away from 0 in the near wall regions. Therefore, a damping function has been suggested :

f d = 1 -tanh[(8r 3 d )]
which is 1 away from the wall an 0 in the near wall regions where rd << 1. Finally, the turbulence length scale becomes :

d = d -f d max(0, d -CDES∆). ( 9
)
If f d = 0, d = d which yields to RANS modelling and if f d = 1, d = min(d, C DES ∆) which yields to the classical DES modelling [START_REF] Spalart | Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[END_REF]. The boundary layer is then shielded from a transition from URANS to LES and this transition is moved farther away from the wall in case of an excessive mesh refinement.

3. Physical analysis of the dynamic stall around a pitching airfoil at 10 5 and 10 6 Reynolds numbers

The term dynamic stall usually refers to unsteady separation and stall phenomena on airfoils that are forced to execute time dependant motion. If the angle of attack oscillates around a mean value of the order of the static stall angle, a large hysteresis cycle develop concerning the aerodynamic forces and moment coefficients. Indeed, during the upstroke motion, the effect of adversure pressure gradient is limited, leading to a dynamic stall angle which largely exceeds the static one. The aerodynamics coefficients are higher than their static counterparts. During the downstroke part of the motion, the effect of adverse pressure gradient is reinforced, leading to a reattachment incidence angle which is lower than in the static case.

Description of the test cases and numerical methodology

In both test cases of this study, the airfoils perform a sinusoidal pitching motion which is described in detail in table 1. Turbulence models used for Berton et two-dimensional computations as well as DDES k -ω SST [START_REF] Menter | Zonal two equation k -ω turbulence models for aerodynamics flows[END_REF] for three-dimensionnal computations. For the McAlister et al. test case, the turbulence models used are k -ω SST, k -ε Chien and k -ε OES for twodimensional computations. Three-dimensional computations have not been carried out in this test case because of the high computational time needed in respect of a lower reduced pitching frequency. The mesh used for the Berton et al test case is a C topology mesh with 256 nodes in the I direction, 80 nodes in the J direction and 40 nodes in the K direction for three dimensional computations. The mesh refinement near-wall allows a y + of an order of 1 in the attached flow region. For the Mc Alister et al test case, an O topology mesh is used with 285 nodes in the I direction and 185 nodes in the J direction. The mesh refinement allows a y + less than one near the wall.

The flow solver is the structured NSMB (Navier Stokes Multi Blocks) code used by the NSMB consortium in which the research group of the authors participates. A detailed description of NSMB, including high-order numerical schemes (implicit and explicit) as well as various turbulence closures (standard and advanced URANS), LES and DES can be found in [START_REF] Vos | Recent advances in aerodynamics inside the nsmb (navier stokes multi blocks) consortium[END_REF]. In the present study, for both test case, the temporal scheme is an implicit backwards second order with dual time stepping. The space discretisation scheme is th upwind-Roe third order with MUSCL fluxes evaluation (van Leer). An adaptive time-step has been used, allowing a CFL number inferior to 1 in the inner loop and inferior to 3 on the outer time-step. The mesh movement in respect of the pitching airfoil is achieved by means of the Arbitrary Lagrangian and Eulerian approach (ALE). The NSMB solver has been used in its full parallelized version for MPI architectures. The super computing centers CINES (Centre Inter-universitaire National d'Enseignement Supérieur), IDRIS (Institut de Recherche en Informatique et Calcul Scientifique), GRID-5000 and CALMIP (Calcul Midi-Pyrennées) have been used for the present computations, providing the CPU allocation.

Simulation of the Berton et al test case

Two dimensional simulations

The flow topology is qualitatively in agreement with the experiment. Fig. 1 shows the isovorticity fields for 10 different incidence values. From 7.2 o to 12 o upstroke, the flow remains attached to the profile. From 14.4 o to 16.8 o the birth of the leading edge vortex can be observed, as well as its progressive growth until it is shed at the maximum incidence. This structure is immediately followed by a trailing edge vortex that also occurs at the maximum incidence. From 18 o to 9.6 o downstroke, the flow is fully stalled and reattachment begins for 7.2 o of incidence downstroke. The above mentioned vortex formation and shedding corresponds to the undulations appearing in the hysteresis loops.

The turbulence models tested display a different behaviour. As seen on Fig. 2, the Spalart and Allmars model over-estimates the maximum lift and drag coefficient, as well as the area under the hysteresis loop. This behaviour is probably due to an excessive dissipative character, which is confirmed by the underestimation of the adverse pressure gradient at high incidence leading to a delay on the prediction of dynamic stall.

The k -ε Chien model also over-estimates the maximum lift and drag coefficient and the area under the loop but in a less critical way than the Spalart-Allmaras model. This can be explained by a less dissipative character than the previous model. However, the dynamic stall is still enhanced by the beginning of the downstroke phase, showing again an underestimation of the adverse pressure gradient.

The two URANS models tested have shown an over-estimation of the lift coefficient and large lift oscillation during the dowstroke phase which are due to an inhibition of the flow in the third dimension because of the two dimensional computation. Therefore, we can observe that the oscillation amplitude diminishes with the decrease of incidence and the Strouhal number of the vortex shedding is found at 0.48 a value that is close to the experimental one for a static stalled airfoil [START_REF] Hoarau | Organized modes and the three-dimensional transition to turbulence in the incompressible flow around a naca0012 wing[END_REF]. This means that the vortices observed are shed on a natural frequency of the flow.

As a conclusion in this part, the k -ε Chien model provides the best results comparing to Spalart and Allmaras model for the case of a Reynolds 10 5 pitching airfoil. Another source of difference with the experiment is the transitional charter of the boundary layer in the experiment. As the three model examined use the hypothesis of a fully turbulent boundary layer, this partly explain the discrepancy between the experiment and the computations.

Three dimensional simulations

The three dimensional computations have shown a practically two-dimensional behaviour during the upstroke part of the motion but a fully three dimensional behaviour during the downstroke phase (Fig. 3). The dynamical increase of incidence is found to enhance two dimensional character of the flow. On the contrary, the decreasing incidence during the downstroke phase of the wing motion reinforces the adverse pressure gradient, producing a fully stalled and highly three-dimensional flow around the airfoil.

Comparing to two-dimensional computations, the hysteresis loops are improved and are able to produce the decrease of the large lift oscillation during the downstroke phase. Therefore, the three dimensional computations show a significant improvement of the results. Considering that for the maximum and minimum incidence of the pitching motion the flow is two dimensional, the extreme lift and drag coefficient are not affected by the threedimensional effects. For this reason, lift and drag coefficient are the same during the upstroke phase in 2-D and 3-D.

Simulation of the McAlister et al test case

To avoid the transitional boundary layer, computations at a Reynolds number of 10 6 are performed according to the experimennt of [START_REF] Mcalister | Dynamic stall experiments on the NACA0012 airfoil[END_REF]. As seen on Fig. 5 the turbulence models, k -ω, k -Chien and k -OES indicate similar behaviour.

The k -ε Chien provides results close to the experiment with a good prediction of lift, a small overprediction of drag but a large delay on the predicted stall angle. This means that the model has an excessive dissipative character that prevents from a correct prediction of the adverse pressure gradient leading to stall. As a result, its behaviour is forced by the movement of the structure.

The k -ω SST gives the most accurate results with a good estimation of the stall angle and good prediction of the lift and drag coefficient during the upstroke motion. However, this model gives the less smooth character which could be due to the addition of the SST limiter that creates abrupt reduction of the eddy viscosity. The k -OES modelling predicts the best angle of dynamic stall but overestimate the lift and drag coefficients.

The three models are over-estimates the lift coefficient during the downstroke motion. The undulations appearing in the lift and drag hysteresis are due to the coherent vortex shedding. This vortex shedding is more pronounced in the 2-D simulations, because of flow confinement induced by the two-dimensional approximation.

Conclusion and prospects

In this study, two-and three-dimensional simulations have been performed comparing the behaviour of various URANS and hybrid turbulence modelling approaches on the flow around a pitching airfoil.

Advanced turbulence models like OES and SST modelling have shown better results than classical URANS models for the occurence of dynamic stall prediction.

The 10 5 Reynolds number computations have shown the need of a transition model. Finally, three dimensional computations have shown an improvement on the downstroke phase allowing the decrease of the coherent structure amplitude comparing to 2-D computations.

However, the two-dimensional study is useful for fast pre-design use, because they are able to capture a significant part of the coherent structures dynamics in a good approximation comparing to the three-dimensional simulations.

A Description of the Spalart and Allmaras one equation turbulence model

The Spalart and Allmaras model [START_REF] Spalart | A one equation turbulence model for aerodynamics flows[END_REF]) is a one equation turbulence model. The transport equation for eddy viscosity is made under dimensional analysis assumptions. The eddy viscositytransport equation is given by

∂ ν ∂t + ũj ν) ∂x j = P + D IF F + D ES , (A.1)
where P , DIF F and DES are respectively production, dissipation and destruction terms of ν expressed as

P = c b1 S ν, (A.2) DIF F = 1 ρ ∂ ∂xj (ν + ν) ∂ ν ∂xj + c b2 ∂ ν ∂xj ∂ ν ∂xj , (A.3) D ES = c ω1 f ω ν d 2 . (A.4)
The eddy viscosity is then defined as follows :

ν t = νf v1 , (A.5)
To ensure that ν is equal to κyu τ in the logarithmic part of the boundary layer, the damping function f v1 is defined as .6) where χ = ν/ν. However, S is modified in S in order to maintain a correct behaviour in the logarithmic part of the boundary layer S = u τ /κy : where g acts as a limiter, avoiding fω to become to high and is expressed as : g = r + cω2 r 6 -r with r = ν Sκ 2 d 2 . r and f ω are both equal to 1 in the logarithmic part of the boundary layer and decrease outside. Constants of the model are :c b1 = 0.1355, c b2 = 0.622, σ = 2 3 , κ = 0.41, c ω1 = c b1 κ 2 + 1+c b2 σ , c ω2 = 0.3, c ω3 = 2 et cv1 = 7.1. The Spalart and Allmaras model is considered as a good compromise between algebraic models and two equation models. This model has become very popular due to its robustness and reasonable predictive efficiency in a wide number of aerodynamic flows.

f v1 = χ 3 χ 3 + c 3 v1 , ( A
S = 2Ω ij Ω ij f v3 + ν κ 2 d 2 f v2 , (A.

  i , f v2 = 1 -χ χ+f v1 and f v3 = 1. Finally, in order to obtain a decrease of the destruction part outside the boundary layer, the function fω i introduced such that :

α

  Fig. 1. Isovorticity field for a pitching cycle using k -εChien modelling

Fig. 4 .Fig. 5 .

 45 Fig. 4. Hysteresis loops obtained with 3-D computations on lift and drag coefficient and compared to experimental results and 2-D computations with Spalart and Allmaras model

  al test case are Spalart and Allmaras and k -ε Chien for Physical parameters of Berton et al test case and Mc Alister et al test case.

	Experiment	Berton et al. (2002) McAlister et al. (1978)
	Reynolds number (based on chord length)	1.10 5	1.10 6
	Mean incidence	12 o	15 o
	Pitch amplitude	6 o	10 o
	Reduced oscillation frequency	0.188	0.1
	Table 1


Ilhem Slama
email: ilhem.slama@imt-atlantique.fr

Oussama Ben-Ammar
email: oussama.ben-ammar@emse.fr

Alexandre Dolgui
email: alexandre.dolgui@imt-atlantique.fr

Faouzi Masmoudi
email: masmoudi.fawzi@gmail.com

Newsboy problem with two-level disassembly system and stochastic lead time

Keywords: Reverse logistics, stochastic disassembly lead-time, Newsboy model

Introduction

Nowadays, there is an increasing interest in disassembly process research. The main objective is to obtain valuable components and separate those that may affect the environment. Actually, the disassembly process has become very vulnerable and far to be deterministic [1]. However, it involves a high degree of variability associated with disassembly lead time. This unpredictability may be caused by technical problems (absenteeism, limited disassembly capacity, etc.) and to economic condition (changes and increase in costs of end of life products, etc). Here, the disassembly lead time implies the total time required to receive the parts after placing an order of disassembling the parent item [2].

Problem description and formulation

In this work, we study the case of single-demand planning in a disassemble-to-order (DTO) environment. In order to fulfill demands for specific components i = 1, ..., n, we disassemble an end-of-life (EoL) product whose real disassembly lead time, denoted L, is a random discrete variable varying between L -and L + and characterized by a cumulative distribution function F (.). We assume that (i) the disassemble capacity is infinite, (ii) for each component-customer demand, the quantity is known, fixed to one and should be delivered at period T i , and (iii) the unit backlogging (b i ) and holding (h i ) costs for components are known. The disassembly process starts at period X and all components are available at period M = X + L. If a given component i is available before T i , we stored it until this period. Then the corresponding holding cost is equal to

h i (T i -T - i ) with T - i = min(M, T i ) . If it is available after T i , a backlogging cost is incurred and equals to b i (T + i -T i ) with T + i = max(M, T i ).
The total cost T C(X, L) is a random discrete variable depending on the variability of L. Its explicit form can be written as follows :

T C(X, L) = n i=1 (b i + h i )T + i -M h i -b i T i ∀X ∈ [[ min ∀i∈[[1,n]] T i -L + ; max ∀i∈[[1,n]] T i -L -]] (1)
The objective is to find the best disassembly date, denoted X * , that minimizes the expected total cost E(T C). In this preliminary work, we give only the explicit forms. Thus, the expression of E(T C) is in Equation 2:

E(T C) = n i=1 (b i + h i ) s≥Ti 1 -F (s -X) + h i T i -X + E(L) n i=1 h i (2)

Optimization approach

The expected total cost expressed in Equation 2 constitutes a non-linear objective function to be minimized. Here, an exact method based on the well-known Newsboy formula is introduced to optimize the studied problem. In order to obtain the optimal disassembly date X * , the computation of X is required to explore the of all possible values of X * and is stopped when we find a solution that checks Equation 3.

Proposition 1 The Newsboy formula gives the optimal disassembly date X * which satisfies the following optimality condition (Equation 3) :

n i=1 (b i + h i ) × F (T i -X * -1) ≤ n i=1 b i ≤ n i=1 (b i + h i ) × F (T i -X * ) (3)
The full proof will be published soon elsewhere. Note that for n = 1, i.e. for a single type of component, our model is equivalent to the one of [3]. Theoretically, the optimum X * is not unique. In order to prove the uniqueness of the optimal solution, it can easily be proved that the objective function is convex.

Conclusions

In this work, we propose a mathematical model to study a single-demand planning and a two-level disassembly system in a DTO environment with stochastic lead time. In order to optimize the expected total cost which is composed of inventory and backlogging costs, we introduce an optimizing approach. In our knowledge, we are the first to propose a Newsboy approach to solve the disassembly scheduling problem under stochastic disassembly lead time.

Several perspectives to this work are identified. Firstly, to develop mathematical models for the case of multi-period planning with dependency between periods. Secondly, it will be interesting to investigate the case where each component is independently disassembled under uncertainty of lead time. Finally, the objective is to generalize the proposed approach for the case of multi-level disassembly systems.

Références

  Lee D-H and Xirouchakis P. Two-phase heuristic for disassembly scheduling with multiple product types and parts commonality International Journal of Production Research, 44 :(1) :195-212, 2006. [3] Ould-Louly Mohamed-Aly and Dolgui Alexandre. A generalized Newsboy model for supply planning of assembly systems with service level constraint Proceedings of the International Conference on Industry, Engineering and Management Systems, 5-7, 2001.

	[1] Kim Hyung-Won, Park Chuljin and Lee Dong-Ho. Selective disassembly sequencing with
	random operation times in parallel disassembly environment International Journal of Pro-
	duction Research, 56 :(24) :7243-7257, 2018.
	[2] Kim H.-J,