Kévin Huguenin
email: huguenin@i3s.unice.fr

M João Rendas
email: rendas@i3s.unice.fr

Distributed Adaptive Sampling Using Bounded-Errors

This paper presents a communication/coordination/ processing architecture for distributed adaptive observation of a spatial field using a fleet of autonomous mobile sensors. One of the key difficulties in this context is to design scalable algorithms for incremental fusion of information across platforms robust to what is known as the "rumor problem". Incremental fusion is in general based on a Bayesian approach, and algorithms (e.g. the Covariance Intersection, CI) which propagate consistent characterizations of the estimation error under this challenging situation have been proposed. In this paper, we propose to base inter-sensor fusion on a deterministic approach which considers that bounds to the observation errors are known, wich is intrinsically robust to the rumor problem. We present the equations that enable the determination of the ellipsoidal domain of uncertainty that covers the intersection of the individual sets describing sensor's uncertainty, and show that they solve some pathologies associated to CI. The results presented corroborate a previous claim of the robustness of our control strategy (the criterion used for adaptively choosing the nodes positions) with respect to the conservativeness of fusion methods able to handle rumor .

I. INTRODUCTION

Over the last few years, sensor networks have been progressively recognized as providing an efficient solution for real-world applications in the context of environmental monitoring, enabling rapid and focused observation of spatially extended fields such as temperature or humidity inside a given region of interest. Recent advances in the areas of wireless communications and robotic sensors have made the design of autonomous mobile distributed sensing systems possible. An example of such networks are fleets of underwater vehicles with embedded sensors, micro-processors and wireless communication devices. They enable more powerful and autonomous applications in the sense that, contrary to fixed sensors, this kind of networks do not rely on preliminary computation for sensor placement, allowing proper observation of dynamic environments, and adapting the network geometry to the characteristics of the measured field. Also, increasing computational power available at each node makes decentralized and distributed path control and data fusion possible. Finally, wireless communications allow collaborative behavior and dynamic data exchange making sampling more efficient.

So far, most of the research on sensor networks has focused on sensor placement and data processing for static sensor networks with one or more data fusion centers [START_REF] Arici | Adaptive Sensing for Environment Monitoring using Wireless Sensor Networks[END_REF], [START_REF] Rabbat | Decentralized Compression and Predistribution via Randomized Gossiping[END_REF]. Dealing with autonomous mobile sensors and designing distributed systems brings new challenges. Due to the dynamic nature of the network, full-connectivity can no longer be guaranteed, making traditional data collection/dissemination schemes and centralized path control infeasible. To date there has been little work on adaptive collaborative sensing using mobile sensors. Rahimini et al. propose an algorithm for adaptive sampling using Networked Info-Mechanical System (NIMS) [START_REF] Rahimi | Adaptive Sampling for Environmental Field Estimation using Robotic Sensors[END_REF] but they restrict to the independent control of each sensor and do not take into account energy and communication constraints since NIMS are power supplied and use a wired communication medium. In the framework of underwater applications, the concept of adaptive sampling using a fleet of autonomous sensors has been advocated for some years by the members of the ASAP (Autonomous Sampling and Prediction) project, see for instance [START_REF] Leonard | Collective Motion, Sensor Networks and Ocean Sampling[END_REF]. However, this project considers a centralized control of the underwater fleet, eliminating most of the problems and difficulties that we address in this paper.

The ultimate goal of an autonomous sensing network is to gather data which collectively allows the estimation of a given spatial field over an area of interest A, this estimate being obtained once data gathering is finished, by globally fusing the data sets collected by the nodes. The key point of such an architecture is the collaboration between nodes to determine the sensing points where the measurements will jointly bring the most information, which requires that nodes share their knowledge about the field. In [START_REF] Huguenin | An Architecture for Distributed Collaborative Adaptive Sensing[END_REF], [START_REF]An Architecture for Distributed Collaborative Adaptive Sensing[END_REF], we presented an architecture integrating communication, control and information processing algorithms that implements a distributed and decentralized strategy for acquiring an estimate of a spatially extended field, based on a probabilistic Bayesian framework, under strong communication constraints. While iterative integration of additional sensor measurements can be easily done in the context of Bayesian estimation theory, the fusion of local estimates is prone to the rumor problem: there is no simple way to distinguish fresh information conveyed by a distinct node from re-circulation, through the network, of a node's information back to itself. In these conditions a naive application of the Bayesian recursion would lead to overconfident error characterizations. This problem is overcome by accepting some degree of information loss, and using fusion methods that propagate upper bounds on uncertainty (conservative error estimates), like the Covariance Intersection (CI) algorithm [START_REF] Julier | A Non-divergent Estimation Algorithm in the Presence of Unknown Correlations[END_REF]. In [START_REF] Huguenin | An Architecture for Distributed Collaborative Adaptive Sensing[END_REF] we presented a new criterion which prevents fresh measurements from being discarded by CI by choosing sensing positions that correspond to different observation structures, and we compared the performance of our distributed sensing architecture using CI and naive Bayes. Our analysis shows the superiority of CI and that our criterion is well adapted to CI.

The entire issue of "rumor propagation" is eliminated if a deterministic, error-bounded [START_REF] Fogel | On the Value of Information in System Estimation-Bounded Noise Case[END_REF] approach to the fusion problem, which propagates upper bounds to the sets of parameters values compatible with the acquired observations, is used. Bounded-error methods have been popular in control and estimation problems, mostly for recursive state estimation problems, and update equations for the set dilation and intersection steps have been published. The problem of distributed fusion requires that two or more uncertainty domains be directly fused. In this paper we present equations that allow the fusion of two ellipsoids that are a special case of the general problem solved in [START_REF] Ros | An ellipsoidal calculus based on propagation and fusion[END_REF]. Besides providing a better solution to the fusion problem, as we show in section VI, the use of a boundederror based methodology for fusion, based on a completely geometric characterization of the uncertainty domain affecting the current estimates at each node, is a natural match to our adaptive sensing criterion, giving a stronger coherency to our overall architecture.

The paper is organized as follows. Section II outlines the problem of distributed collaborative sensing and presents an overview of the global architecture proposed. In Section V we briefly present and justify the adaptive sensing criterion used. Sections III and IV are dedicated to recursive estimation and fusion, focusing in the presentation of the update equations for the bounded-error methods, and pointing out their analogies to Gauss-based Bayesian methodologies. In Sections VI and VII we present simulation results and summarize our conclusions and some ideas for future work.

II. PROBLEM STATEMENT: GENERAL ARCHITECTURE

Our goal is to implement a decentralized strategy for controlling the sensors positions X t = {X tc } t0≤tc≤t such that the information collectively gathered by the network optimizes some aggregated statistic functional of the (error of the) reconstructed field, f (x; X t )

C(X t ) = E A g f (x, y), f (x, y; X t ) dxdy (1)
where E denotes the (ensemble) expectation operator. We consider a set of working assumptions, which reflect constraints typical of ocean-observation problems:

• the observed field f (z) is scalar and assumed static during the time interval required to perform its observation inside a known region, A ⊂ R 2 . • the robots are able to localize themselves (GPS, network of acoustic beacons, distributed location algorithm, . . . ). Their speed is upper-bounded v ≤ v max and they can (and must, to perform a measurement) station-keep.

• error-free fixed radius communications: two nodes can communicate with each other without errors provided that the distance between them is lower than the communication range R com . We present below the major components of our distributed sensing architecture [START_REF] Huguenin | An Architecture for Distributed Collaborative Adaptive Sensing[END_REF]. Our solution involves three main stages that are executed in a cyclic manner: (i) sensing, (ii) communicating and data fusion, and (iii) collaborating to decide on the next network deployment. While sensing is performed independently at each node, choosing the next sensing position can be drastically improved by sharing information across nodes. The second task is thus denoted by fusion insofar as it consists in fusing the estimates stored at each node to compute a more informative one while the third one is called adaptive collaborative observation in the sense that the nodes rely on the exchanged information to determine the set of future sensing points which jointly bring the most additional information about the observed field.

Network partitioning

Obviously, the larger the number of nodes that collaborate to select the future sensing positions, the more their individual decisions will reflect the global state of knowledge of the network. However, under low spatial node density and limited communication range, effective network connectivity may be very limited. In particular, the communication graph may be disconnected, with communication being impossible between some network nodes. Our architecture is based on the dynamic definition of logical groups of collaborating nodes, designed to allow easy information exchange amongst their members. Logical groups (clusters) are defined at the beginning of each cycle, by assigning a virtual group number to each robot and then restricting collaboration to nodes belonging to the same cluster.

