Pierre Monneyron

Marie-Hélène Manero

Jean-Noel Foussard

Measurement and Modeling of Single-and Multi-Component Adsorption Equilibria of VOC on High-Silica Zeolites

published or not. The documents may come

OATAO is an open access repository that collects the work of Toulouse researchers and makes it freely available over the web where possible Any correspondence concerning this service should be sent to the repository administrator: tech-oatao@listes-diff.inp-toulouse.fr evaporation, leading to an initial concentration of about 0.5 mmol L -1 . Stirring times required to reach equilibrium were respectively 2 h for Sil Z and 1 h for Fau Y. Once equilibrium was attained, a gas sample was taken and analyzed by gas chromatography (HP 5890 II). Using the conventional assumption that the clay binderswith a low specific surface area as compared to the zeolitessdoes not take part in the adsorption mechanism (5), adsorbed quantities were normalized for pure zeolite material.

In case of binary mixtures, adsorption procedure was kept unchanged except that co-adsorption was performed using three different initial compositionssbeing equimolar (Z1 ) Z2 ) 0.5) and rich in components 1 and 2 (respectively, Z1 ) 0.85 and Z1 ) 0.15). The selectivity of adsorbent (η1/2) was either directly represented in an X-Y isobaric equilibrium diagram or determined by calculation using eq 1 where Xi and Yi are respectively the adsorbed and gas-phase molar fractions of the ith component and are plotted versus total pressure:

Reproducibility and repeatability of experimental data were checked by three different manipulators. Experimental errors were equal to 3% for equilibrium concentration and ranging from 3 up to 7% for adsorption capacity with increasing pressure.

Theoretical Models

Single-Component Isotherms Correlation. The experimental adsorption isotherms were fitted using Freundlich, Langmuir, Toth (which was intensively used to describe adsorption isotherms on zeolites ( 22)), and Dubinin-Astakhov equations. The models express the adsorbed quantity, Q (mol kg -1 ), as a function of gas-phase pressure, p (Pa), as follows:

Freundlich Equation.

Langmuir Equation Toth Equation.

Dubinin-Astakhov Equation. Dubinin's approach is based on the earlier potential theory of Polanyi and on the concept of micropore filling. The Dubinin-Astakhov (D-A) approach is an extension of Dubinin-Radushkevich equation and represents the isotherm in terms of fractional filling (θ ) W/W0) of the micropore volume with W0 being the maximum adsorption capacity.

In eq 5, Ap is the Polanyi adsorption potential defined by eq 6, and m is the D-A equation adjustment parameter, ranging classically from 1.5 to 3 for activated carbon and from 2 to 4 or 5 in the case of zeolites [START_REF] Garrot | [END_REF]:

Optimization of the models parameters was carried out using the Simplex algorithm, by minimizing the objective function, F(X), where X is the set of adjustable parameters and n is the number of experimental points: Binary Mixtures Modeling with IAST and RAST. Considering the adsorbed phase as a liquid, the IAST is based on the fundamental hypothesis that the multi-component adsorbate phase has the same spreading pressure (π) as that of the ideal single solution of the ith component, the spreading pressure being the difference between the interfacial tensions of the naked solid and that of the surface covered with adsorbate. Algorithm of adsorption capacities determination requires an accurate description of the pure component isotherms, but integration of the chosen model in the Gibbs equation may be more or less direct, depending notably on its number of adjustable parameters. This is clearly highlighted by the spreading pressure expressions obtained either with Langmuir (eq 8) or Toth (eq 9) models in which {0 e θ < 1} and {0 < t < 1}. The improvement of predictions accuracy when changing from 2 to 3 parameters is examined in balance of the very different CPU calculation times induced, which are decisive in dynamic breakthrough curve modeling:

When thermodynamic consistency was satisfied (24), the AST was extended to real AST. Considering deviation from ideality, identified by nonsymmetrical X-Y curves in experimental isobaric equilibrium diagrams, Costa et al. (20) introduced liquid-phase activity coefficients (γi) when establishing chemical potential equality between phases at equilibrium. This extension enables the AST to describe azeotrope-like crossovers. Basic equations of the two models are summarized below using the Wilson model for activity coefficients with two adjustable parameters Λij (i * j). The calculation procedure in both cases was described in the original publications as well as in a previous article for IAST (25).

