Gaëlle Loosli

Stéphane Canu
email: scanu@insa-rouen.fr

Non positive SVM

Keywords:

 have effectively demonstrated SVM classification with a non-positive kernel. However the question of efficiency still applies. In this paper, an efficient direct solver for SVM with non-positive kernel is proposed. The chosen approach is related to existing work on learning with kernel in Krein space. In this framework, it is shown that solving a learning problem is actually a problem of stabilization of the cost function instead of a minimization. We propose to restate SVM with non-positive kernels as a stabilization by using a new formulation of the KKT conditions. This new formulation provides a practical active set algorithm to solve the indefinite SVM problem. We also demonstrate empirically that the proposed algorithm outperforms other existing solvers.

La possibilité d'utiliser les SVM avec des noyaux non positifs a été démontrée récemment avec les Indefinite SVM [START_REF] Luss | Support Vector Machine Classification with Indefinite Kernels[END_REF][START_REF] Ying | Analysis of svm with indefinite kernels[END_REF][START_REF] Chen | Training svm with indefinite kernels[END_REF]. Toutefois la question de l'efficacité de la résolution se pose toujours. Dans cet article, nous proposons un algorithme efficace qui résout directement les SVM avec un noyaux nonpositif. L'approche choisie est liée aux travaux existant sur l'apprentissage avec noyaux dans les espaces de Krein. Dans ce cadre, il est démontré que résoudre un problème d'apprentissage revient en fait à trouver un point stable de la fonction coût plutôt que le minimum. Nous proposons de reformuler les SVM avec noyaux non positifs comme une stabilisation en adaptant les conditions KKT. Cette nouvelle formulation permets de mettre en oeuvre un algorithme de contraintes actives pour résoudre le problème des SVM indéfinis. Nous apportons également une preuve empirique de la supériorité de cet algorithme par rapport aux approches existantes.

Mots clés : Noyaux non positif, résolution de SVM 1 Learning with non positive kernel

From the first stages of SVM [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF], non positive kernels are proposed and used, in particular the tanh kernel. In many application fields, some huge efforts are made to produce true Mercer kernels when the natural kernels turn out to be indefinite (for instance [START_REF] Cortes | Positive definite rational kernels[END_REF][START_REF] Collins | Convolution kernels for natural language[END_REF]). Some author even study some kernels that are definite positive with high probability ( [START_REF] Boughorbel | Nonmercer kernels for svm object recognition[END_REF]). However, until now, there is no adequate solver available. In [START_REF] Luss | Support Vector Machine Classification with Indefinite Kernels[END_REF], the authors propose to solve SVM with indefinite kernel considering that the indefinite kernel is a perturbation of a true Mercer kernel. They end up with a new optimization problem that is able to provide good solutions. The major drawback however would be that the proposed algorithm requires to tune at least one more parameter, plus some different stopping criteria. Independently of the interest of this approach, it makes practical applications difficult. In [START_REF] Ying | Analysis of svm with indefinite kernels[END_REF][START_REF] Chen | Training svm with indefinite kernels[END_REF], the previous method is analyzed and improved. In [START_REF] Haasdonk | Feature space interpretation of svms with indefinite kernels[END_REF], the author states that learning with indefinite symmetric kernels is actually consisting in finding a stationary point, which is not unique but each of those performs correct separation. Moreover, it is shown that the problem is then cannot be seen as a margin maximization although a notion of margin can be defined. We will see that similar notions appear in our formulation.