Distributed clustering is a well-known problem and has been addressed under the single-hop structural constraint [START_REF] Basagni | Distributed Clustering for Ad Hoc Networks[END_REF]. By imposing direct communication between cluster members, the nodes of our network are released from the communication problems and requirements characteristics of Mobile Ad-Hoc Networks (MANET) such as intermittent connectivity and delay tolerant routing. A "Distributed Clustering Algorithm" (DCA) which dynamically creates and destroys these clusters using only gossip messages has been proposed in [START_REF] Huguenin | An Architecture for Distributed Collaborative Adaptive Sensing[END_REF] Distributed data fusion Data fusion consists in updating the estimate of the observed field at each node using either new measurements (we call this operation update) or estimates received from other nodes (that we designate by fusion). Thus, clustering must be immediately followed by data sharing and fusion, that need to occur before the sensing positions are chosen in order to ensure consistency of this last operation.

In [START_REF]An Architecture for Distributed Collaborative Adaptive Sensing[END_REF], [START_REF] Huguenin | An Architecture for Distributed Collaborative Adaptive Sensing[END_REF] we compared two Bayesian methods for update and fusion ("naive Bayes" and "Covariance Intersection") and demonstrated the superiority of CI. The fusion operation is the main new contribution of this paper, which proposes to assess it in a distinct uncertainty modeling framework (deterministic, bounded-error methods), see Section III.

Choice of future positions

Selection of the sensing points is the core problem of adaptive field observation. Its potential impact is evident, since without collaboration nearby nodes will have a similar local "view" of the observed field, and will most likely decide to sense at the same or nearby locations. The criterion presented in [START_REF]An Architecture for Distributed Collaborative Adaptive Sensing[END_REF], [START_REF] Huguenin | An Architecture for Distributed Collaborative Adaptive Sensing[END_REF] to choose the next set of sensing points tries to reach a balance between (i) maximizing the new information jointly brought by the measurements at the target sensing points; and (ii) imposing geographic dispersion of the nodes, to increase the probability they establish communication with nodes from distinct clusters.

Search of the best set of sensing points is performed within a bounded region around the nodes, called action zone using a greedy distributed approach: sensing points are determined one-by-one, taking into account the previously selected ones. Thus, the computational power required to determine the sensing points grows linearly with the number of nodes in the clusters and so does the available computational power. Once the set of the sensing positions has been chosen, each of them is assigned to one node such that the the total amount of energy (this can be replaced by a different criterion) used by the robots to reach them is minimized, see [START_REF] Huguenin | An Architecture for Distributed Collaborative Adaptive Sensing[END_REF] for details.

III. RECURSIVE ESTIMATION

As we stated in the Introduction, scalability constraints forbid storing the entire set of raw data acquired by each node. In this section we briefly present (A) the fundamental equations of Bayesian estimation theory that propagate a sufficient statistic summarizing all relevant information about the estimated field contained in the past observed data, as well as (B) a deterministic estimation method based on feasibleparameter set outer-bounding, both allowing recursive update.

We assume that the sensor observations are noisy versions of the field's values, such that

z k = f (x, y) + k , (x, y) ∈ A,

A. Bayesian estimation

In the context of the probabilistic Bayesian estimation, we assume that a parametric model for the field is known,

f (x, y) = f (x, y; a), a ∼ p a : unknown parameter ,
where p a is the prior distribution of a, such that the task of estimating the field is formally equivalent to the estimation of a. If the noise components affecting distinct measures are statistically independent, then the posterior density p(a|z k ) of a given observed data z k = {z i } k i=1 is a sufficient statistic for the estimation of the field f , and it can be recursively updated [START_REF] Trees | Detection, Estimation, and Modulation Theory, Part I[END_REF] using the Bayes relation. For finite-dimension linear models,

f (x, y; a) = a T f (x, y), f (x, y) = [f 1 (x, y), . . . , f I (x, y)] T ,
where {f i } I i=1 are a set of known functions, and Gaussian distributions (both prior and observation noise), the posterior densities are normal, a|z k ∼ N (µ k , Σ k ), and closed-form expressions can be obtained for their means and variances:

Σ -1 k+1 = Σ -1 k + 1 σ 2 f k+1 f k+1 T (2)
Σ -1 k+1 µ k+1 = Σ -1 k µ k + 1 σ 2 z k+1 f k+1 (3)
In this equation σ 2 is the variance of the observation noise. This version of the Bayes equations are known as the "information form" of the least-squares estimator, and allow asynchronous integration of new measures. The statistics

η k = (Σ -1 k , Σ -1 k µ k
) maintained at each node all have the same form, and their size does not grow with the number of measurements. Since matrix addition is associative and commutative, the fusion equations presented above can be easily extended to several measures.

Note that under the linear-Gauss assumption, the covariance matrix of the posterior density does not depend on the measurement itself (z k ) but only on the local observation operator (f (x k , y k )). Since the amount of information conveyed by z k about a is a function of Σ k , the impact of future measurements on the quality of the estimate can be predicted.

B. Bounded-error

The estimation method presented above relies on knowledge of the statistical properties of the noise, which is a rather strong assumption. Fogel and Huang proposed in their seminal paper [START_REF] Fogel | On the Value of Information in System Estimation-Bounded Noise Case[END_REF] a deterministic approach to parameter estimation under a bounded error assumption. The key idea is that, assuming observation noise is bounded, each measurement defines a set of feasible parameters values (i.e., which are consistent with the measurement, the model and the bound). The actual value of the parameter lies in the intersection of the sets associated to the measures. Several representations for the so-called feasible parameter set have been proposed, including parallelotopes [START_REF] Vicino | Sequential Approximation of Feasible Parameters for Identification with Set Membership Uncertainty[END_REF], polytopes [START_REF] Walter | Exact Recursive Polyhedral Description of the Feasible Parameter Set for Bounded-Error Models[END_REF] and ellipsoids [START_REF] Fogel | On the Value of Information in System Estimation-Bounded Noise Case[END_REF]. Using polytopes gives the optimal solution but the state maintained at each node may potentially grow as the number of measurements, violating the scalability requirement. On the contrary, using ellipsoids leads to a scalable solution (the size of the estimate does not grow with k and is O(I 2 )) and is directly comparable to the the statistical Bayesian approach with Gaussian posteriors. In the following we use ellipsoids to describe the feasible parameter sets:

E k = a : (a -a c k ) T P -1 k (a -a c k ) ≤ 1 and âk = a c k (4)
Assume the measurement errors k are bounded by r:

| k (x, y)| ≤ r.
Then, each new measurement z k defines two

H - k H + k Ek Ek+1 Figure 2.
Minimum-determinant ellipse of the form E k+1 (q) [START_REF] Huguenin | An Architecture for Distributed Collaborative Adaptive Sensing[END_REF] hyperplanes that bound the set of feasible parameters:

H + k = a : 1 r z k -a T f k ≤ 1 and H - k = a : 1 r z k -a T f k ≥ -1
The updated feasible parameter set E k+1 must be of the form (4) and enclose

E k ∩ H + k+1 ∩ H - k+1 .
A simple one-parameter family of sets verifying these properties is [START_REF] Fogel | On the Value of Information in System Estimation-Bounded Noise Case[END_REF]:

E k+1 (q) = a : (a -a c k ) T P -1 k (a -a c k ) + q 1 r 2 z k -a T f k 2 ≤ 1 + q , q ≥ 0 (5)
This ellipsoid can be rewritten in the form (4) using the following recursive updating rules:

P k+1 = ζ k Z k (6) a c k+1 = Z k P -1 k a c k + qr -2 z k f k (7)
Z k = P -1 k + qr -2 f k f k T -1 (8)
ζ k = 1 + q - qr -2 1 + qr -2 G k e 2 k ( 9
)
where

e k = z k -a k T f k , G k = f k T P k f k ,
and q ≥ 0. The value of q is chosen to minimize some optimization criterion. In the traditional case of volume minimization (minimizing det P k+1 ), the optimal value q k of q is the -positive if any, or zero otherwise -root of a second order polynomial and can be expressed in closed form. Finally, E k+1 = E k+1 (q k ).

Contrary to the Bayesian approach, the bounded-error (BE) method enables optimization of the feasible parameter set (i.e, the ellipsoid) according to any criterion -such as the determinant, the trace or the radial spectrum of P k -while integrating new measurements. On the other hand, due to its deterministic nature, a new measurement may not convey any additional information. This case occurs when the hyperplane H + k+1 and H - k+1 do not intersect the ellipsoid E k . In the case where only one hyperplane intersects the ellipsoid, a correction has been proposed so that the one-parameter family of ellipsoids (5) encloses more tightly the intersection [START_REF] Fogel | On the Value of Information in System Estimation-Bounded Noise Case[END_REF].