η 1/2 ) X 1 Y 2 X 2 Y 1 (1) Q ) K f p 1/n (2) Q ) Q m K l p 1 + K l p (3) Q ) Q m K t p [1 + (K t p) t ] 1/t (4) θ ) exp [ -( A p E ) m ] (5) A p ) RT ln [ p s p ] (6)
F(X) ) ∑ k n [Q exp (p k ) -Q calc (p k )] 2 (7) πA RT ) Q m ln(1 + K l C) (8) πA RT ) Q m [ θ -θ t ln(1 -θ t ) -∑ j ∞ θ jt + 1 jt(jt + 1) ] (9

  This is an author's version published in: http://oatao.univ-toulouse.fr/23283 To cite this version: Monneyron, Pierre and Manero, Marie-Hélène and Foussard, Jean-Noel Measurement and Modeling of Single-and Multi-Component Adsorption Equilibria of VOC on High-Silica Zeolites. (2003) Environmental Science & Technology, 37 (11). 2410-2414. ISSN 0013-936X Official URL : https://doi.org/10.1021/es026206c

TABLE 1 . Main Characteristics of Zeolites

 1

		faujasite	ZSM-5
	type (symbol)	(Fau Y)	(Sil Z)
	crystalline framework	R-cages interconnected
			channels
	pore internal diameter (Å)	13	(5.7 × 5.1) and 5.4
	clay binder content (wt %)	25	20
	SiO 2/Al2O3 (mol mol -1 )	360	1880
	BET surface area (m 2 g -1 )	650	300
	active porous volume (cm 3 g -1 ) 0.24	0.18

TABLE 2 . Main Properties of Adsorbates

 2

		methyl ethyl
	property	ketone	toluene 1,4-dioxane
	symbol	MEK	TOL	DIO
	formula molar volume 25 °C (cm 3 mol -1 ) 90.1 C 4H8O	C7H8 107.5 85.8 C4H8O2
	kinetic diameter (Å) vapor pressure at 25 °C (kPa)	5.2 12.60	5.8 3.79	5.3 4.95
	dipolar moment (D)	2.78	0.375 0


Begum Genc
email: begum.genc@insight-centre.org

Mohamed Siala
email: mohamed.siala@laas.fr

Gilles Simonin
email: gilles.simonin@imt-atlantique.fr

Barry O'sullivan
email: barry.osullivan@insight-centre.org

An Approach to Robustness in the Stable Roommates Problem and its Comparison with the Stable Marriage Problem

Recently a robustness notion for matching problems based on the concept of a (a, b)-supermatch is proposed for the Stable Marriage problem (SM). In this paper we extend this notion to another matching problem, namely the Stable Roommates problem (SR). We define a polynomial-time procedure based on the concept of reduced rotation poset to verify if a stable matching is a (1, b)supermatch. Then, we adapt a local search and a genetic local search procedure to find the (1, b)-supermatch that minimises b in a given SR instance. Finally, we compare the two models and also create different SM and SR instances to present empirical results on the robustness of these instances.

Introduction

Robustness to change is an important property that has a variety of definitions in different settings [START_REF] Lussier | On fault tolerance and robustness in autonomous systems[END_REF]. There exist many robustness notions within the context of matching problems. These robustness notions mostly focus on handling uncertainty and erroneous data in the input [START_REF] Drummond | Elicitation and approximately stable matching with partial preferences[END_REF][START_REF] Ashlagi | Dating strategies are not obvious[END_REF][START_REF] Aziz | Stable matching with uncertain linear preferences[END_REF][START_REF] Jacobovic | Perturbation robust stable matching[END_REF]. Genc et al. introduced a novel notion of robustness for the Stable Marriage problem (SM) where the robustness of a solution refers to its capability to be repaired at a small bounded cost in case of an unforeseen event [START_REF] Genc | Finding robust solutions to stable marriage[END_REF]. The notion of (a, b)-supermatches differs from the other robustness notions in this context since it specifies a degree of repairability. This property is often referred as fault-tolerance. The (a, b)-supermatch concept defines the notion of robustness for matching problems by using the fault-tolerance framework [START_REF] Hebrard | Robust solutions for constraint satisfaction and optimisation under uncertainty[END_REF][START_REF] Holland | Super solutions for combinatorial auctions[END_REF].

The SM is defined by a set of men and a set of women, each of which has a set of preferences over people of opposite sex. The task is to find a (monogamous) matching between men and women that is stable. A matching is said to be stable if there are no two pairs that are not matched to each other, but they prefer being together than being with their current partners. The robust variant of the problem is called Robust Stable Marriage (RSM) [START_REF] Genc | Finding robust solutions to stable marriage[END_REF], in which the robustness of a stable matching is measured by the minimum number of changes required to obtain another stable matching in the case of break-up of some pairs. If a pair appears in all the stable matchings, the pair is said to be fixed, otherwise, non-fixed. An (a, b)-supermatch is a stable matching such that if any a non-fixed agents (men/women) break-up, it is possible to find another stable matching by changing the partners of those a agents and also changing the partners of at most b others. The previous work on the RSM includes the proposal of the problem, a complexity study, a polynomial-time verification procedure for a given (1, b)-supermatch, and three different models (constraint programming, genetic algorithm, local search) to find the (1, b)-supermatch that minimises b for a given SM instance [START_REF] Genc | Complexity study for the robust stable marriage problem[END_REF][START_REF] Genc | Finding robust solutions to stable marriage[END_REF]. We investigate in this paper this robustness concept further on a generalised version of the SM, namely the Stable Roommates problem (SR). The Stable Roommates problem is a one-sided generalisation of SM, where any two agents regardless of their gender can be matched. We define the Robust Stable Roommates problem (RSR) analogous to the RSM. To the best of our knowledge, there is no previous research on finding the (a, b)-supermatches of the SR.

The motivation behind studying RSR is due to the large applicability of SR and the importance to handle the dynamism of the real world. Take the example of P2P networks where peers (computers for instance) are connected to each other for file sharing purposes [START_REF] Lebedev | On using matching theory to understand p2p network design[END_REF]. Each peer has a preference list towards the other peers and a matching that respect stability is required. However, as the network evolves during time, peers continuously seek new partners. That is, if a peer that provides the file loses the connection, an alternative peer is needed for downloading a file. In this situation, we have have to maintain stability with possibly the minimum changes to the current solution. An (a, b)-supermatch guarantees finding other peers to the broken ones at a small number of additional changes while preserving stability.

The paper is organised as follows: In Section 2, we give a formal background and introduce the robust stable roommates problem. Then, in Section 3, we show that one can verify in polynomial time if a given stable matching (in the SR context) is a (1, b)supermatch. Next, we adapt a local search procedure and a hybrid (genetic local search) model for finding robust solutions in Section 4. Finally, we present in Section 5 our empirical study.

Background and Notation

The Stable Roommates problem (SR) consists of a set of 2×n agents, where each agent has a preference list in which he/she ranks all other agents in strict order of preference. In this context, given a set of people P , a matching corresponds to a partition of P into disjoint pairs (or partners). A matching is stable if it admits no blocking pairs. A pair {p i , p j } blocks a matching if: p i is unassigned or prefers p j to his/her current partner, or p j is unassigned or prefers p i to his/her current partner. The solution to an SR instance is a stable matching. If such a solution does not exist, then the instance is unsolvable. A pair is stable if it appears in some stable matching. If a pair appears in all stable matchings, it is called a fixed pair. If a person p has at least two different partners among all stable matchings, p is said to be non-fixed. We measure the distance between any two stable matchings M, M by the number of different pairs d(M, M ) =| M \ M |. The stable matching M among all the stable matchings of the instance that has the minimum distance to M is said to be the closest stable matching to M.

Irving defines an O(n 2 ) procedure to find a solution to SR or to report if none exists [START_REF] Irving | An efficient algorithm for the "stable roommates" problem[END_REF]. The procedure consists of two phases. Let us first define some notations to describe these phases. A preference table (denoted by T ) is, for a given problem instance, a set of preference lists for which zero or more entries have been deleted. We use T init to denote the initial preference table. During the two phases, some pairs are removed from T init . We denote the preference list of a person p i in a table T by L T (i). Let f T (p i ), s T (p i ), l T (p i ) denote the first, second and last entries of L T (i). The first phase is based on each person proposing to the first available person in their lists starting from T init until every person has made a proposal that has been accepted, i.e. became semi-engaged. If a person p i becomes semi-engaged to p j , all pairs {p j , p k } such that p j prefers p i to p k are deleted from the table. The table obtained after applying the Phase 1 algorithm is called the Phase-1 table and is denoted by T 0 .

The second phase of the algorithm is based on finding and eliminating rotations starting from T 0 . A rotation ρ is a circular list denoted as ρ = (x 0 , y 0 ), (x 1 , y 1 ), . . . , (x r-1 , y r-1 ), where all x i , y j ∈ P . Each rotation has the property that y i = f T (x i ) and y i+1 = s T (x i ) in a table T for all i, 0 ≤ i ≤ r -1, where i+1 is taken modulo r. The set of people {x 0 , . . . , x r-1 } is called the X-set of ρ, denoted by X(ρ). Similarly, {y 0 , . . . , y r-1 } is called the Y-set of ρ, denoted by Y (ρ). Additionally, given a set of rotations R, X(R) = ∪ ρ∈R X(ρ). Similar for the Y-set. The elimination of a rotation ρ from a table T means for each pair {p i , p j }, where p i = x m and p j = y m and (x m , y m ) ∈ ρ, the deletion of {p i , p j } and all pairs {y m , z} such that y m prefers x m-1 to z from T . In this case, ρ is said to be exposed on T and the table after eliminating ρ is denoted by T /ρ. If after Phase 1 or Phase 2, all lists in T contain exactly one entry, then T represents a stable matching. Note that sometimes we use (p i , p j ) and (x m , y m ) interchangeably. The notation (x m , y m ) is used for denoting the position of the pair (p i , p j ) in ρ. Lemma 4.4.1 from [START_REF] Gusfield | The Stable Marriage Problem: Structure and Algorithms[END_REF] states that {p i , p j } is a stable non-fixed pair if and only if (p i , p j ) or (p j , p i ) is in a non-singular rotation.

There are two types of rotations: singular and non-singular. A rotation ρ = (x 0 , y 0 ), (x 1 , y 1 ), . . . , (x r-1 , y r-1 ) is called a non-singular rotation if ρ = (y 1 , x 0 ), (y 2 , x 1 ), . . . , (y 0 , x r-1 ) is also a rotation. In this case, ρ and ρ are called as duals of each other. If a rotation does not have a dual, then it is a singular rotation. We denote by T S the table where all singular rotations are eliminated from T 0 . A rotation ρ is said to precede another rotation ρ (denoted by ρ ≺ ρ) if ρ is eliminated for ρ to become exposed. In this case, we say ρ is a predecessor of ρ and ρ is a successor of ρ . A rotation ρ is an immediate predecessor of ρ, and ρ is an immediate successor of ρ , if ρ ≺ ρ and there does not exist a ρ * such that ρ ≺ ρ * ≺ ρ. All predecessors and successors of a rotation, not necessarily immediate, are denoted by N -(ρ) and N + (ρ). The set of both singular and non-singular rotations under ≺ defines the roommates rotation poset. The set of non-singular rotations under ≺ defines the reduced rotation poset and is denoted by Π = (V, E). We refer to any two rotations as incomparable if none of them precede the other one, comparable otherwise. Let us illustrate these concepts on an SR instance I. We use a sample instance of 10 people from page 180 in [START_REF] Gusfield | The Stable Marriage Problem: Structure and Algorithms[END_REF]. Fig. 1 represents the T S of I. Fig. 2 represents the reduced rotation poset of I, where the pairs involved in the rotations are given next to their corresponding rotations for convenience. pi LT S (i)

1
} S5 = {ρ3, ρ4, ρ5, ρ6, ρ7} M6 = {(1, 2, 8), (3, 6), (5, 7), (9, 10)
} S6 = {ρ3, ρ4, ρ5, ρ6, ρ7} M7 = {(1, 7), (2, 3), (4, 9), (5, 10), (6, 8)
} S7 = {ρ3, ρ4, ρ5, ρ6, ρ7} 7), (2, 8), (3, 6), (4, 9), (5, 10)
A subset of the rotations in Π, containing one of each dual rotations and all their predecessors, is called a complete closed subset, denoted by S. There exists a 1-1 correspondence between the complete closed subsets of Π and the stable matchings of the underlying instance [START_REF] Gusfield | The Stable Marriage Problem: Structure and Algorithms[END_REF]. Any stable matching can be obtained by eliminating one of each dual rotations starting from T S . A rotation ρ is said to eliminate {p i , p j } if there exists a table T such that {p i , p j } ∈ T and {p i , p j } ∈ T /ρ. On the other hand, a rotation ρ is said to produce {p i , p j } if there exists a table T such that |L T (i)| > 1, |L T (j)| > 1, L T /ρ (i) contains only p j , and L T /ρ (j) contains only p i . We use the term flipping ρ from S as the process of removing ρ ∈ S from S and adding its dual ρ to S. A neighbour rotation ρ is ρ ∈ S and either the rotation has no predecessors (N -(ρ) = ∅) or for all predecessors ρ ∈ N -(S), ρ ∈ S. The set N(S) denotes the set of neighbour rotations. The set of all sink nodes of the graph induced by S is referred as the sink rotations of S, denoted as L(S). Table 1 presents all the 7 stable matchings of I given in Fig. 1 and their corresponding complete closed subsets.

Throughout the paper, we denote by M a given stable matching, and its corresponding complete closed subset by S. If there are any subscripts or superscripts for M such as M * i , then they are applied to the corresponding complete closed subset (i.e. S * i ). Lemma 1 and 2 are included here to be used in our proofs later.

Lemma 1 (Lemma 4.1.1 [START_REF] Gusfield | The Stable Marriage Problem: Structure and Algorithms[END_REF]). Given an instance of the stable marriage problem involving n men and n women, there is an instance (in fact there are many instances) of the stable roommates problem involving those 2n persons such that the stable roommates matchings are precisely the stable matchings for the original SM instance.

Lemma 2 (Lemma 4.3.7 [START_REF] Gusfield | The Stable Marriage Problem: Structure and Algorithms[END_REF]). If ρ, σ are non-singular and π is a singular rotation, then:

(1) ρ ≺ ρ; (2) ρ ≺ σ ⇐⇒ σ ≺ ρ; (3) τ ≺ π =⇒ τ is singular.
Robust Stable Roommates: We refer the problem of finding an (a, b)-supermatch to a given SR instance as the Robust Stable Roommates problem (RSR). A stable matching of an RSR instance is called an (a, b)-supermatch if any a non-fixed pairs do not want to be partners anymore (i.e. leave the stable matching), it is possible to find another stable matching by changing the partners of the people involved in those a pairs and at most b other pairs. The intractability result of the RSM is lifted to the RSR as the SR is a generalisation of the SM.

Theorem 1. RSR is N P-hard.
Proof. The proof is straightforward as it is possible to create an SR instance I SR from any given SM instance I SM with the exact same stable matchings in polynomial-time by padding every other person of the same sex to the preference list of each person (see Lemma 1). Every (a, b)-supermatch in the I SM is also an (a, b)-supermatch in the I SR and vice versa. Hence, RSR is N P-hard because RSM is N P-hard [START_REF] Genc | Complexity study for the robust stable marriage problem[END_REF].

Verification of (1, b)-supermatch in Polynomial Time

We prove in this section that checking if a stable matching M is a (1, b)-supermatch can be done in polynomial time. Indeed, we show in Theorem 3 how to construct the closet matching to M if any non-fixed pair in M wants to leave.

In order to show our main result, we first prove in Theorem 2 that any non-fixed pair can be:(1) produced by a unique rotation and eliminated by another one; or (2) eliminated by two different rotations and produced by two others (see later Example 1). In the first case, we shall denote by ρ e the elimination rotation and by ρ p the production rotation. In the second case, we shall denote by ρ p1 , ρ p2 the two production rotations and by ρ e1 , ρ e2 the two elimination rotations.

We assume w.l.o.g that the input instance admits at least two stable matchings. For any non-fixed stable pair (p i , p j ), there are two possible cases to consider:

Case 1: (A) f T S (i) = p j and l T S (j) = p i , or (B) l T S (i) = p j and f T S (j) = p i ; Case 2: Otherwise.
Case 1 is a special case indicating that if one of the persons in the pair is the other ones' most preferred person in T S (respectively, the other one is the least preferred person in T S ). Note that, in both cases L T S (i) > 1 and L T S (j) > 1, because the pairs are non-fixed. Later, we refer to these cases for identifying scenarios. In Lemma 3, we show how to identify the elimination rotation(s) for a given pair regardless of its case. Lemma 3. A non-fixed stable pair {p i , p j } is eliminated by a rotation ρ if and only if (p i , p j ) ∈ ρ or (p j , p i ) ∈ ρ.

