O Crespo

J E Bergez
email: jbergez@toulouse.inra.fr

F Garcia
email: fgarcia@toulouse.inra.fr

P2 hierarchical decomposition procedure: application to irrigation strategies design

Keywords: Simulation optimization, crop model, irrigation management

Optimization by simulation of agricultural practices can help to improve irrigation water use efficiency. This work introduces an efficient hierarchical decomposition method to design irrigation management strategies that is modelled as a continuous stochastic problem. Various combinations of selection (greedy, Pareto-based), division (middle, pivot, maximization) and evaluation techniques (global, standard deviation) were tested. We present results of an 8-continuous-parameter irrigation strategies design. Two criteria were chosen to evaluate the different combinations: the achieved direct margin, and the number of simulation runs that were needed to reach it. Selection techniques impacted the resolution time, while the evaluation techniques impacted the direct margin efficiency. Based on the two former criteria, the trade-off combination of greedy selection, pivot partition and average value evaluation appeared to be the most efficient to design irrigation strategies.

are already very heavily exploited [START_REF] Smith | The application of climatic data for planning and management of sustainable rainfed and irrigated crop production[END_REF]). Where irrigated agriculture is developed, water used for irrigation can represent more than 90% of the water consumption. Thus we investigate the possibility of increasing the agricultural production by increasing the efficiency of irrigation water use without increasing the quantity of water used.

Agricultural practices have been modelled for some time using decision rules [START_REF] Papy | Farm models and decision support: a summary review[END_REF], [START_REF] Aubry | Modelling decision-making processes for annual crop management[END_REF], [START_REF] Shaffer | Rule-based management for simulation in agricultural decision support systems[END_REF]). Such modeling allows uncertain events such as weather to modify planned technical action on the biophysical system in an adaptive way. For example, the weather influences the sowing date. This type of modeling aims at representing the farmer's behavior: actions are decided with regard to both constraints and goals, and are in addition modified depending on context and local conditions. MODERATO (Bergez et al (2001a)) has been developed using such a modeling framework. The cultural operations are decided with elementary conditional rules. These rules involve condition parameters such as cumulative thermal units, water deficit or irrigation amounts. Its purpose is to evaluate current irrigation strategies for corn and to propose improved strategies.

To optimize management strategies, two main approaches can be applied: either controlbased optimization or simulation-based optimization [START_REF] Bergez | Representing and optimizing management decisions with crop models[END_REF]). The former approach was used with MODERATO, by testing stochastic dynamic programming and reinforcement learning methods for identifying optimal decision rules (Bergez et al (2001b)). Input variables to the optimization problem were the condition parameter of the decision rules. It gave poor results due to the complexity of the problem and the large variable space to optimize. The simulation-based optimization approach deals with large continuous parameters problems and involves explicit techniques to handle with uncertainty. Various methods have been introduced and developed in this field. Reviews from [START_REF] Andradóttir | A review of simulation optimization techniques[END_REF], [START_REF] Ólafsson | Simulation optimization[END_REF] or [START_REF] Fu | Simulation Optimization: A Review, New Developments, and Applications[END_REF] provide a clear insight of this diversity: gradient-based procedures, stochastic approximation, sample path optimization, response surface methodology, ranking & selection, and a few branching approaches dealing with continuous or discrete input variable problems. These methods include stochastic search techniques moving from a current best solution to the next, or deterministic search techniques relying on an approximation of either the response or the objective function.

Among them, a few branching procedures have been introduced in which the feasible region is iteratively broken down before separating promising subsets from unpromising ones. For example, [START_REF] Norkin | On Optimal Allocation of Indivisibles under Uncertainty[END_REF] proposed a version of the branch-and-bound method for discrete parameter optimization where promising and unpromising regions are selected, based on objective function value bounds. This method has recently been introduced for continuous parameter optimization. As another example, PBnB from [START_REF] Prasetio | A probabilistic branch-and-bound scheme for simulation optimization[END_REF] seeks to identify the most promising sub-region containing a desired solution by an iterative pruning-or-branching procedure. Branching action is performed according to the probability of reaching a targeted function measurement, while pruning is performed according to the significant difference of statistical tests between the most-and least-promising regions.

The aim of this paper is to present the P2 algorithm, a development of P2p introduced in [START_REF] Bergez | A hierarchical partitioning method for optimizing irrigation strategies[END_REF]. P2 is a hierarchical decomposition procedure dedicated to continuous parameter optimization, and is one of the simulation-based optimization methods. Promising regions are selected by a heuristic, relying both on the average value and on the standard deviation of simulation responses. The promising region is then divided into two parts, and each of them is evaluated by simulation. Different options for the three steps of selection, division and evaluation were tested on a corn crop irrigation management problem in order to identify the combination of techniques that lead to the most efficient algorithm.

We first describe the principle of the P2 algorithm and the various techniques that are tested for selection, division and evaluation. In section 3, we give some information on the MODERATO bio-decisional simulation model used to simulate the irrigation strategies and we set up the simulation experiments. We present the results of the combinations of techniques and we eventually discuss the best P2 algorithm options and the optimized irrigation strategies.

From P2p to Palgorithm

The P2p algorithm is based on the DIRECT [START_REF] Jones | Lipschitzian optimization without the Lipschitz constant[END_REF]) and the MCS [START_REF] Huyer | Global optimization by multilevel coordinate search[END_REF]) algorithms which have been developed for deterministic optimization. The P2p algorithm is dedicated to large continuous input variable problems. It includes a single objective function, which can be the combination of weighted multiple objectives. P2p belongs to the family of stochastic branching methods, like stochastic branch-andbound or nested partitions methods. It is based on a hierarchical decomposition of the decision space into a binary tree.

This algorithm is dedicated to the optimization problem max

θ ∈Θ J(θ ), ( 1
)
where Θ is the set of possible decisions, and J(θ ) is a performance function that can not be computed analytically. The evaluation of J(θ ) thus relies only on the responses of stochastic simulation runs. The formulation of this kind of optimization involves input continuous D-dimensional variables θ constrained to be contained within a feasible region Θ ⊂ R D . An objective function is defined on these variables, such that J : θ → R, where (θ , ω)] and ω is the uncontrollable input variables vector of the stochastic system. The general formulation of simulation-based optimization is:

J(θ ) = E[L
max θ ∈Θ J(θ ) = E[L(θ , ω)] . (2)
The decision space Θ is a hyper-rectangle, or region, of R D . The P2 optimization aims at finding small hyper-rectangles included in Θ which contain the decision vector that maximizes the expected value of L(θ , ω). We assume a minimal size of these small regions, defined by the user for every dimension as the width p step d of the dimension d ∈ D. Let us call pending regions the regions that are still divisible (or breakable). The principle of the P2 algorithm is described in the Algorithm 1. Initialization allocates the initial decision space as the single pending region of the list pendingRL of all the pending regions. The first step consists of SELECTING out from the pending regions list the region which is potentially optimal: we call it the promising region. The second step DIVIDES this promising region into two parts. These two parts are offspring regions and are collected in the o f f springRL list. During the third step, each of these regions is sampled, simulated and indices are EVALUATED. Eventually, the pending regions list is updated, and the three previous steps are repeated until stopping criteria are reached or the pending regions list is empty. The P2 algorithm main stopping criterion is achieved when there is no more pending region, i.e. there is no more potentially optimal region which could still be divided for further exploration. This division limit is explain in the next section 2.2. For practical reasons, time limits and/or simulation number limits are usually added to the previous stopping criterion. Figure 1 illustrates the three main steps considering a 2-dimensional decision space related to the decision tree produced. The P2p algorithm included the β -selection, middle partition, average value and standard deviation evaluations. The major improvements leading from P2p to P2 consisted in proposing, testing and selecting combinations of new techniques.

Selecting the promising region

The selection step is involved in the process efficiency according to its ability to lead straight to the optimal region. At each iteration, the promising region to be selected is one of the pending regions. Each pending region is sampled, simulated and indices are computed dur-ing the evaluation step. The focus here is to choose one promising region r * within the pending regions set, based on these indices. Three selection techniques have been tested, based on two indices: the average value and the standard deviation of the simulated responses. We extended the consideration of these two indices to score 1 and score 2 which could be either the average value, the standard deviation or the maximum, the minimum or the median value for example. a. The greedy selection consists in choosing the promising region considering only the expected value criterion to be optimized. Considering score 1 as a generalization of the criterion to be optimized and maximized, the promising region r * will be chosen such that r * = arg max r [score 1 (r)]. b. We propose a trade-off between greedy and systematic search. The β -selection takes into account both the average value of the region and the standard deviation of the sampled decision vectors. Let us define score 1 as the average value and score 2 as the standard deviation. The figure 2 represents the set of the convex non dominated1 regions, in which we select the promising one. We use an approach based on multiobjective optimization, which looks for regions with a large standard deviation when regions are large, and for regions with a large average value when regions are small [START_REF] Bergez | A hierarchical partitioning method for optimizing irrigation strategies[END_REF]).

Considering score 1 and score 2 as a generalization of the criteria to be optimized and maximized, let us define β as following.

