Jheisson Lopez
email: jalopez@ungs.edu.ar

Danny Munera
email: danny.munera@udea.edu.co

Daniel Diaz
email: daniel.diaz@univ-paris1.fr

Salvador Abreu

On Integrating Population-Based Metaheuristics with Cooperative Parallelism

Keywords: parallelism, cooperation, metaheuristics, genetic algorithm, hybridization, QAP

Many real-life applications can be formulated as Combinatorial Optimization Problems, the solution of which is often challenging due to their intrinsic difficulty. At present, the most effective methods to address the hardest problems entail the hybridization of metaheuristics and cooperative parallelism. Recently, a framework called CPLS has been proposed, which eases the cooperative parallelization of local search solvers. Being able to run different heuristics in parallel, CPLS has opened a new way to hybridize metaheuristics, thanks to its cooperative parallelism mechanism. However, CPLS is mainly designed for local search methods. In this paper we seek to overcome the current CPLS limitation, extending it to enable populationbased metaheuristics in the hybridization process. We discuss an initial prototype implementation for Quadratic Assignment Problem combining a Genetic Algorithm with two local search procedures. Our experiments on hard instances of QAP show that this hybrid solver performs competitively w.r.t. dedicated QAP parallel solvers.

I. INTRODUCTION

Combinatorial Optimization Problems (COP) are widely used to model and solve real-life problems in several application domains, such as shortest/cheapest round trips, scheduling, planning, time-tabling, resource allocation, network configuration and monitoring, matching, business product line modeling, system management. . . Solving these problems rapidly represents a real challenge due to their inherent difficulty (most of them are NP-hard). The study of COPs as well as the search for efficient algorithms to solve them, has therefore been a very active research area for many years. Due to the computational intractability of most problems, we can find (or generate) problem instances of varying degrees of difficulty.

Easy and medium-sized problems can be solved with exact methods, which guarantee the optimality of the computed solution or prove that a problem does not admit any solution. This is possible because these methods consider the entire search space: either explicitly by exhaustive search or implicitly, by detecting when a portion of the search space can be ignored and get pruned off.

For harder problems, approximation methods are more appropriate. These methods aim at finding a solution of "good" quality (possibly sub-optimal) in a reasonable amount of time.

Metaheuristics have proven to be the most efficient approximation methods to address this class of problems. Metaheuristics are high-level procedures using choices (i.e., heuristics) to limit the part of the search space which actually gets visited, in order to make problems tractable. Metaheuristics can be classified in two main categories: single-solution and populationbased methods. Single-solution metaheuristics maintain, modify and stepwise improve on a single candidate solution (also called a configuration). Examples of single-solutions methods include: Simulated Annealing [START_REF] Kirkpatrick | Optimization by Simulated Annealing[END_REF], Local Search [START_REF] Aarts | Local Search in Combinatorial Optimization[END_REF], Tabu Search [START_REF] Glover | Tabu Search[END_REF], Variable Neighborhood Search [START_REF] Mladenovic | Variable Neighborhood Search[END_REF], Adaptive Search [START_REF] Codognet | Yet Another Local Search Method for Constraint Solving[END_REF], Extremal Optimization [START_REF] Boettcher | Nature's way of optimizing[END_REF]. . . On the other hand, population-based methods, modify and improve a population, i.e. sets of candidate solutions. These methods include: Genetic Algorithms [START_REF] Goldberg | Genetic Algorithms in Search, Optimization, and Machine Learning[END_REF], Ant Colony Optimization [START_REF] Solnon | Ants Can Solve Constraint Satisfaction Problems[END_REF], Particle Swarm Optimization [START_REF] Kennedy | Particle swarm optimization[END_REF]. . . Metaheuristics generally implement two main search strategies: intensification and diversification, also called exploitation and exploration [START_REF] Blum | Metaheuristics in Combinatorial Optimization: Overview and Conceptual Comparison[END_REF]. Intensification guides the solver to deeply explore a promising part of the search space. In contrast, diversification aims at extending the search onto different parts of the search space [START_REF] Hoos | Stochastic Local Search: Foundations and Applications[END_REF]. In order to obtain the best performance, a metaheuristic should provide a useful balance between intensification and diversification. By design, some heuristics are better at intensifying the search while others are so at diversifying it. More generally, each metaheuristic has it own strengths and weaknesses, which may greatly vary according to the problem or even instance being solved.

For yet more difficult problems, the current trend is therefore to design hybrid metaheuristics, by combining different metaheuristics in order to benefit from the individual advantages of each method. An effective approach consists in combining a population-based method with a single-solution method (often a local search procedure). These hybrid methods are called memetic algorithms [START_REF] Moscato | Memetic algorithms[END_REF]. However, hybrid metaheuristics are complex procedures, tricky to design, implement and tune.

An orthogonal approach to address very difficult problems consists in using parallel computation. For instance, several instances of a given metaheuristic can be executed in parallel in order to develop concurrent explorations of the search space, either independently or cooperatively by means of communication between concurrent processes. The independent approach is easiest to implement on parallel computers, since no communication is needed between the processes running a metaheuristics. The whole execution simply stops as soon as any process finds a solution. For some problems this approach provides very good results but in many cases the speedup tends to taper off when increasing the number of processors. A cooperative approach entails adding a communication mechanism in order to share or exchange information among solver instances during the search process. However, designing an efficient cooperative method is a dauntingly complex task [START_REF] Crainic | Cooperative Parallel Variable Neighborhood Search for the p-Median[END_REF], and many issues must be solved: What information is exchanged? Between which processes is it exchanged? When is the information exchanged? How is it exchanged? How is the imported data used? [START_REF] Toulouse | Communication Issues in Designing Cooperative Multi-Thread Parallel Searches[END_REF]. Moreover, most cooperative choices are problem-dependent (and sometimes even instance-dependent). Bad choices result in poor performance, possibly much worse than what could be obtained with independent parallelism. However, a welltuned cooperation scheme may significantly outperform the independent approach. To this end, a framework -called CPLS -which eases the cooperative parallelization of local search solvers has been recently proposed in [START_REF] Munera | Flexible Cooperation in Parallel Local Search[END_REF], [START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF]. This framework, available as an open source library written in the IBM X10 concurrent programming language, allows the programmer to tune the search process through an extensive set of parameters. This framework has been successfully used to tackle stable matching problems [START_REF] Munera | Solving Hard Stable Matching Problems via Local Search and Cooperative Parallelization[END_REF] and very difficult instances of the Quadratic Assignment Problem (QAP) [START_REF] Munera | Solving the Quadratic Assignment Problem with Cooperative Parallel Extremal Optimization[END_REF].

More interestingly, as it is able to run different heuristics in parallel, CPLS has opened a new way to hybridize metaheuristics, by exploiting its solution-sharing cooperative parallelism mechanism. The user only needs to code (in X10) each of the desired metaheuristics, independently, and may rely on CPLS to provide both parallelism and cooperation to get "the best of both worlds". At runtime, the parallel instances of each different metaheuristic communicate their best solutions, and one of them may forgo its current computation and adopt a better solution from the others, hoping it will converge faster. The expected outcome is that a solution which may be stagnating for one solver, has a chance to be improved on by another metaheuristic. CPLS has been successfully used to develop a very efficient hybrid solver for QAP, called ParEOTS, by implicitly combining Extremal Optimization and Tabu Search [START_REF] Munera | Hybridization as Cooperative Parallelism for the Quadratic Assignment Problem[END_REF].

In this work we are interested in hybridizing single-solution metaheuristics with population-based metaheuristics and in particular in memetic algorithms. Theoretically it is possible to plug a population-based metaheuristics into CPLS, we thus tried to hybridize a genetic algorithm with the local search procedures used in ParEOTS (extremal optimization and tabu search) [START_REF] Munera | Hybridization as Cooperative Parallelism for the Quadratic Assignment Problem[END_REF]. It turned out that the resulting solver performed worse than ParEOTS! An explanation for this observation is that the CPLS basic cooperation mechanisms are designed for single-solution metaheuristics, and also that CPLS is better at intensification than diversification. ParEOTS performed well because extremal optimization can be conveniently tuned to provide an appropriate level of diversification. Even though some single-solution metaheuristics are rather good at diversification, when well tuned, population-based metaheuristics offers a wider variety of diversification thanks to the large possibilities of evolution offered by a population. However, to be effective in an hybrid procedure, this population should encompass many solutions (e.g. all solution candidates computed by all solver instances).

In this paper we seek to overcome the current CPLS limitation, extending it to enable population-based metaheuristics in the hybridization process. We experimented with an initial prototype for QAP, which we show to already perform competitively w.r.t. ParEOTS and other state-of-the-art parallel solvers.

The rest of this paper is organized as follows: in section II, we describe the existing CPLS and discuss some of its limitations. In section III we present the adaptation of CPLS used in our prototype implementation in order to hybridize a genetic algorithm with two local search procedures to solve QAP. Section IV provides an initial experimental assessment of our prototype. The paper ends with some considerations for future directions.

II. THE COOPERATIVE PARALLEL LOCAL SEARCH FRAMEWORK