Proof. → Let ρ = (x 0 , y 0 ), (x 1 , y 1 ) . . . , (x |ρ|-1 , y |ρ|-1 ) be a rotation that eliminates {p i , p j }. Observe first that ρ is non-singular (otherwise {p i , p j } is not stable). Recall that the elimination of ρ from a table T means for each pair (x m , y m ) ∈ ρ, the deletion of {x m , y m } and all pairs {y m , z} such that y m prefers x m-1 to z from T . Table 2 gives an illustration of the preferences of x m and y m . The eliminating ρ moves x m from y m to y m+1 and deletes some {y m , z}. In a similar way, eliminating ρ moves y m from x m-1 to x m and deletes some {x m , z }. Since every close complete subset contains either ρ or ρ, then any pair {y m , z} and {x m , z } cannot be part of any solution. Therefore, if ρ eliminates {p i , p j } and {p i , p j } / ∈ ρ then {p i , p j } is not stable. This contradicts the fact that our pair {p i , p j } is a non-fixed stable pair. ← By the definition of eliminating a rotation ρ from a table T , where (p i , p j ) ∈ ρ, the elimination results in the deletion of p j from p i 's list. Similarly, if (p j , p i ) ∈ ρ then it results in the deletion of p i from p j 's list. Lemma 4 identifies the production rotations.

Lemma 4. If a non-fixed stable pair {p i , p j } is eliminated by ρ e , then {p i , p j } is produced by the dual of it, ρ p = ρe .

Proof. A rotation is said to produce {p i , p j } if eliminating it from a table T reduces L T /ρ (i) to a single entry, namely to p j and L T /ρ (j) to p i . We prove the existence of the production rotations over the two cases (Case 1 and Case 2) identified above.

We have two sub-cases in Case 1. First case is when f T S (i) = p j , l T S (j) = p i . In order to reduce p i 's list to only p j , we need a rotation that moves p i from his/her second best choice up to the first choice. We refer to this operation as limiting p i from right. Similarly, to reduce the p j 's list to only p i , we need a rotation that moves p j from his/her second least-preferred person to the least preferred person. We refer to this operation as limiting p j from left. Referring back to Table 2 for notation, the production rotation ρ p of the pair {p i , p j } = (x m , y m ) must contain the pair (y m+1 , x m ) ∈ ρ p to limit x m from right. Additionally, it must contain (y m , x m-1 ) to limit y m from left. To illustrate, the production rotation has the shape: ρ p = . . . , (y m , x m-1 ), (y m+1 , x m ), . . .. Note that, each ordered pair can only appear in exactly one rotation. Observe that, the dual of ρ p contains the pair (x m , y m ) by definition of dual. By Lemma 3, we know that the rotation that contains (x m , y m ) is the elimination rotation of the pair {p i , p j }. Therefore, ρ p = ρe The proof for the second sub-case is similar, where (y m , x m ) ∈ ρ e .

For a pair {p i , p j } of Case 2, each person has both more and less preferred people in their lists. Therefore, in order to produce a pair, their lists must be limited from both left and right. Let ρ p1 denote the rotation that limits p i from left and p j from right, and ρ p2 denote the rotation that limits p i from right and p j from left, respectively. Let the preference lists for the pair {p i , p j } denoted by L T S (i) = [. . . , y m-1 , y m , y m+1 ] and L T S (j) = [. . . , x m-1 , x m , x m+1 ] where {p i , p j } = (x m , y m ). The pair (x m , y m-1 ) must be in ρ p1 to limit p i from left and (x m+1 , y m ) be in ρ p1 to limit p j from right. Additionally, the pair (y m+1 , x m ) must be in ρ p2 to limit p i from right and (y m , x m-1 ) to limit p j from left. Note that, the dual of ρ p1 contains (y m , x m ), the dual of ρ p2 contains (x m , y m ) by the definition of a dual rotation. By Lemma 3, we know these rotations are elimination rotations of the pair {p i , p j }.

Note that the two rotations ρ p1 and ρ p2 do not require one of them to be eliminated from the table first; they are incomparable. Therefore, depending on the order of elimination, both of them are production rotations.

We sum up the findings above for the non-fixed stable pairs. If a pair is of Case 1, then there exists only one elimination rotation for this pair and only one production rotation as the dual of the elimination one. Because the preference list needs to be limited in only one direction. However, for the pairs of Case 2, there exist two elimination rotations for this pair, and also two other production rotations. Observe that, for each non-fixed stable pair {p i , p j } in a stable matching M, the corresponding complete closed subset of M contains all production rotations of {p i , p j }. It is important to note that, especially for the pairs of Case 2, including one production rotation in the complete closed subset and not the other one, results in producing other partners for that pair. Subsequently, Theorem 2 is an immediate result of Lemmas 3 and 4. Theorem 2. Let {p i , p j } be a non-fixed stable pair. If {p i , p j } is of Case 1, then there exists a unique elimination rotation ρ e , where (p i , p j ) ∈ ρ e or (p j , p i ) ∈ ρ e , and a unique production rotation ρ p , where ρ p = ρe . Otherwise (Case 2), there exist two different elimination rotations ρ e1 and ρ e2 , where (p i , p j ) ∈ ρ e1 , (p j , p i ) ∈ ρ e2 and two rotations ρ p1 = ρ e1 , ρ p2 = ρ e2 that produce the pair.

Let S P denote the set of all the complete closed subsets for the underlying SR instance. Lemma 5 gives a characterisation for the complete closed subsets.

Lemma 5. Let S ∈ S P . For each sink rotation ρ of S, the set S \ {ρ} ∪ {ρ} ∈ S P .

Proof. By definition of closed subset, every predecessor ρ ∈ N -(ρ) is in S. Since ρ is a sink rotation, any successor ρ * ∈ N + (ρ) is not in S. Therefore, by definition of the complete closed subset, we have ρ * ∈ S and ρ is not in S. Using Lemma 2, we know that ρ * ≺ ρ. Hence, all predecessors of ρ are already in S, making ρ a neighbour rotation and results in S \ {ρ} ∪ {ρ} ∈ S P .

The distance between two stable matchings d(M, M ) is previously defined in Section 2 as the number of different pairs between M and M . Observe that the distance can be calculated by also using their corresponding complete closed subsets. If S\S = {ρ}, it means ρ ∈ S and ρ ∈ S . We know that, X({ρ}) = Y ({ρ}) and Y ({ρ}) = X({ρ}). Therefore, between M and M , only the people in ρ (or ρ) have different partners. This can also be generalised to a set of rotations. Hence, the distance can also be denoted as d(S, S ) = |X(S \ S ) ∪ Y (S \ S )|/2. Note that d(S , S) = d(S, S ).

Lemma 6 identifies the closest stable matching to a stable matching M, when a rotation from its corresponding complete closed subset is to be removed. Lemma 6. Given a stable matching M and its corresponding complete closed subset S, if ρ ∈ S is a rotation to remove from S, the closest stable matching M to M such that ρ ∈ S is found by the formula 4 :

C(S, ρ) = S = (S \ ({ρ} ∪ N + (ρ))) ∪ {ρ} ∪ ρ * ∈N + (ρ) ρ * (1)
Proof. The proof of the defined set S being a complete closed subset is obvious by using Lemma 2 and Lemma 5 as flipping a sink rotation of S yields in another complete closed subset. However, if ρ is not a sink rotation in S, we must flip all the successors of ρ to obtain a complete closed subset. Let M * denote the stable matching after flipping ρ ∈ S. Then,

d(M, M * ) = d(S, S * ) = |X({ρ}) ∪ Y ({ρ})|/2. Now, let M * denote the stable matching after flip- ping both ρ, σ ∈ S. Then, d(M, M * ) = |X({ρ}) ∪ Y ({ρ}) ∪ X({σ}) ∪ Y ({σ})|/2.
Observe that, flipping more rotations can only increase the distance between matchings. In Formula 1, the required number of flips is minimum. Therefore the function C(S, ρ) returns the closest stable matching to M when ρ ∈ S to be removed from S.

Finally, Theorem 3 concludes how to find the closest stable matching M to M if {p i , p j } ∈ M wants to leave the M. Theorem 3. Given a stable matching M and a pair {p i , p j } to leave M, the closest stable matching M to M is identified by its corresponding S using the Formula 1 as follows:

1. If Case 1, then S = C(S, ρ p ). 2. If Case 2, let M 1 and M 2 be the two stable matchings s.t. S 1 = C(S, ρ p1 ) and

S 2 = C(S, ρ p2 ). Then S = S 1 if d(M, M 1 ) < d(M, M 2 ), otherwise S = M 2 .
Proof. The proof is immediate from Theorem 2 and Lemma 6.

In order to verify if a given M is a (1, b)-supermatch, all closest stable matchings to the given stable matching are found under the assumption that each non-fixed pair wants to leave the stable matching, one at a time. For each pair, its production rotation is identified and then Theorem 3 is applied to find the closest stable matching. Among all the closest stable matchings, the matching that results in the maximum distance to M sets the robustness of M, i.e. b = d(M, M ) -1, where 1 denotes the pair that wants to leave.

Example 1. Computing robustness Let us calculate the closest matching to M 6 given in Table 1. In Table 3, we identify the cases, and the production/elimination rotation(s) for assuming each pair leaves the M 6 at a time, and we apply Theorem 3 to find the robustness. The pair that has the maximum cost to be repaired sets the robustness value of the matching. Therefore, for this case, the robustness of M 6 is 3.

{p i , p j } Case ρp ρe C(S, ρ) S d(M, M ) S b {p1, p7} 1 ρp = ρ6 ρe = ρ6 {ρ3, ρ4, ρ5, ρ6, ρ7} S4 4 S4 3 {p2, p3} 2 ρp1 = ρ7 ρp2 = ρ3 ρe1 = ρ7 ρe2 = ρ3 {ρ3, ρ4, ρ5, ρ6, ρ7} { ρ3, ρ4, ρ5, ρ6, ρ7} S7 S1 2 4 S7 1 {p4, p9} 1 ρp = ρ6 ρe = ρ6 {ρ3, ρ4, ρ5, ρ6, ρ7} S4 4 S4 3
{p5, p10} 2 ρp1 = ρ4 ρp2 = ρ6 ρe1 = ρ4 ρe2 = ρ6 {ρ3, ρ4, ρ5, ρ6, ρ7} { ρ3, ρ4, ρ5, ρ6, ρ7} S2 S1 3 4 S2 2 {p6, p8} 1 ρp = ρ7 ρe = ρ7 {ρ3, ρ4, ρ5, ρ6, ρ7} S7 2 S7 1
The production and elimination rotations of each pair can be identified in a preprocessing step. We show that checking if a stable matching is a (1, b)-supermatch can be performed in O(n × |V|) time after the O(n 3 log(n)) preprocessing step for an instance where 2 × n people are involved. The preprocessing step consists in identifying the rotations and building the reduced rotation poset (O(n 3 logn)) [START_REF] Gusfield | The Stable Marriage Problem: Structure and Algorithms[END_REF]; identifying all the predecessors and successors of each rotation ρ (O(|V| 2 )); and identifying elimination and production rotations for each pair {p i , p j } whenever applicable in (O(n 2 )). Given a stable matching M, its corresponding complete closed subset S is found by finding and adding the production rotation(s) of each pair and their predecessors into S by starting from an empty set (O(n × |V|)). Conversely, given a closet complete S, M can be constructed by eliminating all the rotations in S from T S by respecting their precedence order. The order is found by applying sorting (O(|V| × log|V|)). The main algorithm is to compute for each pair in M, the closest stable matching M by using Theorem 3. Observe that computing the distance between two stable matchings takes O(n) time and flipping a rotation takes a constant time. Moreover, the worst case of finding the closest stable matching is to flip all the non-singular rotations in S, where the number of all non-singular rotations is |V|/2. Therefore, this computation takes O(|V|) time.