It has been shown ( [START_REF] Cheng | Learning with non-positive kernels[END_REF]) that learning with non positive kernel is actually solving the learning problem in a Krein space instead of a Hilbert space. It has also been shown that in this situation, the learning problem is not a minimization anymore but a stabilization problem. This means that the solution is a saddle point of the cost function. Applying this to SVM requires to interpret this stabilization setting. We first give the intuition behind the proposed method. Following [START_REF] Hassibi | Indefinitequadratic estimation and control: a unified approach to H2 and H [infinity[END_REF], we start from the fact that a (unconstraint) quadratic program in a Krein space has a unique solution (if the involved matrix is non singular) which is in general a stationary point. In the case of SVM, we have to apply some box constraints that may exclude this unique solution. Moreover, the optimal constraint solution is not necessarily unique anymore.

In the definite positive case, solving the quadratic program under box constraints is a problem of projection of the unconstrained solution onto the admissible set. To do so, the gradient of the cost function gives the direction towards the minimum admissible point. In the indefinite case, this projection has to be redefined : indeed, since the optimal solution is not the minimum, we can not follow the gradient direction, which leads to the minimum. In our work, we propose another projection, in order to obtain the most stable point within the admissible set. In other words, we want the admissible point that has the lowest gradient.

In the next section, it is shown that the usual SVM problem with positive definite kernels can be directly extended to the indefinite case such that the provided solution is (one of) the most stable admissible point(s). In order to illustrate the fact that the proposed formulation is a stabilization, the original problem is modified such that any critical point becomes a local minimum. We point out that the optimal solution of the modified problem is the same as SVM with non positive kernels. Finally, the projection point of view is provided and it is shown that it also leads to the same problem. Having those three approaches, we conclude that the solution of a non positive SVM is the solution of a stabilization problem, which is also the projection of the unique unconstraint solution onto the admissible set, in the sense of the most stable point.

SVM and KKT conditions of optimality

Let's consider a training set {X , Y} where x i ∈ X ∀i ∈ [1..n] are examples and y i ∈ Y are corresponding labels. In the case of binary classification, possible values for labels are [-1; 1]. We denote K the symmetric kernel matrix, of size n × n. C is the hyper parameter, penalizing the training errors.

We recall in this section the SVM dual quadratic problem ( [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]). This dual is obtained from the margin maximization formulation. A Lagrange multiplier α i is associated to each training example.

   min α 1 2 α ⊤ Gα -α ⊤ 1 subject to α ⊤ y = 0 and 0 ≤ α i ≤ C ∀i ∈ [1..n] (1)
where G(i, j) = y i y j K(i, j).

We give the full KKT conditions for the classical SVM dual system (1). The KKT conditions of optimality are divided into stationarity condition, primal and dual admissibility and complementarity conditions. In the case of SVM dual, the stationarity condition is as follows:

α ⊤ G -1 ⊤ + λy ⊤ -µ ⊤ + η ⊤ = 0 (2)
The primal admissibility is given by

α ⊤ y = 0 α i ≤ C ∀i ∈ [1..n] α i ≥ 0 ∀i ∈ [1..n] (3)
The dual admissibility is given by

µ i ≥ 0 ∀i ∈ [1..n] η i ≥ 0 ∀i ∈ [1..n] (4)
The complementary conditions are

α i µ i = 0 ∀i ∈ [1..n] (α i -C)η i = 0 ∀i ∈ [1..n] (5)
Any solution respecting each of these conditions is a solution of the problem. In the case the kernel matrix is definite positive, the dual SVM problem has a unique solution which is a global minimum.

Point of view 1 : the variational approach of quadratic programming

In the case the kernel matrix is indefinite, the dual SVM problem is not well defined and the solution is not unique. Following [START_REF] Brezinski | Projection methods for linear systems[END_REF] and the variational approach of quadratic programming, we can actually solve the problem using normal residuals (ie. solving Ax = b via A ⊤ Ax = A ⊤ b). Note that normal equations could also be an option (ie. solving

Ax = b via AA ⊤ x ′ = b with x = A ⊤ x ′ ). α ⊤ G = 1 -λy ⊤ + µ ⊤ -η ⊤ α ⊤ GG ⊤ = (1 -λy ⊤ + µ ⊤ -η ⊤ )G ⊤ (6)
This can be seen as least squares. All the other conditions remain identical.

Point of view 2 : a stabilization problem

Optimization technics provide minima and not saddle points. To that purpose, we propose to use a simple trick known as magnitude of the gradient, consisting in changing the problem such that any critical point becomes a local minimum. This is done computing the sum of the squares of the partial derivatives of the function to be stabilized. Let apply this to the following unconstrained function:

J = 1 2 α ⊤ Gα -α ⊤ 1 (7)
Stabilizing J is equivalent to minimizing M:

M(α) = α ⊤ G -1 ⊤ , α ⊤ G -1 ⊤ (8)
This provides the following system:

   min α α ⊤ G -1 ⊤ , α ⊤ G -1 ⊤ subject to α ⊤ y = 0 and 0 ≤ α i ≤ C ∀i ∈ [1..n] (9)
Let now write the KKT condition of optimality for system [START_REF] Haasdonk | Feature space interpretation of svms with indefinite kernels[END_REF]. The stationarity condition is as follows:

(α ⊤ G -1 ⊤ + λy ⊤ -µ ⊤ + η ⊤ )G ⊤ = 0 (10)
The other condition are identical to eq. 3, 4 and 5. Those KKT conditions show that solving the SVM dual problem when the kernel is indefinite is actually a stabilization problem.

Point of view 3 : The projection

As already mentioned, the unconstraint problem has a unique solution. We want to project it onto the admissible set. However, this projection is not obvious : what is the closest point of the admissible set to the unconstraint optimum in the sense of the stabilization? We propose to define it as the most stable point, i.e. the admissible point minimizing the gradient of the cost function (which is α ⊤ G -1 ⊤ ). Solving this minimization with the least squares directly gives the same system as previously (eq. 9).

For illustration purpose, we draw the cost function for a sigmoid kernel of a simple 2 gaussians problem, obtained with two support vectors, for both CSVM cost and NPSVM cost. The plain area shows the admissible solutions. We can easily observe that the usual minimization technics would lead to a non optimal solution (see figure 1).

The solver

The proposed algorithm is derived from active set approach for SVM, similar to [START_REF] Vishwanathan | Simplesvm[END_REF]. The sets of points are defined according to the complementarity conditions (see table [START_REF] Newman | UCI machine learning repository[END_REF].

By default, all training points are in the non support vector set I 0 except for a couple with opposite labels which is in I w . Any other initial situation based on warm-start or a priori does not change the algorithm.

0 0 > 0 I C C > 0 0 I w 0 < α < C 0 0
3.1 Solving linear system in I w

We need to solve the linear system from the stationarity condition (eq. 10) only for unconstrained points, those of I w . This leads to the following equation :

α ⊤ (w) G (w,:) G (:,w) = (1 (:) -λy ⊤ (:) -C1 ⊤ (C) G (C,:) )G (:,w)
This can be solved using QR decomposition of G, for which one can maintain a rank one update at each step of the algorithm. Computing λ can be easily done substituting

α ⊤ (w) in α ⊤ (w) y (w) = -C1 ⊤ (C) y (C) .

Activating constraints in I w

If any α (w) (i) does not lie in [0 C], the current solution is projected on the admissible set such that all α (w) (i) satisfy the primal admissibility and the violating point is transferred towards I 0 or I C according to the violating value.

Relaxing constraints in I 0 or I C

If the current solution is admissible, we check the stationarity conditions for I 0 and I C (using eq. 10 and 5). The most violating point is transferred from its group to I w .

For any point j ∈ I 0 , µ j > 0 and η j = 0. From eq. 10 :

(α ⊤ [w,C] G ([w,C],:) -1 ⊤ + λy ⊤ )G :,j > 0
For any point k ∈ I C , µ k = 0 and η k > 0. From eq. 10 :

(α ⊤ [w,C] G ([w,C],:) -1 ⊤ + λy ⊤ )G :,k < 0
We can observe here that the notion of margin is distorted. Indeed, when using the same active set solver in the positive definite case, the margin clearly appears in the constraint relaxation (for j ∈ I 0 , the condition would be

α ⊤ [w,C] G ([w,C],j) + λy j > 1)
. This means that in the feature space, the solution won't have the same properties as the usual SVM, especially concerning the interpretation of support vectors relatively to the decision boundary.

The algorithm

1: Initialize (one random point for each class in I w , all others in I 0 ) 2: while solution is not optimal do 3:

solve linear system (sec 3.1)

4:

if primal admissibility is not satisfied then 5:

project solution in the admissible domain : remove a support vector from I w (to I 0 or I C ) (sec 3.2)

6:
else if stationarity condition is not satisfied then end if 9: end while Note that the convergence after a finite number of step of the proposed algorithm can always be proved since it can be seen as an active set procedure applied to a convex QP and thus convergence proof in this case applies [START_REF] Nocedal | Numerical Optimization[END_REF].

The code for proposed algorithm is available at ***.

Complexity

Compared to the same active set algorithm in the positive definite case, the proposed formulation increases the complexity. Let's consider that only the original kernel is stored in memory, we list the added operations -sec 3.

Experimental results

In our experiements, we used the well known sigmoid kernel (tanh) : k(x i , x j ) = tanh(scale× x i , x j +bias) and the epanechnikov kernel (epanech) : k(x i , x j ) = max(0, 1-γ * x i , x j ). We also used the positive gaussian kernel (rbf).

UCI datasets

We tested our algorithm NPSVM on some popular UCI datasets [START_REF] Newman | UCI machine learning repository[END_REF] and compared the results to a usual SVM solution (refered as CSVM). For each scaled dataset (heart, sonar, breast cancer, diabetes) we provide the average test accuracy on 10 random split of the database. The statistically better results are bold and underlined (according to chi square independance test at level 2%) (see table 2).

Validation protocol

We describe here the protocol used for each dataset, inspired by [START_REF] Van Gestel | Benchmarking least squares support vector machine classifiers[END_REF]. The procedure is applied 10 times for each and the given results are on average.

-split randomly the dataset, 2/3 for cross validation, 1/3 for test.

-perform 10 fold cross validation on the validation set (C ∈ [0.01, 0.1, 1, 10, 100, 1000], σ ∈ [0.1, 0.5, 1, 5, 10, 15, 25, 50, 100, 250, 500] * √ n] for rbf kernel, scale = [pow2(-5 : 1.5 : 2), -pow2(-5 : 1.5 : 2)] and bias = [pow2(-5 : 1.5 : 2), -pow2(-5 : 1.5 : 2)] for tanh kernel).

-train the svm on the full validation set with the parameters providing the best average performance during cross validation. -test on the separate test set.

Usps digits

We tested the NPSVM on the well known USPS digits database. With kernel, tanh with scale = 0.004 and bias = -1.5, C = 100 we obtained 94.92% of good classification (1 vs 1 multiclass setting). With an epanech kernel (bandwidth = 0.01 and C = 1000) we obtained 94.07% of accuracy. Those experiments convinced us that this approach of non positive SVM could be applied at least to medium sized problem. However the tested non positive definite kernels did not outperform the positive definite ones on that problem.

Monitoring the cost function

We also observed the evolution of the cost function (eq 1) to illustrate the fact that the cost function is stabilized and not minimized (figure 2). However, we can observe that for some particular setting and datasets, the proposed algorithms may converge to a minimum.

Comparison to over solvers libSVM

This part reports an experiment pointing out that the proposed approach can find much better solution than the widely used libSVM. Figure 3 illustrates this on a synthetic problem of checkers. The sigmoid kernel is used, with positive scale and negative bias, as recommended in [START_REF] Hsuan Tien | A study on sigmoid kernels for svm and the training of non-psd kernels by smo-type methods[END_REF]. We can observe that even though libSVM converges, the proposed solution poorly classifies the data while NPSVM converges towards a much better solution. Note that one can find some kernel parameters that works better for libSVM, but the reported result is frequent. Let us remark also that NPSVM is proposed here with an active set method but could also be implemented with an SMO approach.

IndefiniteSVM

We conducted some experiments to compare our results to the Indefi-niteSVM toolbox ( [START_REF] Luss | Support Vector Machine Classification with Indefinite Kernels[END_REF]). We show here one of the best results we manage to obtain (best in the sense of in favor of the indefiniteSVM). Note that we had difficulties to reach some good parameters for IndefiniteSVM. We can observe that even selecting favorable cases, IndefiniteSVM does not attain as good solutions as NPSVM. Moreover, those are usually less sparse.

Global comparison

From results in table 3, we can say that NPSVM clearly outperforms existing solvers handling non positive kernels. In terms of computational time, as expected the proposed implementation is slower than standard SVM solvers. However, it is much faster than IndefiniteSVM and can be used with medium size datasets.

Conclusion

In this article we present a new algorithm to solve SVM with non positive kernels. We base our work on the already stated fact that the solution, in this case, is a stabilized point and not a minimum, and we modified the solver according to that. We illustrate this point of view through three approaches, the variational point of view, the stabilization point of view and the projection point of view. All leads to the same problem for which we provide a practical algorithm. Compared to a usual SVM solver, the NPSVM needs to take into account all training point for intermediate optimizations, thus the computational time is increased. To that point, we think some heuristics can be applied : a subset of points, chosen at random a priori, could be used instead of the whole training set to support the intermediate optimizations. An other lead consists in solving the linear program induced by the projection point of view, if we take norm 1 instead of least squares. This will be studied in the near future. Let us note also recent work on projection in Krein space [START_REF] Ando | Projections in krein spaces[END_REF] that might be applied to the non positive kernels learning methods. Finally, we want to point out the fact that many well known technics, such as regularization paths, could easily be used for NPSVM.

Figure 1 :

 1 Figure 1: SVM cost function with sigmoid kernel, illustrated for 2 support vectors. The plain area shows the admissible solutions. On the left, the CSVM cost function (eq. 1). On the right, the NPSVM cost function (eq. 10). The blue cross on each graph shows the solution computed by solving the usual CSVM problem. The red circle shows the solution proposed by the NPSVM.

7 :

 7 add new support vector to I w (from I 0 or I C ) (sec 3.3) 8:

  1 : 2 matrix by matrix multiplication (O(n|I w | 2 ) and O(n 2 |I w |)) -sec 3.2 : identical -sec 3.3 : 1 matrix by matrix multiplication (O(n 2 (|I 0 | + |I C |)))These can be reduced by various strategies, such as caching, rank-one update, iterative search for violating example, etc.

Figure 2 :

 2 Figure 2: Evolution of the cost function during training. Curves are scaled on vertical axis. On the left, we report the value of the usual CSVM cost function during training. On the right are curves corresponding to the modified cost function that is actually minimized. While the standard CSVM cost function decreases monotonically (rbf CSVM on the left graph), the cost function for NPSVM may increase. Note that test performance are similar for each of the reported experiment on this figure (dataset is synthetic data generated from the apple/banana setting).

Figure 3 :

 3 Figure 3: Results on checkers with NPSVM on the right and libSVM on the left, for an identical sigmoid kernel (scale = 1, bias = -1). Circles are support vectors.

Figure 4 :

 4 Figure 4: Results on checkers with IndefiniteSVM on the left and NPSVM on the right, for an identical epanech kernel. Circles are support vectors.

Table 1 :

 1 Definition of groups for active set depending on the dual variable

	values	Group	α	η	µ
		I 0

Table 2 :

 2 Results on some UCI dataset.

	Solver kernel	Heart	Sonar	Breast
	C-SVM	rbf	82.22% (23.2 sv) 84.78% (90.3 sv)	97.47% (53.8 sv)
	NPSVM	rbf	83.44% (35.9 sv) 86.09% (94.9 sv) 97.37% (56.6 sv)
	NPSVM tanh 82.44% (14.8 sv) 84.06% (70.1 sv) 97.76 % (116 sv)

Table 3 :

 3 Comparison between IndefiniteSVM, libSVM and NPSVM on synthetic data, with the same kernel and C for all solvers . IndefiniteSVM is stopped after 2000 iterations if not converged

	Solver	Accuracy Training time	Solution size	Problem
	IndefiniteSVM	58.25 %	16.34 s	48 sv	Checkers 96 training points
	libSVM	48.56 %	0.003 s	52 sv	tanh kernel
	NPSVM	84.82 %	0.35 s	28 sv	min eig = -25.37
	IndefiniteSVM	68.21%	726 s	496 sv	Checkers 992 training points
	libSVM	51.27%	0.01 s	232 sv	tanh kernel
	NPSVM	90.84%	2.8 s	53 sv	min eig = -136.07
	IndefiniteSVM	-%	-s	-sv	Checkers 4992 training points
	libSVM	53.30%	1.54 s	2369 sv	tanh kernel
	NPSVM	92.41%	53.4 s	51 sv	-

Acknowledgement / Remerciements

This work was supported in part by the IST Programme of the European Community, under the PASCAL2 Network of Excellence, IST-2007-216886. This publication only reflects the authors' views. The authors also thank Ph. Mahey for helpful comments.


B Sahoo

S Poncet

Effets des conditions de glissement sur l'écoulement de Bödewadt d'un fluide de type Reiner-Rivlin

Keywords: écoulement de Bödewadt, fluide de type Reiner-Rivlin, conditions de glissement

L'écoulement stationnaire d'un fluide de Reiner-Rivlin en rotation au-dessus d'un disque fixe est étendu au cas où la surface du disque admet un glissement partiel. Les équations de quantité de mouvement se réduisent à un système d'équations fortement non-linéaires résolues par la méthode des différences finies. Les effets des conditions de glissement et des propriétés du fluide sur les caractéristiques de la couche limite sont discutés en détails. Le paramètre de glissement a un effet dominant sur le champ de vitesse, alors que le caractère non-newtonien influe sur le coefficient de moment.

Introduction

Les écoulements de fluides non-newtoniens engendrés par la rotation d'un disque ont été largement étudiés durant les dernières décennies [START_REF] Attia | Numerical study of flow and heat transfer of a non-Newtonian fluid on a rotating porous disk[END_REF] pour les nombreuses applications industrielles qui en découlent : des disques durs d'ordinateur, au design des générateurs de puissance ou de manière plus évidente dans le domaine de la rhéologie. Dans de précédents travaux [START_REF] Sahoo | Effects of partial slip, viscous dissipation, Joule heating on Von Kármán flow and heat transfer of an electrically conducting non-Newtonian fluid[END_REF], nous nous sommes intéressés à l'écoulement de Von Kármán [START_REF] Von Kármán | Über laminare und turbulente Reibung[END_REF] engendré par la rotation d'un disque infini dans un milieu fluide au repos et aux transferts de chaleur associés dans le cas d'un fluide conducteur de type Reiner-Rivlin avec des conditions de glissement partiel. Ces précédents résultats sont étendus ici au cas de l'écoulement d'un fluide en rotation uniforme au-dessus d'un disque infini au repos.

Ce problème a été largement considéré dans la littérature pour des fluides incompressibles newtoniens. Il a été étudié pour la première fois par Bödewadt [START_REF] Bödewadt | Die Drehströmung über festem Grund[END_REF] par des approximations de couche limite. La solution de Bödewadt montre que les effets de couche limite s'étendent jusqu'à ζ = 8, où ζ est la distance adimensionnée mesurée le long de l'axe de rotation du disque. Cette solution a ensuite été corrigée par Browning (non publié), qui nota une couche limite plus épaisse que pour l'écoulement de Von Kármán. Batchelor [START_REF] Batchelor | Note on a class of solutions of the Navier-Stokes equations representing steady non rotationally symmetric flow[END_REF] généralisa ces analyses au cas où le disque de rayon infini tourne à une vitesse angulaire ω constante dans un fluide tournant à une vitesse Ω différente, ce qui a été étudié plus tard numériquement par Rogers et Lance [START_REF] Rogers | The boundary layer on a disk of finite radius in a rotating fluid[END_REF]. Pour ω = 0, on retrouve le problème considéré par Bödewadt [START_REF] Bödewadt | Die Drehströmung über festem Grund[END_REF]. Nydahl [START_REF] Nydahl | Heat transfer for the Bödewadt problem[END_REF] a étendu le problème de Bödewadt pour introduire des transferts de chaleur. Ses résultats confirment ceux de Bödewadt alors que Rogers et Lance [START_REF] Rogers | The boundary layer on a disk of finite radius in a rotating fluid[END_REF] obtiennent une valeur de H ∞ = 1.3494 (valeur de la vitesse axiale loin du disque) très supérieure à celle précédemment obtenue. Le lecteur peut se référer au livre d'Owen et Rogers [START_REF] Owen | Flow and heat transfer in rotating disk systems[END_REF] pour une revue extensive des travaux sur l'écoulement de Bödewadt jusqu'en 1989. Si la littérature est très fournie sur le problème de Von Kármán dans le cas non-newtonien (voir [START_REF] Sahoo | Effects of partial slip, viscous dissipation, Joule heating on Von Kármán flow and heat transfer of an electrically conducting non-Newtonian fluid[END_REF]), elle l'est beaucoup moins pour l'écoulement de Bödewadt. On peut citer néanmoins les travaux de Kitchens et Chang [START_REF] Kitchens | Newtonian and non-Newtonian liquids rotating adjacent to a stationary surface[END_REF], qui ont étudié l'écoulement de Bödewadt pour un fluide non-newtonien du second ordre. La présente étude est une tentative pour combler ce manque en essayant notamment de quantifier les influences des paramètres non-newtonien et de glissement sur la structure de la couche limite de Bödewadt.

Méthode numérique 2.1 Modèle mathématique

On considère un fluide non-newtonien de type Reiner-Rivlin, pour lequel le tenseur des contraintes τ ij est relié au tenseur des déformations e ij [START_REF] Attia | Numerical study of flow and heat transfer of a non-Newtonian fluid on a rotating porous disk[END_REF] par :

τ ij = -pδ ij + 2µe ij + 2µ c e ik e kj
(1)

e jj = 0 (2)
où p représente la pression, µ le coefficient de viscosité et µ c le coefficient de viscosité croisée.

Le fluide occupe l'espace z > 0 au-dessus d'un disque fixe de dimension infinie situé en z = 0. L'écoulement est dû à la rotation solide du fluide à une vitesse angulaire constante Ω loin du disque. Le problème est écrit en coordonnées cylindriques (r, ϕ, z) avec une hypothèse d'axisymmétrie de l'écoulement : ∂ ∂ϕ ≡ 0. Les conditions aux limites de non-glissement pour le champ de vitesse s'écrivent :

z = 0, u = 0, v = 0, w = 0, z → ∞, u → 0, v → rΩ, p → p ∞ .
(

On introduit alors les transformations de Von Kármán [START_REF] Von Kármán | Über laminare und turbulente Reibung[END_REF] :

u = rΩF (ζ), v = rΩG(ζ), w = √ ΩνH(ζ), z = √ ν Ω ζ, p -p ∞ = -ρνΩP (4)
où u, v, w sont les composantes radiale, tangentielle et axiale de la vitesse et z la distance au disque fixe. Ces transformations permettent de réduire les équations de Navier-Stokes pour un fluide newtonien à un système d'équations différentielles ordinaires. La même démarche est ici adoptée dans le cas nonnewtonien. On définit alors le paramètre non-newtonien L = µ c Ω/µ. En considérant les approximations de couche limite usuelles, les équations de continuité et de quantité de mouvement se réduisent à :

dH dζ + 2F = 0 (5) d 2 F dζ 2 -H dF dζ -F 2 + G 2 -F -1 2 L [( dF dζ ) 2 -3 ( dG dζ ) 2 -2F d 2 F dζ 2 ] = 1 ( 6
)
d 2 G dζ 2 -H dG dζ -2F G -G + L ( dF dζ dG dζ + F d 2 G dζ 2 ) = 0 ( 7
)
d 2 H dζ 2 -H dH dζ -7 2 L dH dζ d 2 H dζ 2 + dP dζ = 0 (8)
Les conditions aux limites (3) deviennent alors :

ζ = 0 : F = 0, G = 0, H = 0, ζ → ∞ : F → 0, G → 1. (9)
Le fluide adhère partiellement à la paroi. Une généralisation des conditions partielles de glissement donne dans les directions radiale et azimutale :

u| z=0 = λ 1 T rz | z=0 , v| z=0 = λ 2 T ϕz | z=0 ( 10
)
où T rz , T ϕz sont les composantes physiques du tenseur des contraintes et λ 1 et λ 2 sont les coefficients de glissement :

λ = λ 1 √ Ω ν µ, η = λ 2 √ Ω ν µ. ( 11
)
En utilisant les transformations (4), les conditions aux limites de glissement partiel (Eq.10) deviennent :

F (0) = λ[F ′ (0) -LF (0)F ′ (0)], G(0) = η[G ′ (0) -2LF (0)G ′ (0)], H(0) = 0 F (∞) → 0, G(∞) → 1 (12)
Les conditions aux limites loin du disque restent ainsi inchangées.

Solution numérique

On résoud le système d'équations différentielles non-linéaires défini par les Équations (5-7) associé aux conditions aux limites (Eq.12) en adoptant le même schéma numérique du second ordre en espace que celui décrit dans [START_REF] Sahoo | Effects of partial slip, viscous dissipation, Joule heating on Von Kármán flow and heat transfer of an electrically conducting non-Newtonian fluid[END_REF]. Les Équations (5-7) sont discrétisées par des opérateurs de différence centrée pour les dérivées :

Le domaine d'intégration semi-infini ζ ∈ [0, ∞) est remplacé par un domaine d'extension fini ζ ∈ [0, ζ ∞ ). En
F i+1 -2F i + F i-1 h 2 -H i ( F i+1 -F i-1 2h ) -F 2 i + G 2 i -F i - 1 2 L [( F i+1 -F i-1 2h ) 2 -3 ( G i+1 -G i-1 2h ) 2 -2F i ( F i+1 -2F i + F i-1 h 2 )] -1 = 0 (13) G i+1 -2G i + G i-1 h 2 -H i ( G i+1 -G i-1 2h ) -2F i G i -G i +L [( F i+1 -F i-1 2h )( G i+1 -G i-1 2h ) + F i ( G i+1 -2G i + G i-1 h 2 )]
= 0 (14)

H i+1 = H i -h(F i + F i+1 ) (15)
où h = 0.05 est le pas d'espace supposé constant sur tout le domaine. Le maillage dans la direction axiale est donc défini par : ζ i = ih pour i = 0, 1, . . . n, avec n le nombre de points. Pour initier la résolution des Équations (13) and ( 14), non seulement les valeurs de F 0 et G 0 sont nécessaires mais il faut également connaître les valeurs de F 1 et G 1 . Celles-ci sont obtenues par des développements en série de Taylor autour de ζ = 0 pour F et G. Le calcul est mené jusqu'à ce que les valeurs de F , G et H soient connues en tout point du maillage. Il est à noter qu'il faut satisfaire les conditions aux limites asymptotiques (Eq.12). Pour satisfaire les conditions en ζ = ζ ∞ (Eq.9), une méthode de shooting associée à une méthode de Runge-Kutta d'ordre 4 est utilisée.

Résultats et discussion

La Figure 1 présente les valeurs de F , G et H pour une condition de non glissement (λ = 0) et un fluide newtonien (L = 0). Près du disque, la composante radiale de la vitesse F est purement centripète. Les comportements des trois composantes de vitesse confirment la représentation de l'écoulement donnée par Schlichting [START_REF] Schlichting | Boundary layer theory (7th Edition)[END_REF]. À partir du profil de G, on peut en déduire l'épaisseur de la couche limite ζ 99 , définie comme étant la position axiale pour laquelle la vitesse tangentielle est égale à la vitesse tangentielle du fluide à l'infini à 1% près. Pour L = 0 et λ = 0, la valeur classique ζ 99 = 8 (voir [START_REF] Owen | Flow and heat transfer in rotating disk systems[END_REF]) est retrouvée par la présente approche.

On considèrera dans la suite une rugosité uniforme, c'est à dire λ = η. Les variations de la composante radiale de la vitesse F avec le paramètre non-newtonien L et le paramètre de glissement λ sont présentées sur les Figures 2a et 2b respectivement. Sur la Figure 2a, il est clair que L a un effet dominant sur F près du disque. L'amplitude de l'écoulement radial centripète diminue près du disque pour des valeurs croissantes de L. Cependant, le profil de la composante radiale de vitesse reste toujours le même quelque soit la valeur de L : F est négative près du disque (écoulement centripète), change de signe loin de celui-ci (l'écoulement est alors dit centrifuge) avant de tendre vers sa valeur asymptotique. La Figure 2b montre la variation de F avec le paramètre de glissement λ(= η), pour des valeurs constantes des autres paramètres. L'effet du glissement sur F est important essentiellement près du disque. Son comportement est ensuite peu modifié pour des valeurs plus grandes de ζ. La valeur asymptotique de F est finalement atteinte sur une distance beaucoup plus courte que dans le cas de base pour L = λ = 0. Les Figures 2c et 2d présentent de même les variations de la composante tangentielle G de la vitesse avec L et λ(= η) respectivement. Le paramètre non-newtonien L a un effet spectaculaire sur G, loin du disque vers ζ = 2 (voir Fig. 2c). Une augmentation de L entraîne une décroissance de la vitesse tangentielle G proche du disque. Par contre, G augmente près du disque lorsque λ(= η) augmente (Fig. 2d). Les deux paramètres L et λ ont donc un effet opposé sur la composante de vitesse G. Les variations de la composante axiale H de la vitesse avec les paramètres L et λ sont données sur les Figures 2e et 2f respectivement. Les résultats montrent que les deux paramètres ont un effet similaire sur H. La composante axiale devient plus faible lorsque l'un de ces paramètres augmente. La valeur de cette composante loin du disque H ∞ est fortement réduite dans les deux cas comparée à la valeur de base 1.3494 obtenue pour L = λ = 0.

On peut déduire des profils de F et G la valeur du coefficient de moment C M , défini sous forme adimensionnée par : 3a), le coefficient de moment C M , en valeur absolue, augmente (resp. diminue) fortement pour des valeurs croissantes du paramètre de glissement λ lorsque λ < 0.6 (resp. λ > 0.6). Il tend rapidement vers 0 pour des très grandes valeurs de λ, ce qui signifie que le couple nécessaire pour maintenir le disque au repos est quasiment nul lorsque le paramètre de glissement est grand. Malgré la différence de conditions aux limites entre les deux problèmes, ces résultats confirment ceux obtenus par Sahoo [START_REF] Sahoo | Effects of partial slip, viscous dissipation, Joule heating on Von Kármán flow and heat transfer of an electrically conducting non-Newtonian fluid[END_REF] pour l'écoulement de Von Kármán. Sur la Figure 3b, le paramètre de glissement λ est fixé à 1 et le paramètre non-newtonien L varie entre 0 et 100 pour Re = 1. C M augmente alors en valeur absolue avec L. Dans ce cas, le couple requis pour maintenir le disque fixe est beaucoup plus grand que celui nécessaire pour maintenir le disque à ω dans le cas de l'écoulement de Von Kármán [START_REF] Sahoo | Effects of partial slip, viscous dissipation, Joule heating on Von Kármán flow and heat transfer of an electrically conducting non-Newtonian fluid[END_REF].