IV. DISTRIBUTED DATA FUSION

The fusion step is the weak link of collaborative sampling architectures. In this section we present Covariance Intersection (CI), the most frequently used fusion algorithm under unknown correlation, and illustrate its performance on simple examples. Also, we highlight pathological behaviors of CI on some special cases. Based on purely geometric considerations, we extend the bounded-error estimation to a fusion method which does not suffer from pathologic cases aforementioned.

A. Covariance Intersection (CI)

For linear-Gaussian problems, fusion is a linear operation. Imposing unbiasedness of the estimator resulting of the fusion of unbiased estimates implies that the linear combination be convex, and finding the fusion operator consists in choosing an element in a I 2 -dimensional space. In their original paper [START_REF] Julier | A Non-divergent Estimation Algorithm in the Presence of Unknown Correlations[END_REF], Julier and Uhlmann proposed CI which reduces the aforementioned problem to a one-dimensional optimization problem, using a convex linear combination:

Σ -1 C = ωΣ -1 A + (1 -ω)Σ -1 B (10) Σ -1 C µ C = ωΣ -1 A µ A + (1 -ω)Σ -1 B µ B ( 11
)
CI has been proved to provide consistent estimates (i.e., Σ C ≥ E (µa)(µa) T ) for any ω ∈ [0, 1] meaning that the free parameter ω can be chosen to optimize some criterion. The consistency of the estimator resulting from CI can be interpreted as a conservativeness property since the fusion method never "creates" virtual information -as the Bayesian method may do -when fusing correlated estimates. However, CI may discard fresh information, considering them as redundant. This is revealed by the fact that the centers of the ellipsoids are not taken into account in the determination of the covariance Σ C : suppose that two nodes have two estimates with the same covariance matrix but different centers as illustrated in Figure 3. Then, CI returns the same covariance matrix for any value of ω leaving thhe center (i.e., the estimate â of a) indetermined.

B. Feasible Parameters Set Intersection

Although CI behaves strangely in some particular cases and ignores the precious information conveyed by the centers of the ellipsoids, its core idea of intersecting uncertainty ellipsoids can be adapted to a deterministic context: let E 1 and E 2 be two estimates of a, then any set enclosing E 1 ∩ E 2 contains the actual value of a, and is thus a valid feasible parameter set. To be directly and fairly comparable to CI, we parametrize our search space with only one parameter. Inspired by both ( 5) and [START_REF] Trees | Detection, Estimation, and Modulation Theory, Part I[END_REF], we naturally define the search space as:

E(ω) = a : (1 -ω) (a -a c 1 ) T P -1 1 (a -a c 1 ) + ω (a -a c 2 ) T P -1 2 (a -a c 2 ) ≤ 1 , 0 ≤ ω ≤ 1 (12)
yielding a set of ellipsoids given by [START_REF] Ros | An ellipsoidal calculus based on propagation and fusion[END_REF] (see [START_REF]An Architecture for Distributed Collaborative Adaptive Sensing[END_REF] for a detailed proof):

P 0 = (1 -ζ 0 )Z 0 (13) a c 0 = Z 0 (1 -ω)P -1 1 a c 1 + ωP -1 2 a c 2 (14)
Z 0 = (1 -ω)P -1 1 + ωP -1 2 -1 (15)
ζ 0 = (1 -ω) (a c 1 -a c 0 ) T P -1 1 (a c 1 -a c 0 ) - (16
)

ω (a c 2 -a c 0 ) T P -1 2 (a c 2 -a c 0 )
Figure 4 shows an example of the minimum-determinant ellipsoid enclosing the intersection of two uncertainty sets (Figure 4(a)) and in the case where the uncertainty sets are the same (Figure 4(b)). One may notice that contrary to the CI method, the estimates â are taken into account yielding a smaller volume of uncertainty and a unique new estimate.

E1 E2 E0 (a)

V. COLLABORATIVE ADAPTIVE SENSING

As stated in [START_REF]An Architecture for Distributed Collaborative Adaptive Sensing[END_REF], the key points of distributed estimation are (i) collaboration between nodes and (ii) adaptation to the sensed field. The first aims at maximizing the joint information conveyed by the measurements performed by the nodes while the latter uses a model for field so as to determine the points that will bring the most information, taking into account the current state of knowledge. In both the deterministic and probabilistic frameworks, the confidence associated to an estimate is described by an ellipsoid. For vector estimates (a is not a scalar) the uncertainty set obtained by integration of additional data sampled does not allow the definition of a complete order relation over A that can be used to select where the next measures should be made. The use of a scalar function of the matrix M describing the uncertainty set (the positive definite matrices P or Σ depending on whether CI or BE is used) -such as trace, determinant or spectral radius -can be used to define a total order relation on which the choice of the next measuring points can be based. These metrics have a direct geometric interpretation, making the analysis of sensing points determination method easier to apprehend.

The most common criterion is det (M ) which is a increasing function of the uncertainty volume of the uncertainty ellipsoid. It can be shown that maximizing the volume is equivalent to maximizing f T M f , which can be interpreted as choosing the sensing point where f is closest to the direction along which âk is the most uncertain. In the deterministic case, choosing the observation direction collinear to the longest axis (eigenvector associated to the largest eigenvalue) of M is the optimal choice, since "cutting" the ellipsoid with hyperplane orthogonal to the observation vector leads to the largest possible reduction of the uncertainty volume and deformation.

Based on a different line of reasoning, a structural algorithm was proposed for choosing the sensing positions. The algorithm uses a orthonormal basis B, which is iteratively updated starting with B = B 0 . At each stage the new sensing point is determined by maximizing the projection of f on Span (B)

⊥ and B is updated by a Gramm-Schmidt iterative orthonormalization algorithm:

x i , y i = argmax x,y∈ G (I -B i-1 B T i-1 )f (x, y) , and
B i = B i-1 ∪ (I -B i-1 B T i-1 )f (x i , y i ) (I -B i-1 B T i-1 )f (x i , y i )
.

The choice of B 0 determines the places where the nodes should not go and is chosen with two goals:

• impose motion to a different position:

B 0 = f (x i k , y i k ) N i=1
• avoiding the directions along which a is best known:

B 0 = {v 1 , . . . , v I-2N }
where v i is the eigenvector associated to the i-th smallest eigenvalue (shortest axis of the uncertainty ellipse) of Σ. The initial basis is obtained by combining B 0 with B 0 .

VI. RESULTS

In this section, we study the performance of the estimation framework presented in this paper using computer simulations. The field to be estimated is composed of 20 Gaussian shaped basis functions with unit amplitude, and random scale and localization. The coefficients α i are drawn from a uniform distribution. We consider a network of 6 robots, initially deployed at random inside a small launching zone in the region of interest A (of size 10,000 × 10,000) -reflecting a realistic deployment scenario. The communication radius and the diameter of the action zone are respectively 2,000 and 2,500. We conducted a Monte-Carlo study, performing two sets of 25 independent simulations of the architecture keeping both the field and the launching zone fixed. One set uses Bayesian estimation together with CI for the fusion step, and the other uses a deterministic approach to both estimation and fusion.

We use two related criteria to evaluate the performance of our architecture: (i) the error on the field estimate f and (ii) the error on the globally reconstructed parameter itself, â (i.e., the fusion of the set of estimates over the whole network): â = fuse (â 1 , . . . , âN ) and f (x, y) = âT f (x, y)

where the fusion of the estimates is performed iteratively.

The first is motivated by the original goal of providing a global picture of the field, as mentioned in the Introduction. We use distortion as a metric to evaluate the quality of the field estimates using equation (1) together with the L 1 norm.

However, this criterion does not directly reflect the precision of the estimate along each component of a, which may -in the underwater framework -be related to interesting natural phenomena. In fact, the strategy that drives the network (the determination of the sensing points) has been designed to reduce the deformation of the uncertainty ellipsoid associated to a (trying to avoid that some of its components be much better known than other). We plot these two criteria as a function of time, observing the rates of convergence and the asymptotic behaviors of the two different estimation/fusion frameworks. The performance improvement of the collaborative architecture with respect to independent node operation has already been demonstrated in [START_REF] Huguenin | An Architecture for Distributed Collaborative Adaptive Sensing[END_REF]. Figure 5 presents the mean square-error on the estimate of a together with the distortion of the field globally reconstructed over time. These plots support our two major conclusions: (i) the BE fusion method is superior to CI leading to better estimates of a (see 5). (ii) the network control strategy based on the geometric criterion is robust to poor error characterization of the fusion methods, as the comparable overall performance of the CI and BE demonstrates. In our architecture, what determines the network dynamics is the subspace generated by the directions of current smallest uncertainty, irrespective of the detailed amount of uncertainty along each one of them. This means that many different matrices will lead to the same orthonormal basis B, and thus to the same set of chosen sampling points which explains the similar performances for the different fusion methods observed in Figure 5(b).