β (r) = score 2 (r) + score 1 (r). tan(γ) with γ = current depth maximum depth + ε . π 2 (3)
current depth is the tree depth of the last region of interest, and maximum depth is the depth of an unbreakable region assuming middle partitioning. ε is an infinite small value excluding π/2 from reachable values. As γ ∈ [0..π/2[ is tree depth dependant, maximizing β is selecting one region among all non dominated ones as shown on figure 2. Then the promising region r * will be chosen among all the pending regions r such that r * = arg max r β (r). c. The previous approach can only reach the convex non dominated front. In order to generalize the method to non convex front, the Concβ -selection combines the previous method with non convex dominance 1 . The selection technique still gives priority to the standard deviation when regions are large, and to average value when regions are small. In that case only the generalized criteria score 1 and score 2 have to be minimized. pro j 1 and pro j 2 are the projections of respectively score 1 and score 2 onto the line angled (π/2) -γ such that pro j 1 (r) = score 1 (r) * cos( π 2 -γ) and pro j 2 (r) = score 2 (r) * sin( π 2 -γ). The promising region r * is still part of the Pareto's set and will be chosen such that

r * = arg max r (min[pro j 1 (r), pro j 2 (r)]) . (4)
The figure 3

Breaking down the promising region

The division step is involved in the resolution time that is needed and in the global optimization efficiency: producing small regions could trap the process in local optima, and producing large regions could be time-consuming. We proposed in P2p to break the D-dimensional promising region into 2 D sub-regions. The depth of the decision tree to reach an unbreakable region is small, but at every nodes one region is taken out from the pending regions list while potentially 2 D are put in. The number of simulations needed to evaluate every new pending regions is then enormously expended at every division. We chose in the P2 method to cut down the D-dimensional promising region into only two parts to save simulations. Hence we first need to choose the parameter p * to be cut. We proposed to choose this parameter with the largest relative range such that

p * = D arg min d=1 ((p ⊕ d -p ⊖ d )/p step d ) (5)
where p ⊕ d (p ⊖ d ) is the upper (lower) bound of dimension d and p step d is its smallest feasible range. Thus three division techniques have been tested. a. The equal partitioning is the easiest way to proceed. Without taking into account the results from samples, the parameter p * is simply cut down right in the middle, producing two equal sub-regions. Instead of blindly producing two sub-regions, two other techniques aim at producing two regions as different as possible. b. The pivot partition chooses the decision vector which divides weak from better responses. The selection of the pivot vector is based on vector indices related to the probability of being the best representation of the threshold dividing the γ% of the 'best' responses from the 1 -γ% of the 'worst' responses. Two user defined indices γr and γs are used. c. The maximization of the difference technique is based on a binary classification and regression tree field [START_REF] Breiman | Technical Note: Some Properties of Splitting Criteria[END_REF]). The best split is taken as the maximizer of a "goodness-of-split" function and is chosen as the one maximizing the difference between potential sub-region indices γr.

Evaluating eligible regions

The system performance is not directly available and we therefore rely on sampled points

L(θ , ω) of the region of interest. N sampled point vectors θ i (i ∈ [1..N]) are simulated M times ( j ∈ [1..M]
) in each new pending region. N is linearly related to the region width and M is a parameter to be chosen. The aim is to compute, based on N ×M sampled performance measures, the indices of the pending regions that are used in the next iteration to select the promising region.

a. As the optimization definition 2 involves the expected value of the sampled performance measures, we first compute the average value Ĵ(θ ) as an approximation of J(θ

) = E[L(θ , ω)]. J(θ ) ≈ Ĵ(θ ) = 1 N 1 M N ∑ i=1 M ∑ j=1 L(θ i , ω j ). (6)
b. The global standard deviation V g can be used in the selection techniques involving two indices.

V g = 1 N 1 M N ∑ i=1 M ∑ j=1 Ĵ(θ ) -L(θ i , ω j ) 2 . (7)
c. If we consider the sampled performance measures as disturbed responses of the objective function J, then the idea of the third proposition is to use the standard deviation of the objective function without considering the standard deviation of the disturbance. We tested two indices V p1 and V p2 . They are computed to focus on the standard deviation due to the controllable input parameters (θ ), excluding the standard deviation due to the uncontrollable input parameters (ω). The aim of this paper is not to establish the formulation of these standard deviations. Nonetheless, V p1 is based on the hypothesis that the response L(θ i , ω j ) depends only on the θ i and ω j effects, although V p2 is based on the hypothesis that the response L(θ i , ω j ) depends on θ i , ω j and the residue due to the combination of θ i with ω j .

Choosing algorithm options by application to irrigation management

We have presented various techniques of selection, division and evaluation of decision hyper-rectangles. The aim of this section is to present the experimental plan which guided our choice of the best combination of technique. First we introduce the irrigation strategy simulator and the irrigation strategies design application. Then, we define the system to be optimized and the experiments we ran.

3.1 MODERATO Simulator MODERATO (Bergez et al (2001a)) is a model aimed at evaluating current irrigation strategies for corn and at proposing improved strategies. It combines a dynamic and biophysical corn crop model with a dynamic decision model. The crop model is described in [START_REF] Wallach | Parameter estimation for crop models a new approach and application to a corn model[END_REF]. The decision model consists of a set of decision rules for different management decisions, and especially for irrigation management decisions. The crop model and the decision model interact every day. The crop model updates the state variables each day and passes their values to the decision model together with the explanatory variables of that day. Within that collection of variables are the indicators of the decision rules. Then, the decision model evaluates the decision rules to decide if a management action is to be taken. If so, this information is passed back to the crop model (for example, amount of water or sowing density).

The timing of irrigation includes for example these rules: Starting irrigation This rule determines the starting day to begin irrigation during the growing season and the water amount for the first irrigation round. Next irrigation round This rule is invoked after a round of irrigation has been terminated. It determines when to start the next round and the irrigation amount for rounds after the first. Stopping irrigation This rule is invoked at the end of an irrigation round. It has one of these three conclusions: either (1) the previous round of irrigation was the last, or (2) another round of irrigation is to be performed and will be the last, or (3) we will re-invoke this rule after another round of irrigation. Granting that the next round is the last, the amount of irrigation is given.

Many of the rules in MODERATO are based on the general form:if (cond 1a ∨ cond 1b) ∧ (cond 2a ∨ cond 2b) then decision; de f ine amount;where cond 1a and 1b concern crop development while cond 2a and 2b refer to water status in the soil. The first condition in each pair (cond 1a and 2a) uses meteorological variables as indicator variables while the second condition in each pair is based on state variables. The user can choose to ignore one of the two conditions in each part of the premise.

Case study

The comparison between the developed options was performed on an eight-parameter strategy as follows.

The main irrigation period starts from T1 as soon as the soil water deficit reaches D1. An amount I1 is applied. Once an irrigation cycle ends, a new cycle starts when the soil water deficit reaches D2. An amount I2 is applied. For the irrigation cycle following T3, if the soil water deficit is greater than D3 before this irrigation cycle starts, a last irrigation cycle is performed; otherwise the irrigation program ends. An amount I3 is applied. The other cultural operations are given in table 2. The irrigation equipment used for the study allows a 3.5 mm/day flow rate. A 180 mm limitation of available water is applied. No flow rate restrictions are imposed during summer except those due to the equipment.

All simulations were performed using a medium clay-silt soil : 0.8m deep, with clay accumulation at depth, locally called "Boulbènes moyennes" (fluvisol). This type of soil is representative of a large area of the Midi-Pyrénées and has a 150 mm cumulative available water capacity. The soil was assumed to be at field capacity at the beginning of the simulation, namely the 1 st of January. Climates used are part of the observed weathers recorded at the Toulouse-Blagnac meteorologic station within 1949 and 1997. On average, July and August receive a total of 92 mm rainfall and the cumulative potential evapotranspiration (ET 0 ) is 290 mm. The average evaporative moisture deficit (ET 0 minus rainfall) for this two-month period is around 200 mm. However, there is a large variation in rainfall during the two summer months as it ranges from 30 to 240 mm, underlining the unpredictable nature of rainfall in the area. Cumulative ET 0 is less variable, ranging from 235 to 372 mm. The objective function to be maximized is the expectation of the direct margin (i.e. the gross margin minus specific costs for a given activity, here irrigation). The direct margin for irrigation can be written as a weighted sum of multiple criteria:

L(θ i , ω j ) = a(θ i , ω j ).B -[C + d(θ i , ω j ).E + f (θ i , ω j ).G] (8)
where L(θ i , ω j ) is the direct margin for climate ω j and the strategy θ i , a(θ i , ω j ) is the grain yield obtained under climate ω j and using the strategy θ i , B is the selling price for corn, C is the operational costs for corn production, d(θ i , ω j ) is the amount of water used under climate ω j and using the strategy θ i , E is the cost of irrigation water, f (θ i , ω j ) is the number of irrigation cycles performed under climate ω j and using the strategy θ i and G is the cost of carrying out a new irrigation cycle. The average selling price for maize (grain) is assumed to be 106.71 ¤/Mg in the Toulouse area. Operational costs (seed, weeding, fertilizer, insurance) are assumed to total 327.77 ¤/ha. The cost of irrigation water is assumed to be 0.76 ¤/mm and the setting up of a new irrigation cycle is assumed to be 7.62 ¤.

Experiments

Simulations were run with a Bi-processor of 3 GHz each, and 2 Go of RAM with Windows XP operating system. Optimization took about 3 hours and 40 minutes for 2 million of simulation runs including less than 9 seconds for the P2 procedure.

We took into account 10 replications for each alternative. We particularly focused on robustness of an alternative through the minimum, maximum and average direct margin values of these 10 replications. The initial feasible region is defined in table 3 as the ranges of the different parameters from the strategy described in section 3.2. Table 3 The eight parameters of the irrigation strategy to be optimized. min and max show the range of each parameter within which the optimum is sought. A step is the minimum feasible range of the parameter. Some procedure parameters were set for all experiments. The maximum number of simulations was set to 2 million. β -and concβ -selections include a 20% probability of randomly choosing the promising region and the sampling follows a uniform distribution. As some of the partitioning alternatives need indices, we called them γr and γs for the pivot partitioning, and γr for the maximum difference partitioning. We tested all possible combinations with γr ∈ {0.2, 0.5, 0.8} and γs ∈ {0.2, 0.5, 0.8}.

To compare alternatives, we defined some criteria in order to measure algorithm efficiency. Each time the direct margin average value increased, we stored it as well as the number of simulation runs used to achieve it. These stored values presented a strictly increasing curve in the simulation runs/direct margin criteria space. The general shape of this curve was made up of two distinct phases: one with a large improvement in direct margin average with a few simulation runs, and another with a tiny improvement with an infinite number of simulation runs. Therefore, we used simulation count and direct margin average corresponding to 95, 97.5, 99 and 99.5% of the best value eventually achieved within 2 million of simulation runs. These indices allowed information about the best performances reached by the alternative to be kept, while focusing on the first increasing phase of interest.

The number of runs criterion ensured that the algorithm efficiency does not depend on the computer's capabilities.

To recognize different techniques and their parameters, we used a XY Z code notation described in table 4. The first item stands for selection alternative: gs for greedy-selection, bs for β -selection or cs for concβ -selection. The second stands for partition alternative: ep for equal-partition, pp for pivot-partition or mp for maximize difference-partition. The last stands for the second evaluation index: ge for global standard deviation V g , p1e for point standard deviation V p1 or p2e for point standard deviation V p2 evaluation. Alternative bs/ep/ge has been used as a test reference. We first compared separately selection, division and evaluation alternatives from bs/ep/ge to assess the impacts of single techniques. Then, we tested a few hybrid alternatives that seemed to be the most interesting according to previous observations.

Results

Observation of P2 alternatives

Selection alternatives

Experimented selection alternatives are β -(bs/ep/ge), concβ -(cs/ep/ge) and greedy-(gs/ep/ge) selections.

Figure 4 displays the alternative averaged indices of 95, 97.5, 99 and 99.5%. Alternative gs/ep/ge was faster and reached a slightly lower direct margin than bs/ep/ge. Alternative cs/ep/ge both achieved a lower direct margin and needed many more simulations to reach the best region. The gs/ep/ge's deviation along the number of simulation is the smallest.

Division alternatives

The tested partitioning alternatives are equal-(bs/ep/ge), pivot-(bs/pp/ge) and max diff-(bs/mp/ge) partitioning. We saw earlier (section 2.2) that bs/pp/ge needs the thresholds γr and γs, and bs/mp/ge needs the threshold γr. Alternatives bs/mp/ge 02 and bs/mp/ge 08 respectively stand for bs/mp/ge with γr = 0.2 and bs/mp/ge with γr = 0.8. Alternative bs/pp/ge 08 02 stands for the alternative bs/pp/ge with γr = 0.8 and γs = 0.2. Figure 5 displays the bs/ep/ge and the non dominated alternatives.

bs/mp/ge 02 was faster but reached a slightly lower direct margin than bs/ep/ge. Both bs/mp/ge 08 and bs/pp/ge 08 02 reached a direct margin larger than bs/ep/ge with less runs. As alternatives are all very close, it is difficult to rank them. Figure 5 highlights the slight improvement due to partitioning methods, though bs/pp/ge 08 02 can be considered as the best. Indeed its maximum value was at least as good as the others, while the minimum and average values were larger.

Evaluation alternatives

The evaluation alternatives include the expected value of the sampled performance measures defined as score 1 . score 2 stands either for the global standard deviation V g (bs/ep/ge), or for the sampling standard deviation V p1 (bs/ep/p1e), or for the sampling standard deviation V p2 (bs/ep/p2e) evaluations. Figure 6 shows that bs/ep/p2e and bs/ep/p1e needed many more simulations than bs/ep/ge. Nevertheless, they both reached a larger averaged direct margin, which amount has never been reached with previous alternatives. The largest number of simulation runs (for bs/ep/p1e) was larger than 1.5 million of simulations. This large spread in simulation runs shows that bs/ep/p1e reached its second phase of a small benefit for an infinite number of simulation runs. Although the direct margin performance was the best, this large spread means that it is difficult to forecast the required time to reach the optimal sub-region.

Hybrid alternatives

There are many hybrid alternatives. We decided to test a few of them, selecting those which could be efficient trade-off between simulation runs number and direct margin achieved. First we hybridized the fast gs/ep/ge alternative with better partition techniques (gs/pp/ge 08 02 and gs/mp/ge 02). Note that it is useless to couple gs/ep/ge with other evaluation techniques since it does not take into account the score 2 value. Then we hybridized the alternatives with larger direct margin achievement with faster ones (bs/pp/p1e 08 02, bs/mp/p1e 02, bs/pp/p2e 08 02 and bs/mp/p2e 02). Eight hybrid alternatives were eventually tested. Figure 7 shows the hybrid alternatives that are non dominated with regard to simulation runs and 99.5% of the final direct margin achieved. Partitioning alternatives maximizing the difference of region indices have completely disappeared as they are always dominated by the division technique cutting the promising region upon the pivot.

Observation of irrigation strategies

In addition to algorithms results, we observed the input variable envelopes translating two computed irrigation alternatives. Figures 8 and9 represent multiple axis charts of the envelopes enclosing the 10 optimal regions reached for the 10 replications run of the gs/ep/ge and bs/pp/p1e alternatives. Couples of dashed lines bound the parameters of the optimal regions that were reached for every replications.

The shape envelopes of figures 8 and 9 are very similar. The main differences involve the T 3 parameter which still allows a wide range of possibilities. Input variables D3 (soil water deficit to stop irrigation) and I3 (amount of water applied during last irrigation round) have been barely broken down. On the other hand, the envelopes of input variables of D1 (soil water deficit to start a new irrigation cycle) and I2 (amount of water applied to a new irrigation round) are thin for every alternative.

The largest averaged direct margin reached over all replications is 547 ¤/ha. It results from the alternative bs/pp/p1e where the envelope denotes the following strategy. The first irrigation is started after observing a small water deficit, and a medium amount of water is applied. The new irrigations are performed as soon as a small water deficit is observed and a large amount of water is applied. The temperature, water deficit and water amount leading to the last irrigation are very variable.

Discussion

P2 alternatives efficiency

With regard to simulation runs and 99.5% of final direct margin achieved, the concβselection is completely dominated by the other alternatives. One can guess that the convex and evolutionary nature of the pending regions set makes the β -selection more efficient. Partitioning alternatives with dominated results have also been removed, so that only the alternatives including γr = 0.8 and γs = 0.2 technique's parameters have been kept. This combination of indices expresses relatively high level of risk acceptance (low sampled decision vector index) and that the pivot was chosen to build large 'good' regions and small 'bad' ones (high level of region index). These selected alternatives are the most efficient at simulating the crop model MODERATO. Concerning the evaluation alternatives, the use of sampled point standard deviations achieved larger direct margins, but needed many more simulation runs. The cutting process does not noticeably improve the P2 procedure, we nevertheless noticed that each of the non dominated hybrid alternatives included the pivot partitioning.

The alternatives shown in figure 10 compete to be the best irrigation strategy designer. The gs/ep/ge and gs/pp/ge alternatives were the most robust with respect to simulation runs. In addition, the gs/pp/ge alternative had a direct margin robustness which was at least as good as those of the others. 50 000 of simulation runs to reach 99.5% of the optimal region. Combinations including β -selection achieve higher direct margins, but it requires a far larger amount of simulation runs. Considering the frontier shape, we regard the gs/pp/ge alternative as an improvement over the gs/ep/ge alternative. Although the three last alternatives improve the direct margin, too many simulation runs are needed to do so. The gs/pp/ge alternative is now the default combination of techniques included in the P2 algorithm. The best combination of techniques is not directly expendable to other stochastic optimization problems. However, the conclusions drawn about the different techniques would allow the user to set it easily, especially when making decision about cultural operations.

Optimal irrigation plans

This section deals with the irrigation plans produced with the gs/pp/ge combination, compared with the irrigation plans produced with our reference combination bs/ep/ge.

Figure 12 shows the optimal decision space envelopes reached for two alternatives gs/pp/ge chosen as the best, and bs/ep/ge which was used as a test reference. The main noticeable difference still concerns the T 3 parameter, and although its range is still wide, its optimal value could be more difficult to find compared to the others. The recurrent thinness of input parameters D1 and D2 point out their sensitivity, while T 3 values variation are strongly related to the prior consumption of the limited amount of water. Nevertheless, the gs/pp/ge alternative improved both the direct margin reached and the number of simulation runs. The optimal decision space was at least no larger than our reference alternative. The agricultural case that was studied was the same as the one we used in [START_REF] Bergez | A hierarchical partitioning method for optimizing irrigation strategies[END_REF]. The direct margin average of 534 ¤/ha was reached within 1 600 000 simulation runs with the P2p algorithm. The direct margin average reached with the gs/pp/ge 08 02 of the P2 algorithm was 545 ¤/ha within 2 000 000 simulation runs, and the direct margin of 534 ¤/ha was achieved within less than 40500 simulation runs.

Conclusion

Although all the combinations are still available, we set the default P2 algorithm as the combination of greedy selection, pivot partition and average evaluation. We reached this conclusion by testing different combinations with an 8-parameter irrigation design problem. These combinations included common and new techniques for the performance of the three main steps of the hierarchical decomposition procedure: selection, division and evaluation. The aim of this paper was to investigate hierarchical decomposition of three main techniques in order to improve an earlier proposal rather than to compare it with other approaches. The optimal direct margin reached was about 10 ¤/ha higher than the one reached with the earlier proposal. About a quarter as many simulation runs were needed to reach the previously optimal direct margin. Hence the results obtained were satisfactory. We are now working on the extension of this approach for a better consideration of uncertainty, and of multiple objectives. In addition, more attention will be paid to the representation of the optimal decision space in order to facilitate its use by advisers.

Fig. 1

 1 Fig.1P2 process illustration: left hand side displays the tree construction related to the decision space exploration, and the right hand side shows iteration 1.

Fig. 2 Fig. 3

 23 Fig. 2 Some of the pending regions are incorporated in the convex non dominated front. The promising region is one of these. N.B. score 1 and score 2 have to be maximized.

Fig. 5

 5 Fig. 5 Evolution of the best partitioning alternatives (10 replications each).

Fig. 6

 6 Fig. 6 Evolution of evaluation alternatives (10 replications each).

Fig. 7

 7 Fig. 7 Evolution of non dominated hybrid alternatives (10 replications each).

Fig. 8 Fig. 9

 89 Fig. 8 Input variable envelope of the 10 optimal regions reached for the 10 replications of alternative gs/ep/ge.

Fig. 12

 12 Fig. 12Input variable envelope of the 10 optimal regions reached for the 10 replications of alternative gs/pp/ge 08 02 (solid line) and Input variable envelope of the 10 optimal regions reached for the 10 replications of alternative bs/ep/ge (dashed line).

Table 2

 2 General description of the strategy simulated.

	Operation	Rules
	Sowing	Sowing is between 20 April and 30 May as soon as the cumulative rainfall during the
		previous 3 days is less than 15 mm. Variety Cécilia is sown at 80 000 plants/ha. Cécilia
		is a late growing variety requiring 1045 accumulated thermal units (ATU) from sowing to
		flowering and 1990 ATU from sowing to maturity (35% grain humidity).
	Fertilization	A single application of 200 kg/ha of nitrogen is made at sowing.
	Harvest	The crop is harvested when grain moisture content reaches 20% or accumulated thermal
		units since sowing reach 2100 ATU and if the cumulative rainfall during the previous 3
		days is less than 15 mm. In any case, the crop must be harvested before 15 October.
	Irrigation
		Sowing Irrigation to facilitate plant emergence (caused either by dryness or crust created
		by heavy rainfall on silty soil) is not taken into account, nor irrigation to dissolve
		fertilizer.
		Starting irrigation Part of the optimization process.
		Next irrigation round Part of the optimization process
		Delay irrigation Precipitation delays irrigation. When the cumulative rainfall over the 5
		previous days is more than 10 mm, one day delay is applied for every 4 mm. The
		delay cannot exceed 7 consecutive days.
		Stopping irrigation Part of the optimization process

Table 4

 4 Triple code notation: definition of X, Y and Z.

	selection β	X bs	division equal partition	Y ep	evaluation average and global standard deviation	Z ge
	Concβ	cs	pivot partition	pp	average and point standard deviation 1	p1e
	Greedy	gs	maximize difference	mp	average and point standard deviation 2	p2e

  Evolution of selection alternatives (10 replications each).

		550
			gs/ep/ge	bs/ep/ge				cs/ep/ge
		545
		540
	direct margin (euros/ha)	525 530 535
		520
							min, average, max of 10 replications
		515							95%
									97.5%
									99%
									99.5%
		510
		0	50000		100000		150000	200000	250000
					nb of simulation runs
	515 520 525 530 535 540 545 550 Fig. 4 510 direct margin (euros/ha)	bs/mp/ge_02 bs/mp/ge_08		bs/ep/ge bs/pp/ge_08_02	min, average, max of 10 replications 95% 97.5% 99% 99.5%
		30000	40000	50000	60000	70000	80000	90000	100000
					nb of simulation runs

  Figure11shows the frontier made of the 99.5% thresholds of the non dominated combinations. Combinations including the greedy selection are the fastest and need less than

		550
		548
									bs/pp/p1e_08_02
		546		bs/pp/ge_08_02	bs/pp/p2e_08_02
	direct margin (euros/ha)	540 542 544	gs/ep/ge gs/pp/ge_08_02
		538
		536					min, average, max of 10 replications
										99.5%
			40000	60000	80000	100000	120000	140000	160000	180000	200000
						nb of simulation runs
	Fig. 11 All non dominated alternatives for 99.5% of the optimal region.
		550	gs/pp/ge_08_02 bs/pp/ge_08_02	bs/pp/p2e_08_02		bs/pp/p1e_08_02
		545
		540	gs/ep/ge
	direct margin (euros/ha)	525 530 535
		520
							min, average, max of 10 replications
		515								95%
										97.5%
										99%
										99.5%
		510
			40000	60000	80000	100000	120000	140000	160000	180000	200000
						nb of simulation runs
	Fig. 10 All non dominated alternatives.

See for example[START_REF] Ehrgott | Multicriteria Optimization[END_REF] for dominance definitions.


P Cathalifaud

J Mauss
email: mauss@cict.fr

J Cousteix
email: jean.cousteix@onecert.fr

Nonlinear aspects of high Reynolds number channel flows

Keywords:

This paper considers the flow in a two-dimensional channel at high Reynolds number, with wall deformations which can lead to flow separation. An asymptotic model is proposed by using the successive complementary expansion method with generalized asymptotic expansions. In particular, the model emphasizes the asymmetry of the channel geometry by introducing a change of variables. It is shown that the model is more general than the models developed with the method of matched asymptotic expansions. Comparisons with Navier-Stokes solutions show that the model is well founded and enables us to treat original problems.

Introduction

We consider a plane Poiseuille flow in a two-dimensional channel perturbed by wall deformations -troughs or bumpswhich are sufficiently severe to induce flow separation.

At high Reynolds number, an asymptotic analysis can be applied. When the Reynolds number tends to infinity, the first step consists of a straightforward simplification of the Navier-Stokes equations leading to the Euler equations. So, it is tempting to call for a hierarchical process. In the vicinity of singular zones, near the walls, the second step consists of trying to correct the first approximation by a boundary layer analysis.

However, it is known that, in many problems involving a strong coupling, this type of hierarchical approach is not possible and the two approximations must be considered simultaneously. In fact, the solution of the inviscid flow equations depends on the Reynolds number and it is better to talk about Euler equations rather than inviscid flow equations as is done very often.

If we decide that the numerical solution of Navier-Stokes equations is not desirable, because either the numerical cost to find such a solution is too big, or a better understanding of physical phenomena is desired, a suitable approach consists of using the asymptotic model called the interactive boundary layer (IBL) model.

Then, we are faced with a paradox. The only constructive way to establish asymptotic approximations is to use a small parameter, here the inverse of the Reynolds number, which goes to zero. This leads to regular asymptotic expansions and to a hierarchical approach.

Excluding a multi-layer approach of triple deck type, which introduces very restrictive hypotheses on the scales, a possibility is to use generalized asymptotic expansions. According to this method, the small parameters of the problem can be included in the functions which form the expansions. This conception is very different because the small parameters are not considered as tending towards zero but are only small. This is easy to conceive if a uniformly valid approximation (UVA) is constructed from the standard method of matched asymptotic expansions.

In fact, the reasoning can be inverted by starting from an assumed form of the UVA. In order to extend the techniques on which the method of matched asymptotic expansions is based, the successive complementary expansion method (SCEM) has been developed [START_REF] Cousteix | Asymptotic Analysis and Boundary Layers[END_REF]. When regular approximations are sufficient, it can be shown that the regular SCEM is equivalent to the method of matched asymptotic expansions associated to a modified Van Dyke matching principle.

Thanks to the generalized expansions, interaction effects between the Eulerian and boundary layer regions are considered simultaneously and the construction of a UVA does not require any matching principle; only the boundary conditions of the problem are used.

In a channel, there is no external flow region and the asymptotic models for the flow perturbations are mainly based on an inviscid rotational core flow region together with boundary layers near the walls; a comprehensive discussion of this structure can be found in Sobey [START_REF] Sobey | Introduction to Interactive Boundary Layer Theory[END_REF]. The asymptotic analysis of these flows has been performed essentially by Smith [START_REF] Smith | On the high Reynolds number theory of laminar flows[END_REF], and a systematic approach has been proposed later on by Cousteix and Mauss [START_REF] Cousteix | Interactive boundary layer models for channel flow[END_REF]. More recently, the modelling of channel flows has been examined by Lagrée et al. [START_REF] Lagrée | Characterization of the pressure drop in a 2D symmetrical pipe: Some asymptotical, numerical and experimental comparisons[END_REF] and by Lagrée and Lorthois [START_REF] Lagrée | The RNS/Prandtl equations and their link with other asymptotic descriptions: Application to the wall shear stress in a constricted pipe[END_REF].

After a presentation of the SCEM (Section 2), the problem is formulated (Section 3), then the flow is analysed with SCEM using standard hypotheses on the wall deformations (Section 4). In Section 5, a new approach is proposed in order to emphasize the role of the asymmetry of the channel geometry. As shown in Section 6, the resulting model contains the standard model obtained by Smith [START_REF] Smith | Flow through constricted or dilated pipes and channels: Part 1[END_REF][START_REF] Smith | Flow through constricted or dilated pipes and channels: Part 2[END_REF][START_REF] Smith | On the high Reynolds number theory of laminar flows[END_REF] but is more general. Finally, the numerical technique is presented in Section 7 and results are given with comparisons with Navier-Stokes solutions (Section 8).

About the SCEM

The purpose of this section is to give a short overview of the SCEM, which is explained in more detail in Cousteix and Mauss [START_REF] Cousteix | Asymptotic Analysis and Boundary Layers[END_REF].

Definition of a generalized asymptotic approximation.

First of all we have to define clearly what is an asymptotic approximation of a function. Two functions Φ(x, ε) and Φ a (x, ε) defined in a domain D are asymptotically identical to order δ(ε) if their difference is asymptotically smaller than δ(ε), where δ is an order function:

Φ(x, ε) -Φ a (x, ε) = o(δ), (1)
where ε is a small parameter arising from the physical problem studied. The function Φ a is called an asymptotic approximation of the function Φ. The converse being also true, it must be stressed that the goal of an asymptotic approximation is to replace Φ by a simpler function Φ a (x, ε). This simplicity comes from the different methods used to construct asymptotic expansions (AE), generally defined by

Φ a (x, ε) = n i=0 δ i (ε)ϕ i (x, ε), (2)
where the asymptotic sequence of order functions δ i is such that δ n δ and such that all functions ϕ i are of strict order 1, ϕ i = O S (1). This is called a generalized asymptotic expansion because ε is included in the developed functions.

Definition of a regular asymptotic approximation.

In many cases, it is possible to construct a regular AE, with an asymptotic sequence noted δ (0)

i , which can be written

Φ a (x, ε) = E 0 Φ = n i=0 δ (0) i (ε)ϕ (0) i (x). (3)
For this reason, when the AE is not regular, it is called a generalized asymptotic expansion. Strictly, it is not necessary to add the qualifying adjective ''generalized'' but, too often, it is thought that an AE is necessarily regular. Such an asymptotic approximation is objectively simpler than the function Φ since the construction is made step by step by limiting processes. This is the reason why, following Eckhaus [START_REF] Eckhaus | Matched Asymptotic Expansions and Singular Perturbations[END_REF], we define, for a regular AE, a special operator E 0 called the ''expansion operator'' at a given order δ; thus, Φ -

E 0 Φ = o(δ).

Method of matched asymptotic expansions (MMAE).

Interesting cases occur when the function Φ is not regular in D, i.e., when AE (2) or (3) of Φ is valid only in a restricted region D 0 ∈ D called the outer region. We are faced with a singular perturbation problem and we must introduce boundary layer domains. Here, in the simplest case, we introduce an inner domain which can be formally denoted D 1 = D -D 0 and located near the origin. The boundary layer variable or inner variable is X = x ζ (ε) , ζ (ε) being the order of thickness of this boundary layer. If a regular expansion can be constructed in D 1 , we can write

Φ a (x, ε) = E 1 Φ = m i=0
δ (1) i (ε)ϕ (1) i (X).

(

) 4
This inner ''expansion operator'' E 1 is defined in D 1 at the same order δ as the outer expansion operator E 0 ; thus, Φ -

E 1 Φ = o(δ).
The links between these expansions E 0 Φ and E 1 Φ are the rules or heuristic principles which define the asymptotic matching.

To our knowledge, the best principle is the modified Van Dyke principle (MVDP), which can be found in [START_REF] Cousteix | Asymptotic Analysis and Boundary Layers[END_REF] and can be written

E 0 E 1 Φ = E 1 E 0 Φ.
As the expansions are defined at the same order δ, we can construct a composite approximation at this order, or an asymptotic approximation, which can be written

Φ a = E 0 Φ + E 1 Φ -E 1 E 0 Φ. (5)
This is clearly a uniformly valid approximation (UVA) defined at order δ.

The successive complementary expansion method (SCEM).

For singular perturbation problems, the successive complementary expansion method (SCEM) is based on the knowledge of generalized asymptotic expansions and upon the idea that the reasoning using the method of matched asymptotic expansions (MMAE) must be inverted. First, the structure of the uniformly valid approximation (UVA) must be assumed and then the method to construct the UVA is deduced. This point of view has been adopted in other methods like WKB or multiple scale methods, but their usefulness is restricted to particular problems.

We assume that two significant subdomains have been identified, classically called the outer and inner domains. According to the SCEM, a UVA Φ a of the function Φ is supposed to have the following structure:

Φ a (x, X , ε) = n i=1 δi (ε) φi (x, ε) + Ψ i (X, ε) . (6)
We say that Φ a is a UVA of Φ defined at order δn if Φ = Φ a + o( δn ). Again, this UVA is called generalized because classical understanding used rather ''regular'' asymptotic expansion where ε is not in the developed functions. For instance, we can write

Φ a = Φ ar + o(δ m ),
where Φ ar is a regular UVA such that δn = O(δ m ), and

Φ ar (x, X , ε) = n i=1 δ i (ε) [ϕ i (x) + Ψ i (X)] . (7)
The sequence of order functions δi may or may not be the same as the sequence δ i since the elements of this latter sequence are gauge functions. Recall that a gauge function is an a priori choice of an order function in its equivalence class.

The SCEM has already been used in its regular form [START_REF] O'malley | Introduction to singular perturbations[END_REF][START_REF] Vasil'eva | The boundary function method for singular perturbation problems[END_REF], but not in its generalized form. Under the name of composite expansions, Nayfeh [START_REF] Nayfeh | Perturbations Methods[END_REF] discusses the related technique developed earlier.

In its regular form, the SCEM is equivalent to the MMAE with the modified form of the Van Dyke matching principle (MVDP). This is a crucial point since the matching is a consequence of the chosen form of the UVA and not the contrary.

The ambiguous writing of the generalized SCEM needs guidelines to enable construction of a UVA. In general, knowledge of the physics of the problem, as will be seen later, provides us with an efficient solution.

One can be guided by the fact that if the solution Φ of integrodifferential equations is formally written L ε Φ = 0, then, for the UVA given by ( 6) we have L ε Φ a = R n (x, ε). Since we have assumed that the boundary conditions are exactly satisfied for Φ a , R n is small in a certain sense. This is ascertained by estimate theorems which can be obtained, in general, only for linear problems. The generalized form proposed here is very well adapted to strong interaction problems leading logically to the models of interactive boundary layers.

Formulation of the problem

The dimensionless Navier-Stokes equations are

∂U ∂x + ∂V ∂y = 0, (8)
L ε U = U ∂U ∂x + V ∂U ∂y + ∂P ∂x - 1 R ∂ 2 U ∂x 2 + ∂ 2 U ∂y 2 = 0, (9)
L ε V = U ∂V ∂x + V ∂V ∂y + ∂P ∂y - 1 R ∂ 2 V ∂x 2 + ∂ 2 V ∂y 2 = 0. ( 10
)
The above equations are written in an orthonormal axis system. Coordinates x and y are reduced by the width W * of the basic channel whose lower wall is y = -1 2 and upper wall is y = 1 2 (Fig. 1). The velocity components are reduced by a reference velocity V * defined below by [START_REF] Nayfeh | Perturbations Methods[END_REF]. The Reynolds number is

R = ̺ * V * W * µ * , (11)
where ̺ * and µ * are the dimensionalized density and viscosity coefficient.

The reference velocity V * is expressed as a function of the basic pressure gradient dp 0 dx or, equivalently, of the dimensionalized mass flow Q * per unit width of the channel:

Q * = ̺ * V * W * 1/2 -1/2 u 0 dy, V * = 6 Q * ̺ * W * , dp 0 dx = - 2µ * ̺ * V * W * , (12)
where the basic plane Poiseuille flow is

u 0 = 1 4 -y 2 , v 0 = 0, p 0 = - 2x R + p c , (13)
and p c is an arbitrary constant. The flow is perturbed by indentations of the lower and upper walls:

y l = - 1 2 + F (x, ε), y u = 1 2 -G(x, ε), ( 14
)
where ε is a small parameter (Fig. 1).

At high Reynolds number, the reduced equations obtained formally by taking their limit when the Reynolds number goes to infinity are of first order. A singular perturbation problem arises.

Direct analysis

To go further, it is usual to consider small wall perturbations leading to hypothesis (H1): From the application of the standard method of matched asymptotic expansions, two zones are considered in the flow. One of them is the core where the flow perturbations are governed by Euler equations. The other zone is the proximity of walls where the flow perturbations develop in a structure of boundary layer type. A perturbation is said to be significant when flow separation is possible. To translate this, it is required that, in the boundary layers of thickness ε, the perturbation of the longitudinal velocity is of the same order as u 0 , i.e., of order O(ε). Thus, according to the SCEM, we are seeking a UVA of the form

(H1) F = εf , G = εg. (15)
U = u 0 (y) + ε û(x, y, ε) + • • • , (16)
V = ε v(x, y, ε) + • • • , (17)
P -p c = - 2x R + ε p(x, y, ε) + • • • . (18)
For large Reynolds numbers, the core flow equations are [START_REF] Cousteix | Asymptotic Analysis and Boundary Layers[END_REF] ∂ û1 ∂x

+ ∂ v1 ∂y = 0, (19)
u 0 ∂ û1 ∂x + v1 du 0 dy + ε û1 ∂ û1 ∂x + v1 ∂ û1 ∂y = - ∂ p1 ∂x , (20)
u 0 ∂ v1 ∂x + ε û1 ∂ v1 ∂x + v1 ∂ v1 ∂y = - ∂ p1 ∂y , (21)
where the index ''1'' denotes the characteristics of the flow perturbation in the core. This means that (19)-( 21) are valid only in the core.

In the whole flow field, the Navier-Stokes equations reduce to [START_REF] Cousteix | Asymptotic Analysis and Boundary Layers[END_REF]

∂ û ∂x + ∂ v ∂y = 0, ( 22
)
u 0 ∂ û ∂x + v du 0 dy + ε û ∂ û ∂x + v ∂ û ∂y = - ∂ p1 ∂x + 1 R ∂ 2 û ∂y 2 . ( 23
)
As shown in [START_REF] Cousteix | Asymptotic Analysis and Boundary Layers[END_REF][START_REF] Cousteix | Interactive boundary layer models for channel flow[END_REF], it must be noted that in the streamwise momentum equation,

∂ p ∂x is replaced by ∂ p1
∂x . This means that the pressure required to solve ( 22)-( 23) satisfies the core flow equations ( 19)-( 21) but ( 22)-( 23) are valid everywhere in the channel and they provide us with a UVA. Since this is the key point of the method, let us recall that the variables of the boundary layer are given by

Y = 1 2 + y ε , Y = 1 2 -y ε . (24)
It can be shown that a uniformly valid approximation of the pressure is given by p(x, y, ε) = p1 (x, y, ε)

+ ε 3 P CL (x, Ŷ , ε) + ε 3 P CL (x, Y , ε), (25)
where P CL and P CL are zero in the core flow. Eqs. ( 22)-( 23) look like the standard boundary layer equations, but they are valid in the whole domain. To calculate the flow, it is required to solve the system comprising the generalized boundary layer equations ( 22)-( 23) and the core flow equations ( 19)-(21). It is also necessary to prescribe the wall conditions y = y ℓ and y = y u : u 0 + ε û = 0, v = 0, (26) and to express the coupling between the core flow equations and the generalized boundary layer equations by identifying û, v and û1 , v1 in the core:

(û, v) → (û 1 , v1 ). ( 27
)
The global interactive boundary layer model described by ( 19)-( 21), ( 22)-( 23), ( 26), ( 27) is the best approximation of the Navier-Stokes model we can propose.

In [START_REF] Cousteix | Asymptotic Analysis and Boundary Layers[END_REF][START_REF] Mauss | Global interactive boundary layer (GIBL) for a channel[END_REF], it is seen that the boundary layer analysis leads to

R = ε -3 . ( 28
)
Then, using the Navier-Stokes operator defined in ( 9)-( 10), we observe that

L ε û = ∂p ∂x - ∂ p1 ∂x - 1 R ∂ 2 û ∂x 2 , L ε v = - 1 R ∂ 2 v ∂x 2 .
From ( 25) and (28), it is easily shown that

L ε û = O(ε 3 ), L ε v = O(ε 3 ), (29)
which shows, as stated before, that this is a very good approximation.

The nonlinear core equations are not easy to solve. Fortunately, as shown below from a deeper analysis, the core flow equations can be linearized, and the solution of the resulting model is much easier.

In the boundary layer whose thickness is of order ε, u 0 = O(ε), û = O(1) and v = O(ε), so the three terms

u 0 ∂ û ∂x + v du 0 dy , ε û ∂ û ∂x + v ∂ û ∂y , ε 3 ∂ 2 û ∂y 2 ,
are of the same order, O(ε).

In the core of the flow, the orders of magnitude are not the same.

The terms û1

, ∂ û1 ∂x , ∂ û1 ∂y , v1 , ∂ v1 ∂x , ∂ v1
∂y are all O(1). Therefore, the field equations are structurally nonlinear, whereas the core flow Euler equations can be linearized rightfully.

If we define

u = u 0 + ε û, (30)
v = ε v, (31)
and

u 1 = u 0 + ε û1 , (32)
v 1 = ε v1 , ( 33
)
p 1 = p 0 + ε p1 , (34)
the field and core flow equations become

∂u ∂x + ∂v ∂y = 0, ( 35
) u ∂u ∂x + v ∂u ∂y = - ∂p 1 ∂x + 1 R ∂ 2 u ∂y 2 , ( 36
)
∂u 1 ∂x + ∂v 1 ∂y = 0, ( 37
)
u 0 ∂u 1 ∂x + v 1 du 0 dy = - ∂ ∂x (p 1 -p 0 ) , (38)
u 0 ∂v 1 ∂x = - ∂ ∂y (p 1 -p 0 ) . (39)
It is interesting to note the singular behaviour of the solution of (37)-(39) in the neighbourhood of the boundaries. For instance, as y → -1/2, we have

u 1 -u 0 = -2p 10 ln 1 2 + y + c 10 + • • • , (40)
where p 10 and c 10 are functions of x and ε.

The logarithmic singularity of this solution near y = -1/2 and y = 1/2 exhibits an essential difference compared with external boundary layers. In the following, for this reason, the evolution of the perturbed flow in the core will be characterized by v 1 . From (37)-(39), it is found that v 1 is a solution of Poisson's equation:

u 0 ∂ 2 v 1 ∂x 2 + ∂ 2 v 1 ∂y 2 = v 1 d 2 u 0 dy 2 , ( 41
)
and the x-component of the pressure gradient required to solve the generalized boundary layer equations is given by (38), in which the continuity equation ( 37) is taken into account:

-u 0 ∂v 1 ∂y + v 1 du 0 dy = - ∂ ∂x (p 1 -p 0 ) . (42)
It can be shown that (41) associated to (42) gives the ymomentum equation (39) if the perturbations vanish at upstream infinity. This establishes the equivalence between (37)-( 39) and ( 41)-(42).

To sum up, the problem to solve comprises (35)-( 36), ( 41) and (42). At the walls, the boundary conditions are

y = y ℓ and y = y u : u = 0, v = 0, (43)
and the coupling between the core flow equations and the generalized boundary layer equations is expressed by identifying u, v and u 1 , v 1 in the core:

(u, v) → (u 1 , v 1 ). ( 44
)
The model presented above belongs to the class of strong coupling models since there is no hierarchy between the boundary layer equations and the core flow equations. The triple deck theory, or more precisely its equivalent for channel flows as developed by Smith [START_REF] Smith | Flow through constricted or dilated pipes and channels: Part 1[END_REF][START_REF] Smith | Flow through constricted or dilated pipes and channels: Part 2[END_REF], also belongs to this class of strong coupling models. In fact, Smith's model is included in the present model since the expansions are regular, whereas in the present model the expansions are generalized. It is interesting to note that the first approximation of Smith's model for v 1 is symmetric with respect to y and corresponds to a geometrically antisymmetric problem. This point will be discussed later when a regular analysis is developed (Section 6). In the core, Smith's model gives

v 1 = -u 0 (y) dA(x) dx , ( 45
)
where A is defined as the displacement function.

It must be noted that ( 45) is an eigensolution of (42) but not of (41). In addition, the approximation (45) does not have the singular behaviour of (40). These two remarks lead us to try to separate as far as possible the symmetric and antisymmetric problems, which leads, as we will see, to a new approach of the asymptotic problem.

Influence of asymmetry

The analysis starts from Navier-Stokes equations ( 8)- [START_REF] O'malley | Introduction to singular perturbations[END_REF] in which we introduce the transformation

X = x, (46) Y = y -H(x, ε), (47) U = U, (48) V = V -U dH dx , (49)
P = P , ( 50
)
where ε is formally the thickness of the boundary layer, i.e., the region where the inertial terms are of the same order as the viscous terms.

This point is also important in transferring the perturbations from the walls into the equations. The asymptotic analysis is also easier to do.

In order to define H, we set

E = F + G 2 , (51)
H = F -G 2 , (52)
so that the problem is geometrically symmetrized in the transformed plane. Obviously, the solution is not symmetrical. Note also that the physical channel is deformed symmetrically when H = 0.

In general, the wall conditions are

Y = Y ℓ = - 1 2 + E and Y = Y u = 1 2 -E : U = 0; V = 0. ( 53
)
We introduce hypothesis (H2):

(H2):

E = O(ε), νµ 2 ε, (54)
with

H = νh(µx),
where ν(ε) and µ(ε) are order functions.

This extension of hypothesis (H1) shows that actually hypothesis (H1) applies only on the symmetric part of the perturbation.

With these hypotheses, the Navier-Stokes equations become

∂U ∂X + ∂V ∂Y = 0, ( 55
) U ∂U ∂X + V ∂U ∂Y = - ∂P ∂X + 1 R ∂ 2 U ∂X 2 + ∂ 2 U ∂Y 2 + O(νµ), (56)
U ∂V ∂X + V ∂V ∂Y = - ∂P ∂Y + 1 R ∂ 2 V ∂X 2 + ∂ 2 V ∂Y 2 + O(νµ). (57)
For small ε, the basic flow corresponding to E = 0 is

U 0 = 1 4 -Y 2 , ( 58
)
V 0 = 0, ( 59
)
P 0 = -2 X R + P c , (60)
where P c is an arbitrary constant.

The Navier-Stokes equations can now be written as

∂U ∂X + ∂V ∂Y = 0, ( 61
) U ∂U ∂X + V ∂U ∂Y = - ∂P ∂X + dH dX ∂P ∂Y + 1 R ∂ 2 U ∂Y 2 + O 1 R , (62)
U ∂V ∂X + V ∂V ∂Y + d 2 H dX 2 U 2 = -1 + dH dX 2 ∂P ∂Y + dH dX ∂P ∂X + 1 R ∂ 2 V ∂Y 2 + O 1 R . ( 63
)
For the core flow whose characteristics are denoted with index ''1'', the Euler equations become

∂U 1 ∂X + ∂V 1 ∂Y = 0, U 1 ∂U 1 ∂X + V 1 ∂U 1 ∂Y = - ∂ ∂X (P 1 -P 0 ) + dH dX ∂ ∂Y (P 1 -P 0 ) , U 1 ∂V 1 ∂X + V 1 ∂V 1 ∂Y + d 2 H dX 2 U 2 1 = -1 + dH dX 2 ∂ ∂Y (P 1 -P 0 ) + dH dX ∂ ∂X (P 1 -P 0 ) . (64)
As in (32)-( 34), the flow characteristics can be expanded about the basic flow characteristics (58)-(60). Then, the following orders of magnitude are found:

U 1 = O S (1), ∂U 1 ∂X = O(ε), ∂U 1 ∂Y = O S (1), V 1 = O(ε), ∂V 1 ∂X = O(ε), ∂V 1 ∂Y = O(ε). ( 65
)
It is recalled that the symbol O S means ''is of strict order of'' whereas the symbol O means ''is at most of order of'' [START_REF] Cousteix | Asymptotic Analysis and Boundary Layers[END_REF].

With hypothesis (H2), the left-hand sides of the second and third equations of (64) are of order ε. It is deduced that

∂ ∂X (P 1 -P 0 ) = O(ε), ∂ ∂Y (P 1 -P 0 ) = O(ε).
As in Section 4, this leads us to the linearized core flow equations:

∂U 1 ∂X + ∂V 1 ∂Y = 0, (66)
U 0 ∂U 1 ∂X + V 1 dU 0 dY = - ∂ ∂X (P 1 -P 0 ) , (67)
U 0 ∂V 1 ∂X + d 2 H dX 2 U 2 0 = - ∂ ∂Y (P 1 -P 0 ) . (68)
Moreover, for large Reynolds numbers, as in [START_REF] Cousteix | Asymptotic Analysis and Boundary Layers[END_REF], it can be shown that, in the boundary layer, ∂P ∂X can be replaced by

∂P 1
∂X in the streamwise momentum equation.

Finally, the complete system to solve comprises the core flow equations (66)-(68) and the field equations ( 61)-(62), simplified as

∂U ∂X + ∂V ∂Y = 0, (69)
U ∂U ∂X + V ∂U ∂Y = - ∂P 1 ∂X + 1 R ∂ 2 U ∂Y 2 . ( 70
)
It must be noted that Eqs. ( 69)-(70) are valid not only in the boundary layers but everywhere in the channel.

The conditions required to solve this system are the boundary conditions (53) and the coupling condition in the core:

V → V 1 . (71) Since d 2 H
dX 2 is an approximation of the curvature along the median line Y = H(X ), the resulting problem obtained using the SCEM is analogous to the problem sets with hypothesis (H1), with an additional term in Eq. (68).

Regular analysis of the transformed problem

The idea is to perform a regular analysis of the model based on (69)-( 70), (66)-( 68), (71) and to show that this model contains in particular the first-order results obtained by Smith. A slow variable X is introduced:

X = X ε α . (72) By defining V such that V = ε α V , (73) the continuity equation is ∂U ∂X + ∂V ∂Y = 0. ( 74
)
It is assumed that α ≥ 0 so that the length of the wall deformations is much larger than the channel width or is of the same strict order. The streamwise momentum equation (70) becomes

U ∂U ∂X + V ∂U ∂Y = - ∂P 1 ∂X + 1 Rε α ∂ 2 U ∂Y 2 . ( 75
)
It is easy to show that, in the vicinity of walls, the regular analysis of (75) leads to boundary layer equations which have the same form as (75). Now, let us examine the core flow equations. By defining V 1 such that

V 1 ε α = V 1 , (76)
the core flow equations (66)-(68) become

-U 0 ∂V 1 ∂Y + V 1 dU 0 dY = - ∂ ∂X (P 1 -P 0 ) , (77)
ε 2α U 0 ∂V 1 ∂X + d 2 H dX 2 U 2 0 = - ∂ ∂Y (P 1 -P 0 ) . (78)
As we know the behaviour of the pressure, we look for regular expansions beginning with

V 1 (X, Y , ε) = ε r V 10 (X, Y ) + • • • , (79) P 1 (X, Y , ε) -P 0 = ε 2 P 10 (X, Y ) + • • • . (80)
With the above expansions, V 1 = O(ε) implies that r ≥ 1. In addition, the regular analysis imposes µ = ε α and H = εh(X).

The analysis of (77)-( 78) with the asymptotic expansions (79)-(80) leads us to consider two cases.

1. α > 1 2 . To first order, the pressure is constant in a cross-section

∂P 10 ∂Y = 0, (81)
and, for a direct comparison with Smith's results, the solution for V 1 can be cast in the form

V 1 = -εU 0 dA dX + dh dX , (82)
where A(X ) is a function of X known as the displacement function.

As the global problem to solve is fully symmetric with respect to Y , V 10 must be antisymmetric, but (82) gives a symmetric solution if A + h = 0. Thus, we have to take A + h = 0. To the considered order, the perturbation is negligible.

2. α ≤ 1 2 . The solution for V 1 is V 1 = -ε dh dX U 0 (Y ) -ε 2-2α dA dX U 0 (Y ) + • • • , (83)
and, for α < 1/2, the pressure variation in a cross-section is given by

∂P 10 ∂Y = d 2 A dX 2 U 2 0 (Y ). (84)
Here, the displacement function A(X ) must be determined from the matching condition between the core and the boundary layer.

For α = 1/2, we have

V 1 = -ε dh dX + dA dX U 0 + • • • , (85) ∂P 10 ∂Y = U 2 0 (Y ) d 2 A dX 2 .
(86)

Finally, the solution for

U 1 , V 1 is 1. α > 1 2 U 1 = U 0 (Y ) + O(ε 2 ), ( 87
) V 1 = O(ε 2 ). ( 88
)
This case is interesting because it gives the limits of a problem where the walls are not exactly symmetric but which can be treated as being symmetric.

2. α ≤ 1

2

U 1 = U 0 (Y ) + εh(X) dU 0 dY + ε 2-2α A(X ) dU 0 dY + • • • , ( 89
) V 1 = -ε 1+α dh dX U 0 (Y ) -ε 2-α dA dX U 0 (Y ) + • • • . ( 90
)
Let us recall that the regular form of the SCEM leads to a uniformly valid approximation to a given order, here ε. In addition, the results are equivalent to those obtained with the standard method of matched asymptotic expansions when the matching is performed to the same order according to the modified Van Dyke principle [START_REF] Cousteix | Asymptotic Analysis and Boundary Layers[END_REF]. For U, we have formally

U = U 1 + εU BL -εU c , (91)
where U is a UVA to order ε, U 1 is the approximation in the core to order ε, and εU BL is the approximation in the boundary layer to order ε. The quantity U c results from the matching condition to order ε. For example, in the vicinity of Y = 1/2, the boundary layer variable is

Y = 1 2 -Y ε . (92)
With

U 0 = 1 4 -Y 2 ,
we obtain the following matching conditions:

1. α < 1 2 : lim Y →∞ (U BL -Y ) = -h(X), 2. α = 1 2 : lim Y →∞ (U BL -Y ) = -h(X) -A(X ), 3. α > 1 2 : lim Y →∞ (U BL -Y ) = 0.
These results can be interpreted in terms of the original variables.

Taking into account the transformation (47), and defining the boundary layer variable ỹ = 1 2 -y /ε, we have Y = ỹ + h(X ), (93) so the matching conditions become

1. α < 1 2 : lim ỹ→∞ (u BL -ỹ) = 0, 2. α = 1 2 : lim ỹ→∞ (u BL -ỹ) = -A(X), 3. α > 1 2 : lim ỹ→∞ (u BL -ỹ) = h(X ).
If all the results obtained in this section with the transformed variables are interpreted in terms of the original variables, they are identical to those obtained in [START_REF] Cousteix | Asymptotic Analysis and Boundary Layers[END_REF][START_REF] Smith | Flow through constricted or dilated pipes and channels: Part 1[END_REF][START_REF] Smith | Flow through constricted or dilated pipes and channels: Part 2[END_REF]. This means that the model based on (69)-( 70), ( 66)-( 68), (71) contains in particular the model developed by Smith to first order, but is more general.

According to the analysis proposed in Section 5, the problem is to solve the generalized boundary layer equations ( 69)-( 70) associated to the core flow equations (66)-(68). The velocity components U and V vanish at the walls, and the coupling condition between the boundary layer solution and the core flow solution is given by (71).

First, the core flow equations are rewritten by introducing a change of variable:

V * 1 = V 1 + dH dX U 0 . (94)
Then, it is shown that V * 1 obeys Poisson's equation:

U 0 ∂ 2 V * 1 ∂X 2 + U 0 ∂ 2 V * 1 ∂Y 2 = V * 1 d 2 U 0 dY 2 . (95)
The coupling condition becomes

V = V * 1 - dH dX U 0 . (96)
To solve the equations, a change of coordinates is introduced,

ξ = X , (97) η
= Y ζ (X) with ζ (X) = 1 -2E(X ), (98)
so that the domain is limited by the lines η = -1/2 and η = 1/2 which define the channel walls.

In agreement with the analysis of Section 5, the core flow equations are linearized in the (ξ , η)-plane, so that Poisson's equation for V * 1 keeps the same form as (95). In this equation, the velocity of the basic flow is

U * 0 = 1 4 -η 2 . ( 99
)
The generalized boundary layer equations are rewritten by introducing the following velocity transformation:

U * = U, (100) V * = V -Uη dζ dX , (101)
so the generalized boundary layer equations are only slightly modified.

The numerical procedure used to solve the whole set of equations is discussed below.

Poisson's equation for V * 1 is discretized using a five-point stencil. A successive line overrelaxation method is applied to solve the resulting equations. According to this technique, an iterative procedure is implemented in which the calculation domain is swept from upstream to downstream. When the solution is calculated at a given ξ -station, the value of V * 1 at the upstream ξ -station is the updated value of V * 1 calculated at the current iteration, whereas the value of V * 1 at the downstream station is taken at the previous iteration. The sweeping of the calculation domain is repeated until convergence is achieved.

At a given ξ -station, the longitudinal component of the pressure gradient is calculated from the core flow streamwise momentum equation and is used to solve the generalized boundary layer equations. These equations are discretized using a standard finite difference method with three points in the η-direction and two points in the ξ -direction.

In fact, Poisson's equation and the generalized boundary layer equations are strongly coupled, and the whole set of equations is solved as simultaneously as possible. The procedure is described below.

The calculation of the values of U * , V * and V * 1 at the different grid points for a given ξ -station consists of solving a system of algebraic equations in which the coefficients of the unknowns form tridiagonal matrices. The main problem is that the values of V * 1 are not known a priori along the lines η = -1/2 and η = 1/2. At a given ξ -station, the solution of Poisson's equation can be written as

V * 1i = A i V * 1 - 1 2 + B i V * 1 1 2 + C i , (102)
where V * 1i is the value of V * 1 at a certain grid point i along a line ξ = const. The coefficients A i , B i and C i do not depend on the values

V * 1 (-1/2) and V * 1 (1/
2), and are calculated explicitly.

In the same way, after replacing the streamwise component of the pressure gradient by its expression as a function of V * 1 and of its transverse derivative, the solution of the generalized boundary layer equation for U * has the form

U * i = a i V * 1 - 1 2 + b i V * 1 1 2 + c i , (103)
and for V * , we have

V * i = α i V * 1 - 1 2 + β i V * 1 1 2 + γ i .
(104) Solution (103) for U * is obtained by applying the wall conditions on U * . This means that the coefficients a i , b i and c i are calculated by taking into account these boundary conditions. For V * , the above solution (104) is obtained by using the condition V * = 0 along one wall only, for example the lower wall. Indeed, only one boundary condition can be taken into account because V * is deduced from U * by integrating the continuity equation which is of first order with respect to the transverse coordinate.

At this stage, two conditions are not yet prescribed. The first condition is V * = 0 along the other wall (the upper wall); this condition is equivalent to the prescription of a constant mass flow in the whole channel. The second condition is the coupling between the boundary layer solution and the core flow solution. These two missing conditions are written from the expressions of V * 1i , U * i and V * i . For example, the coupling condition is

A c V * 1 - 1 2 + B c V * 1 1 2 + C c - dH dξ U * 0 (η c ) = α c V 1 * - 1 2 + β c V * 1 1 2 + γ c + U * c η c dζ dξ , (105)
where the index c denotes the given core line along which the coupling conditions hold. Eq. ( 105) gives the first relation to calculate the values of V * the calculation can proceed to the next ξ -station. In this way, the whole domain is swept from upstream to downstream. The sweep is repeated until convergence is achieved. This procedure ensures a quasi-simultaneous resolution of the whole system of equations and makes the algorithm very efficient [START_REF] Veldman | Matched asymptotic expansions and the numerical treatment of viscous-inviscid interaction[END_REF].

Results

To assess the validity of the interactive boundary layer (IBL) method, we examined the evolution of the skin-friction coefficient which is a very sensitive flow feature,

C f = 2 R τ w ,
where τ w is the reduced wall shear stress. For all cases, the Reynolds number is R = 1000.

First, comparisons between the IBL and Navier-Stokes results are given for the purely symmetric and antisymmetric wall deformations in Fig. 2. The lower wall is deformed by a trough located in the domain -2 ≤ x ≤ 2 (L = 4) with h l = -0.3. The upper wall is deformed in the same domain -2 ≤ x ≤ 2 (L = 4) with h u = 0.3 for the symmetric case and h u = -0.3 for the antisymmetric case. Even though the amplitude of the wall deformation is not really small as required by the formal theory, an excellent agreement with the Navier-Stokes results is observed.

It can be noted that there is a very slight upstream influence in the symmetric case. In fact, the pressure is also nearly constant in a cross-section. This is not the case for antisymmetric geometries where the upstream influence is clearly evidenced.

The IBL model enables us also to treat original problems. In the case of a bend, when the channel does not recover its initial axis location at the downstream end, the usual techniques of small perturbations do not work any longer. As an example, the walls are deformed in the domain -2 ≤ x ≤ 0 (L = 4) with h l = 0.5 and h u = 0.5; for x > 0, we have y l = 0, y u = 1, so the channel axis is displaced from y = 0 upstream to y = 0.5 downstream.

In this case again, a good agreement with Navier-Stokes results is observed (Fig. 3).

In this latter case, if the flow is far from being separated, a slight constriction can be added in the bend to induce separation deliberately. Thus, the walls are deformed by an additional symmetric constriction given by (106) with L = 2 in the domain -3 ≤ x ≤ -1 with h l = 0.1 and h u = -0.1. The same configuration with h l = 0.2 and h u = -0.2 has also been calculated. The results are given in Figs. 4 and5.

Then, we can observe that an antisymmetric deformation is much less disturbing than a symmetric deformation of same amplitude except if the flow is excited by a small upstream perturbation.

Along the same lines, it should be possible to consider channels whose downstream cross-section is not equal to the upstream cross-section. However, it is possible to give an idea of the quality of the results obtained with the IBL model presented here. To show this, the calculations are performed with the transformation (97)-(98) and with the basic flow (99). For instance, Figs. 6 and7 show the results obtained for a constricted channel and a dilated channel. The shape of these channels is defined from (106). The walls are deformed in the domain -2 ≤ x ≤ 0 with L = 4. In the case of the constricted channel, we have h l = 0.2 and h u = -0.2; in the downstream part (x > 0), we have y l = -0.3 and y u = 0.3. In the case of the dilated channel, we have h l = -0.4 and h u = 0.4; in the downstream part (x > 0), we have y l = -0.9 and y u = 0.9.

A very good agreement with the Navier-Stokes results is obtained again, even if a significant separated zone is present in the case of the dilated channel.

Conclusion

The method presented in this paper enables us to consider significant deformations of the channel walls. A few representative cases have been selected and the skin-friction coefficients have been compared to Navier-Stokes solutions.

The skin-friction coefficient is a very sensitive flow feature. A very good agreement is observed in difficult cases which cannot be treated easily with other methods such as the method developed by Smith which is the counterpart of the triple deck theory for internal flows.

The idea of strong coupling between the linear core flow equations and the nonlinear uniformly valid equations is called here the interacting boundary layer, IBL, for simplicity. We use also GIBL, G for global, since the equations are valid in the whole domain. This approach is systematically used so that a non-hierarchical process is employed even if it is not necessary. Moreover, there is no matching and the boundary conditions are exactly satisfied.

The main goal of this study was to separate, as far as possible, the symmetrical problem and the antisymmetrical one. The final problem is geometrically symmetrized and the interaction with the antisymmetrical effects appears in the transverse Euler equation (68) through the second derivative of H, knowing that the problem is fully symmetrical for H = 0.

Separation is more easily produced with a symmetrical dilatation. In that case, for H = 0, numerical calculations show that V 1 = 0, and the transverse pressure variation is null which is a classical hypothesis and a good approximation. The drawback of taking a constant pressure in a cross-section is that this assumption prevents us from observing any upstream influence. However, in practice, this upstream influence is very slight.

In contrast, for the antisymmetrical case, the upstream influence is important as the transverse pressure gradient is not negligible. Moreover, when E = 0, as H can be of order 1, Poiseuille flow is only displaced so that the approximation (45) is only good locally. The analysis of the upstream interaction has been done in this framework by Smith [START_REF] Smith | Upstream interactions in channel flows[END_REF]. According to hypothesis (H2), within our degree of approximation, the second derivative of H is the curvature of the curve y = H(x). This suggests that it should be possible to consider channels whose downstream direction is no longer the direction of the x-axis, and thus is different from the upstream direction. In fact, the hypothesis that the first derivative of H is small is not necessary. In the same way, in the numerical applications, we considered channels whose downstream cross-section is not equal to the upstream cross-section. Thus, it is possible to change the fact that E is small. Obviously, the study of such cases will be considered in the frame of an extended theory with curvilinear coordinates, but interesting results have been obtained with the present approach.

This study confirms the soundness of the SCEM which, we recall, is based on the construction of a uniformly valid approximation. This is particularly fruitful in the case of strong coupling when, for example, boundary layer separation can occur, which implies a significant modification of the flow. The numerical calculations are faster for these simplified models than for Navier-Stokes solutions, but this is not the primary goal of the development of these models. In fact, this new asymptotic analysis, the so-called SCEM, helps us to understand the flow structure; this step was necessary for further important problems. It is clear that future developments are not lacking. If channel flows proved to be an interesting test to approach other problems such as flow ducts, it seems valuable to look at external flows and wakes with flow control, in laminar or turbulent regimes, which are a fundamental target and for which the SCEM seems to be particularly well adapted.

Fig. 1 .

 1 Fig. 1. Flow in a two-dimensional channel with deformed walls. In this figure, all quantities are dimensionless.

Fig. 2 .

 2 Fig. 2. Flow produced by symmetric and antisymmetric wall deformations. R = 1000. NS = Navier-Stokes results, IBL = interactive boundary layer results. The Navier-Stokes calculations were performed with the commercial CFD package FLUENT. The test-case geometry is a channel whose walls are deformed in a domain x 1 ≤ x ≤ x 2 : F = h l 2 1 + cos 2π x L ; G = -h u 2 1 + cos 2π x L . (106)

Fig. 3 .

 3 Fig. 3. Flow produced by a channel bend. R = 1000. NS = Navier-Stokes results, IBL = interactive boundary layer results.

Fig. 4 .

 4 Fig. 4. Flow produced by a the superposition of a bend and of a symmetric constriction. R = 1000. NS = Navier-Stokes results, IBL = interactive boundary layer results.

Fig. 5 .

 5 Fig. 5. Flow produced by the superposition of a bend and of a symmetric constriction. Streamlines calculated in the separated region with FLUENT.

Fig. 6 .

 6 Fig. 6. Flow produced by a symmetric constriction. R = 1000. NS = Navier-Stokes results, IBL = interactive boundary layer results.

Fig. 7 .

 7 Fig. 7. Flow produced by a symmetric dilatation. R = 1000. NS = Navier-Stokes results, IBL = interactive boundary layer results.