Finding Robust Solutions to the SR

We consider in this section two meta-heuristic approaches to solve the problem of finding a (1, b)-supermatch to a given Stable Roommates instance that minimizes the value of b.

Local Search

Considering the structural similarities between the RSM and the RSR, we tailored the local search model (LS) for the RSM, as it is shown that the LS model produces near optimal solutions for RSM and is better than the proposed genetic algorithm [START_REF] Genc | Finding robust solutions to stable marriage[END_REF]. In the generic LS model, there exists a neighbourhood N for the current solution. The algorithm works by searching the neighbourhood of the current solution, finding the best neighbour M n in the neighbourhood and then proceeding the search by checking the neighbourhood of M n . The aim is to find the stable matching that has the minimum b value. The search is restarted by a random stable matching at every few iterations to avoid getting stuck at a local optimum. The search continues until a termination criterion is met.

In our model, we have four termination criteria. The first one is a cut-off limit lim cutof f , which "counts" the number of steps since the last best solution is found. The second one is the depth limit lim depth , which indicates the depth of the neighbourhood search starting from a random stable matching. Another criterion is the optimality opt, which indicates if the algorithm has already found a solution with b = 1. Finally, we use a time limit lim time for each instance.

The procedure starts by creating a random stable matching M c as follows. We first mark all the non-singular rotations as available. Let A denote the set of rotations that are available. Then, we randomly select a rotation ρ from A and add it to the initially empty S c . Subsequently, we remove ρ and ρ from A. We also add all predecessors ρ of ρ that are not in S c to S c and remove ρ and ρ from A. This operation operates in a loop until |S c | = |V|/2. Once the complete closed subset S c is found, its corresponding stable matching is computed by eliminating all rotations in S c from T S by respecting their precedence order.

After creating a random stable matching M c , the neighbourhood N of M c is found by checking all the sink rotations in S c . By using Lemma 5, we know that flipping any sink rotation in M c creates another stable matching M n , which we refer as a neighbour of the M c . The general procedure is the same as the one developed for the RSM [START_REF] Genc | Finding robust solutions to stable marriage[END_REF]. In brief, the process starts by descending from the M c by finding N of M c . The next iteration descends from the neighbour of M c that has the lowest b value. This loop is restarted every lim depth iteration by a random M c . The stable matching that has the minimum value of b as found during the search is returned as the solution.

The complexity of the LS procedure depends on the computation of the b values. Finding neighbours is based on the identification of the sink rotations of S c , where there can be at most |V|/2 sink rotations and then a constant cost for flipping each sink rotation. The best neighbour is identified after computing b values of |N | stable matchings. This procedure takes

O(k × n × |V| × |N |),
where k is the number of iterations and n is the number of non-fixed people.

Genetic Local Search (Hybrid)

Combining different search techniques to enhance the performance of a single model is proven to improve solution quality and the models [START_REF] Kolen | Genetic local search in combinatorial optimization[END_REF][START_REF] Nico | Genetic local search algorithms for the traveling salesman problem[END_REF]. Genc et al. propose three different models (constraint programming, local search and genetic algorithm) for finding (1, b)-supermatches to the RSM in [START_REF] Genc | Finding robust solutions to stable marriage[END_REF]. The results indicate that genetic algorithm (GA) procedure has poor performance when compared to the LS. In this work, we consider combining the two metaheuristics: the genetic algorithm and the local search to provide a hybrid procedure. We denote this hybrid model as HB. The overview of the GA procedure we use in the HB model is the same as the one used for RSM (details can be found in [START_REF] Genc | Finding robust solutions to stable marriage[END_REF]).

The procedure begins by initialising a population of random stable matchings. Then, the population is evolved by randomly selecting individuals from the population, applying crossover, searching for neighbours of the products of crossover, applying mutation. This process is repeated until some termination criteria is met (no improvement, timelimit exceeded, optimal solution found). The procedure below gives a pseudo-code of the evolution phase of HB.

1: procedure EVOLUTION() 2: M1 ← SELECTION() 3: M2 ← SELECTION() 4: if M1 = M2 then 5: (Mc1, Mc2) ← CROSSOVER(M1, M2) 6: N ←FINDNEIGHBOURS(Mc1) 7: Mc1 ← BEST(N ) 8: N ←FINDNEIGHBOURS(Mc2) 9: Mc2 ← BEST(N ) 10: REFINE(Mc1, Mc2 ) 11: EVALUATION() 12:
M f it ← GETFITTEST(P ) 13:

Mm ← SELECTION() 14:

rand ← RANDOM(0, 1) 15:

if Mm = M f it and rand < pm then 16:

MUTATION(Mm)

As can be seen from the procedure, the only LS enhancement to the GA algorithm is the search for the neighbours of the stable matchings after crossover (see Lines 6-9). Let M c1 , M c2 be the two stable matchings produced by the crossover. We update M c1 by its best neighbour after the neighbour search (same applies to the M c2 ). Creating a random stable matching and finding neighbours are already discussed in Section 4.1.

If the original methods from LS and GA as described in [START_REF] Genc | Finding robust solutions to stable marriage[END_REF], where the evolution phase is updated with the one here are used, we obtain the HB model for the RSM. In the RSR model, only the crossover and mutation operations are different than the original GA model defined for RSM. Instead of defining the crossover by adding rotations to the closed subset or removing them as we did for the RSM, we use the terminology flip for the RSR. Considering the Lemma 6, we define the crossover procedure for two stable matchings M 1 , M 2 as follows. First, we find a random rotation ρ 1 ∈ S 1 , and a random rotation ρ 2 ∈ S 2 . If ρ 1 / ∈ S 2 , then ρ1 ∈ S 2 due to the completeness property of the closed subsets in SR. Therefore, we flip ρ1 in S 2 and the duals of all of its predecessors ρ ∈ N -(ρ) if ρ is not included in S 2 . We apply the same procedure to the other stable matching as well. Moreover, for the mutation operation, we select a random rotation ρ from the reduced rotation poset of the underlying instance and also a stable matching M. If ρ ∈ S, we flip ρ and all the required predecessors. If its dual ρ ∈ S, then we flip ρ and the predecessors.

Experiments

In this section, we first compare the performances of the HB and the LS proposed for the RSR. Then, we investigate the robustness of different sets of RSM and RSR instances 5 . The code is implemented in Java, reusing the RSM experiments from [START_REF] Genc | Finding robust solutions to stable marriage[END_REF]. All experiments are performed on a Dell M600 with 2.66 Ghz processors under Linux, using three different randomisation seeds and fixing time limit lim time = 20 mins, number of iterations without improvement lim cutof f = 10000, the number of stable matchings in LS that descend from a random stable matching lim depth = 50. We use the population size for the HB as |P | = 30 and the mutation probability as p m = 0.7. We use a high p m as GA is suffering from getting stuck at local minima and randomisation helps with it. We discuss the size of the population for HB later. HB v. LS for the RSR. Our first experiment is about the comparison of LS and HB models. Random SR instances have only a small number of stable matchings as we verify on the dataset RANDOM later [START_REF] Pittel | The" stable roommates" problem with random preferences[END_REF]. For the comparison of HB and LS models, we look for instances that are likely to contain many stable matchings to gain more insight on their performances. For this purpose, we first created a dataset of purely random SM instances as each SM instance contains at least one stable matching, then we converted these instances to the SR. This conversion tackles the problem of random SR instances having only a few stable matchings, while preserving the randomness. Our SM dataset consists of 30 random instances of each size n ∈ {100 × k | k ∈ {1, . . . , 10}}. Note that, the resulting SR instances have size 2 × n.

Fig. 3 and Fig. 4 provide detail on the comparison between the LS and the HB.6 Fig. 3 compares the average minimum b value found by the two models for each instance in the set. In the x-axis, the range shows the size of the instances such that all the instances that have x-values between [0, 200] is of size 200, [201 -400] is of size 400, etc. We confirm by our experiments and also observe in Fig. 4 that for each instance that has size 200 ≤ n ≤ 600, both models complete the search within the given time limit. Additionally, they either produce similar results (b values) or HB performs slightly better as can be observed in Fig. 3. The reason for exceeding the time limit in Fig. 4 is due to us not interrupting the construction of a stable matching. The construction of a stable matching consists of exposing all rotations in its complete closed subset in order starting from T S . Then, the b value is computed. For large stable matchings, this computation becomes very costly. We can conclude that for small instances, both HB and LS perform well in terms of finding solutions with low b values. If the time is essential, HB model can be preferred over LS as it converges faster. Additionally, HB is able to find better solutions for larger instances.

Random RSM v. Random RSR. Next, we perform some tests for SM-SR comparison on our dataset RANDOM. Our dataset RANDOM consists of 30 purely randomly created SM and SR instances for each size n ∈ {100 × k | k ∈ {1, . . . , 10}}. Note that, for an SM instance of size n, there exists n men and n women in the problem. We have 2 × n people in the corresponding SR instance. However, both have n pairs. All SR instances in RANDOM have at least two stable matchings. Considering the good performance of LS for small instances, we used the LS models for both RSM and RSR.

Table 4 and Table 5 present a summary of the robustness of random RSM and RSR instances. The columns report for each size the average value of: the total number of pairs in the instance (n), the number of rotations in the rotation poset or the reduced rotation poset (|V|), the number of different stable matchings found during the search of LS (sm), the number of non-fixed pairs (np), the b value of the solution found (b), the ratio b np (ratio), and the time spent until finding the best solution by LS in seconds (t best ).