VII. CONCLUSIONS AND FUTURE WORK

In this paper we extended the overall distributed algorithm for collaborative adaptive sensing proposed in [START_REF] Huguenin | An Architecture for Distributed Collaborative Adaptive Sensing[END_REF] proposing a deterministic approach to parameter estimation and fusion. Assuming bounded-errors, and bounding the uncertainty domains by ellipsoids, provided a comprehensive geometric view of the information processing aspects of the architecture (estimation, fusion and adaptive determination of the sensing points). All blocks of the architecture now rely on the same formalism, becoming matched to each other, and improving the architecture consistency and coherency. Besides this overall coherency, the deterministic setting provided an inherent robustness with respect to one of the difficult problems affecting distributed fusion (the rumor problem). It also improved estimation when compared to the traditional Bayesian estimation with CI as fusion method, as the computer simulations demonstrated.

More generally, design of efficient distributed fusion methods, together with robust criteria for sampling site selection, are critical steps for being able to fully exploit the potential advantages of the flexible and powerful framework of mobile sensor networks, providing effective answers to the difficulties resulting from the stringent constraints of these systems.

Since our architecture is built at the application level, the communication aspects are idealized (fixed radius and error-free communications with unlimited bandwidth). For real-world implementation, these aspects must be taken into account, especially in underwater environment where the communication conditions are harsh.

Our work can be improved in several aspects. The final fusion is done off-line with no computational power limitation or time restrictions, and thus may be performed more efficiently than simply fusing iteratively (pairwise) the node estimates. The same problem affects fusion inside each cluster, but under stronger constraints in terms of memory and computational power. Designing a fusion method that directly assesses the problem of fusing multiple ellipsoids is a point for future study. Also, looking for the subspace that results in the largest observation rank instead of iteratively optimizing a geometric criterion should lead to a more efficient set of sampling points.

Figure 1 .

 1 Figure 1. Coarse Block-based overview of execution cycle of the Distributed Algorithm for Collaborative Adaptive Sensing

Figure 3 .

 3 Figure 3. Fusion of two estimates with the same covariance using (a) Naïve Bayes (NB) (b) Covariance Intersection (CI)

Figure 4 .

 4 Figure 4. (a) Sample minimum-determinant ellipse of the form E(ω) (12), (b) fusion of two estimates with the same covariance using ellipsoid intersection

Figure 5 .

 5 Figure 5. Metrics on fusion of the collection of network estimates over time: (a) a -â 2 ; (b) L 1 norm.


Maria Adela Puscas
email: maria-adela.puscas@cea.fr

Romain Lagrange
email: romain.lagrange@cea.fr

Interaction of two cylinders immersed in a viscous fluid. On the effect of moderate Keulegan-Carpenter numbers on the fluid forces

Keywords: Fluid-structure interaction, Keulegan-Carpenter number, Nonlinear effects, Stokes number, Viscous effects, Self-added coefficients, ALE method

This work deals with the hydrodynamic interaction of two parallel circular cylinders, with identical radii, immersed in a viscous fluid initially at rest. One cylinder is stationary while the other one is imposed a harmonic motion with a moderate amplitude of vibration. The direction of motion is parallel to the line joining the centers of the two cylinders. The two dimensional fluid-structure problem is numerically handled by the Arbitrary Lagrangian-Eulerian method implemented in the open-source CFD code TrioCFD. First, we show that the fluid forces on the two cylinders are aligned with the direction of the imposed motion. Second, we show that the moderate oscillations of the moving cylinder create nonlinear effects in the fluid that strongly affect the characteristics (Fourier harmonics) of the hydrodynamic force acting on the stationary cylinder. The fluid force on the moving cylinder is shown to be poorly affected by the nonlinear effects, which makes it possible to extend the linear concept of self-added mass and damping coefficients. First, we show that the self-added coefficients decrease as Sk -1/2 , with Sk the Stokes number (dimensionless number constructed from the imposed vibration frequency). Second, we show that the self-added mass (resp. damping) decreases (resp. increases) as -KC 3 (resp. +KC 3 ), with KC the Keulegan-Carpenter number (ratio between the imposed amplitude vibration and the separation distance between the cylinders). These variations are included in new power laws derived from nonlinear regressions of the numerical results. These new power laws for the self-added coefficients combine the effect of both Sk and KC, covering the viscous (Sk ≥ 500) and weakly nonlinear (KC ≤ 0.3) regimes.

Introduction

The analysis of forced vibrations of structures immersed in flowing or quiescent fluids is of fundamental importance in many scientific and engineering fields, spanning a wide range of applications, from the nuclear industry and the flow-induced vibrations of heat exchanger tubes [START_REF] Chen | Vibration of nuclear fuel bundles[END_REF][START_REF] Chen | Dynamics of heat exchanger tube banks[END_REF][START_REF] Axisa | Random excitation of heat exchangertubes by cross-flows[END_REF], the dynamics of trees and plants [START_REF] De Langre | Effects of wind on plants[END_REF], to the energy harvesting of flexible structures [START_REF] Michelin | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF][START_REF] Virot | Coupling between a flag and a spring-mass oscillator[END_REF][START_REF] Eloy | Aeroelastic instability of cantilevered flexible plates in uniform flow[END_REF]. Since the seminal works of [START_REF] Poisson | Sur les mouvements simultanés d'un pendule et de l'air environnant[END_REF][START_REF] Green | Researches on the vibration of pendulums in fluid media[END_REF][START_REF] Stokes | On some cases of fluid motion[END_REF], a considerable amount of theoretical, numerical and experimental studies have been conducted, considering single or multiple immersed structures, rigid or flexible, with various geometries and kinematics (imposed or fluid coupled), and different fluid conditions or assumptions, see [START_REF] Chen | Flow-induced vibration of circular cylindrical structures[END_REF][START_REF] Blevins | Flow-induced vibration[END_REF][START_REF] Paidoussis | Fluid-Structure Interactions: Cross-Flow-Induced Instabilities[END_REF][START_REF] Nakamura | Flow-induced vibrations: classifications and lessons from practical experiences[END_REF][START_REF] Bazilevs | Computational fluid-structure interaction: methods and applications[END_REF][START_REF] Sigrist | Fluid-structure interaction: an introduction to finite element coupling[END_REF] for an extensive review.

For small oscillations of an immersed body, the nonlinear convective acceleration of the Navier-Stokes equation is small compared with the unsteady acceleration, see [START_REF] Lagrange | A new analytical approach for modeling the added mass and hydrodynamic interaction of two cylinders subjected to large motions in a potential stagnant fluid[END_REF]. The fluid forces are linear combinations of the velocity and the acceleration of the immersed body. The coefficients entering in these linear combinations are referred to as the added-mass term and the added-damping term. These coefficients are sensitive to the geometry of the moving body, the number of immersed structures, the degree of confinement of the problem, the viscosity of the fluid, and the characteristics of the motion (for e.g., the frequency of the motion, or in the case of a flexible structure the shape of its vibration modes, see [START_REF] Lagrange | Hydrodynamic interaction between two flexible finite length coaxial cylinders: new theoretical formulation and numerical validation[END_REF]). In many practical cases, the Stokes number, defined as the ratio of the unsteady acceleration of the Navier-Stokes equations versus the viscous term, is large enough to consider the fluid forces to be purely inertial. In such cases, the added-damping coefficient is disregarded, whereas the added-mass coefficient can be obtained from a potential theory, solving the Laplace equation for a fluid potential, with given boundary conditions on the surface of the body and on the limits of the fluid domain. The reader is referred to the books of [START_REF] Patton | Tables of hydrodynamic mass factors for translational motion[END_REF]20] for a detailed presentation on the concept of the added-mass coefficient and the related calculation in different configurations.