C M = -2πG ′ (0)[1 -2LF (0)] √ Re ( 16

Conclusion

L'écoulement dû à la rotation d'un fluide non-newtonien de type Reiner-Rivlin près d'un disque fixe a été étudié pour la première fois pour diverses conditions de glissement et différentes valeurs du paramètre non-newtonien. Les équations différentielles fortement non linéaires qui en résultent sont résolues par un schéma aux différences finies du second ordre en espace. Les effets combinés du glissement λ = η et du paramètre non-newtonien L sur le champ de vitesse ont été regardés en détails. Il est intéressant de noter que λ et L ont des effets similaires sur les composantes radiale et axiale de la vitesse. Pour des valeurs croissantes de λ ou L, un léger épaississement de la couche limite du disque fixe est observé. Par contre, ces deux paramètres ont des effets opposés sur la composante tangentielle de la vitesse et sur le coefficient de moment. Le couple ainsi requis pour maintenir le disque fixe tend vers 0 pour des valeurs croissantes de λ et augmente en valeur absolue avec L, montrant un effet dominant du paramètre non-newtonien.

Figure 1 -Figure 2 -

 12 Figure 1 -Profils de vitesse pour un fluide newtonien avec λ(= η) = 0.

Figure 3 -

 3 Figure 3 -(a) Variations du coefficient de moment C M (a) avec λ pour L = 1 et Re = 1 ; (b) avec L pour λ = 1 et Re = 1.

  pratique, ζ ∞ doit être choisi suffisamment grand pour que la solution numérique tende vers les conditions aux limites à l'infini. Ici, on fixe ζ

∞ = 27, valeur qui est à comparer avec celle utilisée pour le problème de Von Kármán : ζ ∞ = 10 [2].


Guillaume Fontaine

Sébastien Poncet

Eric Serre

Approche

Approche multidomaine pseudospectrale pour la résolution des équations de Navier-Stokes en géométrie cylindrique

come

Les écoulements confinés en rotation sont présents dans un grand nombre d'applications industrielles liées aux machines tournantes et mettent en jeu des comportements physiques présents également dans des systèmes ouverts de plus grande échelle (circulations océanique et atmosphérique) qui demeurent pertinents même au-delà des aspects spécifiques de rotation [START_REF] Launder | Laminar, transitional and turbulent rotor-stator flows[END_REF]. Les sujets majeurs abordés actuellement concernent les scénarios de transition à la turbulence et le contrôle ainsi que la modélisation des transferts de masse et de chaleur dans les écoulements turbulents. La précision des solutions requiert l'utilisation de méthodes d'ordre élevé qui assurent une très grande précision de la solution sur des temps de simulation longs avec une très faible dissipation et dispersion numérique. Ces études nécessitent des résolutions spatiales et temporelles fines afin de prendre en compte avec précision les structures de petite échelle très anisotropes qui se développent au voisinage des parois et qui sont responsables de la production de turbulence. Ces résolutions ne sont accessibles qu'aux travers d'algorithmes parallèles optimisés sur des supercalculateurs ou des fermes de PC. Différentes méthodes existent aujourd'hui dans la litérature depuis les éléments finis hp jusqu'aux éléments spectraux et qui satisfont les critères mentionnés ci-dessus (voir dans [START_REF] Kardianakis | Spectral/hp elements methods for computational fluid dynamics[END_REF]). Néanmoins, ces méthodes qui permettent de traiter des géométries extrêmement complexes sont difficiles à mettre en oeuvre et relativement couteuses en temps de calcul, ce qui ne les rend pas optimales pour l'étude d'écoulements en géométries simples ou semi-complexes comme celles qui nous intéressent dans ce travail. Dans la méthode proposée, les composantes de la vitesse aux interfaces artificielles entre sous-domaines sont calculées par une technique de matrice d'influence, construite à partir de bases respectant la contrainte d'incompressibilité. La matrice d'influence a été appliquée depuis plus de dix ans avec succés en géométrie cartésienne [START_REF] Sabbah | A Divergence-Free Multidomain Spectral Solver of the Navier-Stokes Equations in Geometries of High Aspect Ratio[END_REF] et cylindrique uniquement pour des écoulements axisymétriques décrits en formulation fonction de courant vorticité [START_REF] Raspo | A spectral multidomain technique for the computation of the Czochralski melt configuration[END_REF]. La méthode développée ici est de type collocation spectrale Chebychev (Chebychev-Fourier en 3D) et utilise des solveurs directs de l'équation de Helmholtz et du problème de Stokes généralisé.

Modélisation géométrique et mathématique

Une des difficultés dans la modélisation des écoulements en géométrie cylindrique est liée aux effets de courbure en 1/r. Une cavité annulaire composée de deux disques parallèles est définie par ses rayons intérieur R int et extérieur R ext ainsi que par son espacement interdisque H. Ces caractéristiques permettent de définir deux paramètres géométriques globaux caractérisant le rapport d'aspect et la courbure, soient respectivement :

R m = R ext + R int R ext -R int , L R = R ext -R int H (1)
Le mouvement du fluide incompressible est gouvernée par les équations de Navier-Stokes écrites en variables primitives et en coordonnées cylindriques (r, θ, z) :

   ∂V ∂t + V.∇V = -∇p + ν∆V + F dans Ω V = W sur Γ = ∂Ω ∇.V = 0 sur Ω = Γ ∪ Ω (2)
où V est le vecteur vitesse de composantes (u, v, w) respectivement dans les directions radiale, azimutale et verticale, p est la pression, ν la viscosité cinématique, F le terme force et ∆ l'opérateur Laplacien.

Modélisation numérique

La discrétisation des équations est basée sur le travail de Raspo et al. [START_REF] Raspo | Spectral projection methods for the simulation of complex three-dimensional rotating flows[END_REF]. La discrétisation spatiale est basée sur une méthode collocation-Chebychev (dans les directions non-homogènes (r, z)) et Galerkin-Fourier en azimut. La discrétisation temporelle est basée sur un schéma temporel semi-implicite d'ordre 2 avec une évaluation Adams-Bashforth explicite des termes convectifs, et d'un schéma d'Euler retardé pour les termes visqueux. Soit,

   3V n+1 -4V n +V n-1 2δt + 2V.∇V n -V.∇V n-1 = -∇p n+1 + ν∆V n+1 + F n+1 dans Ω V n+1 = W n+1 sur Γ ∇.V n+1 = 0 sur Ω (3)
où δt est le pas de temps.

Le couplage vitesse-pression est résolu au moyen de l'algorithme de prédiction-correction. Un champ de vitesse prédit non-solénoïdal est tout d'abord calculé en utilisant un prédicteur de pression estimé à partir de l'équation de Poisson obtenue en prenant la divergence de l'équation de quantité de mouvement. Par la suite ce champ prédit est projeté sur un espace de champs à divergence nulle au moyen d'un correcteur de pression.

Cet algorithme se ramène à la résolution successive de problèmes d'Helmholtz bidimensionnels pour chaque mode de Fourier.

Matrice d'Influence de Continuité (MICO)

La technique est présentée ici pour une découpage en domaines dans la direction radiale uniquement (Figure 1) mais pourrait s'étendre à la direction verticale sans aucune difficulté conceptuelle. La géométrie cylindrique nécessite la définition de rapports de courbure (R mi ) et d'allongements (L Ri ) locaux dans chacun des sous-domaines compatibles avec les paramètres globaux définis précédemment. Ces paramètres permettent de transformer localement le domaine géométrique couvert par le domaine dans l'intervalle [-1, 1] de l'approximation Chebyshev.

Résolution d'un problème d'Helmholtz

Dans un souci de simplicité de notation, considérons un domaine en 2 dimensions (r et z) divisé en 2 sous-domaines afin d'établir les étapes de ce schéma de projection multidomaine. Le raisonnement serait parfaitement identique pour un nombre plus élevé de domaines.

   ∆ i φi -σ i φi = Ŝi dans Ω i A i φi = b i sur Γ i ∩ Γ φ1 = φ2 = φ sur Γ i ∩ Ω = ξ (4)
avec A i = Id pour les composantes de vitesse et A i = ∂ ∂n pour la pression. (n représente la dérivée normale à la frontière où la condition limite est définie).

Bien entendu, on ne connaît pas la valeur de φ à imposer sur l'interface ξ. Cette valeur est déterminée au moyen de la matrice d'influence en recherchant la solution φi sous la somme de la solution φi d'un problème à conditions de Dirichlet homogène à l'interface et de la solution d'un problème stationnaire φi calculé en préprocessing. A chaque itération temporelle on doit résoudre le problème "homogène" suivant :

   ∆ i φi -σ i φi = S i dans Ω i A i φi = b i sur Γ i ∩ Γ φ1 = φ2 = 0 sur ξ (5)
dont la solution φ1 ∪ φ2 sur le domaine global est continue à travers ξ (car φ1 (z ∈ ξ) = φ2 (z ∈ ξ) = 0) mais qui présente un saut de dérivée normale à travers la frontière noté H(z ∈ ξ). Ce saut de dérivé est alors corrigé au moyen de la matrice d'influence M calculée à partir de dérivées des solutions φik des problèmes d'Helmholtz élémentaires stationnaires définis sur le sous-domaine i par :

∀k ∈ [2, N P Z -1]    ∆ i φik -σ i φik = 0 dans Ω i A i φik = 0 sur Γ i ∩ Γ φik (z l ∈ ξ) = δ kl sur ξ (6)
où N P Z est le nombre de points sur la frontière (Figure 2). La matrice d'influence M est l'endomorphisme qui à la valeur de la solution recherchée φ = λ sur ξ associe le saut de dérivée H de la solution homogène. La contrainte imposée pour développer cet endomorphisme est qu'on veut ajouter à la solution homogène (C 0 mais pas C 1 ) une combinaison linéaire de coefficients λ des solutions élémentaires définies plus haut qui conduit à une solution finale C 0 et C 1 à travers la frontière. Par inversion de cette matrice d'influence on pourra alors retrouver par multiplication matricielle la valeur de λ à imposer sur la frontière pour obtenir φ(ξ), à partir du saut de dérivée H de la solution homogène φ à travers la frontière. ∀θ ∈ [1, N P A], les solutions sont obtenues au moyen de la résolution du système linéaire suivant :

λ = M -1 H (7)
Ainsi, pour un problème à 3 dimensions, on résoudra le problème aux frontières précédent autant de fois qu'on a de modes de Fourier dans le problème. Dans un calcul à valeurs propres réelles, cela se traduit par l'utilisation de N P A/2 matrices d'influence. Dans le cas où les valeurs propres sont complexes, on utilisera N P A matrices d'influence.

Remarque 1 -Pour la pression et la pseudo-pression, les solutions sont définies à une constante additive près, ce qui se traduit par un opérateur Laplacien non inversible. En monodomaine on filtre ce mode. Cependant on ne retrouve pas cette singularité en multidomaine, car on implémente des conditions-limites de type Dirichlet sur certaines frontières. Cependant ce problème se retrouve sur la matrice d'influence du mode 0 qui n'est à son tour pas inversible. L'équivalent du filtrage effectué en monodomaine est alors assuré par une technique de diagonalisation sur la matrice d'influence du mode 0 : on identifie l'indice de la valeur propre nulle de la matrice d'influence, et on annule la composante correspondante (dans la base de diagonalisation) dans le saut de dérivée.

Besançon, 29 août au 2 septembre 2011

Remarque 2 -Dans le cas considéré précédemment (2 sous-domaines, problème plan), M est diagonale par bloc. Dans le cas d'un nombre plus élevé de sous-domaines, elle devient tridiagonale par blocs. En effet, la correction des sauts de dérivées sera aussi influencée par les solutions élémentaires définies sur les frontières situées de l'autre côté du sous-domaine considéré.

Précision de la méthode

La précision du nouveau schéma multidomaine a été vérifiée à partir de tests de consistence utilisant les champs de vitesse et de pression stationnaires suivants :

      
U ana (r, z, θ) = 1 2π sin(πr) 2 sin(2πz)cos(θ) V ana (r, z, θ) = -1 2π sin(πr) 2 sin(2πz)sin(θ) W ana (r, z, θ) = 1 2πG sin(πz) 2 sin(2πr)cos(θ) P ana (r, z, θ) = [cos(πz) + cos(πr)]cos(θ) L'erreur en norme L 2 est tracée pour différents maillages dans la direction radiale et nombres de sousdomaines. Une cavité de rapport de forme L R = 1 et de rayon de courbure R M = 2 est considérée. Le nombre de Reynolds est fixé à 500 et le pas de temps à 10 -1 Ω -1 . Le nombre de points dans la direction azimutale est fixé à 8. Dans une premier temps, le nombre de points dans les directions radiale et axiale varie de de 11 à 35 dans chacun des deux sous domaines. Puis, le nombre de points par sous domaine est fixé à 25 dans la direction radiale et à 9 dans la direction axiale et c'est le nombre de sous domaine NSD qui varie de 2 à 9. Dans la première situation, on retrouve une convergence exponentielle de l'erreur caractéristique de la précision spectrale. Lorsque le nombre de degré de liberté est accru par le nombre de sous-domaine la convergence de l'erreur de troncature est seulement algébrique. Ce résultat est attendu dans le sens où la régularité de la solution globale étant uniquement C 1 aux interfaces, lorsque celles-ci deviennent trop nombreuses la convergence spectrale de la solution globale est perdue.

Finalement, une simulation 3D en cavité rotor-stator à Re = 25000 montre la capacité de la méthode pour obtenir des solutions non laminaires (Figure 4). Un mode 3D d'instabilité est observé sur la couche du stator comme attendu par l'expérience et les analyses de stabilité. Ce qui est intéressant est de voir que ces pertubations sous forme de rouleaux traversent la frontière entre les domaines sans être altérées.

Remarques de conclusion

Nous avons présenté une méthode multidomaine pseudospectrale qui résoud les équations de Navier-Stokes en géométrie cylindrique pour un fluide incompressible. Cette méthode est "divergence-free" et conserve la précision spectrale de l'algorithme, à savoir la convergence exponenetielle de l'erreur de troncature lorsque le nombre de sous-domaines restent raisonnable par rapport à l'extension géoémétrique. L'utilisation d'une technique de matrice d'influence conserve le caractère direct (non itératif) de l'algorithme de projection ce qui le rend particulièrement efficace en terme de coût de calcul. Des calculs préliminaires en configuration rotor-stator ont montré son aptitude à simuler des écoulements 3D instationnaires en accord avec les résultats de la littérature. La parallèlisation MPI de notre algorithme est en cours ce qui permettra très prochainement d'accroître la résolution afin de simuler des nombres de Reynolds pertinents avec l'application industrielle.

  Mots clefs : Méthode pseudospectrale, Ecoulement tournant, Matrice d'influence de continuité 1 Introduction

Figure 1 -

 1 Figure 1 -Maillage et cavité annulaire (à gauche). Exemple de décomposition multidomaine dans la direction radiale (à droite).

Figure 2 -

 2 Figure 2 -Illustration d'un cas à 2 sous domaines.

Figure 3 -

 3 Figure 3 -Evolution de l'erreur de troncature L 2 avec le nombre de points (gauche) et le nombre de de sous-domaines (droite).

Figure 4 -

 4 Figure 4 -Solution 3D en cavité rotor-stator à Re = 25000, R M = 1.8, L R = 6.26. Rouleaux dans la couche limite du stator (à gauche) se développant en structures spiralées (à droite).