Observe from the tables that the random RSM instances contain many more stable matchings than the random RSR instances of similar sizes. Recall that, the value of sm denotes only the number of a subset of the stable matchings found during the search. However, we can confirm the RSR instances not containing many stable matchings by looking at the number of rotations in their rotation posets. Note that, for the RSR instances, when 1000 pairs are included, the corresponding rotation posets, on the average, contain |V| ≈ 5 rotations. This is mainly caused by the large numbers of fixed-pairs in the random RSR instances. For instance, the average number of non-fixed pairs in the RSM instances of size 100 is np = 75.12. However, we observe in the large RSR instances that contain 1000 pairs that there are only 55. [START_REF] Pittel | The" stable roommates" problem with random preferences[END_REF] non-fixed pairs on the average, which is less than the smallest sized RSM instances that we tested. Note that, we measure the robustness ratio over the non-fixed pairs of the instances. It is desirable to obtain a smaller value for the ratio to indicate a better robustness for the instance. Because a smaller ratio indicates that a smaller proportion of the people that have alternative partners need to change their partners for a repair. Observe that, the ratio of the RSR instances is lower when compared to RSM. The ratio shows that the breakage of the pairs in the RSR instances are less costly to be repaired. Thus, we conclude that purely random RSR instances require a smaller proportion of the people to change their partners in the case of a breakage, when compared to the RSM.

Large RSM and RSR instances. In this experiment, we search for instances with potentially many number of stable matchings and low b values. Therefore, we generated a dataset called MANY consisting of 100 SM instances for each size n = {16, 32, 64, 128} using the family described by Irving and Leather [START_REF] Irving | The complexity of counting stable marriages[END_REF], and then used in [START_REF] Siala | Rotation-based formulation for stable matching[END_REF]. Note that, each SM instance in this set has a corresponding SR instance (see Lemma 1), where the corresponding SR instance has a reduced rotation poset of twice the size of the rotation poset of the SM instance. First, let us introduce this family of instances described by Irving and Leather. They prove that any instance in the original family contains at least 2 n-1 stable matchings for an instance of size n = 2 i . They define this family over two matrices for the preferences of each gender, and the preference lists of these large instances are obtained recursively by appending the following matrices until the desired instance size is found. In our dataset MANY, we slightly modify each instance of this original family by first randomly selecting two random men m i , m j . Then, we modify m i 's preference list by swapping the positions of two randomly selected women within the list. We repeat the same for m j . We also modify the preference lists of two random women in the same way. In other words, the original preference set between the original and the modified instances have a Hamming Distance of 8.

Table 6 reports for each size the average value of: the number of all men or women (n), the number of non-fixed men(np), the number of rotations in the rotation poset (|V|). Additionally, it reports the average minimum b found by the model LS, HB where population size |P | = 10, and HB where population size |P | = 60 (b); followed by the total time spent in minutes for each of the three models (t (min)).

This dataset shows that the robustness of instances that have many stable matchings is very high (i.e the value of b is low). For each instance size, our best model for that size is able to find solutions whose average b values evaluate to a b value that is opt = 1 < b < 2. For instance, for size n = 16, the LS model finds solutions such that for the breakage of any man on the solution, on the average, 1.12 other men need to breakup from their current partners. Similarly, for size n = 128, HB models find that the solution is guaranteed to be repaired by only 1.02 additional men's break-up.

As one can observe from Table 6, we ran the HB model by using different sizes of population. Observe that, reducing the number of individuals in the population of HB (60 to 10) causes the algorithm to find slightly worse solutions (i.e. larger b). For instance, for size n = 64 , the average minimum b is found as 1.74 by a population of size 10, and 1.28 by a population of size 60. This is due to having an increased chance of getting stuck at local minima for a smaller population. On the other hand, LS finds competitive values for b for sizes 16 ≤ n ≤ 64. However, as we can see for n = 128, LS finds solutions that are far away from the optimal solution. We conclude that, an improvement for HB by changing population size is possible in exchange of obtaining slightly worse solutions. LS performs well for smaller instances.

Recall that, each SM instance in MANY has a corresponding SR instance that has exactly the same stable matchings. We do not run the RSR models on this dataset as they are much slower. However, this test provides an insight to some RSM and RSR instances that are repairable at low additional costs.

Conclusions

We study the notion of (a, b)-supermatch in the context of Stable Roommates problem. We propose a polynomial-time algorithm based on the reduced rotation poset to verify if a stable matching is a (1, b)-supermatch. Next, we use this procedure to design local search (LS) and hybrid genetic local search (HB) models to find robust solutions for the (1, b) case (i.e., (1,b)-supermatch with (possibly) the minimum b). We empirically show that the HB model usually performs better than LS for RSR. Furthermore, we perform an RSM/RSR comparison and identify a family of instances that are rich in stable matchings and very robust.

Acknowledgement

This material is based upon works supported by the Science Foundation Ireland under Grant No. 12/RC/2289 which is co-funded under the European Regional Development Fund.

Fig. 1 .

 1 Fig. 1. The TS for an SR instance I of size 10. pi LT S (i) 1 3 4 7 2 4 3 8 9 3 5 6 2 1 4 9 1 6 2 5 7 10 8 3 6 8 3 4 10 7 1 5 8 10 2 5 6 9 2 10 4 10 6 5 9 8

Definition 1 (

 1 (a, b)-supermatch). Given an SR instance I, and two positive integers a, b ∈ N, a stable matching M of Iis said to be an (a, b)-supermatch if for any set Ψ ⊆ M of non-fixed stable pairs, where |Ψ | = a, there exists a stable matching M such that M ∩ Ψ = ∅ and d(M, M ) ≤ b + a.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Avg min b value found by LS and HB.