In this paper, we consider the two dimensional problem of two circular cylinders immersed in a quiescent viscous fluid. One cylinder is stationary while the other is imposed a harmonic motion in the direction of the line connecting the center of the two cylinders. We expand on the findings of [START_REF] Lagrange | A new analytical approach for modeling the added mass and hydrodynamic interaction of two cylinders subjected to large motions in a potential stagnant fluid[END_REF][START_REF] Lagrange | New estimations of the added mass and damping of two cylinders vibrating in a viscous fluid, from theoretical and numerical approaches[END_REF] to explore the effect of finite amplitude vibrations on hydrodynamic loading. Specifically, we depart from small structural oscillations and linear fluid-structure interaction to consider finite-amplitude structural vibrations with Keulegan-Carpenter number up to 30%. The Keulegan-Carpenter number is defined as the ratio of the vibration amplitude versus the distance between the two cylinders and is used to measure the relative effect of convective acceleration with respect to local acceleration. Here, we consider Stokes numbers from 500 to 2000 and Keulegan-Carpenter numbers from 0.005 to 0.3. In this regime, viscous effects and convection-driven nonlinearities in the fluidstructure interaction are expected to be not negligible, as demonstrated in the literature on rigid cylinders in oscillatory flows, see [START_REF] Gabbai | An overview of modeling and experiments of vortex-induced vibration of circular cylinders[END_REF][START_REF] Sarpkaya | Force on a circular cylinder in viscous oscillatory flow at low Keulegan-Carpenter numbers[END_REF][START_REF] Sarpkaya | Hydrodynamic damping, flow-induced oscillations, and biharmonic response[END_REF], and recent works on vibrations of slender beams in unbounded viscous fluid, see [START_REF] Aureli | Low frequency and large amplitude oscillations of cantilevers in viscous fluids[END_REF][START_REF] Aureli | Nonlinear finite amplitude vibrations of sharp-edged beams in viscous fluids[END_REF]. The fluid-structure problem is numerically handled by the Arbitrary Lagrangian-Eulerian (ALE) method [START_REF] Donea | Arbitrary Lagrangian-Eulerian Methods[END_REF] implemented in the open-source CFD code, TrioCFD, see [START_REF] Angeli | Overview of the TrioCFD code: Main features, V&V procedures and typical applications to nuclear engineering[END_REF][START_REF] Angeli | FVCA8 benchmark for the Stokes and Navier-Stokes equations with the TrioCFD code -benchmark session[END_REF][START_REF] Panunzio | FSI-Vibrations of immersed cylinders. Simulations with the engineering open-source code TrioCFD. Test cases and experimental comparisons[END_REF].

The paper is organized as follows. § 2 presents the problem and the governing equations for two circular cylinders immersed in a viscous fluid initially at rest. The numerical approach used to solve the fluid-structure problem is briefly described in this section. § 3 presents the results of the numerical simulations for the fluid forces when one cylinder is stationary, and the other is imposed harmonic oscillations with different frequencies and increasing amplitudes. In § 4, nonlinear regressions of the numerical results are performed to extend the linear concept of self-added coefficients to the weakly nonlinear regime of moderate vibration amplitudes. Finally, some conclusions are drawn in § 5.

Definition of the problem and numerical approach

Presentation of the problem

We consider the two-dimensional problem of two parallel circular cylinders C j , radius R, boundaries ∂C j , immersed in a viscous stagnant fluid, characterized by its volume mass density ρ and its kinematic viscosity ν. The cylinders are separated by a distance E, as illustrated in Fig. 1. The cylinder C 1 remains stationary while C 2 is imposed a harmonic displacement U = Q sin (Ωτ ) e x , with Q the amplitude, τ the time, Ω the angular frequency and e x the direction of motion, parallel to the line joining the cylinder centers. The oscillations of C 2 generate an incompressible fluid flow (V,P) governed by the Navier-Stokes equations

∇ • V = 0, (1a) ∂V ∂τ + (V • ∇) V + 1 ρ ∇P -ν∆V = 0, (1b)
V = 0 on ∂C 1 , (1c)
V - dU dτ = 0 on ∂C 2 . ( 1d
)
The linear fluid force acting on C j is the sum of a pressure and a viscous term and writes

F j = - ∂Cj P ndL + ρν ∂Cj ∇V + (∇V) T • ndL, (2)
with n the outward normal unit vector to ∂C j , (∇V) T the transposate tensor of ∇V and dL an infinitesimal line element of ∂C j .

The dimensionless parameters of this study are

KC = Q E , Sk = R 2 Ω ν , ε = E R , (3)
known as the Keulegan-Carpenter number, the Stokes number, and the dimensionless separation distance, respectively. In what follows, we note t = Ωτ the dimensionless time, f j = F j / ρQR 2 Ω 2 the dimensionless fluid force and f i its x-component.

Numerical approach

To solve Eq. (1), numerical simulations are performed with the open-source code TrioCFD, coupled with an Arbitrary Lagrangian-Eulerian (ALE) module. The reader is referred to [START_REF] Angeli | Overview of the TrioCFD code: Main features, V&V procedures and typical applications to nuclear engineering[END_REF][START_REF] Angeli | FVCA8 benchmark for the Stokes and Navier-Stokes equations with the TrioCFD code -benchmark session[END_REF][START_REF] Panunzio | FSI-Vibrations of immersed cylinders. Simulations with the engineering open-source code TrioCFD. Test cases and experimental comparisons[END_REF][START_REF] Donea | Arbitrary Lagrangian-Eulerian Methods[END_REF] for extensive details on the CFD code and the ALE module.

In the ALE approach, the fluid flow is computed in a domain that is deformed in order to follow the movement of the fluid-solid interface. The approach treats the mesh as a frame that moves with the arbitrary velocity V ALE . For moderate deformations, V ALE is usually defined as the solution of an auxiliary Laplace problem

C 1 • C 2 • E R R U L x L y Ω f luid
∆V ALE = 0, (4a) V ALE = 0 on ∂Ω f luid \ ∂C 2 , (4b)
V ALE - dU dτ = 0 on ∂C 2 . ( 4c
)
from which the kinematics of the mesh grid is updated. In this new frame of reference, the Navier-Stokes equation (1b) writes

∂JV ∂τ + J (V • ∇) V -(V ALE • ∇) V + 1 ρ ∇P -ν∆V = 0, (5a)
with J the Jacobian of the transformation between the ALE and the Lagrange descriptions.

In our simulations, a first-order backward Euler scheme is used to time discretize the Navier-Stokes equations and the space discretization is based on the hybrid Finite Element-Volume method [START_REF] Angeli | FVCA8 benchmark for the Stokes and Navier-Stokes equations with the TrioCFD code -benchmark session[END_REF][START_REF] Fiorini | Sensitivity equation method for the navier-stokes equations applied to uncertainty propagation[END_REF] for unstructured grids, leading to the following discrete system

DV n+1 h = 0, (6a)
M J n+1 V n+1 h -J n V n h ∆t -J n+1 AV n+1 h -L(V n h )V n+1 h + L(V n )V n+1 h,ALE -GP n+1 h = 0, (6b)
where V h is the discrete fluid velocity, P h is the discrete fluid pressure, V h,ALE the discrete mesh velocity, ∆t the time step, M is the mass matrix, A is the discrete diffusion operator, L(V h ) is the non-linear discrete convection operator, G is the discrete gradient operator, and D is the discrete divergence operator. The superscripts n and n + 1 indicate the time step at which the variable is computed. The nonlinear convective term,

L(V n h )V n+1 h
, is approximate with the upwind MUSCL (Monotonic Upstream-centered Scheme for Conservation Laws, see [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF]) scheme.