The Cooperative Parallel Local Search framework (CPLS) [START_REF] Munera | Flexible Cooperation in Parallel Local Search[END_REF], [START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF] was designed to facilitate the implementation of cooperative parallel local search solvers. CPLS is available as an open source library written in the IBM X10 concurrent programming language. From a user point of view, one needs only code the desired metaheuristic, in X10, using the CPLS Application Programming Interface (API). This API provides object-oriented abstractions to hide all parallel management details from the user, who writes a purely sequential program, without having to worry about parallelism and communication. The resulting solver can be executed on a wide variety of machine (multicores, manycores, clusters,. . . ), implicitly using the parallelism provided by the X10 programming system. Furthermore, CPLS augments the independent parallelism strategy with a tunable communication mechanism, which enables cooperation between the multiple metaheuristic instances. It also offers several parameters which control this mechanism. To activate cooperative parallelism, the user simply annotates the existing code with calls to CPLS functions to periodically execute 2 actions1 :

• report on its current configuration (Report action).

• retrieve a solution 2 and, if it is good enough, the solver adopts it, i.e. it abandons its current solution and replaces it with the received one (Update action). At runtime, each local search instance is encapsulated in an explorer node. The point is to use all available processing units by mapping each explorer node to a physical core. CPLS organizes the explorers into teams, where a team consists of NPT explorers which aim at intensifying the search in a particular region of the search space using intra-team communication. NPT can range from 1 to the maximum number of cores. This parameter is directly related to the trade-off between intensification and diversification since it is expected that different teams explore different regions of the search space. When NPT is 1, the framework coincides with the independent parallelism strategy, it is expected that each 1node team be working on a different region of the search space, without any effort to seek parallel intensification. When NPT is equal to the maximum number of nodes (creating only 1 team in the execution), the framework hits the maximum level of intensification (note that a certain amount of diversification is inherently provided by parallelism, between 2 cooperation actions, due to the stochastic nature of metaheuristics.) Tuning the value of NPT allows the user to adjust the trade-off between intensification and diversification.

Head Node

Explorer n

Fig. 1: CPLS team structure By design, a team seeks to intensify the search in the most promising neighborhood found by any of its members. The parameters which guide the intensification are the Report Interval R and Update Interval U -every R iterations, each explorer node sends its current configuration and the associated cost metric to its head node. The head node is a particular explorer inside the team which periodically collects and processes this information, retaining the best configurations in the Elite Pool EP of size |EP|. Every U iterations, explorer nodes randomly retrieve a configuration from EP, in the head node. An explorer node may adopt the configuration from the EP, if it is "better" than its own current configuration with a probability pAdopt. Figure 1 depicts the structure of a team.

In addition to the notion of teams of parametric size, CPLS provides an additional mechanism for diversification through inter-team communications. For this, a particular node, called the master node, periodically compares each team's best configuration (taken from its elite pool) using a distance function (e.g., using a Hamming distance). If two or more teams are too close, the master node selects a team to perform a corrective action (this can be the "worst" team, i.e. the one whose best configuration cost is highest). The corrective action -which is indeed a diversification action -aims at forcing all explorers (or a fraction thereof) in the selected team, to explore a different area of the search space.

It is worth noting that CPLS abstracts the notion of metaheuristic and thus accepts any metaheuristic, including others besides local search. Several configurations are possible; the most common one consists in executing in parallel several instances of exactly the same metaheuristic (clones). One may also execute the same metaheuristic but configured differently, i.e. with different parameter tuning. For instance, the same tabu search procedure could be executed with different tabu tenure values. A portfolio approach can be obtained with the independent parallel execution of different metaheuristics. Finally, it is also possible to execute different metaheuristics, within the cooperative scheme just described. The resulting solver behaves like a hybrid solver, just that it's achieved without the complexity required by the design and coding of the hybrid solver.

ParEOTS [START_REF] Munera | Hybridization as Cooperative Parallelism for the Quadratic Assignment Problem[END_REF] is a powerful hybrid solver for QAP combining Extremal Optimization and Tabu Search. Being interested in memetic algorithms, we tried to add a Genetic Algorithm (GA) as third component to this hybrid solver. The integration was simple but we the resulting solver was not any better than ParEOTS. We sought the CPLS characteristics which prevent the efficient use of population-based metaheuristics for hybridization. It turns out that the CPLS basic cooperation mechanisms are designed for single-solution metaheuristics; with the Elite Pool, CPLS is intrinsically better at intensification than diversification. The diversification mechanism described above is very difficult to tune (we could never obtain any gain from it). Additionally, ParEOTS does not use it. This mechanism has two major drawbacks:

1) It is a centralized process: the master node has the big responsibility of making the right decision (it risks breaking a winning team). It acts with a very limited view of the computation: basically the best solutions computed so far by a team. In a sense, diversification is not sufficiently well informed to be effective when the likely next moves are taken from larger sets, as happens with population-based metaheuristics such as GA. 2) CPLS is biased towards execution on distributed architectures based on clusters of multicore computers: teams are mapped to nodes and explorers run on individual cores. In this setting, inter-team communications tend to be relatively expensive because they require node-tonode communication. 3 It can be argued that diversification actions are costly and should therefore be used with parsimony. ParEOTS performs well because extremal optimization can be conveniently tuned to provide an appropriate level of diversification. Even though some single-solution metaheuristics are rather good at diversification, population-based metaheuristics generally offer a wider spectrum of diversification, when well tuned, because of the large evolutionary possibilities afforded by a population. However, to be effective in a hybrid procedure, the population should encompass many solutions (e.g. all solutions computed inside a team). This was not the case with our previous GA experiment: the GA was treated like any other metaheuristic, it has its own private data structures and only sends and receives single solutions from time to time. This drove us to conjecture that the notion of population should be promoted so as to share a global vision of what happens across the team, thereby contributing to a better diversification decision. This is the idea underlying our development, which we now flesh out.

III. HYBRIDIZING WITH GA

In order to provide an efficient interaction, within the CPLS framework, between population-based and local search metaheuristics, we conducted several experiments seeking to enhance the intra-team communication strategy. We tested several architectures, parameters and components, and finally we developed a prototype, dedicated to solving QAP, that actually improves on the performance of the original CPLS.

Our proposal, as shown in Figure 2, considers the inclusion of a genetic algorithm (GA) in the head node interacting with three explorer nodes: two implementing the robust tabu search (RoTS) method and one implementing the extremal optimization (EO) method. These methods are detailed in [START_REF] Munera | Solving the Quadratic Assignment Problem with Cooperative Parallel Extremal Optimization[END_REF], [START_REF] Munera | Hybridization as Cooperative Parallelism for the Quadratic Assignment Problem[END_REF]. We extend the CPLS model by having two independent pools of configurations inside the head node: the diversification pool (DivPool) and the intensification pool (IntPool).

Head Node

Update individual from IntPool

Report individual to DivPool

Explorer Nodes: EO -RoTS

Head Node: GA Fig. 2: CPLS new team structure LS nodes start with a random configuration and the head node with a random population composed of individuals (i.e. configurations). They all interact with the pools as follows:

1) LS Explorers to IntPool: every R S iterations, explorer nodes implementing local search methods (EO and RoTS) report their current configuration and corresponding cost to the IntPool. 2) IntPool to Population: every U P iterations, the GA replaces a random individual in its population with a random configuration from the IntPool. 3) Population to DivPool: every R P iterations, the GA reports an individual to the DivPool. 4) DivPool to LS Explorer: every U S iterations, the explorers nodes implementing LS methods randomly retrieve a configuration from the DivPool and may adopt this configuration, if it is "better" than its current configuration. During the execution of the solving process, the head node's IntPool collects elite configurations from LS explorer nodes. Then, the GA generates offspring configurations (possibly stemming from some those directly generated by the LS metaheuristics) and feeds the DivPool with some individuals. LS explorer nodes consume configurations from the DivPool and possibly adopt them, closing the interaction loop inside the team. Note that this process emulates the basic idea of memetic algorithms, while keeping the main functionality of each metaheuristic unaltered and providing hybridization through cooperative parallelism.

Finally, we propose a mechanism to ensure entropy in the GA's population, mainly because GA has the role of a "diversification agent" inside the team. If the GA's population becomes uniform (not diverse), it could mean that enough intensification has been done on this region, and it would be more interesting to explore other regions of the search space. The diversity of the population may be measured every diversity check interval dCkI with a diversity check function dCkF. If the dCkF value is lower than a minimum allowed diversity value mdV, then a percentage of the population (pPr) gets restarted. Nevertheless, a 1 -pPr portion of the population is kept as memory of the search process, and it only gets replaced if the explorations in the new region give rise to better configurations.

A. Parameter selection

Implementing these changes requires finding suitable values for a new set of parameters. We now proceed with discussing the design and implementation decisions which we followed to select particular values for these parameters.

In our proposal, IntPool and DivPool behavior is controlled by the CPLS pool parameters, e.g., pool size, entry policy and request policy [START_REF] Munera | Flexible Cooperation in Parallel Local Search[END_REF], [START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF]. We use the CPLS default values for both pools' parameters (the size of the pool is 4 and stores the best solutions; on demand a random solution is picked form the pool). This turns out to be convenient for our model, as we already ensure that configurations reported by a given metaheuristic only compete with other configurations reported by metaheuristics of the same nature.

To simplify the values selection for the report/update parameters, we use global R and U values, as in the original CPLS. Based on these global values, we computed the individual values for R P , U P , R S and U S , as follows: R p = R/N P T and U p = U/N P T ; whereas in the explorer nodes, R s = R/(problem size) and U s = U .