Table 1 .

 1 Fig.2. Reduced rotation poset of I given in Table1. A list of all the stable matchings and their corresponding complete closed subsets of I.

		(1, 3), (2, 4)	(8, 10), (9, 2)	(3, 5), (10, 6)
			ρ 3	ρ 5		ρ 4
	3 4 7
	2 4 3 8 9
	3 5 6 2 1	(2, 3), (6, 8)				(4, 9), (7, 1), (10, 5)
	4 9 1 6 2 5 7 10 8 3 6 8 3 4 10	ρ 7	ρ 6 (9, 10), (1, 4), (5, 7)	ρ 7 (3, 6), (8, 2)	ρ 6
	7 1 5
	8 10 2 5 6
	9 2 10 4		ρ 4	ρ 5		ρ 3
	10 6 5 9 8	(6, 3), (5, 10)	(2, 8), (10, 9)	(4, 1), (3, 2)
	M1 = {(1, 3), (2, 4), (5, 7), (6, 8), (9, 10)} S1 = { ρ3, ρ4, ρ5, ρ6, ρ7}
	M2 = {(1, 7), (2, 8), (3, 5), (4, 9), (6, 10)} S2 = {ρ3, ρ4, ρ5, ρ6, ρ7}
	M3 = {(1, 4), (2, 9), (3, 6), (5, 7), (8, 10)} S3 = {ρ3, ρ4, ρ5, ρ6, ρ7}
	M4 = {(1, 4), (2, 3), (5, 7), (6, 8), (9, 10)} S4 = {ρ3, ρ4, ρ5, ρ6, ρ7}
	M5 = {(1, 4), (

Table 2 .

 2 An illustration of the preferences p Preference lists . . . . . . xm . . . , ym, z , ym+1 . . . . . . . . . ym . . . , xm-1, z, xm, . . . . . . . . .

Table 3 .

 3 Computing the closest matching to M6

Table 4 .

 4 Results on uniformly random instances for RSM.

	n	|V| sm	np	b ratio t best
	100 22.02 47.9 75.12 48.27 0.64 0.02
	200 41.59 116.9 166.19 115.34 0.69 0.10
	300 60.22 182.4 263.94 193.08 0.73 0.37
	400 74.51 244.1 356.58 265.98 0.74 0.77
	500 91.47 322.5 456.00 350.16 0.76 2.18
	600 103.82 394.9 551.10 425.51 0.77 3.67
	700 117.08 449.6 646.69 505.61 0.78 5.89
	800 131.81 527.6 749.98 595.64 0.79 9.09
	900 146.34 585.5 848.32 679.82 0.80 14.60
	1000 156.00 632.4 943.23 758.82 0.80 21.16

Table 5 .

 5 Results on uniformly random instances for RSR.

	n |V| sm np	b ratio t best
	100 3.91 3.78 17.91 5.31 0.3 0.003
	200 3.87 3.94 26.76 8.52 0.32 0.003
	300 4.36 4.56 35.53 11.22 0.32 0.017
	400 4.71 5.92 37.64 10.93 0.29 0.048
	500 4.29 4.81 37.62 11.70 0.31 0.066
	600 4.16 4.48 42.44 14.47 0.34 0.130
	700 4.58 5.50 48.71 16.02 0.33 0.312
	800 4.93 5.99 55.02 17.39 0.32 0.498
	900 4.82 7.07 57.64 18.50 0.32 0.662
	1000 4.60 5.19 55.16 18.36 0.33 0.557

Table 6 .

 6 Summary of the results on large instances for RSM.

	instance	LS	HB, |P | = 10 HB, |P | = 60
	n np	|V|	b t (min) b	t (min)	b	t (min)
	16 15.99 100.43 1.12 0.003 1.21 0.003	1.1 0.004
	32 31.99 447.26 1.03 0.054 1.30 0.024 1.04 0.029
	64 64 1889.95 1.685 3.158 1.74 0.824 1.28 0.916
	128 128 7788.02 14.055 8.367 1.02 13.989 1.01 17.609

The parentheses are used to indicate priority.

Our datasets are publicly available at: github.com/begumgenc/rsmData

The reader is referred to the online version for coloured version.


Simon Gsell

Rémi Bourguet
email: remi.bourguet@imft.fr

Marianna Braza

One-versus two-degree-of-freedom vortex-induced vibrations of a circular cylinder at Re = 3900

The one-versus two-degree-of-freedom vortex-induced vibrations of a circular cylinder are investigated on the basis of direct numerical simulation results. The Reynolds number, based on the oncoming flow velocity and cylinder diameter, is set to 3900. Three cases are examined: the elastically mounted body is free to oscillate either in the direction aligned with the current (in-line direction; IL case), in the direction normal to the current (crossflow direction; CF case), or in both directions (IL+CF case). In each case, the behavior of the flow-structure system is studied over a range of values of the reduced velocity (inverse of the oscillator natural frequency). The in-line and cross-flow responses observed in the IL+CF case substantially differ from their one-degree-of-freedom counterparts, especially in the intermediate reduced velocity region. In this region, no vibrations develop in the IL case and in-line oscillations only occur if cross-flow motion is allowed. These in-line oscillations are accompanied by a major increase of the cross-flow responses, compared to the CF case. The two-degree-of-freedom vibrations are associated with the emergence of large-amplitude higher harmonics in the fluid force spectra. These aspects and more specifically the impact of the existence of a degree-of-freedom and oscillations in a given direction, on the fluid force and structural response in the perpendicular direction, do not seem to be systematically connected to changes in wake topology. Here, they are discussed in light of the orientation and magnitude of the instantaneous flow velocity seen by the moving body.

Introduction

Vortex shedding downstream of a bluff body immersed in a cross-flow is accompanied by unsteady fluid forces exerted on the body. If the body is flexible or flexibly mounted, these forces may lead to structural vibrations, called vortex-induced vibrations (VIV). VIV occur when the body oscillation and the unsteady wake synchronize, a mechanism referred to as lock-in. Many natural and industrial systems are subjected to VIV. Their physical analysis and prediction have motivated a number of research works, as reviewed by [START_REF] Bearman | Vortex shedding from oscillating bluff bodies[END_REF], [START_REF] Sarpkaya | A critical review of the intrinsic nature of vortex-induced vibrations[END_REF], [START_REF] Williamson | Vortex-induced vibrations[END_REF] and [START_REF] Païdoussis | Fluid-Structure Interactions: Cross-Flow-Induced Instabilities[END_REF].

Even though most of real systems subjected to VIV involve slender flexible bodies (e.g. chimneys, marine risers), this phenomenon has been extensively studied through the canonical problem of a rigid circular cylinder mounted on an elastic support allowing oscillations in the cross-flow direction (i.e. perpendicular to the oncoming flow), the direction along which large oscillation amplitudes are generally expected (e.g. Huera-Huarte and [START_REF] Huera-Huarte | Wake structures and vortex-induced vibrations of a long flexible cylinder -part 1: dynamic response[END_REF]. Through this simplified configuration, VIV can be analyzed with a limited number of structural parameters and a single structural mode (e.g. [START_REF] Feng | The measurement of vortex-induced effects in flow past stationary and oscillating circular and D-section cylinders[END_REF][START_REF] Hover | Forces on oscillating uniform and tapered cylinders in crossflow[END_REF][START_REF] Khalak | Motions, forces and mode transitions in vortex-induced vibrations at low mass-damping[END_REF][START_REF] Govardhan | Modes of vortex formation and frequency response of a freely vibrating cylinder[END_REF][START_REF] Blackburn | A complementary numerical and physical investigation of vortex-induced vibration[END_REF][START_REF] Shiels | Flow-induced vibration of a circular cylinder at limiting structural parameters[END_REF][START_REF] Leontini | The beginning of branching behaviour of vortex-induced vibration during two-dimensional flow[END_REF]. Significant body oscillations are known to occur over a well-defined range of values of the reduced velocity, defined as the inverse of the oscillator natural frequency non-dimensionalized by the inflow velocity and body diameter. In this range called the lock-in range, the oscillation frequency coincides with wake frequency. The body response amplitude exhibits a bell-shaped evolution as a function of the reduced velocity. Peak amplitudes of the order of one body diameter can be observed, depending on the structural properties [START_REF] Khalak | Fluid forces and dynamics of a hydroelastic structure with very low mass and damping[END_REF] and Reynolds number (Re), based on the oncoming flow velocity and body diameter [START_REF] Govardhan | Modes of vortex formation and frequency response of a freely vibrating cylinder[END_REF].

The response of an elastically mounted cylinder restricted to move in the in-line direction (i.e. parallel to the oncoming flow) has also been addressed in prior studies (Aguirre, 1977;[START_REF] Naudascher | Flow-induced streamwise vibrations of structures[END_REF][START_REF] Okajima | Flow-induced in-line oscillation of a circular cylinder[END_REF][START_REF] Cagney | Wake modes of a cylinder undergoing free streamwise vortex-induced vibrations[END_REF]. In-line oscillations are known to occur over a range of lower reduced velocities, compared to the typical lock-in range associated with cross-flow VIV. The evolution of the response amplitude as a function of the reduced velocity is characterized by two bell-shaped branches, separated by a region of low-amplitude oscillations. The peak amplitudes are comparable in both branches; they are typically of the order of 0.1 body diameters.

Recently, more attention has been paid to the case where the body is allowed to oscillate in both the in-line and crossflow directions. [START_REF] Cagney | Streamwise vortex-induced vibrations of cylinders with one and two degrees of freedom[END_REF] studied the two-degree-of-freedom VIV of a cylinder in the low reduced velocity region, where in-line oscillations occur for a single-degree-of-freedom system. No cross-flow oscillations appear in the first branch of in-line oscillations. In the second branch, cross-flow oscillations develop with amplitudes similar to those noted in the in-line direction. In both branches, the in-line response is almost unaltered compared to that observed for a singledegree-of-freedom system. The effect of adding a degree of freedom in the in-line direction on the system behavior for higher reduced velocities, i.e. where cross-flow oscillations are expected to occur, has been examined by [START_REF] Jauvtis | The effect of two degrees of freedom on vortex-induced vibration at low mass and damping[END_REF]. When the structural mass and damping are small, large-amplitude in-line oscillations are superimposed to the crossflow responses, and the body typically exhibits figure-eight-shaped trajectories, as also confirmed by several studies [START_REF] Dahl | Dual resonance in vortex-induced vibrations at subcritical and supercritical reynolds numbers[END_REF][START_REF] Navrose | Free vibrations of a cylinder: 3-d computations at Re= 1000[END_REF][START_REF] Gsell | Two-degree-of-freedom vortex-induced vibrations of a circular cylinder at re = 3900[END_REF]. The in-line oscillations are accompanied by an increase of the cross-flow oscillation amplitudes, compared to the one-degree-of-freedom case.

The behaviors of the one-and two-degree-of-freedom systems have been widely studied but most of prior works have considered each system separately. The present study aims at a joint analysis of these system responses, over a range of reduced velocities. Among other aspects, the alteration of the in-line vibrations in the intermediate reduced velocity range, and the evolutions of fluid force statistics and frequency content, when adding/removing a degree of freedom to the system, still need to be investigated. More generally, several connections between the in-line and cross-flow responses remain to be clarified.

In the present work, the behaviors of three distinct systems where the cylinder is free to oscillate either in the in-line direction, in the cross-flow direction, or in both directions, are examined over an interval of reduced velocities encompassing the region of cross-flow VIV. The Reynolds number is set to 3900 as a typical case of the early turbulent regime. The comparison between the structural responses, fluid forces and wake patterns observed for the one-and two-degree-offreedom systems is based on direct numerical simulation results. Particular attention is paid to the relative impact of the in-line and cross-flow oscillations on fluid forces and simple mechanisms are proposed to shed some light on the responses of the two-degree-of-freedom system.

The paper is organized as follows. The methodology employed in this study is presented in Section 2. The behaviors of the three systems are described Section 3. Some elements concerning in-line/cross-flow motion interaction and forcing are discussed in Section 4. The principal findings of this work are summarized in Section 5.

Method

The physical systems are described in Section 2.1. The numerical method and data processing approach are presented in Section 2.2. The tools employed to analyze fluid forces are introduced in Section 2.3.

Physical systems

A sketch of the physical configuration is presented in Fig. 1(a). An elastically mounted, rigid circular cylinder of diameter D and mass per unit length ρ c is immersed in a cross-flow. The cylinder axis is parallel to the z axis. The flow, parallel to the x axis, is characterized by its velocity U, density ρ f and dynamic viscosity µ. The Reynolds number based on U and D, Re = ρ f UD/µ, is set to 3900. The flow dynamics is governed by the three-dimensional incompressible Navier-Stokes equations. The cylinder is elastically mounted and free to oscillate either in the in-line direction (x axis), in the cross-flow direction (y axis), or in both directions, as schematized in Fig. 1(b); the three distinct cases are referred to as IL, CF and IL+CF, respectively. The structural stiffness and damping ratio in the i (x or y) direction are designated by k i and γ i . All the physical quantities are made non-dimensional by D, U and ρ f . The non-dimensional cylinder displacement, velocity and acceleration in the i direction are denoted by ζ i , ζi and ζi . The force coefficient in the i direction is defined as C i = 2F i /ρ f DU 2 , where F i Fig. 1. Sketch of the physical configuration: (a) general two-degree-of-freedom system, and (b) three cases addressed in this paper, where the cylinder is allowed to oscillate in the in-line direction (IL), in the cross-flow direction (CF ), or in both directions (IL+CF ).

denotes the span-averaged force in the i direction. The non-dimensional mass is defined as m = ρ c /ρ f D 2 ; it is set to 2. The body dynamics in the i direction is governed by a forced second-order oscillator equation:

ζi + 4π γ i U * i ζi + ( 2π U * i ) 2 ζ i = C i 2m
.

(1)

The reduced velocity in the i direction is defined as

U * i = 1/f nat,i , where f nat,i is the natural frequency in vacuum, f nat,i = D/2π U √ k i /ρ c .
In the one-degree-of-freedom cases (IL and CF ), the reduced velocity U * i is denoted by U * , and the natural frequency f nat,i by f nat . In the two-degree-of-freedom case (IL+CF ), the structural stiffnesses are the same in both directions; the reduced velocity and natural frequency of the oscillator are referred to as

U * = U * x = U * y and f nat = f nat,x = f nat,y
. The damping ratio is set equal to zero in both directions to allow maximum amplitude oscillations (γ i = 0).

Numerical method and data processing

The behavior of the coupled flow-structure system is predicted by direct numerical simulation of the three-dimensional Navier-Stokes equations. The numerical procedure is identical to that employed in [START_REF] Gsell | Two-degree-of-freedom vortex-induced vibrations of a circular cylinder at re = 3900[END_REF]. The computations are performed using the finite-volume code Numeca Fine/Open (www.numeca.com), which is based on second-order spatial schemes and a second-order dual-time-stepping time integration. The Navier-Stokes equations are expressed in the cylinder frame which avoids any grid deformation. The frame motion is taken into account by adding inertial terms in the Navier-Stokes equations. At each physical time step, the body motion equations (1) are solved implicitly through a dual-timestepping scheme. Flow and body solutions are therefore updated simultaneously.

A detailed convergence study has been performed in order to set the numerical parameters, as reported in [START_REF] Gsell | Two-degree-of-freedom vortex-induced vibrations of a circular cylinder at re = 3900[END_REF]. The flow is discretized on a non-structured grid in a rectangular computational domain. The streamwise (x), crossflow (y) and spanwise (z) lengths of the domain, non-dimensionalized by D, are 120, 60 and 3. Periodic boundary conditions are used in the spanwise and cross-flow directions. The grid size in the wall-normal direction at the cylinder surface is ∆n = 1.5 × 10 -3 . In the spanwise direction, 80 cells are considered. The total number of cells is equal to 11.5 × 10 6 .

All the computations are initialized with a static body. The reliability of the simulation approach was assessed in [START_REF] Gsell | Two-degree-of-freedom vortex-induced vibrations of a circular cylinder at re = 3900[END_REF], where the two-degree-of-freedom system responses and associated fluid forces were found to be close to the experimental results of [START_REF] Jauvtis | The effect of two degrees of freedom on vortex-induced vibration at low mass and damping[END_REF]. Additional elements of validation can be found in Section 3.1 for the CF case: the present simulation results match the experimental data reported by [START_REF] Hover | Forces on oscillating uniform and tapered cylinders in crossflow[END_REF] for a similar physical system.

The physical quantities are analyzed over time series of more than 20 oscillation cycles, collected after convergence of the structural response. The maximum amplitude of a time-dependent signal s, denoted by s m , is defined as the average of the highest 10% of its amplitudes. The time-averaged value of s is denoted by s, ands = ss designates the fluctuating part of s. The root-mean-square (RMS) value of s is denoted by s ′ . The dominant frequency of the body response in the i direction, based on the Fourier transform of ζ i time series, is denoted by f i . The frequency ratio is defined as

f * i = f i /f nat . As
in [START_REF] Gsell | Two-degree-of-freedom vortex-induced vibrations of a circular cylinder at re = 3900[END_REF], the signals were high-pass filtered in order to avoid low-frequency fluctuations which are not occurring through lock-in. A span and phase averaging procedure of the flow quantities is employed to determine the wake patterns in Section 3.3. The phase averaging is performed over 4 oscillation cycles. For each cycle, a series of 5 snapshots close to the targeted phase are selected. The phase-averaged fields are thus computed with 20 snapshots.

Neglecting higher harmonic terms, the square of the instantaneous flow velocity magnitude can be expressed as follows:

Ψ ≈ Ψ + Ψ 2 sin(4π f 1 t + φ Ψ 2 ), (7)
where Ψ 2 and φ Ψ 2 denote the amplitude and phase of the second harmonic of Ψ . Combining expressions (4), ( 6) and ( 7), the fluid force coefficients in the η and ξ directions can be modeled as follows:

C η ≈ [ 1 + κ η (Ψ -1) ] C f x + 1 2 κ η Ψ 2 C f x,2 cos(φ Ψ 2 -φ C f x,2 ) +κ η C f x Ψ 2 sin(4π f 1 t + φ Ψ 2 ) + [ 1 + κ η (Ψ -1) ] C f x,2 sin(4π f 1 t + φ C f x,2 ) - 1 2 κ η Ψ 2 C f x,2 cos(8π f 1 t + φ Ψ 2 + φ C f x,2 ), (8)
and

C ξ ≈ [ 1 + κ ξ (Ψ -1) ] C f y,1 sin(2π f 1 t + φ C f y,1 ) + 1 2 κ ξ Ψ 2 C f y,1 cos(2π f 1 t + φ Ψ 2 -φ C f y,1 ) - 1 2 κ ξ Ψ 2 C f y,1 cos(6π f 1 t + φ Ψ 2 + φ C f y,1
).

(9)

The connections between fluid forces and the two kinematic quantities α and Ψ will be explored in Section 4.

Results

The results issued from the numerical simulations for the one-and two-degree-of-freedom systems are reported in this section. The body responses are quantified in Section 3.1, the fluid forces are examined in Section 3.2 and some elements regarding wake topology are presented in Section 3.3.

Structural responses

The structural responses obtained in the three studied cases (IL, CF and IL+CF ) are depicted in Fig. 3. In the CF case, the oscillation amplitudes match the experimental data of [START_REF] Hover | Forces on oscillating uniform and tapered cylinders in crossflow[END_REF] for U * < 9 (Fig. 3(a)). The peak amplitude, observed for U * = 5, is approximately equal to 0.8 diameters. For U * > 9, the oscillation amplitudes reported by Hover et al.

(1998) decrease down to 0.2 diameters, while relatively large vibrations are still observed in the present numerical results. As shown in [START_REF] Gsell | Two-degree-of-freedom vortex-induced vibrations of a circular cylinder at re = 3900[END_REF], the body response is particularly sensitive to structural damping in the high reduced velocity region: a small structural damping can substantially decrease the response amplitude in this region. In the experiments of [START_REF] Hover | Forces on oscillating uniform and tapered cylinders in crossflow[END_REF], the structural damping is small but still not negligible; this may explain the differences observed with the present results. The cross-flow oscillation amplitudes obtained in the IL+CF case are also reported in Fig. 3(a). The cross-flow response is globally amplified when the in-line degree of freedom is added to the system, as previously noted by [START_REF] Jauvtis | The effect of two degrees of freedom on vortex-induced vibration at low mass and damping[END_REF]. The peak oscillation amplitude is lower than the maximum amplitude reported by [START_REF] Jauvtis | The effect of two degrees of freedom on vortex-induced vibration at low mass and damping[END_REF]. This may relate to the value of the Reynolds number, which is higher in the experiments as discussed in [START_REF] Gsell | Two-degree-of-freedom vortex-induced vibrations of a circular cylinder at re = 3900[END_REF]. Otherwise, the cross-flow amplitudes remain globally close to the experimental data. The differences noted in the high reduced velocity region can be explained by the above mentioned effect of the structural damping.

The in-line oscillation amplitudes obtained in the IL and IL+CF cases are plotted in Fig. 3(b). In the IL+CF case, the amplitudes match those reported by [START_REF] Jauvtis | The effect of two degrees of freedom on vortex-induced vibration at low mass and damping[END_REF]. Small oscillations occur in the low reduced velocity region (U * = 3). This region has been thoroughly analyzed by [START_REF] Cagney | Streamwise vortex-induced vibrations of cylinders with one and two degrees of freedom[END_REF] for a slightly different system (pivoted cylinder). They found that the in-line vibrations are not impacted by the addition of the cross-flow degree of freedom. This observation is confirmed by the present results. In contrast, at higher reduced velocities, the IL and IL+CF cases exhibit distinct behaviors: negligible vibrations are noted in the IL case, while significant oscillations appear in the IL+CF case, especially around U * = 6. In this region, which coincides with the region of peak cross-flow oscillation amplitudes, in-line response amplitudes up to 0.3 diameters are encountered. The in-line vibration amplitude rapidly decreases beyond this region where only residual oscillations persist. The present results suggest that in the IL+CF case, the in-line oscillations occurring in the intermediate reduced velocity range are closely connected to the presence of cross-flow motion. This aspect is discussed in Section 4.

Body responses are generally close to harmonic. Their frequencies are quantified in Fig. 3(c) which represents the evolutions of the frequency ratios in both directions, as functions of the reduced velocity. In the IL case, the oscillations appearing for U * = 3 occur close to the oscillator natural frequency. The frequency of the oscillations of negligible amplitude noted at higher U * matches the frequency of the in-line fluid force in the fixed body case, equal to twice the Strouhal frequency (i.e. vortex shedding frequency in the fixed body case, denoted by f st in the plot). In the CF case, the oscillation frequency remains close to the Strouhal frequency for U * ∈ [3,6]. At higher reduced velocities, f * y reaches a plateau, and the oscillation frequency significantly departs from the Strouhal frequency. In the IL+CF case, the in-line response frequency is equal to twice the cross-flow response frequency over the range of U * under study. For U * = 3, the in-line The RMS value of the fluctuating in-line force coefficient is plotted in Fig. 4(b). The CF case results indicate a major amplification of Cx when the body oscillates in the cross-flow direction, compared to the fixed body case. Such amplification of the in-line force fluctuation is expected to impact the in-line responses when the body is allowed to oscillate in both directions. An increase of C ′

x , even though less pronounced than in the CF case, is also noted in the IL case (for

U * = 3). C ′
x is generally larger in the two-degree-of-freedom case in the intermediate reduced velocity region. In particular, the peak of C ′ It is recalled that f 1 denotes the dominant frequency of the wake unsteadiness, synchronized with body oscillations under the lock-in condition. The frequency ratio between the in-line and cross-flow forces is expected due to the symmetry of the system. Higher harmonic components may also emerge. In particular, a third harmonic (3f 1 ) appears in C y spectrum and a fourth harmonic (4f 1 ) is sometimes noted in C x spectrum. The contributions of the force higher harmonic components

The amplitudes of the principal spectral components of the in-line (C x,2 and C x,4 ) and cross-flow (C y,1 and C y,3 ) forces are plotted in Fig. 6, as functions of the reduced velocity, in the three studied cases. In each direction, the evolution of the dominant spectral component amplitude (C x,2 and C y,1 ) is comparable to the evolution of the RMS value of the entire signal (C ′

x and C ′ y , Fig. 4(b,c)). In the in-line direction (Fig. 6(a)), the magnitude of the fourth harmonic component remains small compared to C x,2 . An increase of C x,4 can however be noted for U * = 5 and U * = 6 in the IL+CF case. In the cross-flow direction (Fig. 6(b)), a large-amplitude third harmonic component occurs in the IL+CF case, in the intermediate reduced velocity region: for U * = 6, the third harmonic amplitude reaches 40% of the first harmonic amplitude. The IL case results

for U * = 3 also indicate a significant amplification of the third harmonic when the body oscillates in the in-line direction. On the other hand, C y,3 remains negligible in the CF case. These observations suggest that the large-amplitude third harmonic noted in the IL+CF case may be essentially related to the in-line oscillations of the cylinder; this phenomenon is examined in Section 4.

Wake patterns

Visualizations of wake patterns encountered in the three studied cases are presented in Fig. 7. They consist of iso-contours of the span-and phase-averaged, spanwise vorticity, for U * = 6, i.e. in the region where body responses and fluid forces significantly differ from one case to the other (Figs. 3, 4 and6). Two phases are considered for each case. In each plot, the phase-averaged trajectory of the body is indicated as well as its actual position. In the IL case (Fig. 7(a,b)), the body is almost stationary and the flow exhibits a typical vortex street pattern, similar to that observed when the body is fixed, where single counter-rotating vortices are alternatively shed. Following the nomenclature introduced by [START_REF] Williamson | Vortex formation in the wake of an oscillating cylinder[END_REF], this flow topology can be referred to as a 2S pattern. Wake structure is altered when the body oscillates in the cross-flow direction (Fig. 7(c,d)). A 2P pattern can be identified in this case: two pairs of counter-rotating vortices form per oscillation cycle. This is consistent with prior experimental results reported by [START_REF] Brika | Vortex-induced vibrations of a long flexible circular cylinder[END_REF] and [START_REF] Govardhan | Modes of vortex formation and frequency response of a freely vibrating cylinder[END_REF], in the region of peak oscillation amplitudes. In contrast, a 2S pattern, globally comparable to that noted in the IL case, develops in the IL+CF case, which corroborates the observations of [START_REF] Navrose | Free vibrations of a cylinder: 3-d computations at Re= 1000[END_REF], but departs from the visualizations of [START_REF] Jauvtis | The effect of two degrees of freedom on vortex-induced vibration at low mass and damping[END_REF], due to the different Reynolds numbers, as discussed in [START_REF] Gsell | Two-degree-of-freedom vortex-induced vibrations of a circular cylinder at re = 3900[END_REF].

The topology of the flow may thus vary from one case to the other. However, cases with distinct behaviors, as those observed in the intermediate reduced velocity region, may also exhibit comparable wake patterns: for example large higher harmonic components appear for U * = 6 in the IL+CF case but not in the IL case, while both cases exhibit a 2S wake pattern. The differences noted between the IL+CF case and its one-degree-of-freedom counterparts, and more generally the interactions between in-line and cross-flow forces/responses, do not seem to be systematically connected to wake topology. Some possible mechanisms driving the contrasted behaviors of the one-and two-degree-of-freedom systems are investigated in the next section.

Discussion

The results reported in Section 3 show that the fluid force and body response in a given direction (x or y) may be substantially impacted by the existence of a degree-of-freedom and oscillations in the perpendicular direction (y or x). Three striking features observed in the intermediate reduced velocity region are particularly addressed in this section: (i) in this range of U * , in-line oscillations only occur if cross-flow motion is allowed (Fig. 3(b)); (ii) cross-flow oscillations are also influenced by the in-line degree of freedom in this region, since an increase of the oscillation peak amplitude is noted between cases CF and IL+CF (Fig. 3(a)); (iii) concerning fluid forces, large higher harmonic contributions are observed in the cross-flow direction when the body is subjected to both in-line and cross-flow oscillations, while in the CF case they remain negligible (Fig. 6(b)).

In the following, these features are discussed in light of the evolutions of the two kinematic quantities introduced in Section 2.3, α and Ψ , which relate to body motion. The angle α is the angle between the x axis and the instantaneous oncoming flow velocity (expression (2)). The evolution of the maximum angle α (α m ) as a function of the reduced velocity is plotted in Fig. 8(a), in the three studied cases. As expected, α vanishes in the IL case. Large values of α m are encountered when cross-flow oscillations occur. In both the CF and IL+CF cases, the region of peak α m matches the region of peak oscillation amplitudes (Fig. 3). The values of α m are generally larger in the IL+CF case, where angles up to 0.9 radians (≈ 50 • ) can be noted.

The second kinematic quantity, Ψ , relates to the magnitude of the instantaneous oncoming flow velocity (Ψ = |V in | 2 ). When the cylinder is fixed then Ψ is equal to 1 and it may substantially depart from this value when the body oscillates. A spectral analysis shows that the fluctuating part of Ψ is dominated by its second harmonic component (2f 1 ) so that expression (7) provides a reasonable approximation. The values of Ψ and Ψ 2 are plotted in Fig. 8(b) as functions of U * . In the IL case, Ψ is almost unaltered when the body oscillates (U * = 3). In contrast, it significantly increases in the CF and IL+CF cases, especially in the intermediate reduced velocity region. In the IL+CF case, a maximum increase of 50% is noted compared the fixed cylinder case. Ψ 2 is also altered when the body oscillates. Both in-line and cross-flow oscillations are associated with an increase of Ψ 2 , as shown by the results obtained in the IL and CF cases. The largest values of Ψ 2 are observed in the two-degree-of-freedom case; in the region of peak oscillation amplitudes (U * ≈ 6), Ψ 2 is close to the value of Ψ in the fixed body case (Ψ 2 ≈ 1).

The role of α and Ψ in the alteration of fluid forcing is examined hereafter.

or y) can be written as follows:

C i = D i + L i , (10)
where D i denotes the drag-like contribution to C i (D x = C η cos(α) and D y = C η sin(α)), in reference to its alignment with the instantaneous flow velocity, and L i denotes the lift-like contribution (L x = -C ξ sin(α) and L y = C ξ cos(α)), i.e. perpendicular to the instantaneous flow velocity. In the fixed body case and IL case, α vanishes, therefore C x = D x and C y = L y . Assuming that C i , D i and L i are periodic functions of time, they can be expressed as Fourier series,

C i = ∞ ∑ n=0 C i,n sin(2π nf 1 t + φ C i,n ), ( 11a
)
D i = ∞ ∑ n=0 D i,n sin(2π nf 1 t + φ D i,n ), (11b)
L i = ∞ ∑ n=0 L i,n sin(2π nf 1 t + φ L i,n ), (11c)
where C i,n , D i,n and L i,n are the spectral amplitudes of the nth harmonics, and φ C i,n , φ D i,n and φ L i,n the corresponding phases.

The amplitude of the nth harmonic of C i relates to the nth harmonics of D i and L i as follows:

C i,n = D i,n cos(φ D i,n -φ C i,n ) + L i,n cos(φ L i,n -φ C i,n ). ( 12
)
The results reported in Fig. 4(b) reveal a major amplification of the fluctuating in-line force when the body oscillates in the cross-flow direction. According to (10), this amplification may relate to two effects: an amplification of D x related to body motion, and the appearance of a lift-like contribution (L x ) induced by the angle α. A spectral analysis of C x shows that its fluctuating part is dominated by the second harmonic component. According to (12), the magnitude of this harmonic can be expressed as follows:

C x,2 = D x,2 cos(φ D x,2 -φ C x,2 ) + L x,2 cos(φ L x,2 -φ C x,2 ). (13)
The relative weights of the different terms in (13), issued from the present simulation results, are plotted in Fig. 9, in the CF and IL+CF cases. In the CF case (Fig. 9(a)), the lift-like contribution is clearly dominant, especially in the region of peak amplitudes of C x,2 . For U * = 4 and U * = 5, the drag-like contribution is even negative: it tends to decrease the amplitude of the in-line force fluctuation. At higher reduced velocities, both contributions are positive but the relative contribution of L x to C x,2 remains larger than 70%. This result suggests that the large amplitude of the fluctuating in-line force noted in Fig. 4(b) in the CF case, is mainly related to the emergence of a lift-like contribution associated with the angle α. A comparable behavior is observed in the IL+CF case (Fig. 9(b)). A simple mechanism of interaction between the in-line and cross-flow motions can be proposed: as the body moves in the cross-flow direction, a fluctuating lift-like component emerges and considerably alters the amplitude of Cx , and the resulting in-line vibration, if the body is allowed to move in this direction. This may explain why, in the intermediate range of U * , in-line oscillations appear when cross-flow oscillations occur. Similarly, the cross-flow force coefficient C y can be altered by the emergence of a drag-like contribution related to the angle α. According to expression (12), the amplitude of the first harmonic of C y (which dominates its spectrum) writes

C y,1 = D y,1 cos(φ D y,1 -φ C y,1 ) + L y,1 cos(φ L y,1 -φ C y,1 ). ( 14
)
An analysis of the relative weights of each term in ( 14) shows that the drag-like contribution to C y,1 is negligible.

Influence of Ψ

A simple model has been introduced in Section 2.3 to shed some light on the modulations of the fluid forces expressed in the moving frame, due to the variation of the instantaneous flow velocity magnitude. Following this model, the timeaveraged value of the force coefficient aligned with the instantaneous flow velocity can be expressed as

C η ≈ (1 -κ η )C f x + κ η C f x Ψ + 1 2 κ η C f x,2 cos(φ Ψ 2 -φ C f x,2
)Ψ 2 .

(15)

In ( 15), C η depends on the time-averaged and fluctuating parts of Ψ (Ψ and Ψ 2 ). However, as C f x,2 is generally small compared to C f x , the contribution of the term related to Ψ is expected to dominate, especially when Ψ ≫ Ψ 2 . In this case, the model suggests a linear trend of C η versus Ψ :

C η ≈ C f x + κ η C f x (Ψ -1). ( 16
)
The evolution of C η as a function of Ψ -1, issued from the simulations, is plotted in Fig. 10(a). In this plot, the color of the symbols indicates the relative weights of Ψ and Ψ 2 : high/low values of the ratio Ψ /Ψ 2 are colored in dark blue/white. The The phenomenological model also suggests a trend for the fluid force normal to the instantaneous oncoming flow velocity, C ξ . According to (9), the magnitude of the first harmonic can be expressed as follows:

C ξ ,1 ≈(1 -κ ξ )C f y,1 cos(φ C f y,1 -φ C ξ ,1 ) +κ ξ C f y,1 Ψ cos(φ C f y,1 -φ C ξ ,1 ) + 1 2 κ ξ C f y,1 Ψ 2 sin(φ C ξ ,1 -φ Ψ 2 + φ C f y,1
),

where φ C ξ ,1 denotes the phase of the first harmonic of C ξ . The third term of (17) becomes negligible when Ψ ≫ Ψ 2 . In this case, C ξ ,1 is expected to follow a linear evolution as a function of Ψ ,

C ξ ,1 ≈ C f y,1 cos(φ C f y,1 -φ C ξ ,1 ) + κ ξ C f y,1 cos(φ C f y,1 -φ C ξ ,1 )(Ψ -1). ( 18
)
The results reported in Fig. 11 (19b) where ζ i,n are the spectral amplitudes and φ ζ i,n the corresponding phases, the time-averaged value of Ψ is equal to

Ψ = 1 + 2π 2 f 2 1 ( 4ζ 2 x,2 + ζ 2 y,1
) .

(20)

Fig. 4 .

 4 Fig. 4. Fluid force statistics as functions of the reduced velocity, in the three studied cases: (a) time-averaged in-line force coefficient and RMS values of the fluctuating (b) in-line and (c) cross-flow force coefficients.

x

  in this region is approximately twice larger in the IL+CF case than in the CF case.The evolution of C ′ y versus U * is depicted in Fig.4(c). As shown by the results obtained in the one-degree-of-freedom cases, both in-line and cross-flow oscillations are accompanied by an amplification of Cy . It appears that small-amplitude in-line oscillations can induce a substantial increase of C ′ y (U * = 3 in the IL case). The cross-flow force fluctuation is often amplified in the IL+CF case compared to one-degree-of-freedom cases. Typical time evolutions of the fluid force coefficients are presented in Fig. 5. Selected time series of the force coefficients and their spectral amplitudes based on Fourier transform are plotted for U * = 3 in the IL case and U * = 6 in the CF and IL+CF cases, i.e. in the regions of peak oscillation amplitudes. The spectra show that the fluctuating in-line and cross-flow forces are dominated by frequencies equal to 2f 1 and f 1 , respectively; this is generally the case over the parameter space investigated.

Fig. 9 .

 9 Fig. 9. Drag-and lift-like contributions to the fluctuating in-line force as functions of the reduced velocity: evolutions of the three terms in (13) in the (a) CF and (b) IL+CF cases.

  (a), which are issued from the present simulations, globally confirm this trend for Ψ ≫ Ψ 2 (dark-blue symbols).As shown inFig. 11(b), the amplification of C ξ ,1 is generally accompanied by an increase of C y,1 . Therefore, the modulation of C ξ ,1 by Ψ may play a role in the alteration of the cross-flow response when in-line oscillations occur. The structural responses are often close to harmonic. For harmonic oscillations defined asζ x =ζ x,2 sin(4π f 1 t + φ ζ x,2 ),(19a)ζ y =ζ y,1 sin(2π f 1 t + φ ζ y,1 ),