A projection-correction technique [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Temam | Une méthode d'approximation de la solution des équations de Navier-Stokes[END_REF] is employed to solve the velocity-pressure coupling. A predicted velocity V * h is computed

J n+1 1 ∆t M -A + L(V n h ) V * h = 1 ∆t M J n V n h + J n+1 L(V n h )V n+1 h,ALE -J n+1 GP n h , (7)
and the mass conservation is then enforced by solving a Poisson equation for pressure

DM -1 G P n+1 h -P n h = 1 ∆t DV * h . ( 8
)
The velocity is corrected using the predicted velocity

V * h V n+1 h = V * h -∆tM -1 G P n+1 h -P n h . (9)
Iterative solvers from the PETSc library (Portable, Extensible Toolkit for Scientific Computation, see [START_REF] Balay | Petsc users manual[END_REF]) solve the discrete equations ( 7) and ( 8).

Presentation of a case study

We now present the results of our two dimensional numerical simulations, considering the case in which the cylinder C 1 is fixed while C 2 is imposed a dimensionless displacement in the e x direction, of the form u = sin (t). To investigate the weakly nonlinear effects generated by the oscillations of C 2 on the fluid forces, the Keulegan-Carpenter number is varied in the range KC ∈ [0.005, 0.3]. Four representative values were chosen for the Stokes number, Sk ∈ {500, 1000, 1500, 2000}, thereby covering a viscous regime for Sk = 500 and an almost inviscid regime for Sk = 2000. Finally, in all our numerical simulations, the dimensionless separation distance is ε = 2. Concerning the numerical setup, the fluid properties are ν = 1.007 • 10 -6 m 2 /s, ρ = 1000 kg/m 3 and the two cylinders have the same radius, R = 1 m. The computational domain is rectangular and its size (L x = 11.5 m and L y = 10.5 m) is considered sufficiently large to minimize the end effects. To discretize the fluid domain, a locally refined grid of triangles generated by the open-source Gmsh mesh generator, see [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in preand post-processing facilities[END_REF], is used. An analysis of mesh sensitivity yields convergent results with a mesh of 379610 triangles, see Table B.5 in Appendix B. This mesh refinement is used in all our numerical simulations. In Fig. 2, we show the time evolution of the x and y-components of the dimensionless fluid force f i , i ∈ {1, 2}, for KC ∈ {0.005, 0.3} and Sk = 500. First, we observe that f i has no component along the y-direction so that the fluid force f i = f i e x is aligned with the direction of motion of C 2 . Second, we show that f i is maximum on C 2 , and is a periodic function whose amplitude, angle with respect to the imposed displacement, and more generally the overall shape, are affected by the value of the Keulegan-Carpenter number. For small values of KC, i.e. small vibration amplitudes of C 2 , the magnitude of the convective term of the Navier-Stokes equations remains negligible. It follows that f i is poorly affected by the nonlinear effects from the fluid and has a sinusoidal shape, as the imposed displacement. For larger values of KC, such as KC = 0.3, the vibrations of C 2 create nonlinear effects in the fluid that strongly affect the shape of f i (especially f 1 ).

Fluid forces

To quantitatively determine the effects of KC on f i , a nonlinear least-squares formulation is used to fit a model to numerical data. As f i is a periodic function, a Fourier model of the form

f i (t) = a 0 2 + N j=1 [a j cos (jωt) + b j sin (jωt)] = a 0 2 + N j=1 h j sin (jωt + ϕ j ) , (10)
is used with {a j , b j , ω} a set of 2(N + 1) fitting parameters. These fitting parameters, namely the real Fourier coefficients a j , b j and the angular frequency ω, are functions of the dimensionless separation distance, ε, the Stokes number, Sk, and the Keulegan-Carpenter number, KC. The magnitude h j and the phase ϕ j of the j-th Fourier harmonics of f i , see second equality of [START_REF] Stokes | On some cases of fluid motion[END_REF], are defined as h j = a 2 j + b 2 j and ϕ j = arg (b j + ia j ). With such a definition of ϕ j , the j-th Fourier harmonics of f i is in phase (resp. in phase opposition) with the imposed displacement if ϕ j = 0 (resp. ϕ j = π). From the fitting parameters {a j , b j }, we also define

T HD = 1 - a 2 1 + b 2 1 N j=1 a 2 j + b 2 j , M = a 0 2 , A = a 2 0 4 + 1 2 N j=1 a 2 j + b 2 j , (11)
as the total harmonic distorsion, the mean value, and the Root Mean Square (RMS) amplitude of f i , respectively. The total harmonic distorsion is an indicator of the importance of the nonlinear effects present in the fluid force signal f i . The root mean square amplitude A indicates the deviation of f i from its mean value a 0 /2. Finally, we note A 0 and ω 0 the limit of A and ω, respectively, as KC → 0. In Fig. 3, we show that the magnitude of all the harmonics of the fluid force f i , i ∈ {1, 2}, decreases as the Stokes number Sk increases. Thus f i is all the more intense as the viscous effects are important, as already observed by [START_REF] Lagrange | New estimations of the added mass and damping of two cylinders vibrating in a viscous fluid, from theoretical and numerical approaches[END_REF] in the case of small oscillations of C 2 , i.e. KC → 0. For moderate oscillations of C 2 , the magnitude of all the harmonics of f 1 also increases as the Keulegan-Carpenter number KC increases. For KC = 0.3, we observe that the magnitude of the second harmonic of f 1 is about 37% of the magnitude of the fundamental harmonic, thus explaining the nonlinear effects shown in the time evolution of f 1 , depicted in Fig. 2. As expected, the strengthening of all the harmonics of f 1 as KC increases eventually yields an increase of its total harmonics distortion. On the other side, the magnitude of the fundamental harmonic (resp. sub-harmonics) of f 2 decreases (resp. increases) as KC increases. This behavior indicates that the kinetic energy of C 2 is transferred to the fluid, and distributed from the fundamental harmonic of the pressure field to the sub-harmonics. For KC ≤ 0.3, the magnitude of the sub-harmonics of f 2 remains very small compared to the magnitude of the fundamental harmonic. Still, this effect sufficiently sensitive to generate a slight increase of the total harmonic distortion of f 2 as KC increases. Finally, the variations of h j shown in Fig. 3 confirm that the moderate oscillations of C 2 create some nonlinear effects in the fluid that tend to intensify (resp. weaken) the hydrodynamic force acting on the stationary (resp. moving) cylinder, as expected from the observation of the time signals in Fig. 2. Having analyzed the variations of the magnitude of the harmonics of f i , we now proceed with the discussion on the variations of the phase ϕ 1 of the fundamental harmonic of f i . In Fig. 4 we show that the two forces have fundamental harmonics that are nearly in phase opposition, as expected from the time evolutions of the fluid forces shown in Fig. 2. At first order, the fundamental harmonic of f 1 (resp. f 2 ) is in phase opposition (resp. in phase) with the imposed sinusoidal displacement of C 2 . The phase of the fundamental harmonic of f i , i ∈ {1, 2}, is sensitive to the Stokes number and is shown to increase as Sk decreases, as already observed by [START_REF] Lagrange | New estimations of the added mass and damping of two cylinders vibrating in a viscous fluid, from theoretical and numerical approaches[END_REF] in the case of small oscillations of C 2 , i.e. KC → 0. For moderate oscillations of C 2 , the phase of the fundamental harmonic of f 1 is poorly sensitive to the Keulegan-Carpenter number whereas the phase of the fundamental harmonic of f 2 exhibits a slight decrease as KC increases. The variations of ϕ j with the Stokes and the Keulegan-Carpenter numbers imply that the direction of the fluid force f i depends on both Sk and KC. For one oscillation of the moving cylinder C 2 , the fluid force f i vanishes and reverses its direction when f i (t) = 0, i.e. when t = t

Effect of the Keulegan-Carpenter number

(1) i and t = t

(2) i . At these two specific times, the dimensionless displacement of the moving cylinder is u

(1) i = sin t (1) i and u (2) i = sin t (2) i > u (1)
i . In Fig. 5, we show that f 1 is positive (resp. negative) when

-1 ≤ u ≤ u (1) 1 (resp. u (2) 1 ≤ u ≤ 1). As f 2 is at first order in phase opposition with f 1 , it is negative (resp. positive) when -1 ≤ u ≤ u (1) 2 (resp. u (2) 2
≤ u ≤ 1). It follows that the fluid forces cause the cylinders to attract (resp. repel) each other when -1 ≤ u ≤ min u

(1) 1 , u (1) 2 (resp. max u (2) 1 , u (2) 2 ≤ u ≤ 1).
In the narrow range u

(1) i ≤ u ≤ u (2) i
the sign of f i depends on the direction of motion of the cylinder C 2 . When C 2 approaches C 1 , f 1 (resp. f 2 ) is negative (resp. positive) so that the fluid forces cause the cylinders to repel each other when max u

(1) 1 , u (1) 2 ≤ u ≤ min u (2) 1 , u (2) 2
. On the contrary, when C 2 moves away from C 1 , f 1 (resp. f 2 ) is positive (resp. negative). In that case, the fluid forces cause the cylinders to attract each other when max u

(1) 1 , u (1) 2 ≤ u ≤ min u (2) 1 , u (2) 2 .
To assess the validity of our numerical predictions, we refer to the work of [START_REF] Stokes | On some cases of fluid motion[END_REF] who considered the small oscillations, i.e. KC → 0, of an isolated cylinder, i.e. ε → ∞, immersed in a viscous fluid initially at rest. In such a case the fluid problem is fully linear and can be solved by introducing a stream-function from which the fluid force f 2 (t) = ℜ{e it f2 } is derived, ℜ being the real part operator. For a harmonic displacement, i.e. u (t) = ℜ{e it ū} , the complex representation of the force writes f2 = ūπ

  1 + 4 √ iSk K1 √ iSk K0 √ iSk   , (12)
with K 0 and K 1 the modified Bessel functions of second kind. The change of sign of the fluid force function occurs when t = -arg f2 /ū + kπ, k ∈ Z, leading to

(j) 2 = (-1) j | sin arg f2 /ū |. ( 13
)
In Fig. 5 , we show that the graph of Eq. ( 13) is in line with the one obtained numerically for KC = 0.005, thereby corroborating our predictions. Having analyzed the properties of the Fourier harmonics of the fluid forces, we now proceed with the determination of the variations of the angular frequency ω, mean value M and RMS amplitude A of f i , as KC and Sk are varied. The quantities ω 0 , M 0 = 0 and A 0 , corresponding to KC → 0, are used as references. In Fig. 6, we show the variations of ω/ω 0 , M and A/A 0 for KC ≤ 0.3. For all Stokes numbers Sk ∈ {500, 1000, 1500, 2000}, we observe that ω/ω 0 is a constant function of KC. It follows that the moderate oscillations of C 2 do not significantly affect the period of f i , so that ω = ω 0 = 2π. These moderate vibration amplitudes mostly affect the dimensionless mean value M and the RMS amplitude ratio A/A 0 of f i . The mean value of f 1 (resp. f 2 ) is shown to linearly decrease (resp. increase) as KC increases. On the other hand, the RMS amplitude ratio of f 1 (resp. f 2 ) increases as KC 2 (resp. decrease as -KC 3 ) as KC increases.

Self-added coefficients and power laws

In the previous section, we have shown that the amplitude of the sub-harmonics of f 2 remain small compared with the amplitude of the fundamental harmonic, and that even in the case of moderate vibration amplitudes of C 2 . Thus, in first approximation, f 2 is a linear combination of cos (t) and sin (t). Introducing m self = b 1 /π and c self = -a 1 /π as the self-added mass and damping coefficients, see [START_REF] Lagrange | New estimations of the added mass and damping of two cylinders vibrating in a viscous fluid, from theoretical and numerical approaches[END_REF], f 2 writes

f 2 (t) ≈ π [m self sin (t) -c self cos (t)] . (14)
The self-added coefficients relate the fluid force acting on the cylinder C 2 due to its own motion. As a 1 and b 1 , these coefficients depend on both the Keulegan-Carpenter number, the Stokes number and the dimensionless separation distance. As shown in § 3.2, the fluid force on C 1 is affected by the nonlinear effects and so does not express as a linear combination between the velocity of C 2 and its Therefore, the linear concept of cross-added mass and damping coefficients cannot be extended to the case of moderate oscillations of C 2 . In Fig. 7, we show that m self → m 0 as KC → 0 and decreases monotonically as KC increases in the range [0, 0.3]. The self-added damping coefficient first slightly decreases from its asymptotic limit c self → c 0 as KC → 0, and eventually increases as KC increases. From these variations and in first approximation, it is assumed that the fluid added-coefficients follow power laws of the form m self = m 0 (1 + m 1 Sk qm KC pm ) and c self = c 0 (1 + c 1 Sk qc KC pc ), with m 0 , c 0 some functions of (Sk, ε) and m 1 , c 1 some functions of ε only. The self-similarity of the lines in Fig. 7, obtained for different Stokes numbers Sk ∈ {500, 1000, 1500, 2000}, also suggests that p m and p c do not depend on Sk. In our analysis, we assume that q m , p m , q c and p c actually do not depend on ε, and therefore are constant power coefficients. To precisely determine these constants, we first study the variations of m 0 and c 0 with the Stokes number, by carrying out numerical simulations for small vibrations amplitudes of C 2 , i.e. KC = 0.005. The evolutions of m 0 and c 0 are depicted in Fig. 8, for Sk ∈ [100, 900]. In agreement with the theoretical predictions of [START_REF] Lagrange | New estimations of the added mass and damping of two cylinders vibrating in a viscous fluid, from theoretical and numerical approaches[END_REF], we show that m 0 and c 0 decreases as Sk -1/2 . Noting m 00 the limit of m self as Sk → ∞ (inviscid limit, i.e. c 00 = 0), we assume that m 0 and c 0 follow power laws of the form

m 0 ≈ m 00 (ε) 1 + m 01 (ε) Sk 1/2 , c 0 ≈ c 01 (ε) Sk 1/2 , (15)
with m 00 , m 01 and c 01 , some functions of ε. Having obtained power laws for the small oscillation terms, m 0 and c 0 , a nonlinear least-squares formulation is used to fit the numerical data shown in Fig. 7 (and given in Appendix A) to the models for moderate oscillations, i.e. m self = m 0 (1 + m 1 Sk qm KC pm ) and c self = c 0 (1 + c 1 Sk qc KC pc ). This formulation yields power constants q m = -1/4, p m = 3, q c = 1/8 and p c = 3 so that the added-mass and damping coefficients eventually write

m self ≈ m 00 (ε) 1 + m 01 (ε) Sk 1/2 1 + m 1 (ε) Sk 1/4 KC 3 , c self ≈ c 01 (ε) Sk 1/2 1 + c 1 (ε) Sk 1/8 KC 3 , (16)
with m 1 and c 1 some functions of ε. For ε = 2, the nonlinear least-squares formulation yields m 00 = 1.009, m 01 = 2.8, m 1 = -6.5, c 01 = 3.14 and c 1 = 4.81. In Fig. 9, we show that the predictions from the power laws [START_REF] Sigrist | Fluid-structure interaction: an introduction to finite element coupling[END_REF] are in very good agreement with the numerical results, and that for all the Stokes and Keulegan-Carpenter numbers simulated. Indeed, the relative deviation δ between the numerical results and the predictions from the power laws is less than 0.25% for the self-added mass coefficient and less than 7% for the self-added damping coefficient, see Tables A.1-A.4 in Appendix A. The power laws [START_REF] Sigrist | Fluid-structure interaction: an introduction to finite element coupling[END_REF] have been obtained for Sk ≥ 500 and in the limit of moderate Keulegan-Carpenter numbers, KC ≤ 0.3. These new laws extend the concept of the self-added coefficients, initially introduced for linear problems involving small oscillations of an immersed structure, to the case of weakly nonlinear problems involving moderate oscillations.

Conclusion

In the present work, we have considered the interaction of two parallel circular cylinders immersed in a fluid at rest. Two dimensional numerical simulations have been performed to determine the fluid forces acting on the when one of them is imposed a harmonic motion with moderate vibration amplitudes.The direction of the imposed motion is parallel to the line joining the cylinder centers. First, we showed that the fluid forces on the two cylinders have only one component, along the direction of motion. Second, we showed that moderate vibration amplitudes enhanced nonlinear effects in the fluid, leading to increasing sub-harmonic Fourier components in the time signal of the fluid force acting on the stationary cylinder. The force on the moving cylinder has a strong fundamental harmonic, negligible sub-harmonics, so that it is poorly affected by the fluid nonlinear effects. From this observation, we have extended the concept of the self-added coefficients, initially introduced for linear problems involving small oscillations of the cylinder, to the case of weakly nonlinear problems involving moderate oscillations. Nonlinear regressions of our numerical results yielded new power laws for the self added-coefficients, expressed in terms of the Stokes and the Keulegan-Carpenter numbers. We showed that the two coefficients scale as Sk -1/2 , whereas the self-added mass (resp. damping) decreases (resp. increases) as -KC 3 (resp. +KC 3 ).

In future work, we will perform numerical simulations considering various values of the dimensionless separation distance, in order to derive new laws for the functions of ε appearing in [START_REF] Sigrist | Fluid-structure interaction: an introduction to finite element coupling[END_REF]. We will then investigate the highly viscous regimes, i.e. Sk ≪ 500, which require adequate meshes to accurately reproduce the fluid flow in the boundary layer of the moving cylinder. Finally, we shall investigate the case of large oscillations, i.e. KC ≫ 0.3, to reveal the limitations of the linear concept associated to the self-added coefficients.

KC

Appendix B. Mesh convergence

In this appendix, we present the results of a mesh sensitivity analysis conducted for Sk = 900, KC = 0.005 and a dimensionless separation distance ε = 2. The Table B.5 shows that the two dimensional mesh with 379610 triangles is sufficiently refined to accurately estimate the self-added coefficients. This mesh uses two different local sizes for the elements: a small local size of 0.02 m for the elements close to the cylinders, and a large global size of 0.08 m for the other elements. m 0 c 0 Theoretical predictions of [START_REF] Lagrange | New estimations of the added mass and damping of two cylinders vibrating in a viscous fluid, from theoretical and numerical approaches[END_REF] 1

Figure 1 -

 1 Figure 1 -Schematic diagram of two parallel circular cylinders with radius R, separated by a distance E, immersed in a viscous stagnant fluid domain Ω f luid . The cylinder C 1 is fixed while C 2 is imposed a displacement vector U.

Figure 2 -

 2 Figure 2 -Time evolution of the dimensionless fluid force, for two different values of the Keulegan-Carpenter number, KC ∈ {0.005, 0.3}. Top figures : x-component, f i , of the fluid force. Bottom figures : y-component of the fluid force. On these figures the red (KC = 0.005) and blue (KC = 0.3) lines are superimposed. Left figure : fluid force acting on the fixed cylinder. Right figure : fluid force acting on the moving cylinder. The Stokes number is Sk = 500.

Figure 3 -

 3 Figure 3 -Evolution of the magnitude, h 1 of the fundamental harmonic of the fluid force, the normalized magnitude, h j /h 1 of the j-th harmonic, and the total harmonic distortion, T HD, as a function of the Keulegan-Carpenter number, KC, for different values of the Stokes number, Sk ∈ {500, 1000, 1500, 2000}. Left figure : fluid force f 1 acting on the fixed cylinder. Right figure : fluid force f 2 acting on the moving cylinder.

Figure 4 -

 4 Figure 4 -Phase, ϕ 1 of the fundamental harmonic of the fluid force, as a function of the Keulegan-Carpenter number, KC, for different values of the Stokes number, Sk ∈ {500, 1000, 1500, 2000}. Left figure : fluid force f 1 acting on the fixed cylinder. Right figure : fluid force f 2 acting on the moving cylinder.

Figure 5 -

 5 Figure 5 -Direction (+e x or -e x ) of the fluid forces, depending on Sk ∈ {500, 1000, 1500, 2000}, KC ∈ {0.005, 0.1, 0.2, 0.3} and the dimensionless displacement u of C 2 . Left figure : fluid force f 1 acting on the fixed cylinder. Right figure : fluid force f 2 acting on the moving cylinder. On this figure the arrows show the direction of motion of C 2 . The black line shows the graph of Eq. (13), derived in the case of an isolated cylinder with KC → 0. In both figures, the open symbols represent numerical results. The solid (resp. dashed) lines correspond to u = u(1) i (resp. u = u (2) i ). .

Figure 6 -

 6 Figure 6 -Evolution of the angular frequency ratio, ω/ω 0 , dimensionless mean value, M , and RMS amplitude ratio A/A 0 of the dimensionless fluid force f i , as a function of the Keulegan-Carpenter number, KC, for different values of the Stokes number, Sk ∈ {500, 1000, 1500, 2000}. Left figures : fluid force f 1 acting on the fixed cylinder. figures : fluid force f 2 acting on the moving cylinder.

Figure 7 -

 7 Figure 7 -Evolution of the dimensionless self-added coefficients, m self and c self , as a function of the Keulegan-Carpenter number, KC.

Figure 8 -

 8 Figure 8 -Evolution of the dimensionless self-added coefficients, m 0 and c 0 , as a function of the Stokes number, Sk. The solid lines are theoretical predictions from [21] and the symbols are numerical results. The horizontal dashed lines correspond to the asymptotic limits Sk → +∞. The Keulegan-Carpenter number is KC = 0.005.

Figure 9 -

 9 Figure 9 -Comparison between the dimensionless self-added coefficients, m self and c self , obtained numerically, and those predicted by the fitted power laws (16). The shaded gray areas show the relative deviations δ ≤ 1% and δ ≤ 5%. The dashed lines correspond to δ = 0%.

Table A .

 A 2 -Comparison between the numerical results and the fitted power laws[START_REF] Sigrist | Fluid-structure interaction: an introduction to finite element coupling[END_REF] for the dimensionless self-added coefficients. The maximum relative deviation is highlighted in red color. The Stokes number is Sk = 1000.

		Numerics TrioCFD	Fitted power laws	Relative deviation (%)
		m self	c self	m self	c self	m self	c self
	0.005	1.1529	0.1383	1.1529	0.1404	0	1.5459
	0.025	1.1526	0.1382	1.1529	0.1404	0.0260	1.6198
	0.05	1.1520	0.1386	1.1527	0.1406	0.0642	1.4084
	0.1	1.1499	0.1412	1.1514	0.1419	0.1279	0.5003
	0.15	1.1463	0.1452	1.1476	0.1454	0.1108	0.1218
	0.2	1.1403	0.1522	1.1403	0.1522	0.0028	0.0367
	0.25	1.1300	0.1631	1.1282	0.1634	0.1597	0.1717
	0.3	1.1115	0.1802	1.1102	0.1801	0.1178	0.0356
	Table A.1 -Comparison between the numerical results and the fitted power laws (16) for the dimensionless
	self-added coefficients. The maximum relative deviation is highlighted in red color. The Stokes number is
	Sk = 500.
	KC	Numerics TrioCFD	Fitted power laws	Relative deviation (%)
		m self	c self	m self	c self	m self	c self
	0.005	1.1143	0.0997	1.1143	0.0993	0	0.3987
	0.025	1.1139	0.0992	1.1142	0.0993	0.0267	0.1498
	0.05	1.1134	0.0991	1.1141	0.0994	0.0612	0.3325
	0.1	1.1117	0.1001	1.1130	0.1004	0.1114	0.3693
	0.15	1.1088	0.1030	1.1099	0.1031	0.1015	0.1049
	0.2	1.1036	0.1081	1.1040	0.10842	0.0296	0.2123
	0.25	1.0945	0.1160	1.0941	0.1170	0.0331	0.8416
	0.3	1.0769	0.1304	1.0795	0.1299	0.2274	0.3564
	KC	Numerics TrioCFD	Fitted power laws	Relative deviation (%)
		m self	c self	m self	c self	m self	c self
	0.005	1.0970	0.0839	1.0970	0.0811	0	3.3651
	0.025	1.0968	0.0828	1.0970	0.0811	0.0191	2.0319
	0.05	1.0964	0.0822	1.0969	0.0812	0.0436	1.1914
	0.1	1.0951	0.0823	1.0959	0.0820	0.0729	0.3622
	0.15	1.0926	0.0843	1.0931	0.0844	0.0509	0.0211
	0.2	1.0881	0.0883	1.0878	0.0889	0.0212	0.6583
	0.25	1.0800	0.0950	1.0791	0.0963	0.0828	1.4013
	0.3	1.0641	0.1096	1.0661	0.1074	0.1796	2.0887

Table A .

 A 3 -Comparison between the numerical results and the fitted power laws[START_REF] Sigrist | Fluid-structure interaction: an introduction to finite element coupling[END_REF] for the dimensionless self-added coefficients. The maximum relative deviation is highlighted in red color. The Stokes number is Sk = 1500.

	KC	Numerics TrioCFD	Fitted power laws	Relative deviation (%)
		m self	c self	m self	c self	m self	c self
	0.005	1.0867	0.0752	1.0867	0.0702	0	6.6559
	0.025	1.0866	0.0735	1.0866	0.0702	0.0064	4.4836
	0.05	1.0863	0.0723	1.0865	0.0703	0.0173	2.7282
	0.1	1.0853	0.0721	1.0856	0.0711	0.0243	1.3541
	0.15	1.0833	0.0735	1.0831	0.0732	0.0147	0.4382
	0.2	1.0792	0.0769	1.0782	0.0772	0.0940	0.4347
	0.25	1.0720	0.0822	1.0702	0.0839	0.1738	2.0658
	0.3	1.0580	0.0938	1.0581	0.0938	0.0139	0.0288

Table A .

 A 4 -Comparison between the numerical results and the fitted power laws[START_REF] Sigrist | Fluid-structure interaction: an introduction to finite element coupling[END_REF] for the dimensionless self-added coefficients. The maximum relative deviation is highlighted in red color. The Stokes number is Sk = 2000.

Table B .

 B 5 -Mesh convergence. The Stokes number is Sk = 900 and the Keulegan-Carpenter number is KC = 0.005.

		.11	0.105
	Dimension of the mesh	Numerics Numerics
	29620	1.09	0.134
	110116	1.11	0.112
	242916	1.12	0.108
	379610	1.12	0.106
	678 088	1.12	0.106

Appendix A. Tables of results

In this appendix, we provide the numerical and fitted values, see Eq. ( 16), of the self-added coefficients, depicted in