We decided to use the normalized entropy, proposed in [START_REF] Fleurent | Genetic hybrids for the quadratic assignment problem[END_REF], to measure the diversity of the GA's population in the range from 0 to 1. An entropy value equal to 0 indicates a population consisting of multiple repetitions of the same individual, while an entropy value equal to 1 specify a population composed of individuals that are all different from each other.

B. The Genetic Algorithm Design

Genetic algorithms (GA) are inspired by Darwin's law of natural selection. The basic idea is that by combining different bits of the best adapted individuals you can get even better individuals. Hybridization of GAs is one of the most successful strategies that have been proposed to solve QAP [START_REF] Drezner | Extensive experiments with hybrid genetic algorithms for the solution of the quadratic assignment problem[END_REF]. For this reason, we have designed a GA and implemented it in CPLS, to test the new intra-team communication scheme that we described above. On one hand, a genetic algorithm basically consists of five stages that are repeated cyclically (each of these cycles is named generation): initialization of its population, selection of individuals for crossing, crossing, mutation and selection of the population for the next generation. On the other hand, a GA must define a way to encode individuals as well as two main parameters: the population size and the mutation rate. We now discuss the decisions we made for the stages and the parameters mentioned above. The pseudo-code for our GA is presented as Algorithm 1. We encode the individuals as permutations, as this allows for direct evaluation of the cost (it does not require a decoding process) and therefore takes less processing time. Besides, as suggested by Ahuja, et al. [START_REF] Ahuja | A greedy genetic algorithm for the quadratic assignment problem[END_REF], the GA is largely independent from the initialization mechanism for its population, as long as sufficiently diverse populations are generated. For this reason, we chose a random initialization of the population as it produces individuals uniformly distributed throughout the search space.

Algorithm 1 GA for hybridization

We decided to use a population size (pz) equal to the instance size. As the size of the population determines the size of the search space region to be explored by each team, it makes sense to size it proportionately to the instance. Moreover, experimentation shows that populations greater than the instance size have a negative impact on the performance of the algorithm: the population must be sorted and the offspring should be compared with it in each generation. In a probatory experiment, we observed that population sizes less than the instance size n have a very low initial entropy, while population sizes greater that n have an initial entropy barely greater than one of size n.

Regarding the selection for crossing, in [START_REF] Ahuja | A greedy genetic algorithm for the quadratic assignment problem[END_REF] the authors conclude that the use of a mechanism that favors the best individuals (i.e. an elitist one) does not necessarily lead to better results, and can even negatively affect them by hastening the convergence towards local optima. In order to establish an intermediate point of elitism, we adopt the strategy proposed by Lalla-Ruiz et al. [START_REF] Lalla-Ruiz | A hybrid biased random key genetic algorithm for the quadratic assignment problem[END_REF]: the total population is divided into a group of elite individuals and another of non-elite ones; an element from each group is randomly selected and they are used to carry out the crossing process. The individuals generated through this process are called the offspring. We have used Uniform Crossover (UX) because, according to Benlic and Hao [START_REF] Benlic | Memetic Search for the Quadratic Assignment Problem[END_REF], it offers the best results for QAP. The crossing process generates two offspring in each generation.

Classically, for QAP, the mutation operator consists in swapping two of the individual's positions (2-exchanges are widely used for permutation problems). The mutation rate is a parameter which determines the probability of applying the mutation operator on the offspring. However, in our prototype, the mutation rate (mr) is the percentage of individuals that will be mutated (taking into account that a swap modifies 2 individuals). In this setting, we use a mutation rate around 40%. This value is very high for classical genetic algorithms but in our contect, the genetic algorithm is mainly used to ensure a high degree of diversification for other explorers.

Finally, we decided to use an elitist criterion in the evaluation of the offspring for deciding whether they are to be included in the population for the next generation: an offpring individual replaces the worst individual of the population if it has better cost and is not already present in the population.

IV. EXPERIMENTATION

We performed an experimental evaluation of our prototype implementation 4 . All experiments have been carried out on server with 4 × 16-core AMD Opteron 6376 CPUs, running at 2.3 GHz and 128 GB of RAM. This is the same machine as used for ParEOTS evaluation in [START_REF] Munera | Hybridization as Cooperative Parallelism for the Quadratic Assignment Problem[END_REF] but we could only use 32 cores (while the reference paper presents results for 128 cores).

For this evaluation we selected the 33 hardest instances of the QAPLIB benchmarks highlighted in [START_REF] Munera | Solving the Quadratic Assignment Problem with Cooperative Parallel Extremal Optimization[END_REF], [START_REF] Munera | Hybridization as Cooperative Parallelism for the Quadratic Assignment Problem[END_REF]. Each

TABLE I: Comparison of GA-CPLS with ParEOTS

We first compare our prototype implementation (called GA-CPLS in the rest of this section) with ParEOTS, our reference. Table I reports for each benchmark:

• the number of times the BKS is reached out of 10 runs (#BKS). • the Average Percentage Deviation (APD) which is the average of the 10 relative deviation percentages, computed as follows: 100 × F (sol)-BKS

BKS

, where F (sol) is the cost of the obtained solution. This information is useful for comparing the average quality of solutions (across the 10 executions) found by the solvers.

• the Best Percentage Deviation (BPD) which corresponds to the relative deviation percentage of the best solution found among the 10 executions. • the average execution time given in second which corresponds to the actual elapsed times, and include the time to install all solver instances, solve the problem, communications and the time to detect and propagate the termination.

To compare the two solvers, we first compare the number of BKS found, then (in case of a tie), the APD and finally the execution time. For each benchmark, the best-performing solver row is highlighted and the discriminant field is enhanced in bold font. GA-CPLS clearly outperforms ParEOTS. Several instances are now systematically solved at each replication (e.g. sko90 or tai100b or wil100). Moreover, for the very difficult tho150 which was never solved by ParEOTS, GA-CPLS can reach the BKS 3 times. ParEOTS is better on tai50a-tai100a. These problems require a very strong intensification to reach the BKS. It would be interesting to experiment with different parameters for GA-CPLS to increase its level of intensification on these problems. The "Summary" row gives relevant information: GA-CPLS reaches the BKS 261 times (over the 330 possible) while ParEOTS only reaches 221. The APD are on average better than with ParEOTS; similarly for the best reached solution. Finally, regarding execution time, our prototype solver is, on average, 1.5 times fastet than ParEOTS. The difference is particularly impressive on instances like: tai35b (ParEOTS needs 18.784sec while GA-CPLS only needs 0.068sec. Clearly, the hybridization with a GA now works really well and definitely improves performance.

It is also interesting to compare GA-CPLS with state-ofarts parallel solvers for QAP. We have selected the following challengers:

• COSEARCH [START_REF] Talbi | COSEARCH: A parallel cooperative metaheuristic[END_REF], which is a cooperative parallel method based on tabu search, genetic algorithms and a kick operator. The (published) results are obtained from a heterogeneous grid computing platform using around 150 processors.

• CPTS [START_REF] James | A Cooperative Parallel Tabu Search Algorithm for the Quadratic Assignment Problem[END_REF]: a cooperative parallel tabu search algorithm in which processors exchange information to intensify or diversify the search. The results are obtained from the cited paper using 10 Intel Itanium processors (1.3 GHz).

• PHA [START_REF] Tosun | On the Performance of Parallel Hybrid Algorithms for the Solution of the Quadratic Assignment Problem[END_REF]: a hybrid approach based on parallel island genetic algorithm and robust tabu search. Experimentation was performed on a cluster of 46 nodes, each with 8 cores. • TLBO-RTS [START_REF] Dokeroglu | Hybrid Teaching-Learning-Based Optimization Algorithms for the Quadratic Assignment Problem[END_REF]: a hybrid teaching-learning based optimization approach. This algorithm combines the teaching-learning based optimization to optimize a population, then, the population is given to a robust tabu search technique. TLBO-RTS supports a parallel execution using a parallel island technique. The results are obtained using 40 to 50 processors.

Table II reports all available information from the cited papers. Execution times are given either in seconds for small values (as a decimal number) or in a human readable form (as mm:ss or hh:mm:ss). These times are mainly given for information purpose. Unfortunately, some methods do not provide the number of obtained BKS (obviously an APD = 0.000 corresponds to #BKS = 10). For this reason TABLE II: Comparison of GA-CPLS with other parallel solvers we consider the APD, and highlight the solver with the best value. GA-CPLS has a better, or at least equal, APD than COSEARCH for all common instances. Comparing the performance with COSEARCH, out of 11 instances, it is apparent that COSEARCH performs better than GA-CPLS in only one instance: tai256c; note that no execution time data is available for COSEARCH. Regarding CPTS, out of 31 common instances, it appears that it only performs better than GA-CPLS in two cases: tai150b and tai256c, but taking a very long time: more than 122 hours! Here, we obtain a slightly worse value but in just 5 minutes. A similar situation appears with TLBO-RTS which only performs better than GA-CPLS on two instances: tai40a and tai150b, for which it requires 7 hours. The comparison with PHA is more balanced; both systems are able to reach the BKS of most problems at each replication. PHA is better in 6 cases (sko100f, tai40a-tai60a, tai150b and tai256c). These instances require strong intensification and/or are large size (greater than 100). Clearly, a 5 minutes time limit is too short for these problems (PHA needs around 3 hours to obtain its performance on tai150b and more than 11 hours for tai256c). However, it is worth mentioning, our method returns better solutions for 3 large instances (tai80a, tai100a and the very difficult tho150). The "summary" row confirms that PHA provides slightly better solutions than GA-CPLS while our prototype is clearly faster. In conclusion, in terms of solution quality, GA-CPLS is one of the best performing parallel methods. Moreover, in terms of execution time, our method has the best run times while maintaining a high quality of the solution and making efficient use of parallel resources. We are aware that COSEARCH and CPTS were reported on several years ago, on slower machines, and that both TLBO-RTS and PHA run with a different (higher) core count, but the results are significant nonetheless.

V. CONCLUSION AND FURTHER WORK

Hybrid cooperative parallel metaheuristics is currently one of the most efficient approaches to address very hard Combinatorial Optimization Problems. Designing such a solver is a very complex task which requires very strong skills and considerable effort in parallel programming, modeling, metaheuristics, hybridization and runtime tuning. Recently, the CPLS framework has been proposed to ease the design and the implementation of such hybrid parallel solvers. However, CPLS is oriented towards single-solution metaheuristics, a bias which we have experimentally confirmed. It is all the more regrettable, considering that population-based metaheuristics have shown to be very efficient when combined with other local search procedures, as testified by the success of memetic algorithms.

In this paper we analyzed the base mechanisms of CPLS which prevent efficient hybridization with population-based metaheuristics. We then proposed a modification of CPLS, for which we did an initial prototype for the very difficult Quadratic Assignment Problem (QAP). Our experiments with this prototype indicate very good performance on the hardest instances of QAPLib, making it already competitive with other state-of-the-art parallel solvers.

Our medium-term goal is to propose a new framework, based on CPLS, to fully support hybridization with populationbased metaheuristics. For this we first need to experiment more extensively, in order to analyze the impact of parameters such as the size of the population, the reset percentages of the population, different values for parameters which control the communication. Regarding QAP, we will try to further improve our Genetic Algorithm (e.g. testing with the Insert Path Crossover as operator and crossover).

We will then test this improved solver on even more difficult QAP instances (e.g. instances proposed by Drezner, Palubeckis and Carvalho & Rahmann). We also plan to intensify the experimental evaluation, to cover a wider base of test cases and runtime configurations, namely the progression in speedup attained by scaling the computational resources (cores and nodes).

These calls must be inserted inside the main loop of the metaheuristic.

The received configuration is most likely reported by another local search instance.

This aspect tends to be mitigated by high-throughput and low-latency network interconnects, such as Infiniband.

The source code and QAP instances are available from https://github.com/ jlopezrf/COPSolver-V 2.0

problem is executed 10 times, stopping as soon as the Best Known Solution (BKS) is reached. This execution is done with a time cap of 5 minutes (in case the BKS is not reached). Such experiments give useful information regarding the quality of quickly obtainable solutions. GA


Andrei G Shvarts

Vladislav A Yastrebov

Fluid flow across a wavy channel brought in contact

Keywords: wavy channel, contact, fluid flow, sealing, fluid-solid coupling, monolithic approach, finite element analysis

A pressure driven flow in contact interface between elastic solids with wavy surfaces is studied. We consider a strong coupling between the solid and the fluid problems, which is relevant when the fluid pressure is comparable with the contact pressure. An approximate analytical solution is obtained for this coupled problem. A finite-element monolithically coupled framework is used to solve the problem numerically. A good agreement is obtained between the two solutions within the region of the validity of the analytical one. A powerlaw interface transmissivity decay is observed near the percolation. Finally, we showed that the external pressure needed to seal the channel is an affine function of the inlet pressure and does not depend on the outlet pressure.

Introduction

The problem of a thin fluid flow in narrow interfaces between contacting or slightly separated surfaces occurs in different applications. The first example is the sealing problem: seals are used to minimize or prevent leakage of fluids from and into internal chambers of numerous engineering systems, such as gas cylinders, water circuits, lubricated bearings and gears, heat engines and others. Dynamic and static seals are distinguished, the former seal interfaces between surfaces with no relative motion, the latter deal with relatively moving surfaces. Contact and non-contact seals are also distinguished: the former possess contacting parts in the sealing interface, the latter do not.

Hydraulic fracturing is another application which involves interaction of fluid and solid with possible contact between crack faces or/and with a third body, like sand particles [START_REF] Bažant | Why fracking works[END_REF]. The fluid extraction of shale gas and oil from rocks represents an antipodal problematic to sealing applications, but the physics and fluid-solid coupling remain identical. A slightly different problem involving fluid, solid and contact appears in fatigue-crack growth in lubricated rolling or cyclically sliding contacts [START_REF] Bower | The influence of crack face friction and trapped fluid on surface initiated rolling contact fatigue cracks[END_REF]. Such an interaction between fluids and solids in contact can be also found in poromechanics [START_REF] Coussy | Poromechanics[END_REF] and at larger scales in landslides [START_REF] Viesca | Nucleation of slip-weakening rupture instability in landslides by localized increase of pore pressure[END_REF], slip in pressurized faults [START_REF] Garagash | Nucleation and arrest of dynamic slip on a pressurized fault[END_REF], bazal sliding of glaciers [START_REF] Fischer | Stick slip sliding behaviour at the base of a glacier[END_REF], and in other applications.

Depending on the application, the fluid can be considered compressible or incompressible (both in gas and liquid states), and it might flow under capillary effect or pressure difference. Fluids in liquid state at high pressure may evaporate due to a pressure decrease along the fluid path or due to temperature increase induced by frictional heating, which results in a mixed gas-liquid-solid problem: a notable example is cavitation in lubrication 1 problems. Compressibility and viscosity of fluid can depend significantly on pressure and temperature within a certain range of loading parameters. Polymers are used in most sealing applications, however, due to fluid-uptake, chemical and thermal degradation, and also the glass transition, their usage is limited; several applications require usage of metal-to-metal contacts in seals. In nature, the problems of interfacial fluids in contact interfaces are relevant for rock materials in terms of hydrogeology, shale gas and oil extraction as well as fracking, and also magma rise in volcanology [START_REF] Lister | Fluid mechanical models of crack propagation and their application to magma transport in dykes[END_REF]. In biological systems, the relevant materials are soft tissue and the applications include circulation of blood and other fluids in organisms.

Inevitable roughness, sometimes complemented by on-purpose patterning of engineering and natural surfaces, affects their sealing properties. Inversely, the presence of a fluid in the contact interface may affect the mechanical properties of seals by adding extra loadcarrying capacity, changing the interface stiffness and the friction coefficient. For soft materials and high fluid pressures, interface fluids can considerably deform the solids in elastohydrodynamic lubrication, however it may be also relevant for static seals [START_REF] Nau | Fluid Sealing Technology: Principles and Applictions[END_REF]. In contact seals this nonlinear fluid-solid interaction is intensified by nonlinear contact constraints. This coupling presents the topic of the current study.

The roughness of contacting surfaces [START_REF] Whitehouse | Handbook of surface and nanometrology[END_REF][START_REF] Thomas | Rough Surfaces[END_REF] has strong implications in mechanics and physics of contact: adhesion, friction and wear in dry and lubricated contacts are controlled to a great extent by parameters of the roughness of contacting solids. Mass and energy transport through and across contact interfaces strongly depend on the surface roughness, for instance: electric contact resistance, heat conduction between contacting solids and the sealing problem -the topic of the interest of the present paper. The roughness, or more generally the surface geometry, may contain some deterministic features (turned surfaces, patterned surfaces [START_REF] Prodanov | Contact mechanics of laser-textured surfaces[END_REF]) or be purely random, self-affine [START_REF] Nayak | Random process model of rough surfaces[END_REF] down to atomistic scale [START_REF] Krim | Experimental observations of self-affine scaling and kinetic roughening at sub-micron lengthscales[END_REF][START_REF] Ponson | Anisotropic self-affine properties of experimental fracture surfaces[END_REF]. Surface morphology may be determined by surface processing (polishing, work-hardening), underlying microstructure and its deformation (grain boundaries, plasticity induced roughness [START_REF] Šiška | Finite element simulations of the cyclic elastoplastic behaviour of copper thin films[END_REF], persistent slip marks [START_REF] Bao-Tong | Overview of fatigue behavior in copper single crystals -I. Surface morphology and stage i crack initiation sites for tests at constant strain amplitude[END_REF], rumpling [START_REF] Tolpygo | Surface rumpling of a (ni, pt) al bond coat induced by cyclic oxidation[END_REF]), corrosion and oxidation; for coatings the roughness is determined by the deposition method (physical vapor deposition, gas dynamic cold spray deposition, electroplating and others), in biology the surface is determined by the tissue growth processes and related instabilities or assigned functionalities [START_REF] Amar | Swelling instability of surface-attached gels as a model of soft tissue growth under geometric constraints[END_REF][START_REF] Li | Surface wrinkling of mucosa induced by volumetric growth: Theory, simulation and experiment[END_REF][START_REF] Wang | A three-dimensional phase diagram of growth-induced surface instabilities[END_REF]. The resulting surface morphology may be rather complex and span over many scales from atomistic to structural ones, it can be characterized by numerous parameters, such as standard deviation of heights and height gradient, height distribution, in particular its kurtosis and skewness, power spectral density, spectral moments, etc. For mechanical contact problems, since in most applications only the highest asperities come in contact, an approximation of the roughness by a number of isolated spherical or elliptic asperities results in a rather accurate and helpful model [START_REF] Archard | Elastic deformation and the laws of friction[END_REF][START_REF] Greenwood | Contact of nominally flat surfaces[END_REF][START_REF] Bush | The elastic contact of a rough surface[END_REF]. On the other hand, the fluid flow through the free volume1 is mainly affected by "deeper" surface features: grooves, valleys and dimples.

Numerical analysis of the fluid flow through contact interfaces was carried out between real [START_REF] Amyot | Critical point network for drainage between rough surfaces[END_REF][START_REF] Vallet | Real versus synthesized fractal surfaces: Contact mechanics and transport properties[END_REF][START_REF] Vallet | Sampling effect on contact and transport properties between fractal surfaces[END_REF][START_REF] Durand | Multiscale Modelling of Contact of Rough Surfaces[END_REF] or model rough topographies [START_REF] Plouraboué | Conductances between confined rough walls[END_REF][START_REF] Dapp | Self-affine elastic contacts: Percolation and leakage[END_REF][START_REF] Paggi | Evolution of the free volume between rough surfaces in contact[END_REF]. In contrast to the complexity and lack of scale separation of nominally flat realistic surface roughness, surface patterning allows to use the concept of scale separation to a certain extent and to limit the analysis to major geometrical features of the surface [START_REF] Sahlin | A mixed lubrication model incorporating measured surface topography. part 1: Theory of flow factors[END_REF][START_REF] Sahlin | A mixed lubrication model incorporating measured surface topography. part 2: Roughness treatment, model validation, and simulation[END_REF][START_REF] Pérez-Ràfols | Modelling of leakage on metal-to-metal seals[END_REF]. On the other hand, analysis of simple models of surface geometry, for example wavy and bi-wavy models [START_REF] Westergaard | Bearing pressures and cracks[END_REF][START_REF] Kuznetsov | Effect of fluid lubricant on the contact characteristics of rough elastic bodies in compression[END_REF][START_REF] Johnson | The contact of elastic regular wavy surfaces[END_REF][START_REF] Johnson | The adhesion of two elastic bodies with slightly wavy surfaces[END_REF][START_REF] Carbone | Adhesion and friction of an elastic half-space in contact with a slightly wavy rigid surface[END_REF][START_REF] Yastrebov | The contact of elastic regular wavy surfaces revisited[END_REF][START_REF] Dapp | Contact mechanics of and reynolds flow through saddle points: On the coalescence of contact patches and the leakage rate through near-critical constrictions[END_REF][START_REF] Xu | Elastic contact between a geometrically anisotropic bisinusoidal surface and a rigid base[END_REF][START_REF] Matsuda | Real contact area and friction property of rubber with two-dimensional regular wavy surface[END_REF], helps to understand better local deformation mechanisms in rough contact and the role of patterning for macroscopic behavior. A wavy channel also serves as an important test model in fluid flow analysis [START_REF] Nishimura | Flow characteristics in a channel with symmetric wavy wall for steady flow[END_REF][START_REF] Brown | Applicability of the reynolds equation for modeling fluid flow between rough surfaces[END_REF][START_REF] Sui | Fluid flow and heat transfer in wavy microchannels[END_REF][START_REF] Ničeno | Numerical analysis of fluid flow and heat transfer in periodic wavy channels[END_REF][START_REF] Dapp | Contact mechanics of and reynolds flow through saddle points: On the coalescence of contact patches and the leakage rate through near-critical constrictions[END_REF]. Here we also consider a periodic wavy channel, but contrary to other flow studies, it is brought in mechanical contact with a rigid flat and the fluid flows across the wavy section in channels delimited by mechanical contact zones. In the first approximation this model represents a "rough" surface with parallel grooves.

To analyze the fluid flow in the contact interface, two problems should be solved: a mechanical problem of the contact between solids with rough surfaces and the problem of the fluid flow through the resulting free volume. If during the loading process the fluid pressure is negligible in comparison to the contact pressure, it is possible to assume that both physics are weakly coupled. This implies that they can be solved separately: firstly, the solid contact problem is solved to obtain the deformed geometry of the channels. This geometry is then passed to the fluid solver, which resolves the fluid flow under the assumption of rigid boundaries. This solution strategy leads to a one-way coupling, whereby the contact problem is independent of the fluid pressure, while the fluid problem depends on the geometry computed by the solid solver. However, it is not rare that in the load interval of interest a stronger coupling between fluid and solid equations is required. For example, it is the case when the local contact pressure is comparable with the hydrostatic fluid pressure. Note that it is always the case near edges of contact zones at which the contact pressure (in the uncoupled case) decreases to zero as ∼ √ ξ, when the distance from the contact edge ξ decreases ξ → 0 [START_REF] Manners | Some observations on persson's diffusion theory of elastic contact[END_REF]. The so-called two-way coupling strategy allows to take into account the fluid pressure distribution and its effect on the deformation of the solid and vice versa: the effect of the elastic deformation on the fluid flow.

In many industrial applications the thickness of the free volume interface between contacting surfaces is quite small, and the flow is often laminar even for gas [START_REF] Nau | Fluid Sealing Technology: Principles and Applictions[END_REF]. Nevertheless, considering turbulent flow could be essential for more demanding applications, like drag delivery, microfluidic chemical reactors, etc. Moreover, we postulate that the variation of the mechanical loading conditions is slow enough compared to characteristic time of the fluid flow, so the flow is assumed to be stationary; capillary actions are neglected, only pressure driven flow is considered. This two-dimensional flow through the contact interface is accurately described by the Reynolds equation [30, 21, e.g.], which is used in this study.

The paper is structured as follows: in Section 2 we formulate the problem to be solved; in Section 3 we recall classical solutions for a wavy profile with a pressurized fluid present in the interface; in Section 4 we obtain an approximate analytical solution; in Section 5 a monolithic numerical scheme which couples the fluid and solid equations is briefly outlined; finally, in Section 6 the numerical results are presented and discussed, and in Section 7 we make conclusions.

Problem set-up

We consider an array of wavy channels of length L along OY direction (see Fig. 1(a)) with a sine-shape section z(x

) = ∆[1 -cos (2x )], (1)
where x = πx/λ, brought in contact with a rigid flat2 , the pressure driven flow in this channel (∆p f = p ip o , where p i and p o are the inlet and outlet pressures, respectively) of incompressible viscous fluid is governed by the stationary Reynolds equation for the Poiseuille flow.

We assume an isothermal fluid flow at a temperature at which it does not evaporate under the pressure drop on its way from the inlet to the outlet. The system of equations to

u z V Ω f Ω c Γ c Γ c (a) inlet, p i outlet, p o L flow (b) contact contact x y λ 0 a(y)
Figure 1: (a) Problem set-up: an elastic wavy surface comes in contact with a rigid flat, and an incompressible fluid flows under pressure difference from the inlet to the outlet; the resulting interface pressure distribution p(x, y) is sketched. (b) Sketch of the contact interface, note that due to fluid pressure acting on the surface of the solid, the contact patches in the sections parallel to OX are wider at the outlet, than at the inlet, which is associated with the direction of the fluid pressure drop. be solved takes the following form:

∇ • g 3 ∇p f = 0 in Ω f (2)
p f y=0 = p i , p f y=L = p o , [∇p f • e x ] x=0,λ/2 = 0 (3) ∇ • σ = 0 in V (4) u x | x=0,λ = 0, u y y=0,L = 0, σ zz | z=-∞ = -p ext (5) g ≥ 0, p -p f ≥ 0, (p -p f )g = 0 in ∂V, (6)
where Eq. ( 2) is the Reynolds equation for pressure driven stationary incompressible viscous Poisseuille flow, the distance between immobile walls is given by at least C 1 -smooth gap distribution g = g(x, y) in the domain Ω f , which is the closure of the solid volume ∂V projected on the rigid flat, and p f denotes the fluid pressure. Eq. ( 3) summarizes boundary conditions for the fluid problem: the fixed inlet p i and outlet fluid pressure p o and zero flux at crests of the surface resulting from the problem symmetry. Expression ∇p f • e x is a quantity proportional to the fluid flux in the OX direction (orthogonal to the main flow direction), and e x is a nonzero in-plane vector collinear with the axis OX. The fluid flux is given by q = -g 3 (∇ • p f )/12µ, where µ is the dynamic viscosity. Prescribing the fluid flux at the inlet/outlet of the fluid domain would only slightly change the numerical treatment, and since for contact static seals the fluid flow with prescribed hydrostatic inlet and outlet fluid pressures presents a more common situation, for the rest of the paper we will stick to this particular boundary condition. Eq. ( 4) is the momentum balance equation for the quasi-static solid mechanical problem in absence of volumetric forces, while (5) summarizes boundary conditions for the solid problem, where p ext is the squeezing pressure applied at infinity. Due to the symmetry, horizontal displacements are zero at lateral walls orthogonal to OX, which corresponds to the infinite periodic set-up. Vertical walls on the inlet and outlet sides are assumed to remain flat. Conditions (6) will be explained in the following section. Note that the mechanical contact pressure p is defined in contact regions Ω c , whereas the fluid hydrostatic pressure is a complementary field which is defined in the fluid zone Ω f , thus for every point p • p f = 0 with an exception of the contact-fluid interface (contact edge), ∂Ω f ∩ ∂Ω c , at which p = p f . Here, the condition on the contact pressure pp f ≥ 0 is valid only on the contact edge shared with the fluid boundary ∂Ω f ∩ ∂Ω c , as p f = 0 in the interior of the contact zone Ω c \ ∂Ω c . Note that such formulation makes possible existence of contact zones with contact pressure smaller than the pressure of surrounding fluid if and only if these zones are separated from the contact-fluid interface by zones with contact pressure higher than the fluid pressure; a suction cup could be given as an example here. Linear isotropic elasticity is considered for the elastic half-space V, so the stress-strain relation is given by the Hooke's law. Finally, we have one unknown vector field in three dimensions, displacements u(x, y, z) in V, and one unknown scalar field in two dimensions, which is the hydrostatic fluid pressure p f (x, y) in Ω f . We assume that full contact is not reached in all sections parallel to OX, so the fluid can always circulate.

In Fig. 1(a) we also sketch the resulting interface pressure, note that in the contact zones it is not uniform along the OY axis, but rather increasing towards the outlet side, while the width of contact patches in sections orthogonal to OY is also increasing towards the outlet, in accordance to the direction of the fluid pressure drop, see also Fig. 1(b).

We study the evolution of the two fields u and p f with the increasing external pressure p ext . In particular, we are interested to know how the contact profile a(y) delimits the contour of the fluid channel and how its depth g(x, y) changes. In addition, we can determine under which conditions the flow is possible through the interface, i.e. what are the critical values of p i , p o , p ext resulting in channel closure.

Wavy profile with pressurized fluid in the interface

Before making an attempt to solve the three-dimensional problem formulated in the previous section, we focus our attention on a simpler, planar contact problem with a pressurized fluid in the interface. Understanding of this problem will be helpful for the derivation of the approximated solution for the full problem, which is presented in the following section.

An elastic solid with a wavy surface (Fig. 2, Eq. ( 1)) is brought in contact with a rigid flat in a fluid environment, which is retained under a constant pressure p f . A plane strain problem is considered. The solid mechanical contact problem was solved for this configuration by Westergaard without fluid pressure (p f = 0) [START_REF] Westergaard | Bearing pressures and cracks[END_REF][START_REF] Johnson | The contact of elastic regular wavy surfaces[END_REF] under the assumption of the infinitesimally small slope of the wavy profile; the pressure distribution in contact region was found to be

p W (x , a ) = 2 p(a ) cos(x )
sin 2 (a ) sin 2 (a )sin 2 (x ), [START_REF] Bažant | Why fracking works[END_REF] where x = πx/λ, a = πa/λ, and a is the half-width of the contact zone; the mean contact pressure reads as

p(a ) = πE * ∆ λ sin 2 (a ), (8)
i.e. p(a ) = a 0 p W (x , a ) dx , and E * is the effective elastic modulus defined by

1 E * = 1 -ν 2 1 E 1 + 1 -ν 2 2 E 2 ,
where E i , ν i are the Young's moduli and Poisson's ratios of the two contacting solids i = 1, 2, respectively. Associated effective displacements (taken with a negative sign) in the contact interface are given by

u z (x , a ) =            -∆ cos(2x ) + C, in contact cos(x ) > cos(a ) -∆ [cos(2x ) + 2 sin(x )h(x , a ) -2 sin 2 (a ) × × ln sin(x )+h(x ,a ) sin(a ) + C, out of contact cos(x ) ≤ cos(a ), (9)
where h(x , a ) = sin 2 (x )sin 2 (a ). The solution of the contact problem for the same configuration in the pressurized environment was found by Kuznetsov in [START_REF] Kuznetsov | Effect of fluid lubricant on the contact characteristics of rough elastic bodies in compression[END_REF]. If we assume that the fluid pressure acts only vertically 3 and that the profile slope is infinitesimal, the stress state in the contact interface in the presence of the additional fluid pressure, applied outside the contact patches, can be considered as the superposition of the stress state corresponding to the same contact area, but without influence of the fluid, i.e. the Westergaards solution [START_REF] Bažant | Why fracking works[END_REF], and a uniform field of the fluid pressure p f :

p(x , a ) = p f + p W (x , a ). (10)
A detailed rigorous analysis of the trapped fluid in the contact interface without the assumption of infinitesimal slopes and with the fluid pressure acting normally to the surface was carried out in [START_REF] Shvarts | Trapped fluid in contact interface[END_REF], but here the simplified vision ( 10) is sufficient to investigate the strongly coupled fluid flow along the wavy channel brought in contact.

In the classical Hertzian contact the pressure decreases to zero towards contact edges, but in a pressurized environment such a situation is impossible as the contact would be open by the environmental pressure. So the pressure p f represents an offset which can be complemented by the mechanical pressure rising in contact. Since a constant p f does not change the shape of the surface in infinitesimal slope assumption, displacements obtained for pressures [START_REF] Bažant | Why fracking works[END_REF] and [START_REF] Brown | Applicability of the reynolds equation for modeling fluid flow between rough surfaces[END_REF] differ only by a constant, thus Eq. ( 10) satisfies the contact conditions, which can be formulated for the coupled problem in the following way:

p -p f ≥ 0, g ≥ 0, (p -p f )g = 0, ( 11
)
where p is the contact pressure, g is the gap, and p f is the pressure of the environment. That is the form which appears in the main system of equations to be solved [START_REF] Bao-Tong | Overview of fatigue behavior in copper single crystals -I. Surface morphology and stage i crack initiation sites for tests at constant strain amplitude[END_REF].

The relation between the mean gap and the contact force almost does not change compared to the unpressurized case [START_REF] Johnson | The contact of elastic regular wavy surfaces[END_REF][START_REF] Kuznetsov | Effect of fluid lubricant on the contact characteristics of rough elastic bodies in compression[END_REF], only a force offset p f λ is added for every period. The mean gap can be computed by integrating g(x) = z(x) + u z (x) in the non-contact region giving

ḡ = ∆ 1 - p p * 1 -ln p p * , for p ∈ [0, p * ] ( 12
)
where p * = πE * ∆/λ is the pressure needed to bring the wavy surface in full contact in the absence of fluid pressure and the mean contact pressure in absence of fluid p is given by Eq. ( 8). In Fig. 2(b) the full normalized pressure p/p * = ( p + p f )/p * is plotted with respect to the mean gap ḡ evolution, where p + p f is the full pressure applied to the system. It is naturally assumed that the fluid can be squeezed-out from the contact interface.

Based on this planar solution the following preliminary conclusion can be drawn for the three-dimensional problem. If the hydrostatic pressure changes only weakly along the channel, i.e. p i ≈ p o , then the contact lines would remain almost parallel to the axis OY and the derivative of the gap with respect to y may be neglected. Then the hydrostatic pressure will be an affine function of the coordinate y, i.e. p f = p i + (p op i )y/L. The flux, which would have a non-zero component only along y axis, can be readily found as q y (x) = -g 3 (x)(p o -p i )/(12µL); note that it depends only on the x-coordinate. Thus, naturally for the situation p i ≈ p o the channel would be sealed at p ext ≈ p * + p o .

Approximate solution

To solve approximately the coupled problem formulated in Eqs. ( 2)-( 6) we need to make several strong assumptions: we assume (i) that in every section y = const, the pressurized Westargaard-Kuznetsov solution (10) is satisfied for a = a(y) and p f = p f (y). However, it is clear that it should imply that ∂p f /∂x = 0, which could seem to exclude the channel narrowing. But since in the following, the problem will be reduced to a one-dimensional flow along OY axis, this assumption (i) is not contradictory: the fluid pressure can be considered as the mean pressure in the section

p f (y) = 1 λ -2a λ-a a p f (x, y) dx.
We also assume (ii), which is the strongest and the least realistic assumption, that in every section

p ext = p + p f = const, (13)
for that we require that p f ≤ p ext in Ω, which is equivalent to require that p i ≤ p ext . Another simplification would be (iii) to reduce the two dimensional Reynolds equation to a onedimensional equation for the average gap [START_REF] Carbone | Adhesion and friction of an elastic half-space in contact with a slightly wavy rigid surface[END_REF], which implies that the hydrostatic pressure is independent of the x-coordinate: p f = p f (y). Under these assumptions Eq. ( 2)-(3) can be rewritten as ḡ3

p f = C 1 , p f (0) = p i , p f (L) = p o , (14)
where the prime sign denotes partial derivative with respect to y, and C 1 is the integration constant. The condition of the zero flow at x = 0, λ (see Eq. ( 2)) is automatically satisfied as p f is assumed not to depend on x.

Now we can use the relation between the mean gap and the pressure (12) through the relation [START_REF] Coussy | Poromechanics[END_REF], which being substituted in (14) yields:

-p * ∆ 3 C 1 ρ(y) ρ i 1 -ρ 1 -ln ρ 3 dρ = y, (15)
where ρ(y) = (p extp f (y))/p * , and ρ i = ρ(0) = (p extp i )/p * . From [START_REF] Dapp | Contact mechanics of and reynolds flow through saddle points: On the coalescence of contact patches and the leakage rate through near-critical constrictions[END_REF] it follows that ρ(y) should be a monotonically increasing function of y in the range y ∈ [0; L]. The boundary conditions now read as:

ρ(0) = ρ i = (p ext -p i )/p * , ρ(L) = ρ o = (p ext -p o )/p * .
The integral I(ρ) = 1 -ρ 1ln ρ 3 dρ from [START_REF] Dapp | Contact mechanics of and reynolds flow through saddle points: On the coalescence of contact patches and the leakage rate through near-critical constrictions[END_REF] with zero integration constant is evaluated as:

I(ρ) =ρ -α 1 ρ 2 + α 2 ρ 3 -α 3 ρ 4 + β 0 ρ 2 (1 -β 1 ρ + β 2 ρ 2 ) ln(ρ) + ρ 3 (1 -γρ) ln 2 (ρ) + ρ 4 4 ln 3 (ρ),
where α 1 = 9/4, α 2 = 17/9, α 3 = 71/128, β 0 = 3/2, β 1 = 16/9, β 2 = 13/16, γ = 15/16. The solution cannot be provided in the form p f (y), but rather y(p f ), which reads as

y = -p * ∆ 3 C 1 (I(ρ) + C 2 ),
where the integration constants can be found through boundary conditions: C 2 = -I(ρ i ) and C 1 = -p * ∆ 3 I(ρ o ) -I(ρ i ) /L. The final approximate solution, which determines the average fluid pressure distribution along the channel coordinate y is given below:

y L = I(ρ) -I(ρ i ) I(ρ o ) -I(ρ i ) . ( 16
)
Resulting curves for the variation of hydrostatic pressure, mean gap and contact half-width along the channel are depicted in Fig. 3, 4(a), 4(b), respectively. This approximate result is able to capture the non-linear hydrostatic pressure distribution along the channel, to account for the induced deformation of the solid and thus to obtain the narrowing of the channel. Note that the obtained fluid pressure distribution becomes a concave function and resembles the pressure distribution of a compressible fluid.

With the derived analytical solution we may calculate the fluid flux in the y-direction as

q y (x, y) = - g 3 (x, y) 12µ dp(y) dy , (17)
where the gap g(x, y) is obtained in each section y = const using the Westergaard's solution [START_REF] Bower | The influence of crack face friction and trapped fluid on surface initiated rolling contact fatigue cracks[END_REF] corresponding to the mean pressure p = p extp f (y), and the derivative dp/dy is calculated using [START_REF] Dapp | Fluid leakage near the percolation threshold[END_REF] as

dp dy = - 1 p * L I(ρ o ) -I(ρ i ) dI(ρ) dρ -1 . ( 18
)
This result will be used in the following for computation of the mean flux and the effective transmissivity of the interface and comparison with the numerical solution.

Numerical coupling schemes for fluid-solid interaction

Two different approaches for numerical solution of fluid-structure interaction problems (including the fluid transport through the contact interface) are generally distinguished: the first is the partitioned approach, under which solvers for solid mechanics problem and for the fluid transport are separated and work sequentially, thus the data exchange between them must be established; the second is the monolithic approach, which operates with a single combined solver for both physics, so that solutions of the solid and fluid problems are obtained simultaneously.

In this study we use the monolithic approach, some details of its implementation in the finite element framework are presented below. To satisfy the contact constraints, we add contact elements on the potential contact zone, here we utilize the mortar approach and the augmented Lagrangian method to fulfill the constraints [START_REF] Alart | A mixed formulation for frictional contact problem prone to newton like solution methods[END_REF][START_REF] Yastrebov | Numerical methods in contact mechanics[END_REF][START_REF] Wriggers | [END_REF][START_REF] Puso | A mortar segment-to-segment contact method for large deformation solid mechanics[END_REF]. To solve the Reynolds problem, fluid pressure values are added as degrees of freedom to the surface nodes of the structural mesh, and the surface elements for solving the fluid problem are appended to this mesh. Finally, to enforce the Neumann boundary conditions of the fluid pressure in the out-of contact zones of the surface, we use the mortar method again [START_REF] Yang | A mortar-finite element approach to lubricated contact problems[END_REF]. Therefore, the equations governing contact constraints, fluid flow and fluid/structure interaction are embedded into a single global system of equations, which is solved at each iteration of the Newton-Raphson method, required due to intrinsic nonlinearity of the coupled problem. When convergence criteria are satisfied for displacement and pressure degrees of freedom, the solution of the coupled problem is obtained for given boundary conditions.

Note that contrary to many multiphysical coupled problems, in which different physics require significantly different temporal and spatial discretizations, the particular coupling considered here deals with stationary equations for the fluid and the solid, which can be resolved numerically on almost coinciding surface meshes. The mesh used for the fluid is obtained by projecting the surface of the solid mesh on the rigid flat, with which the solid comes in contact.

The coupling algorithm was implemented in the in-house finite element software Z-set [START_REF] Besson | Large scale object-oriented finite element code design[END_REF][START_REF]Finite element suite Z-set[END_REF]. Due to the reflection symmetry of the geometry and loads, only the half wavelength is simulated using structured finite element mesh of hexahedral linear elements, with 128×128 square-shaped faces on the surface, the mesh gradually coarsens with the depth. The mesh is depicted in Fig. 5, it contains approximately 109 000 nodes and approximately 98 000 elements with 8 integration points per element. At the bottom surface z = -H, the displacement vector is prescribed as u x = u y = 0, u z = kt, where k is a load factor, t is the time. Geometrical parameters of the problem are the following L = 1 cm, λ = 2 cm, ∆ = 0.2 mm, H = 1.4 cm, Young's modulus E = 2 GPa and Poisson ratio ν = 0.3, thus E * ≈ 2.2 GPa.

Results

We carried out several simulations for three different fixed inlet fluid pressures p i = 2, 10, 50 MPa, respectively, and for outlet pressures p o /p i = {0; 0.25; 0.5; 0.75}. The solid is gradually brought in contact by applying vertical displacement on its bottom side. We present the detailed results of numerical simulations in Fig. 6 and Fig. 7: the distribution of the fluid pressure and the contact pressure, as well as the gap and the fluid flux for the case p i = 50 MPa, p o = 0 and for two particular load steps: p ext /p * = 0.48, 0.8, respectively. On the initial stage of loading no contact occurs and the load is supported by the fluid solely, which flows along the entire channel (this classical situation is not presented in the figures). For higher loads, the solid comes in contact with a rigid flat, the contact starts from the outlet zone (see Fig. 6). With the increasing load the contact zone spreads out and at a certain load reaches the inlet zone (see Fig. 7); starting from this moment the approximate solution ( 16) becomes applicable. The gap, being bigger at the inlet due to higher inlet pressure, narrows towards the outlet region. The flow is localized within the trough of the wavy profile, and the flux intensifies towards the outlet. Due to the narrowing of the channel along the flow direction, the current lines converge towards the outlet, it results in a small but still distinguishable fluid pressure gradient in the OX direction, which was not explicitly taken into account in the derivation of the approximate solution.

In Fig. 8 numerical results for the variation along the channel of mean fluid pressure, mean gap and the contact width, as well as the contact and fluid pressure profiles in the section y = L/2 are compared with approximate solution from Section 4. A rather good agreement is obtained in the range of validity of the approximate result: the inlet pressure and external pressure satisfy the following conditions p i ≤ p ext , p ext < p * , and they are chosen such that the contact zone reaches the inlet, see [START_REF] Coussy | Poromechanics[END_REF] and the discussion in Section 4. These limitations are quite strong and in reality ensure only a limited range of validity of the approximate solution.

In Fig. 9(a) the numerically computed evolution of the length of the contact zone on the inlet and the outlet sides with the increasing pressure is depicted for three different inlet pressures p i /p * = 0.03, 0.14, 0.72 and p o = 0. These results are compared with the analytical solution, which is valid if the external pressure is in the interval p i < p ext < p * . In accordance with the assumptions of the approximate solution, the evolution of the length of the contact zone on the outlet side is independent of the inlet fluid pressure. However, in the numerical results curves both for inlet and outlet sides shift into the region of higher external pressures with the increasing inlet pressure. Note also that in the numerical results the outlet contact length grows faster than the inlet one. Finally, the strongly coupled numerical model shows that the higher is the inlet pressure p i , the higher external load is needed to completely seal the channel.

In Fig. 9(b) the effective transmissivity K eff of the wavy channel brought in contact is plotted: it is defined as

K eff = - 12µQL ∆ 3 (p o -p i ), ( 19
)
where Q is the mean flux over the area λ × L, i.e.

Q = 1 λL λ 0 L 0 q y dxdy. ( 20
)
For computation of the local flux q y (x, y) we used the results of numerical simulations as well as the approximate solution ( 16)- [START_REF] Fischer | Stick slip sliding behaviour at the base of a glacier[END_REF]. We considered the same three cases as before with different inlet pressures. For each case we highlight the corresponding external pressure necessary for the contact to appear at the inlet and outlet sides. As soon as the contact appears on both sides, the evolution of the transmissivity becomes exponential with respect to the external pressure normalized by p * with the exponent ≈ -8 in all three cases. We note that this coefficient is lower than the one observed in simulations of the interface transmissivity for the surfaces with representative random roughness, where it was reported to be of order ≈ -12 [START_REF] Dapp | Fluid leakage near the percolation threshold[END_REF]. Closer to the complete sealing (the percolation limit) the transmissivity decays faster and can be very accurately described by a power law with respect to the difference between the critical external pressure p c , necessary to seal the channel, and p ext , i.e. K eff ∼ (p cp ext ) γ . Note that for our results γ was estimated as 6 ± 0.5, while in accurate but one-way coupled studies of the percolation limit of bi-sinusoidal surfaces it was found to be equal to 3.45 [START_REF] Dapp | Contact mechanics of and reynolds flow through saddle points: On the coalescence of contact patches and the leakage rate through near-critical constrictions[END_REF][START_REF] Dapp | Fluid leakage near the percolation threshold[END_REF].

The transmissivity for the lowest inlet pressure p i /p * = 0.03 almost coincides with that of the one-way coupling analysis and is well described by Kuznetsov's analytical solution. For higher inlet pressures we obtain significantly bigger transmissivity. In the region of the exponential decay in case of two-way coupling we have ≈ 1.6 times higher transmissivity than in the case of the one-way coupling for p i /p * = 0.14 and ≈ 32 times higher for p i /p * = 0.72. Note that the transmissivity curves based on the analytical approximation ( 17)-( 18) are in a very good agreement with the numerical ones in the range of the validity of the former. However, the analytical result cannot be used to study the flow near the percolation since the pressure needed to seal the channel given by the analytical solution (i.e. p ext = p * + p o ) strongly underestimates the real one, which can be accurately studied using the numerical approach.

In the inset of Fig. 9(b) we plot the effective transmissivity with respect to the real contact area fraction, curves coincide for three different cases in the beginning of loading, while the complete sealing occurs at different values of the real area fraction A/A 0 = 80% -90%. Note that in one-way coupled studies of the percolation limit of the randomly rough surfaces in contact, the complete sealing was found to correspond to ≈ 42% of the real contact area [START_REF] Dapp | Self-affine elastic contacts: Percolation and leakage[END_REF].

Finally, we evaluated the critical external pressure necessary to seal the channel, i.e. to prevent the fluid flow across the interface. The results are presented in Fig. 10(a). The relationship between the inlet pressure and the critical external pressure is found to be close to linear, the results of the least squares fitting are presented in Fig. 10(a external pressure p c needed to seal the channel can be approximately found as:

p c ≈ p * + 0.8p i , (21)
where 0.8 is a fit parameter. We performed additional simulations with different outlet pressures for a given inlet one. The results are presented in Fig. 10(b), where we plot evolution of the mean flux Q (20) with the increasing external pressure. In the beginning of loading the mean flux is lower in case of a smaller pressure drop, however, under the increasing external load curves converge, and complete sealing of the channel occurs at the same value of the critical external pressure, which therefore is determined only by the inlet pressure. The explanation comes from the fact that close to the complete sealing of the channel, the fluid pressure drop occurs in the vicinity of the outlet, while in the remaining channel the pressure is close to the inlet one. Therefore, remarkably, the load-carrying capacity of the fluid in the interface close to percolation is defined only by the inlet pressure. In addition, in Fig. 10(a) we highlight the external pressure at which the contact appears on the inlet side, this corresponds to the onset of validity of the approximate solution. The end of its validity corresponds to p ext /p * = 1. The narrow validity range of our approximate solution can be improved by dropping the assumption (ii) in Section 4 and by including strong elastic interaction between y-sections, which would lead to a much more complicated analysis and is not addressed here.

Concluding remarks

We presented a theoretical study of the pressure driven creeping flow in contact interface formed between solids with wavy surfaces. We considered the case when the fluid flows across the wavy section in channels delimited by mechanical contact zones. This problem is relevant for certain applications of thin fluid flow in contact interfaces, including sealing, hydrogeology and biological systems. A strong (two-way) coupling between fluid flow and deformation of the solid was assumed, which is crucial for applications in which the fluid pressure is comparable with the mean contact pressure, for example, for the soft matter or biological tissue.

We derived an approximate analytical solution based on the Westergaard-Kuznetsov solution and a one-dimensional formulation of Reynolds equation. This solution describes both the solid deformation and the fluid pressure distribution in the strongly coupled case. A finite-element monolithically coupled framework for fluid and solid equations was also used to solve this non-linear multi-field problem and to prove the validity of the approximate solution. Despite a rather limited interval of loads within which the latter is applicable, it can provide a useful first-order approximation for the analysis of transmissivity of contact interfaces. At the same time, numerical results showed that in a wide range of the external loads up to the complete sealing of the channel, the transmissivity of the interface can be described by a power law, which has already been reported in the studies of contact interfaces having representative and model roughness.

Both numerical and analytical results, which take into account two-way coupling, showed that the interface transmissivity is significantly higher than this predicted by the one-way coupling if the fluid pressure is high enough. An additional result of this study is the affine dependence of the external critical pressure which seals the channel on the inlet fluid pressure: this relation may be shown useful in sealing applications as well as in soft porous or cracked media, in zones where the flow can be described by Reynolds equation. Remarkably, this critical pressure was found to be independent of the outlet pressure.

Figure 2 :

 2 Figure 2: (a) -Contact between an elastic wavy profile and a rigid flat with a pressurized fluid in the interface, (b) -variation of the mean pressure on the surface with the normalized gap.

Figure 3 :Figure 4 :

 34 Figure 3: Approximate solution (16) for the mean hydrostatic fluid pressure distribution along the wavy channel in contact with a rigid flat for p o = p i /2 and p o = 0 (a) p ext /p * = 0.05, (b) p ext /p * = 0.1, (c) p ext /p * = 0.2; the fluid pressure is normalized by the inlet pressure p f /p i and the coordinate is normalized by the channel length y/L.

Figure 5 :

 5 Figure 5: Finite element mesh with 128 × 128 elements on the contact surface, which was used to solve the coupled solid-fluid problem.

Figure 6 :

 6 Figure 6: Results of the numerical simulation: (a) normalized fluid pressure p f /p i , (b) normalized contact pressure |σ n |/|σ n | max , (c) normalized fluid flux intensity |q|/|q max |, (d) normalized gap g/∆; p ext /p * = 0.48, p i /p * = 0.72, which corresponds to p i /p ext = 1.5.

Figure 7 :

 7 Figure 7: Results of the numerical simulation: (a) normalized fluid pressure p f /p i , (b) normalized contact pressure |σ n |/|σ n | max , (c) normalized fluid flux intensity |q|/|q max |, (d) normalized gap g/∆; p ext /p * = 0.8, p i /p * = 0.72, which corresponds to p i /p ext = 0.9.

Figure 8 :

 8 Figure 8: Comparison of results of numerical simulations with the approximate analytical solutions for (a) normalized contact region half-width a/λ, (b) normalized mean gap ḡ/∆ and (c) normalized fluid pressure p f /p i evolution along the channel, (d) normalized contact and fluid pressure profiles p/E * in the section y = L/2: p = |σ n | in Ω c , p = p f in Ω f and p = |σ n | = p f at ∂Ω c ∩ ∂Ω f .

  ). So the critical

	(a)		(b)
			One-way coupling
			Two-way coupling:
			p i /p ext Numerical	Analytical
			0.04
			0.18
			0.90
		One-way coupling
		Two-way coupling:
		p i /p ext Numerical	Analytical
		0.04
		0.18
		0.90
	(c)		(d)
			One-way coupling
			Two-way coupling:
			p i /p ext Numerical	Analytical
			0.04
			0.18
			0.90
	One-way coupling
	Two-way coupling:
	p i /p ext Numerical	Analytical
	0.04
	0.18
	0.90

By free volume here we mean the separation field between contacting surfaces.

Note that all the discussions are valid not only for an elastic solid with a wavy surface in contact with a rigid flat, but for two elastic solids with the effective wavy roughness given by z = z 1z 2 + c, where z 1 , z 2 determine surface geometries of the two contacting solids, and c is an arbitrary constant. However, for simplicity hereinafter we will assume that an elastic wavy surface is brought in contact against a rigid flat.

In[START_REF] Shvarts | Trapped fluid in contact interface[END_REF], this assumption was shown to be too prohibitive for certain applications even if the surface slope is assumed infinitesimal.

Acknowledgments

The authors are grateful to Julien Vignollet for useful comments to the paper and to Andreas Almqvist and Francesc Pérez Ràfols for enriching discussions. The financial support of Safran Group and MINES ParisTech (Thèse-Open) is acknowledged.

 17)-( 18) are shown as thick dashed curves, numerical results are presented using markers, while full curves are fittings of power law K eff ∼ (p cp ext ) γ , where γ = 6 ± 0.5. "Cross" and "star" markers are used to highlight the external pressure necessary for the contact to appear at the outlet and inlet sides, respectively. Inset in (b) shows the evolution of the effective transmissivity with respect to the ratio of the real contact area to the apparent one. zone III.a corresponds to the combination of loads for which the approximate analytical solution remains valid; IV -the contact is completely sealed, no flow passes through, the fluid is under the inlet hydrostatic pressure and ensures some load-bearing capacity in the non-contact region. (b) Evolution of the mean flux across the interface with the increasing external pressure for p i /p * = 0.14 and four cases of the outlet pressure p o /p i = 0, 0.25, 0.5, 0.75, showing that the critical pressure, necessary to seal the channel, is independent of the outlet pressure. Markers represent numerical results, while full curves are fittings to the power law discussed above.