Distributed Source Seeking without Global Position Information

Ruggero Fabbiano, Federica Garin, Carlos Canudas de Wit

Abstract-We present a distributed control law to steer a group of autonomous communicating sensors towards the source of a diffusion process. The graph describing the communication links between sensors has a time-invariant topology, and each sensor is able to measure (in addition to the quantity of interest) only the relative bearing angle with respect to its neighbour, but has no absolute position information and does not know any relative distance. Using multiple sensors is useful in wide environments (e.g., under the sea), or when the function describing the diffusion process is slowly changing in space, so that a single sensor may have to travel long distances before having a good gradient estimation. Our approach is based on a twofold control law, which is able to bring and keep the set of sensors on a circular equispaced formation, and to steer the circular formation towards the source via a gradient-ascent technique. The effectiveness of the proposed algorithm is both theoretically proven and supported by simulation results.

Index Terms-Source seeking, Distributed control, Formation control, Gradient search, Poisson integral, Sensor networks.

I. INTRODUCTION

T HE source seeking problem consists in localizing a source emitting from an unknown position, by exploiting the fact that the concentration of the emitted quantity is higher near the source and decays with distance. This problem arises in various applications, e.g., detection of pollutants, explosives, drugs, or localization of sound, heat or vent sources.

In recent years, leveraging the richness of results on collaborative sensing, distributed control and multi-agents formations, a new line of research has emerged, focusing on multi-agent source seeking, where multiple mobile sensors collaborate in order to find the source, improving over single-sensor techniques, e.g., by allowing faster convergence, exploration of a wider area, robustness to errors and failures.

One of the first contributions is [START_REF] Moreau | Coordinated gradient descent: a case study of Lagrangian dynamics with projected gradient information[END_REF], where it is assumed that each vehicle can measure the full gradient, and the authors develop a twofold algorithm with a gradient-descent term and inter-vehicle forcing terms. Other early contributions are [START_REF] Ögren | Cooperative control of mobile sensor networks: adaptive gradient climbing in a distributed environment[END_REF] and [START_REF] Rabbat | Distributed optimization in sensor networks[END_REF]. In the former, formation control uses potential functions to drive agents in a desired formation moving along a virtual body, whose gradient-climbing trajectory depends on a computation of the gradient obtained filtering all sensor measurements. In [START_REF] Rabbat | Distributed optimization in sensor networks[END_REF], the authors propose to solve the sourcelocalization problem with distributed optimization algorithms making use of measurements from sensors spread all over the R. Fabbiano was and F. Garin is with INRIA (NeCS team, Inria Grenoble Rhône-Alpes) and GIPSA-lab, Grenoble, France; e-mails: ruggero.fabbiano@sfr.fr, federica.garin@inria.fr.

C. Canudas-de-Wit is with CNRS, GIPSA-lab (NeCS team), Grenoble, France; e-mail: carlos.canudas-de-wit@gipsa-lab.fr. region of interest, and not moving towards the source. Papers [START_REF] Moore | Source-seeking via collaborative measurements by a circular formation of agents[END_REF], [START_REF] Briñón Arranz | Collaborative estimation of gradient direction by a formation of AUVs under communication constraints[END_REF] have introduced a collaborative control law to steer a fleet of AUVs (autonomous underwater vehicles) to the source of a signal distribution using only direct signal measurements, by arranging the AUVs in a circular formation and using measurements along the circle to get a gradient approximation. The formation control needed to achieve this goal has been further studied in [START_REF]Cooperative control design for time-varying formations of multiagent systems[END_REF]. The gradient approximation formula has been studied in [START_REF] Fabbiano | Source localization by gradient estimation based on Poisson integral[END_REF], giving a mathematical framework for its validity and an extension to the 3-dimensional case; in this paper, though, gradient ascent was performed in a centralized way, with a single agent using measurements from all sensors. The recent paper [START_REF] Briñón Arranz | Distributed sourceseeking via a circular formation of agents under communication constraints[END_REF] builds upon ideas from [START_REF] Briñón Arranz | Collaborative estimation of gradient direction by a formation of AUVs under communication constraints[END_REF] to perform an improved collaborative search, with a Newton-Raphson method.

A gradient-based method is also considered in [START_REF] Atanasov | Distributed algorithms for stochastic source seeking with mobile robot networks[END_REF], where the authors propose two algorithms, based on the estimation of the gradient at the formation centroid, with different techniques for the two cases where the function describing the spatial distribution of the measured quantity is known or unknown; then, stochastic approximation is used to ensure convergence.

In [START_REF] Wu | Robust cooperative exploration with a switching strategy[END_REF], the gradient estimation makes use also of time series of measurements, instead of instantaneous measurements only, so that agents can also have an individual local gradient direction, and alternate phases of individual exploration with a collaborative exploration in formation.

Another class of algorithms are extremum seeking ones, where the gradient direction is indirectly estimated from the measurements, by probing the environment with suitable movements. In the early work [START_REF] Biyik | Gradient climbing in formation via extremum seeking and passivity-based coordination rules[END_REF], a leader agent performs extremum seeking, while all other agents are controlled to behave as followers. In [START_REF] Khong | Multi-agent source seeking via discrete-time extremum seeking control[END_REF], distributed formation control is combined with a centralized extremum seeking controller, which gives a gradient-ascent reference to the formation.

Most of these algorithms make explicit use of the absolute position of the sensors. A notable exception is [START_REF] Atanasov | Distributed algorithms for stochastic source seeking with mobile robot networks[END_REF], where the formation control and gradient estimation rather rely on relative localization, with respect to other sensors and to the formation centre.

Paper Contribution

We propose a distributed algorithm to estimate the gradient of a function describing a diffusion process, and use this information to drive a set of moving sensors towards the source. The effectiveness of such an algorithm is both theoretically proven and supported by simulation results.

The framework is close to the one of previous works [START_REF] Moore | Source-seeking via collaborative measurements by a circular formation of agents[END_REF], [START_REF] Briñón Arranz | Collaborative estimation of gradient direction by a formation of AUVs under communication constraints[END_REF], [START_REF] Briñón Arranz | Distributed sourceseeking via a circular formation of agents under communication constraints[END_REF], but our contribution differs as we avoid using full position information, which may be unavailable in various operating environments such as underwater or indoor vehicles, or in applications where inertial navigation systems are too expensive or not sufficiently accurate. We only suppose the sensors to be able to measure the relative angle with respect to their neighbours, differently from all above-mentioned distributed source seeking methods, which require either the knowledge of sensors absolute position, or at least relative localization. The main motivation for the restriction to relative angles is that in vision-based navigation angles are much easier to obtain than distances, see e.g. [START_REF] Scaramuzza | Precise bearing angle measurement based on omnidirectional conic sensor and defocusing[END_REF] for an example of bearing-angle sensor. We have already considered the same assumption in our preliminary work [START_REF] Fabbiano | Distributed source localisation with no position information[END_REF]; in that paper, though, our twofold objective (formation control and motion towards the source) was expressed in an antagonistic control law, which, together with a formation control algorithm without theoretical stability guarantees [START_REF] Moshtagh | Bearingonly control laws for balanced circular formations of ground robots[END_REF], resulted in an algorithm with slower convergence in simulation and no theoretical convergence proof. In this paper we introduce and study a new algorithm, where the formation control is based on techniques from [START_REF] Marshall | Formations of vehicles in cyclic pursuit[END_REF], together with a suitable change of variables and a distributed gradient-approximation algorithm, so that we can express the global objective of driving the vehicles to a circular formation whose centre moves along the gradient direction, by making use only of the limited local information.

Notation

We will use the notation z = ρe iθ to denote the point (x, y) = ρ(cos θ, sin θ) ∈ R 2 , using the usual bijection z = (x, y) between points in R 2 and elements of C. Multiplication e iα z will denote rotation of vector z by an angle α (counterclockwise for α > 0), as in usual multiplication of complex numbers. However, all other operations will be intended as operations on reals and on vectors in R 2 ; in particular, all derivatives are intended as (partial) derivatives of real-valued functions, and integrals are intended as entry-wise integrals of the real-valued entries of the vectors. ∇f is the gradient of the function f , and ∇ 2 is the Laplacian operator. ∂Ω denotes the border of an open set Ω, and its closure is denoted by Ω = Ω ∪ ∂Ω. B r (z c ) indicates the circle of radius r centered in z c . Sums of sensor indexes 1, . . . , n are always intended mod n, and sums of angles are mod 2π.

II. PROBLEM FORMULATION A. Autonomous Sensors

We consider a group of n communicating sensors moving in a region Ω ⊂ R 2 , where a diffusion process is taking place, with the nonholonomic kinematics żi (t) = v i (t)e iθi(t)

v i (t) θi (t) = u i (t), (1)
where z(t) is the rotational point of the sensor in the plane, and θ(t) is its heading angle. The sensors, which can take pointwise measurements of the emitted quantity, have no position information, but know the value of n and are able to measure the bearing angle with respect to their following neighbour (i.e., the angle between the orientation of an agent and the vector from itself to its neighbour, as explained in Fig. 1), defined by

i i + 1 θ i α i
α i (t) = arg (z i+1 (t) -z i (t)) -θ i (t). (2)
We consider a ring communication network, where sensor i can receive messages sent by sensor i -1, irrespective of the sensors position and of time (see Fig. 2). This assumption is not essential (other communication strategies are possible), but is natural when the sensors are not too far from a circular formation, as it will be the case in this paper. Indeed, we will show that a circular formation with sensors equispaced along the circle is desirable in order to compute a gradient approximation from their measurements, and we will describe a suitable control law which ensures that such a formation is obtained (from initial conditions not too far from this equilibrium) and maintained while the formation centroid is driven towards the source.

B. Diffusion Process

We consider steady-state behaviors of homogeneous diffusion processes caused by an isotropic source emitting at a constant rate. Such a process is governed by the well-known diffusion equation

∂f (z, t) ∂t -c∇ 2 f (z, t) = 0, ∀z ∈ Ω, t ≥ 0, (3)
where f is the concentration variable, c is a diffusion coefficient, and Ω ⊂ R 2 . In particular, as depicted in Fig. 3, we consider the region of interest Ω ⊂ R 2 as a connected bounded set Ω = Ω ∪ Ω s , where Ω s is the connected bounded set that identifies the source, and therefore we have that ∂Ω = ∂ Ω ∪ ∂Ω s . The values of f on the inner boundary ∂Ω s are imposed by the source, so we can assume that values of f on ∂Ω s are higher than the ones on ∂ Ω. As we will see below, we have that max f (z) ∈ ∂Ω, that means, for our previous consideration, that it lies on ∂Ω s ; our source localisation problem is then mathematically equivalent to the problem of finding the maximum of f . Once the steady-state has been reached, supposing that the source is still emitting at a constant rate (this happens in many cases of practical interest, e.g., in a heating process, or in the dispersion of a chemical substance), or that possible source variations are slow in the time-scale of interest, the diffusion equation ( 3) can be approximated by the Laplace equation, i.e.,

∇ 2 f (z) = 0, ∀z ∈ Ω, (4)
whose solutions are called "harmonic functions" and have many mathematical properties which can be exploited to devise a gradient-ascent source seeking.

III. SOURCE-SEEKING STRATEGY

A. Gradient ascent

The source-localisation problem, formulated in previous section as the problem of finding the maximum of a harmonic function f , can be solved by a gradient ascent algorithm. We refer the reader to our previos work [START_REF] Fabbiano | Source localization by gradient estimation based on Poisson integral[END_REF] for a detailed discussion, while we briefly recall here the main results, which will be used in the remainder of this paper as a tool for constructing a distributed source-seeking algorithm.

Harmonic functions satisfy the so-called local maximum principle, which ensures that (except for the trivial case of a constant function), f can have maxima only on the border ∂Ω, not in the interior of Ω. This justifies the use of a gradient ascent technique, despite the fact that harmonic functions might not be convex.

Then, harmonicity of f ensures the possibility to determine f and all derivatives of f at a point z c , using measurements along a curve encircling z c . In particular, in the case where the curve is a circle, explicit integral formulae exist, called Poisson integrals for the circle. We will make use of a Poisson integral to find the gradient f (z c ), based on the values of f along a circle centered in z c , as follows:

∇f (z c ) = 1 πr 2π 0 e iθ f z c + re iθ dθ. (5)
In the following, we will approximate the above integral with a finite sum of measurements along the circle, taken by sensors, while a suitably designed control law will steer the sensors in a circular formation whose centre moves in the direction of the computed gradient.

As discussed in [START_REF] Briñón Arranz | Collaborative estimation of gradient direction by a formation of AUVs under communication constraints[END_REF], [START_REF] Fabbiano | Source localization by gradient estimation based on Poisson integral[END_REF], the Poisson formula [START_REF] Briñón Arranz | Collaborative estimation of gradient direction by a formation of AUVs under communication constraints[END_REF], which holds exactly for harmonic functions, can also be used to approximate the gradient of more general functions, under the only requirement of C 1 regularity. Hence, the distributed source-seeking method developed in this paper can be applied to some more general sources, provided they have enough regularity to justify the gradient approximation and the gradient ascent technique.

B. Control Law

We need to develop an algorithm which is able to both steer the sensors towards the gradient direction, and bring them on a circular formation, such that the gradient itself is correctly estimated.

We introduce a control law assuming that each sensor i is able to compute the difference θi (t) = θ ref (t)θ i (t), where θ ref is a common reference (we will discuss in Section IV how to obtain it). Let v c > 0 and v r > v c , with v = vc vr , be fixed parameters; we propose the following control law:

v i (t) = v c cos θi (t) + v 2 r -v 2 c sin 2 θi (t) (6A) θi (t) = b i (t) θref (t) + 1 -b i (t) ω i (t), (6B)
where

b i (t) = v c v i (t) cos θi (t) (6C) ω i (t) = k α i (t) + arcsin v sin θi (t) , (6D)
k > 0 is a positive constant gain, and α i is defined as in [START_REF] Ögren | Cooperative control of mobile sensor networks: adaptive gradient climbing in a distributed environment[END_REF].

We note that such a control law only involves local variables, as it is function of only the measured relative bearing angle and the angle θi , which is an angle difference and therefore independent of the coordinate system. We will show that this control law brings the autonomous sensors to settle on an equispaced circular formation around a centre moving along the reference direction. As a first step, the following proposition shows an equivalent description of system (1), [START_REF]Cooperative control design for time-varying formations of multiagent systems[END_REF].

Proposition 1: The closed-loop system given by dynamics (1), under the control law (6A)-(6C), is equivalent to the system

żi (t) = v c e iθ ref (t) + v r e iγi(t) (7A) γi (t) = ω i (t).
(7B)

Proof: From (7A) we compute

| żi | 2 = v 2 c + 2v c v r (cos θ ref cos γ i + sin θ ref sin γ i ) + v 2 r , ( 8A
)
arg żi = arctan v c sin θ ref + v r sin γ i v c cos θ ref + v r cos γ i , (8B)
and moreover,

v r cos γ i = | żi | cos(arg żi ) -v c cos θ ref v r sin γ i = | żi | sin(arg żi ) -v c sin θ ref . (9)
Using ( 9) we get

v r (cos θ ref cos γ i + sin θ ref sin γ i ) = | żi | cos(θ ref -arg żi ) -v c (10
) and, by substituting it back into (8A), it yields

| żi | 2 = 2v c cos (θ ref -arg żi ) | żi | -v 2 c + v 2 r . (11)
Taking the positive solution of this quadratic equation, we obtain the expression for the velocity:

| żi | = v c cos (θ ref -arg żi ) + cos 2 (θ ref -arg żi ) -1 v 2 c + v 2 r = v c cos (θ ref -arg żi ) + v 2 r -v 2 c sin 2 (θ ref -arg żi ). ( 12
)
Computing the time derivative of (8B) we get

d dt (arg żi ) = v c θref [v c + v r (cos θ ref cos γ i + sin θ ref sin γ i )] v 2 c + 2v c v r (cos θ ref cos γ i + sin θ ref sin γ i ) + v 2 r + v r γi [v r + v c (cos θ ref cos γ i + sin θ ref sin γ i )] v 2 c + 2v c v r (cos θ ref cos γ i + sin θ ref sin γ i ) + v 2 r ;
substituting again [START_REF] Wu | Robust cooperative exploration with a switching strategy[END_REF]:

d dt (arg żi ) = v c | żi | cos (θ ref -arg żi ) θref 2v c | żi | cos (θ ref -arg żi ) -v 2 c + v 2 r + v 2 r + v c | żi | cos (θ ref -arg żi ) -v 2 c γi 2v c | żi | cos (θ ref -arg żi ) -v 2 c + v 2 r
and, finally, by [START_REF] Biyik | Gradient climbing in formation via extremum seeking and passivity-based coordination rules[END_REF],

d dt (arg żi ) = v c cos (θ ref -arg żi ) | żi | θref + 1 - v c cos (θ ref -arg żi ) | żi | γi . ( 13
)
We can notice that, with notation v i = | żi | and θ i = arg żi as in [START_REF] Moreau | Coordinated gradient descent: a case study of Lagrangian dynamics with projected gradient information[END_REF], the final expressions of equations ( 12) and ( 13) correspond to our choice of control input u i = [v i , θi ] T in (6A) and (6B), thus showing the equivalence of systems (1), (6A)-(6C) and [START_REF] Fabbiano | Source localization by gradient estimation based on Poisson integral[END_REF].

Remark 1: We can see from this equivalence that the dynamics of each sensor is given by the superposition of a motion with speed v c in the common direction θ ref and a motion with speed v r and heading γ i such that γi = ω i . In particular, since the direction θ ref is common to each autonomous sensor, the motion with speed v c results in a translation of the whole formation in such direction. We need to show that our law ensures that the sensors settle on a circular equispaced formation (thus allowing for a gradient estimation based on the Poisson integral formula), i.e., they asymptotically tend to a formation in which each sensor has position and velocity

z i (t) = z c (t) + re iζi(t) żi (t) = v c e iθ ref (t) + v r e iγi(t) , with ζ i (t) = ζ 0 (t) + 2π n i and γ i (t) = ζ i (t) + π 2 .

C. Convergence to formation

Here, we show how the control law ( 6) can bring the sensors to the desired circular formation. We start by recalling some results by Marshall et al. [START_REF] Marshall | Formations of vehicles in cyclic pursuit[END_REF].

Let us consider a fleet of autonomous underactuated vehicles described by the dynamics [START_REF] Moreau | Coordinated gradient descent: a case study of Lagrangian dynamics with projected gradient information[END_REF]. By a suitable change of coordinates [START_REF] Marshall | Formations of vehicles in cyclic pursuit[END_REF], the system (1) can be described by the new state variables

d i (t) = |z i+1 (t) -z i (t)| (14A) α i (t) = arg(z i+1 (t) -z i (t)) -θ i (t) (14B) β i (t) = θ i (t) -θ i+1 (t) -π. (14C)
With such state variables we lose any information about the absolute position of each vehicle; nonetheless, they are more suitable to describe the formation reached, as they define the position of each vehicle in terms of relative variables with respect to each other's neighbours, such as the relative distance d i and the relative bearing angle α i . For the subsequent analysis, it is important to notice that the above-defined variables satisfy the following constraints1 :

n i=1 β i = -nπ (15)
n i=1 d i e i[αi-(i-1)π-i-1 j=1 βj ] = 0. (16)
Eq. ( 15) is immediate from the definition of β i , while Eq. ( 16) is obtained from n i=1 (z i+1 -z i ) = 0, where the vector z i+1z i can be expressed with its modulus d i and its argument arg(z i+1 -

z i ) = α i + θ i = α i + θ 1 - i-1 j=1 (β j + π).
The authors in [START_REF] Marshall | Formations of vehicles in cyclic pursuit[END_REF] propose the following control law

v i (t) = ṽ θi (t) = kα i (t) (17)
with ṽ, k > 0 constant parameters, and show that the closedloop dynamics, in the state variables [START_REF] Fabbiano | Distributed source localisation with no position information[END_REF], becomes ḋi (t) = -ṽ cos α i (t) + cos(α i (t)

+ β i (t)) αi (t) = ṽ sin α i (t) + sin(α i (t) + β i (t)) d i (t) -kα i (t) βi (t) = k α i (t) -α i+1 (t) . (18)
Their study shows that the equilibrium points of such dynamics correspond to the vehicles disposed as the vertices of an ordinary regular polygon, and thus being equispatially disposed on a circle (whose radius is given by r = nṽ kπ ); moreover, they analyze stability of such equilibria. The following proposition summarizes some of their most important findings, while in the remainder of this section we will show how such results can be applied to our problem formulation.

Proposition 2 ([16, Lemma 2, Thm. 5 and Coroll. 4]): The system [START_REF] Garin | A Survey on Distributed Estimation and Control Applications Using Linear Consensus Algorithms[END_REF], restricted to the invariant manifold M ⊂ R 3n described by the constraints in ( 15)-( 16), has 2(n -1) equilibria, among which the two following ones:

di = d = 2ṽ k ᾱ sin ᾱ (19A) ᾱi =ᾱ = ± π n (19B) βi = β = π -2ᾱ. (19C)
The two above-described equilibria are locally asymptotically stable.

The interested reader can find the proof in [START_REF] Marshall | Formations of vehicles in cyclic pursuit[END_REF]; it is based on geometrical considerations to characterize equilibria, on linearization around equilibria of the system restricted to the invariant manifold, and on a careful study of the eigenvalues of the resulting linear system. Since our control law is different from the one proposed in [START_REF] Marshall | Formations of vehicles in cyclic pursuit[END_REF], also our closed-loop dynamics will differ from [START_REF] Garin | A Survey on Distributed Estimation and Control Applications Using Linear Consensus Algorithms[END_REF], and the results of Theorem 2 do not apply directly. Nevertheless, defining the variables

αi (t) = arg(z i+1 (t) -z i (t)) -γ i (t) (20A) βi (t) = γ i (t) -γ i+1 (t) -π, ( 20B
)
with γ i as in ( 9), we will show that in the new state d i (t), αi (t), βi (t) our closed-loop system converges to the same equilibrium points as [START_REF] Hendrickx | Finite-time consensus using stochastic matrices with positive diagonals[END_REF]. First, we notice that the above-defined variables satisfy the same constraints as the variables in ( 14), namely

n i=1 βi = -nπ (21)
n i=1 d i e i[ αi-(i-1)π-i-1 j=1 βj ] = 0. (22)
Then, in the following, we prove that, under the control law [START_REF]Cooperative control design for time-varying formations of multiagent systems[END_REF], the variables d i , αi , βi defined in (14A) and [START_REF] Dimakis | Gossip algorithms for distributed signal processing[END_REF] have the same dynamics as [START_REF] Garin | A Survey on Distributed Estimation and Control Applications Using Linear Consensus Algorithms[END_REF], so that the same equilibrium analysis applies.

Proposition 3: For the system [START_REF] Fabbiano | Source localization by gradient estimation based on Poisson integral[END_REF], with ω i as in (6D), the dynamics in the state variables (14A), [START_REF] Dimakis | Gossip algorithms for distributed signal processing[END_REF] is

ḋi (t) = -v r cos αi (t) + cos αi (t) + βi (t) αi (t) = v r sin αi (t) + sin αi (t) + βi (t) d i (t) -k αi (t) βi (t) = k αi (t) -αi+1 (t) . (23)
Proof: We start by showing that the control law (7B) is equivalent to the control law γi = k αi . We define

δ i (t) = θ i (t) -γ i (t).
Taking the sine and using (9), we obtain

sin δ i = sin θ i v i cos θ i -v c cos θ ref v r -cos θ i v i sin θ i -v c sin θ ref v r = v c v r (cos θ i sin θ ref -sin θ i cos θ ref ) = v sin θi , (24)
from which it follows that δ i (t) = arcsin v sin θi (t) . Therefore, we have that

γi = k(α i + δ i ) = k arg(z i+1 -z i ) -θ i + θ i -γ i = k αi . i -1 i i + 1 αi-1 + γ i-1 αi-1 + γ i-1 + π αi + γ i η i
Fig. 4: The internal angle η i between consecutive sensors.

Now we have

żi (t) = v c e iθ ref (t) + v r e iγi(t)
γi (t) = k αi (t).

From this system, the dynamics (23) are obtained, similarly to the fact that system (1), [START_REF] Fabbiano | Collaborative source-seeking control[END_REF] gives dynamics [START_REF] Garin | A Survey on Distributed Estimation and Control Applications Using Linear Consensus Algorithms[END_REF] in [START_REF] Marshall | Formations of vehicles in cyclic pursuit[END_REF] (full details of the calculations can be found in the proof of Proposition 5.3 in [START_REF] Fabbiano | Collaborative source-seeking control[END_REF]). Proposition 4: The system ( 14), restricted to the invariant manifold M ⊂ R 3n described by the constraints in ( 21)-( 22), has 2(n -1) equilibria, among which the two following ones:

di = d = 2v r k ᾱ sin ᾱ ᾱi =ᾱ = ± π n βi = β = π -2 ᾱ. ( 25
)
The two above-described equilibria are locally asymptotically stable.

Proof: By Prop. 3 the system d i (t), αi (t), βi (t) has dynamics equivalent to [START_REF] Garin | A Survey on Distributed Estimation and Control Applications Using Linear Consensus Algorithms[END_REF]. Hence, applying Prop. 2 yields to the equilibria (25).

Proposition 5: At the equilibria (25), the relative positions z i+1z i draw an ordinary regular polygon.

Proof: The internal angle between sensors z i-1 , z i , z i+1 (see Fig. 4) can be expressed as

η i (t) = arg(z i (t) -z i-1 (t)) + π -arg(z i+1 (t) -z i (t)) = αi-1 (t) + γ i-1 (t) + π -αi (t) -γ i (t) = αi-1 (t) -βi-1 (t) -αi (t);
(26) at the equilibrium, αi = ᾱ and βi = β, ∀i, gives

η i (t) = -β, ∀i.
Hence, the relative positions at the equilibrium are both equilateral (because di = d ∀i) and equiangular, and thus represent an ordinary regular polygon.

D. Gradient-Ascent Motion of the Formation Centroid

Our goal is to use the sensors measurements to compute a discrete approximation of the gradient formula [START_REF] Briñón Arranz | Collaborative estimation of gradient direction by a formation of AUVs under communication constraints[END_REF], and then to perform a gradient ascent using the (approximated) gradient argument as reference θ ref for the control law [START_REF]Cooperative control design for time-varying formations of multiagent systems[END_REF]. The computation of the gradient approximation will be discussed in Section IV, while here we discuss how the control law (6) indeed drives the centroid of the formation to follow the reference direction and perform gradient ascent towards the source. We will make use of the result of Prop.s 4 and 5, which ensure local asymptotic convergence to an equispaced circular formation, and we will study the motion of the formation centroid.

Proposition 6: If the initial conditions of the agents are such that the corresponding system (23) is in the attractive basin of one of the two locally asymptotically stable equilibria (25), then arg żc (t)θ ref (t) → 0.

Proof: The dynamics of the formation centroid (i.e., the barycentre of positions z i 's) is given by

żc (t) = 1 n n i=1 żi (t) = v c e iθ ref (t) + v r n n i=1
e iγi(t) .

We can easily show that, at any of the equilibria (25),

żc (t) = v c e iθ ref (t) .
Indeed, the result of Prop. 5 implies that, at equilibria (25), they occupy positions z i which are equispaced along a circle.

In particular, Eq. (20B) implies that γ i+1γ i = -βiπ ∀i which, at the equilibria (25), gives γ i = γ 0 + 2π n i for all i, and hence n i=1 e iγi(t) = 0. Then, the local asymptotic stability of such equilibria ensures the asymptotic convergence of żc (t)v c e iθ ref (t) = vr n n i=1 e iγi(t) to zero, and hence also ensures that arg żc (t)θ ref (t) → 0. This proposition means that, if the agents have initial positions such that they will converge to circular formation, then, while this convergence happens, their centroid will follow with increasing precision the reference direction. In particular, given any ε > 0, there exists a time t such that, for all t ≥ t, |arg żc (t)θ ref (t)| ≤ ε; moreover, velocities of individual agents and of the formation centroid are all bounded by v c +v r , so that there is no risk of divergence in finite time.

In the next section we will discuss how to compute a reference θ ref (t) which is an approximation of the argument of the gradient ∇f (z c (t)) (at least for sufficiently large times, when circular formation is nearly attained). Now, let us assume that such a reference is given, and that for all t ≥ t, |θ ref (t)arg ∇f (z c (t))| ≤ ε for some time t and for some error bound ε < π 2 . Then, we can clearly find ε such that ε + ε < π 2 , and a time t ≥ t such that, for all t ≥ t,

|arg żc (t) -arg ∇f (z c (t))| ≤ ε + ε < π 2 .
This inequality ensures that the formation centroid performs gradient ascent, namely that f (z c (t)) increases along trajectories; indeed, the derivative of f along trajectories of z c (t) is given by the scalar product ∇f (z c (t))• żc (t), which is positive since |arg żc (t)arg ∇f (z c (t))| < π 2 . Gradient ascent is usually applied to convex functions, while here we are considering harmonic functions, which might not be convex. However, the local maximum principle ensures that f does not have any local maximum inside Ω: hence, gradient ascent search is ensured not to get trapped in any local maximum, except possibly on the outer boundary ∂ Ω; some simple rule can be introduced that allows to distinguish the external boundary from the source, e.g., some rough bound on the source intensity.

IV. IMPLEMENTATION DETAILS

We discuss here how to implement the control law proposed in Sect. III. We assume to have no absolute position information, and very little relative one: this means that the sensors neither know their absolute position, nor their relative one with respect to the formation centroid or to each other, and they do not know their absolute orientation θ i either. They are only able to measure the pointwise value of the diffusing quantity at their position, i.e., f (z i (t)), and the relative bearing angle with respect to their neighbour α i ; in addition, they know the common parameters v c , v r , k and their number n.

To drive the sensors to the source location, we make use of the control input defined in [START_REF]Cooperative control design for time-varying formations of multiagent systems[END_REF], which involves the difference θi between the reference input θ ref and the current heading of the sensor θ i , and where we want the reference θ ref to be the angle of the approximation, computed by each sensor i, of the gradient at the formation centroid. To obtain such an approximation we propose a discrete-time algorithm, involving an exchange of messages among the agents, which is run within every sampling time; the value of the gradient approximation g i (t), and consequently of its argument, is then held constant until the next value is available. The algorithm is designed to use the information collected by each sensor and produce, in a distributed way, a discrete approximation of the integral formula (5), i.e., to compute the Riemann sum

∇f (z c (t)) = 2 nr 2 n i=1 f (z i (t))(z i (t) -z c (t)), (27)
and then obtain the difference θi = θ refθ i , where ideally θ ref = arg ∇f . The difficulty lies in the fact that the vector z i -z c is not directly available, as we assumed not to have any position information. We propose the following technique.

During every time interval (t, t + ∆t), with ∆t sufficiently long, each sensor i performs n -1 iterations which involve an exchange of messages with its preceding and following neighbours: each sensor i receives the current value computed by its predecessors, rotates it of an angle κ = ± 2π n (where the sign depends on the direction of the rotation of the agents at the equilibrium), adds its own measurement f (z i ) and sends the result to its successor.

In the sequel we will suppose for simplicity to be in the attractive region of the equilibrium with ᾱ > 0 (i.e., a counterclockwise rotation), hence, we will use a rotation angle of κ = -2π n . A possible more general implementation can use a rotation angle depending on the current relative bearing angle measured, i.e., κ i (t) =sign α i (t) 2π

n . The argument of the vector obtained in this way is the ith reference, which is then used to compute the difference θi and then, together with the relative bearing angle α i which each sensor i measure with respect to its follower, the new control u i . The detailed algorithm is reported below.

Algorithm 1 g i (0) = f (z i )e i0 {initialisation} for l = 1, . . . , n -1 do g i (l) ← g i-1 (l -1){node i receives from i -1} g i (l) ← g i (l)e -i 2π n + f (z i )e i0 {update} end for θ (i) ref = arg g i (n -1) δ i = -arctan vc cos θ (i) ref vr+vc sin θ (i) ref θi = θ (i) ref -δ i -π 2
The following proposition shows the effectiveness of the proposed algorithm if the sensors are already at the equilibrium.

Proposition 7: If the n sensors described by the dynamics (23) are at the equilibrium (25) with ᾱ > 0, then θi computed with Algorithm 1 corresponds to the difference θi = arg ∇f (z c )θ i .

Proof: When the sensors are at the equilibrium (25) with ᾱ > 0, they are turning equispaced in circle around their centroid z c ; hence we have that

γ i = arg(z i -z c ) + π 2 γ 0 + 2π n i , ∀i; (28)
in particular, the second one implies that γ j = γ i + 2π n (ji) for all i, j. Using these properties, we can rewrite the argument of the Riemann sum in (27) as

arg ∇f (z c ) = arg n j=1 f (z j )e i(γj -π 2 ) = arg n j=1 f (z j )e i[γi+ 2π n (j-i)-π 2 ] , ∀i.
Node i, at the end of the n -1 iterations in Algorithm 1, has computed

g i (n -1) = n j=1 f (z j )e i 2π n (j-i) , (29) so that θ
(i) ref = arg g i (n -1) = arg ∇f (z c ) -γ i + π 2 . ( 30
)
We want to show that the solution of system

δ i = -arctan v c cos θ (i) ref v r + v c sin θ (i) ref (31A) θi = θ (i) ref -δ i - π 2 (31B)
and the one of

δ i = θ i -γ i (32A) θi = arg ∇f (z c ) -θ i (32B)
coincide, thus showing that system (31) used in the algorithm correctly computes the quantities of interest defined by (32).

We already know from (24) that δ i = θ iγ i implies sin δ i = v sin θi ; thus, combining it with (30), system (32) becomes

sin δ i = v sin θi θi = θ (i) ref + γ i - π 2 -θ i = θ (i) ref -δ i - π 2 .
Solving for sin θi we get

sin θi = sin δi v -cos θ (i) ref cos δ i -sin θ (i) ref sin δ i ;
equating the solutions and multiplying by v cos δi yields

tan δ i = - v cos θ (i) ref 1 + v sin θ (i) ref ,
and assuming that δ i ∈ -π 2 , π 2 and θi ∈ [-π, π] we obtain

δ i = -arctan v c cos θ (i) ref v r + v c sin θ (i) ref
.

Guaranteeing a common reference

As one may notice, Algorithm 1 gives the argument of the correct gradient approximation only when the sensors are in formation. Otherwise, not only θ (i) ref is not ensured to be a good approximation of the gradient argument, but it is not even a common reference, i.e., θi (t) + θ i (t) is not the same for all i (while it should be equal to a common reference in order to apply Prop. 1). Indeed, Eq. ( 28) holds only at the formation. This contradicts the assumption we made about the reference θ ref being the same for every sensor, which guarantees the convergence to the equilibria (25) and thus to a circular equispaced formation; nevertheless, we proposed this implementation for its simplicity, since it shows a good behaviour in simulations, as we will show in Section V.

Below, we propose an alternative algorithm which ensures that all sensors have a common reference θ ref , regardless their relative positions and orientations.

To do so, we need to designate a leader sensor, say node 1. Moreover, we need the measure of the relative bearing angle with respect to the predecessor; i.e., in addition to the alreadydefined angle α i , we assume that each sensor i can also measure the angle

α - i (t) = arg(z i-1 (t) -z i (t)) -θ i (t).
The idea is to perform, after a first phase as in Algorithm 1, a second communication round enforcing θi = θ refθ i for every i, with

θ ref = θ1 + θ 1 .
Algorithm 2 g i (0) = f (z i )e i0 {initialisation, ∀i} for l = 1, . . . , n -1 do g i (l) ← g i-1 (l -1){node i receives from i -1} g i (l) ← g i (l)e -i 2π n + f (z i )e i0 {update} end for θ Proof: It is easy to see that, with Algorithm 2, θ1 is the same as in Algorithm 1, and

θi+1 = θi + α - i+1 -α i -π = θi + arg(z i -z i+1 ) -θ i+1 -arg(z i+1 -z i ) + θ i -π = θi + θ i -θ i+1 ,
from which, recursively, we get θi = θ1 + θ 1θ i .

Alternative communication strategies

Algorithms 1 and 2 are based on a ring communication topology, namely messages are sent from sensor i to i + 1, but they can be easily adapted to a variety of other communication setups, by making use of well-known linear average-consensus algorithms. Consider a given communication graph describing which sensors are allowed to send messages to one another, and a linear consensus algorithm, described by n × n matrices P (0), P (1), . . . , where P ij is 0 if there is no edge from j to i. At each iteration of the consensus algorithm, each node i sends its current state, say ξ i (l), to its neighbors, and then computes the state update ξ i (l + 1) = j P ij (l)ξ j (l), using the received neighbors states and its own state. Under mild assumptions on the communication graph, there are well-known constructions of asymptotic average consensus algorithms, such that lim l→∞ ξ(l) = 1 n j ξ j (0) (e.g., with time-invariant doubly-stochastic P [START_REF] Garin | A Survey on Distributed Estimation and Control Applications Using Linear Consensus Algorithms[END_REF]), of finite-time average consensus algorithms, such that ξ(L) = 1 n j ξ j (0) for some L [START_REF] Hendrickx | Finite-time consensus using stochastic matrices with positive diagonals[END_REF], and gossip consensus, where P (k)'s are random matrices consistent with random asynchronous communications while ensuring some stochastic convergence near the average of initial values [START_REF] Dimakis | Gossip algorithms for distributed signal processing[END_REF]. Algorithm 1 can be modified as follows, where we use a given linear consensus algorithm, but we rotate state vectors at each iteration in a suitable way. Algorithm 3 g i (0) = f (z i )e i0 {initialisation} for l = 1, . . . , L do i receives g j (l-1) from all neighbors j{communication} g i (l) ← j P ij (l -1)e i 2π n (j-i) g j (l -1){update} end for θ

(i) ref = arg g i (L) δ i = -arctan vc cos θ (i) ref vr+vc sin θ (i) ref θi = θ (i) ref -δ i -π 2
At the end of the cycle of communications and updates, we get

g i (L) = j L-1 • • • j0 P ij L-1 (L-1)•. . .•P j1j0 (0)e i(j0-i) f (z j0 ).
If we define ξ j (0) = f (z j )e i 2π n j , the standard linear consensus (with no rotation) would give

ξ i (L) = j L-1 • • • j0 P ij L-1 (L -1) • . . . • P j1j0 (0)f (z j0 )e i 2π n j0 .
Hence, g i (L) = e i(-2π n i) ξ i (L), and inherits the convergence properties of the chosen consensus algorithm: if this is a finitetime average consensus, we get

g i (L) = e i(-2π n i) 1 n j ξ j (0) = 1 n j f (z j )e i 2π n (j-i) .
Comparing this with (29), we can see that the only difference with Algorithm 1 is a factor 1 n , which is irrelevant when computing the argument θ (i) ref = arg g i (L), so that the result of Algorithm 3 is exactly the same as the one of Algorithm 1. When consensus only has asymptotic convergence, the same applies to Algorithm 3, which then gives an approximation of the value obtained in Algorithm 1.

In a similar way, also Algorithm 2 can be adapted, provided that each node can measure the relative angle α ij = arg(z jz i )-θ i with respect to all of its in-and out-neighbors. The first part, to find θ1 , follows the same technique as in Algorithm 3. Then, node 1 initializes its state to χ 1 (0) = e i θ1 , while all other nodes set χ i (0) = 0. Again, a modified consensus is run, where messages are suitably rotated versions of the local states: at iteration l, node i receives e i(-αji-π 2 ) χ j (l -1) from each neighbor j, and then computes

χ i (l) = j P ij (l -1)e i(αij -π
2 ) e i(-αji-π 2 ) χ j (l -1).

After L iterations, each node i computes θi = arg χ i (L), thus ending the algorithm.

It is easy to show that, with finite-time average consensus, χ i (L) = 1 n e i( θ1+θ1-θi) , so that its argument gives exactly the same θi as in Algorithm 2. In case of asymptotic convergence to consensus, this method only allows an approximation of a common reference, without guaranteeing it exactly. To illustrate our algorithm, we consider a set of n sensors with dynamics (1) and with the control law described in [START_REF]Cooperative control design for time-varying formations of multiagent systems[END_REF]. The reference direction approximating the gradient is computed with Algorithm 1, and its time derivative θref is obtained with discrete backward approximation.

Figure 5 shows the result of three simulations with n = 5 sensors, starting from random initial positions and with random initial orientations, with velocities of v r = 0.5 and v c = 0.1 m/s; the values for the control gain k is chosen such to have a formation radius of r = 0.5 m, i.e., k = nvr πr ≈ 1.5915. We can see the trajectory of the sensors as well as the one of their centroid (in black); the black-bordered triangles represent the initial conditions of each sensor, while the whitebordered ones show their position and orientation when the formation reached the heater (the circle in the middle). The sensors reach a circular configuration soon, and then keep moving towards the source following the gradient of the field, approximately. Even if the theoretical convergence is local (for initial conditions near circle formation), the algorithm demonstrated to be effective also for randomly generated initial conditions. Moreover, despite the theoretical request to have the same common reference for each sensor in order to guarantee convergence, in practice Algorithm 1, which does not fulfill this assumption, gives good results.

Figures 6 and7 show simulations similar to the ones in Figure 5, but illustrating the effect of noise and delay, respectively. In Figure 6, both temperature and bearing angle measurements are corrupted by white Gaussian noise of standard deviation σ = 0.5. In Figure 7 the updates of the reference angle θ ref happen with a delay of 1s, so as to show that the behaviour is not very degraded even in a case where distributed computation of θ ref is heavily affected by delays in transmissions and computations.

Our algorithm can be applied with a different number of sensors, the minimum number being three, because two measurements are not enough to gather information about a 2dimensional gradient. In our simulations (not reported for lack of space), we have seen that three sensors are enough to have a nice behaviour, qualitatively similar to the one illustrated here with five sensors. A larger number of sensors ensures a better gradient approximation (at least when all sensors are in circular formation), but on the other hand the delay in the distributed computation of θ ref grows linearly with n, due to the round of communications and updates in Algorithm 1. In Figure 8, we show the trajectories obtained with different values of n, taking into account the corresponding delay. Sensors are initially placed in equispaced formation around a same circle, and we plot only the trajectory of the formation centroid. The parameters used in these simulations are v r = 0.5 and v c = 0.1 m/s; the gain k = nvr πr whith radius r = 1 m. Due to the use of Algorithm 1, θ ref is computed with a delay which grows linearly with n. However, the delay duration strongly depends on the software and hardware used for computation and communication. In Figure 8, we consider an example of application using MATLAB © on a general-purpose computer with a 2GHz clock, and an Arduino board. In this case, temperature acquisition and wireless IrDA transmission take negligible time, while computations (which need to be done sequentially: node i waits for the result of the computation at node i -1), results in a total delay of approximately 5n ms. Such a delay does not have a major impact on the trajectories, and with growing n the trajectory of the centroid becomes a better2 approximation of the ideal gradient ascent, which is plotted in black line for comparison. We have simulated a case where delay grows more dramatically, as 25 n ms instead of 5 n ms. For n ≤ 9 this delay results in trajectories similar to the ones with smaller delay depicted in Figure 8, while for n = 10 the centroid trajectory is visibly degraded, as shown by the dashed line added to the figure. This suggests to use a small number of sensors in applications with high delays. If a large number of sensors is necessary, e.g., because a large formation radius is required and communication range is small, it is still possible to mitigate the effect of delays, by suitably tuning the speed parameters v c and v r , although at the price of a slower motion towards the source.

VI. CONCLUSIONS

We considered the problem of localising a source by a set of communicating mobile sensors providing pointwise concentration measures and able to measure the relative bearing angle with respect to their neighbour. We have focused our attention on steady-state homogeneous isotropic diffusion processes, and developed a control law which allows the sensors to compute, in a distributed way, an estimation of the gradient argument of the function describing the diffusion of the measured quantity; such a control law did not use any explicit expression for the solution of the PDE, nor any global position information. The convergence of the source-seeking has been theoretically proven and, as an illustrative example, we have considered the search of a heat source in a 2-dimensional room.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Relative bearing angle α i

Fig. 3 :

 3 Fig. 3: Sets as defined in Section II: Ω = Ω \ Ω s , and ∂Ω = ∂ Ω ∪ ∂Ω s .

( 1 )Proposition 8 :

 18 ref = arg g 1 (n -1){only node 1} δ =arctan vc cos θ δ -π 2 {only node 1} for i = 1, . . . , n -1 do ι = θiα i -π 2 {node i} node i + 1 receives ι from i θi+1 = ι + α - i+1 -π 2 {node i + 1} end for If each sensor i computes θi as in Algorithm 2, then θi (t) = θ ref (t)θ i (t) ∀i, where θ ref = θ1 + θ 1 ,and θ1 is the same as in Algorithm 1.

  algorithm we consider the search of a heat source in a 2-dimensional room 20 × 12 m large; the walls are considered perfectly isolating, and we suppose to have an open window with an external temperature of T ext = 5 °C, and a circular heater of T s = 50 °C in the middle. Our scenario is described by the problem with mixed Dirichlet-Neumann boundary conditions f (z) = T ext on the window border f (z) = T s on the source border ∇f (z) • n = 0 on the walls ∇ 2 f (z) = 0 inside the room, with n the outward unit normal to the wall.

Fig. 5 :

 5 Fig. 5: Three simulations with five sensors in random initial positions.

Fig. 6 :

 6 Fig. 6: Three simulations with five sensors in random initial positions, with noisy measurements.

Fig. 7 :

 7 Fig. 7: Three simulations with five sensors in random initial positions, with delay in the computation of θ ref .

Fig. 8 :

 8 Fig. 8: Trajectories of the formation centroid, from initial equilibrium formation around a same circle, for different values of n. The computation of θ r is affected by delay growing linearly with n, as 5 n ms. Ideal gradient ascent is shown for comparison. Dashed line has n = 10 and higher delay 25 n = 250 ms.

We use the notational convention that a sum over an empty set of indices is equal to zero, namely, 0 j=1 β j = 0

In our simulations with same parameters as in Figure8, we have always obtained that larger n brings trajectories nearer to the ideal gradient ascent curve. In many cases, as in Figure8, this also gives better trajectories, both in terms of time and distance necessary to reach the source. However, this is not true in general, in some cases a deviation from gradient ascent due to lower n can also lead to a shorter trajectory.

ACKNOWLEDGMENTS

The authors would like to thank Antoine Rousseau for the FreeFem++ code used for the finite elements simulation.


Karim Hitti
email: khitti@ndu.edu.lb

Stephanie Feghali
email: stephanie.elfeghali@gmail.com

Marc Bernacki
email: marc.bernacki@mines-paristech.fr

PERMEABILITY COMPUTATION ON A REPRESENTATIVE VOLUME ELEMENT (RVE) OF UNIDIRECTIONAL DISORDERED FIBER ARRAYS

Keywords: Mathematics subject classification: 65N30, 76S05 Permeability, Fibrous Media, Finite Element Method, Level-set

An efficient method to compute the permeability of disordered fibrous arrays is proposed. A stabilized mixed finite element method is used with an immersed domain approach to represent the porous material at its microscopic scale. Therefore, the Stokes equations are solved in the whole domain (including solid part) using a penalization method. The accuracy is controlled by refining the mesh around the fluid-solid interface defined by a level-set function. Using homogenization techniques, the permeability of an RVE is obtained. Furthermore, a new method to generate disordered fibers in function of the porosity, φ, and other microstructural parameters is proposed and a study of the effect of inter-fiber spacing on K, the permeability tensor, is performed. This task was achieved using parallel computation and over 460 simulations were carried out in two-dimensional RVEs consisting of over 555 fibers.

Introduction

Liquid Composite Moulding (LCM) processes are injection processes used for manufacturing large and complex composite materials with fiber reinforcements. Numerical simulation at the macroscopic scale is based on the resolution of Darcy's law [START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] Darcy | Les fontaines publiques de la ville de Dijon[END_REF] to predict flow front progression, filling time and injection pressure and to improve the design of tools and molds.

In the case of ordered fiber arrays, several analytical relations have been established to predict the permeability of fibrous media [START_REF] Carman | L'écoulement des gaz à travers les milieux poreux[END_REF][START_REF] Gebart | Permeability of unidirectional reinforcements for RTM[END_REF][START_REF] Happel | Viscous flow relative to arrays of cylinders[END_REF][START_REF] Hoareau | Injection sur renfort : Etude du remplissage de moule et détermination théorique de la perméabilité des tissus[END_REF][START_REF] Keller | Viscous flow through a grating or lattice of cylinders[END_REF][START_REF] Tamayol | Analytical determination of viscous permeability of fibrous porous media[END_REF]. All of these studies consider simple geometries, such as square or triangular packing of unidirectional arrays of cylinders. And the analytical relations are only a function of φ. By considering some hypothesis, different authors provide analytical solutions of the Stokes equations to obtain the permeability. These different laws give a valid solution for different fiber volume fractions. In what concerns numerical analysis, several studies had been done at the microscopic and mesoscopic scales through the finite element or finite volume methods [START_REF] Gebart | Permeability of unidirectional reinforcements for RTM[END_REF][START_REF] Hoareau | Injection sur renfort : Etude du remplissage de moule et détermination théorique de la perméabilité des tissus[END_REF][START_REF] Nordlund | Numerical study of the local permeability of noncrimp fabrics[END_REF][START_REF] Song | Asymptotic expansion homogenization of permeability tensor for plain woven fabrics[END_REF]. In these studies, fibers or yarns were considered as impermeable solids. Then, only the fluid domain was meshed and no-slip boundary conditions were imposed on solid boundaries. Velocity and pressure fields were computed on the fluid domain and permeability of the volume was then obtained by a homogenization method.

In the case of disordered fiber arrays, which represent real fiber performs, the use of porosity alone cannot define their permeability. Other microstructural parameters should be taken into account. Chen and Papathanasiou [START_REF] Chen | Micro-scale modeling of axial flow through unidirectional disordered fiber arrays[END_REF][START_REF] Chen | The transverse permeability of disordered fiber arrays: a statistical correlation in terms of the mean nearest interfiber spacing[END_REF] studied the effect of the mean nearest interfiber spacing, denoted δ 1 , or the degree of disorder, on both the transverse and longitudinal permeabilities. They found that the latter decreases on all porosity levels and the former increases on porosity levels ranging from 0.45 to 0.7 and decreases above these levels when δ 1 increases. In these studies a Monte-Carlo procedure [START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF] was used to generate the fiber distributions which are governed by the choice of the porosity and the minimum allowable inter-fiber distance δ min .

In this paper, we propose an effective method to generate fiber distributions with specific porosities. It is based on advancing front methods [START_REF] Bagi | An algorithm to generate random dense arrangements for discrete element simulations of granular assemblies[END_REF][START_REF] Benabbou | Numerical modeling of nanostructured materials[END_REF] and dropping and rolling techniques [START_REF] Hitti | Optimized Dropping and Rolling (ODR) method for packing of polydisperse spheres[END_REF]. Also, the inter-fiber distances are chosen either by imposing a δ min or generated by a Gaussian distribution law. The influence of the standard deviation of this law, which is proportional to the fiber array's degree of disorder and inversely proportional to Chen and Papathanasiou's δ 1 , on the permeability is studied. Moreover, a new method to compute permeability of a fiber array is proposed where a monolithic approach to solve fluid flow equations is used, followed by a homogenization method to compute permeability on the whole domain. Using a monolithic approach, a unique equation is solved on a mesh containing both fluid and solid domains, and, in our case, using a mixed finite element approximation. The interface is implicitly represented by the zero iso-surface of a level-set function and a penalization method is used to take into account the motion of the solid part [START_REF] Coupez | Multidomain finite element computations: Application to multiphasic problems, Arbitrary Lagrangian-Eulerian and Fluid-Structure Interaction[END_REF].

All numerical calculations mentioned in this paper were performed with the CimLib finite element C++ library [START_REF] Digonnet | Cimlib: a fully parallel application for numerical simulations based on components assembly[END_REF].

Model Equations

Governing equations

The injected fluid is considered incompressible Newtonian. Due to the fluid's low injection pressure and its high viscosity, inertia and gravity terms can be neglected. Consequently, Stokes equations, describing fluid flow, are written as:

η f ∆v -∇p = 0 ∇ • v = 0 (2.1)
with v the fluid velocity, p the pressure and η f the dynamic viscosity. Darcy's law is traditionally used to model flow motion in porous media at the macroscopic scale. These macroscopic equations are obtained from a volume average of the Navier-Stokes ones [START_REF] Whitaker | The Method of Volume Averaging, Theory and Applications of Transport in Porous Media[END_REF][START_REF] Pillai | Governing equations for unsaturated flow through woven fiber mats. Part 1. Isothermal flows[END_REF][START_REF] Tucker | Flow and Rheology in Polymer Composites Manufacturing[END_REF], describing the flow motion at the microscopic scale.

We will use the definitions of the averages given in [START_REF] Tucker | Flow and Rheology in Polymer Composites Manufacturing[END_REF], where the spatial average of a quantity B is defined in the whole volume:

B = 1 V Ω B dΩ . (2.2)
The intrinsic phase average of a quantity B f is defined in the fluid phase:

B f f = 1 V f Ω f B f dΩ , (2.3)
with Ω f the fluid domain.

The porosity φ and the fiber volume fraction φ s are defined by the ratio between the volume occupied by the fluid V f and total volume V :

φ = V f V and φ s = 1 -φ . (2.4)
We suppose that the solid skeleton is static and non-deformable, and that the porous media is saturated. Neglecting effects of viscosity on the flow with respect to the resistance of the skeleton, Darcy's equation is obtained:

v = - 1 η f K • ∇p f . (2.5)
The permeability K is a tensor for anisotropic porous media, like fibrous media for example. In the case of unidirectional fiber packings, if the z axis of the coordinate system is taken in the same direction as the axis of the fibers, the permeability tensor can be then written as follows:

K = K ⊥ 0 0 K , (2.6)
with K ⊥ a transverse permeability tensor and K the permeability along the fiber axis. For regular packings (rectangular, triangular or hexagonal) the matrix K ⊥ is proportional to the unit matrix and then only a scalar K ⊥ has to be computed. However for disordered arrays, we need to look for a more suitable form:

K ⊥ = K xx K xy K xy K yy or K ⊥ = K xx 0 0 K yy .
(2.7)

Immersed Domain Method for the Computation of Flow Motion

First, the Stokes equations (2.1) are solved on the RVE by an immersed domain method [START_REF] Hitti | Direct numerical simulation of complex Representative Volume Elements (RVEs): Generation, resolution and homogenization[END_REF]. The approach proposed here can rather be seen as an extension of the Stokes problem with two fluids, for which the solid behavior of one phase is imposed by using the viscosity as a penalty factor. All the phases of the multi domain problem are implicitly represented by a level-set function α and the interface between the two phases is defined by its iso-zero [START_REF] Coupez | Multidomain finite element computations: Application to multiphasic problems, Arbitrary Lagrangian-Eulerian and Fluid-Structure Interaction[END_REF].

To obtain an accurate description of the interface and small errors concerning the computed flow, efficient anisotropic adaptive meshing tools are used [START_REF] Bernacki | Finite element model of primary recrystallization in polycrystalline aggregates using a level set framework[END_REF]. The viscosity η(α) is defined on the whole computational domain using a mixture law [START_REF] Hachem | Stabilized finite element solution to handle complex heat and fluid flows in industrial furnaces using the immersed volume method[END_REF] (see Figure 3.1). The rigidity of the solid part is taken into account with a high value of the viscosity (namely η s = 10 3 η f [START_REF] Laure | Numerical methods for solid particles in particulate flow simulations[END_REF]), this viscosity η s acts as a penalty coefficient. Using this approach, a zero velocity boundary condition imposed on the solid's interface is propagated into the entire solid domain.

The FE formulation

First, the function spaces used in this paper are defined. The scalar function space for the pressure P = L 2 (Ω), the function space for the velocity V = (H 1 (Ω)) d and the weighting function space V 0 = (H 1 0 (Ω)) d , where d is the space dimension. L 2 (Ω) is the Lebesgue space of square assumable functions on Ω and H 1 (Ω) is the Sobolev space included in L 2 (Ω), are respectively defined by:

L 2 (Ω) = {q/ Ω q 2 dΩ < ∞} H 1 (Ω) = {q ∈ L 2 (Ω)/∇q ∈ L 2 (Ω)} (3.1)
H 1 0 (Ω) = {q ∈ H 1 (Ω)/q = 0 on ∂Ω)}
The FE formulation begins by writing the weak form of the Stokes equations (2.1). The problem consists in finding (v, p) ∈ (H 1 (Ω)) d × L 2 (Ω) such that:

2 Ω η(α)ε(v) : ε(w)dΩ -Ω p∇.wdΩ = 0 Ω q∇.vdΩ = 0 (3.2)
where the trial function q for the pressure is defined in L 2 (Ω) and the trial function w for the velocity is defined in (H 1 0 (Ω)) d . The domain Ω is decomposed into N el triangles T such that they cover the domain and there are either disjoint or share a complete edge. The triangulation will be denoted T h . Using this partition, the above-defined functional spaces are approached by a finite dimensional spaces spanned by continuous piecewise polynomials such that:

Q h = {q h ∈ C 0 (Ω)/q h|T ∈ P 1 (T ), ∀T ∈ T h } V h = {v h ∈ (C 0 (Ω)) d /v h|T ∈ (P 1 (T )) d , ∀T ∈ T h } (3.3) V h,0 = {v h ∈ V h /v h|Γ = 0}
where P 1 (T ) is the vector space of first degree polynomials defined on the simplex T . The Galerkin discrete problem consists now in solving the mixed problem by finding the pair

(v h , p h ) ∈ V h × Q h such that ∀(w h , q h ) ∈ V h,0 × Q h : 2 Ω η(α)ε(v h ) : ε(w h )dΩ -Ω p h ∇.w h dΩ = 0 Ω q h ∇.v h dΩ = 0 (3.4)
It is well known that the P 1 -P 1 approximation for the discrete Stokes problem (3.4) gives an unstable formulation, since the Brezzi-Babuska condition is not respected [START_REF] Franca | On the stability of residual-free bubbles for convection diffusion problems and their approximation by a two-level finite element method[END_REF]. For this reason, we use the MINI-element (see Figure 3.2), in which the functional space V h for the velocity is enriched with the space of the bubble function V ′ :

V ′ = {v ′ , v ′ | Ti ∈ P 1 (T i ) ∩ H 1 0 (T i ), ∀T ∈ T h , i = 1, ..., D}, (3.5)
where D is the topological dimension and T i is a decomposition of T in D sub-triangles, that have as a common vertex the barycenter, G T , of T . In other words, the choice of this bubble function is continuous inside the element, considered as linear on each sub-triangle and vanishes at the boundary of T . The velocity field is now an element of the function space generated by the following direct sum:

V h = V h ⊕ V ′ . (3.6)
Hence, on each element T , the solution takes the form:

v h | T = D i v i T N i + v ′ T b T , (3.7)
where N i is the interpolation function associated with node i. The bubble function b T is defined on each element T as follows: b T = 0 on ∂T ; b T (G T ) = 1, Furthermore, b T satisfies the orthogonality condition: Using the orthogonality properties of the bubble space, the discrete Stokes equation (3.4) is then written:

T ∂ x k N i ∂ x l b T dΩ = 0 ∀k, l; ∀N i ∈ P 1 (T ) . (3.8)
find (v h , v ′ , p h ) ∈ V h × Q h such that        2 Ω η(α)ε(v h ) : ε(w h )dΩ -Ω p h ∇.w h dΩ = 0 2 Ω η(α)ε(v ′ ) : ε(w ′ )dΩ -Ω p h ∇.w ′ dΩ = 0 Ω q h ∇.(v h + v ′ )dΩ = 0 (3.9)
This system can be put naturally under the following matrix form:

  A vv 0 A t vp 0 A bb A t bp A vp A bp 0     v h v ′ p h   =   0 0 0   , (3.10)
where:

A vv = 2 Ω η(α)ε(v h ) : ε(w h )dΩ , A bb = 2 Ω η(α)ε(v ′ ) : ε(w ′ )dΩ , A vp = Ω p h ∇.w h dΩ , A bp = Ω p h ∇.v ′ dΩ . (3.11)
Finally, the static condensation process [START_REF] Arnold | A stable finite element for the stokes equations[END_REF] consisting in solving the second line for the bubble function v ′ which by inserting into the third line of (3.10) results the final matrix system to solve which is the condensed matrix scheme for large-scale unknowns v h and p h reading:

A vv A t vp A vp C v h p h = 0 0 , (3.12)
with C = -A bp A -1 bb A t bp .

The Representative Volume Element (RVE)

In order to get suitable averaged quantities describing the macroscopic behavior of the material, it is necessary to perform computations on an RVE. The RVE is a two dimensional rectangle of width, L and height H on which three different boundary conditions are applied:

1. the pressure gradient is in the horizontal direction x, v y is null on the external boundary and

∂ y v x (x, 0) = ∂ y v x (x, H) = 0 ; v z (x, y) = 0, ∂ z p(x, y) = 0 ; p(0, y) = P, p(L, y) = 0. (4.1)
2. the pressure gradient is in the vertical direction y, v x is null on the external boundary and

∂ y v y (0, y) = ∂ y v y (L, y) = 0 ; v z (x, y) = 0, ∂ z p(x, y) = 0 ; p(x, 0) = P, p(x, H) = 0. (4.2)
3. the pressure gradient is in the longitudinal direction z, v x and v y are null,

∂ x v z (0, y) = ∂ x v z (L, y) = 0 ; ∂ y v z (x, 0) = ∂ y v z (x, H) = 0 ; ∂ z p(x, y) = P/L . (4.3)

Disordered packing

In order to generate disordered fiber packing with pre-chosen porosity, radius and inter-fiber spacing, a sphere packing method called the "Inwards Packing Method (IPM)" [START_REF] Bagi | An algorithm to generate random dense arrangements for discrete element simulations of granular assemblies[END_REF][START_REF] Benabbou | Numerical modeling of nanostructured materials[END_REF], coupled with a "dropping and rolling" algorithm [START_REF] Hitti | Optimized Dropping and Rolling (ODR) method for packing of polydisperse spheres[END_REF][START_REF] Hitti | Precise generation of complex statistical Representative Volume Elements (RVEs) in a finite element context[END_REF] and a disc shrinking technique, were chosen to generate the initial medium. At first, tangent discs are generated with radii bigger than the wanted radius and equal to R + δ i /2 where R is the wanted radius and δ i is the i th inter-fiber spacing. And then, the radii are decreased by removing the δ i /2 leaving for two fibers f i and f j , which were in contact after the IPM and "dropping and rolling" algorithms, with an inter-fiber spacing of δi+δj 2 . Furthermore, the choice of δ i has a great influence on the porosity and on the degree of disorder of the fiber arrays. Hence, we propose the use of a Gaussian distribution law for the choice of the inter-fiber spacings. The reason for choosing a Gaussian distribution law comes from the fact that, contrary to other classical probability distributions, the sum of independent Gaussian distributions N (µ i , σ 2 i ) is the Gaussian distribution N i µ i , i σ 2 i . So when the inter-fiber distance is wanted according to N (δ, σ 2 ), it is actually generated following N (δ/2, σ 2 /2) and added to the wanted radius then when subtracted, the inter-fiber distance between two fibers will follow N (δ, σ 2 ). Furthermore, a simple, yet efficient, choice of δ mean , which is the mean of the Gaussian distribution law, leads to obtaining the wanted porosity. This choice is influenced by the fact that the optimal density, ρ opt , of equally spaced discs (i.e. the hexagonally packed) is given by (see Figure 4.1):

ρ opt = 2πR 2 √ 3(2R + δ) 2 . (4.4)
But since the IPM does not reach the optimal density, this equation is multiplied by ρ(R) which is the disc density obtained by the IPM using the radius R. Finally, the wanted porosity, equal to 1ρ opt .ρ(R), yields that:

δ mean = R 2πρ(R) (1 -φ) √ 3 -2 . (4.5)
The standard deviation of the Gaussian distribution is written σ = δmean-δmin

3

, where δ min is the minimal allowable inter-fiber spacing. Typical fiber distributions are shown in Figure 4.2 for different values of φ and σ.

It is noticed that σ has a greater effect on the fiber distribution when the porosity is large and that it is proportional to the fiber array's degree of disorder. Moreover, by varying σ a spectrum of fiber distributions can be generated at the same φ which means that this standard deviation is linked to the microstructure heterogeneity and gives a good insight into the fiber microstructure.

Algorithm of fiber packing generation

1. The wanted fiber radius R and inter-fiber spacing are chosen. And ρ opt is computed by (4.4).

2. The IPM is performed using the radius R to get ρ(R).

3. δ mean and σ of the Gaussian distribution are computed so that the δ i will follow N (δ mean ; σ).

4.

The IPM coupled with dropping and rolling is used to generate discs, representing the fibers, with radii R + δ i /2.

5.

The discs are shrinked by removing δ i /2 and thus obtaining the wanted fiber packing.

Ripley's K r

As one of the objectives of this study is to compute the permeability of different random fiber arrangements as a function of other parameters beyond porosity, the randomness of the microstructure needs to be properly quantified. The Ripley's K r function [START_REF] Ripley | Modelling spatial patterns[END_REF][START_REF] Ripley | Spatial Statistics[END_REF] can be used to differentiate between regular, Completely Spatially Random (CSR), and clustered point patterns. This function is defined as:

K r (r) = A N 2 N k=1 I k (r) w k (r) . (4.6)
where I k (r) is the number of points found within a distance r of the point k, N is the total number of points in the area of interest A, and w k (r) is a correction factor taking account of the fact that it is possible that only a part of the observation area πr 2 falls within the area of interest A. This correction factor is computed as the proportion of the perimeter of the circle centered in k with radius r which is inside the area of interest [START_REF] Goreaud | On explicit formulas of edge effect correction for ripley's k-function[END_REF]. K r (r) describes characteristics of point patterns at many length scales. The K r -function of a Poisson distribution is πr 2 and draws a dividing line between a regular and a clustered pattern. Estimates of K r (r) are expected to be smaller than πr 2 if the points form a regular pattern, and to be larger than πr 2 in the presence of clustering [START_REF] Diggle | Statistical Analysis of Spatial Point Patterns[END_REF]. Furthermore, the extent of the deviation of K r (r) from πr 2 and the length scale at which such deviations occur give some additional insight into a microstructure. The L-plot, a linearized plot of K r (r) defined as: L(r) = K r (r)/π, is frequently used to show these deviations and the length scales at which they occur. The L-plot of a Poisson distribution is simply a straight line of 45-degree slope through the origin. Figure 4.3 shows the resulting L(r)-function applied to the fiber distributions in Figure 4.2. It is evident that at large distances all fiber distributions approach the CSR pattern. This concludes that our method can generate random distributions usable in this type of studies.

Algorithm of Ripley's K r (r)

1. for r = 0 to 20 * f iber ′ s radius K r (r) = 0

for i = 0 to Number of fibers (N f ) a) I(r) = 0 b) Compute the distance, d, between f iber i and f iber j (i = j) if d ≤ r I(r) = I(r) + 1 end if c)
Compute the edge correction factor w(r); 28 cases exist (see [START_REF] Goreaud | On explicit formulas of edge effect correction for ripley's k-function[END_REF])

end for i K r (r) = Area of RV E N 2 f N f i=1 I(r) w(r)

Permeability Computation

Computation of the transverse permeability tensor

We have chosen to use the least squares method to compute the transverse permeability. The computations are performed with the first two sets of boundary conditions for 2D Stokes equations (2.1).

First assuming that the permeability tensor K ⊥ has the form

K ⊥ = K xx K xy K xy K yy , (5.1)
and second, that 2 computations were performed which gives v i and ∇p i f with i = 1, 2.

The least square method consists in minimizing the function J

J(K xx , K yy , K xy ) = 2 i=1 K xx ∂ x p i f + K xy ∂ y p i f + η f v x,i 2 + K xy ∂ x p i f + K yy ∂ y p i f + η f v y,i 2
.
(5.

2)

The minimum of Eq.(5.2) is obtained when the partial derivatives of the function J with respect to the permeability components are zero i.e.

∂J ∂K xx = 0 , ∂J ∂K yy = 0 , ∂J ∂K xy = 0 . (5.3)
This leads to the following system of equations:

AK xx + CK xy = E , BK yy + CK xy = F , CK xx + CK yy + (A + B)K xy = G , (5.4)
with

A = 2 i=1 ∂ x p i f 2 , B = 2 i=1 ∂ y p i f 2 , C = 2 i=1 ∂ x p i f ∂ y p i f , E = -η f 2 i=1 ∂ x p i f v x,i , F = -η f 2 i=1 ∂ y p i f v y,i , G = -η f 2 i=1 ∂ y p i f v x,i + ∂ x p i f v y,i ,
which gives

K xx = -BCG + C 2 F -C 2 E + EAB + EB 2 d , K yy = A 2 F -CAG + C 2 E -C 2 F + ABF d , K xy = BAG -BCE -F AC d ,
where d = (A + B)(AB -C 2 ).
If K xy is assumed to be zero, then,

K xx = E A , K yy = F B ,
If K xy = 0 and K xx = K yy , a scalar transverse permeability K can be computed as

K = E + F A + B = -η f 2 i=1 ∂ x p i f v x,i + ∂ y p i f v y,i 2 i=1 ∂ x p i f 2 + ∂ y p i f 2 = -η f 2 i=1 v i • ∇p i f 2 i=1 ∇p i f • ∇p i f .
(5.5)

Computation of longitudinal permeability

The longitudinal permeability is obtained by computing the scalar v z solution of the reduced Stokes equation

η f ∂ 2 v z ∂y 2 + ∂ 2 v z ∂x 2 = ∂p ∂z on Ω f , (5.6)
and

K = -η f v z ∂ z p .
(5.7)

Numerical Results

The RVE's size

The randomness of the generated fiber distributions leads to a scatter in the computed permeability. Hence, an average permeability should be calculated along with a data variation in order to be statistically representative. In this work, 20 random microstructures were generated for each class of fiber distributions characterized by the same φ and σ. The average permeability and its standard deviation, σ p , are given by:

K = 1 20 20 i=1 K i , (6.1)
σ p = 1 19 20 i=1 (K i -K ) 2 . (6.2)
Moreover, the permeability of fiber distributions is not only a function of microstructural parameters but also a function of the RVE's size which is the same as the number of fibers N f in the domain. An extremely large computational domain is surely representative of the fiber array but this will cause a drastic increase in the computational cost which makes this type of studies quasi-impossible. In order to determine the appropriate size of the RVE a study of the size effect on K , K and φ was performed and computations were carried out on elementary volumes of different size extracted from microstructures generated with σ = 0.06R and φ = 0.5 as illustrated in Figure 6.1. The results are shown in Figure 6.2 where the permeability values are adimensionalised. Fluctuations are observed in the values of K /R 2 , K /R 2 and φ when a relatively small number of fibers is considered then the values reach a plateau after N f > 385. Also, after N f > 555 the associated standard deviations are reduced. As a conclusion, an RVE's size of N f > 555 is used in the remainder of this paper.

Isotropic properties

We computed the values of the components of the transverse permeability tensor K ⊥ , K xx , K yy and K xy , in order to see if the scalar transverse permeability K can be computed. In Figure 6.3, the averaged components are computed at different porosity levels and plotted against σ/R. The values at different porosities are illustrated in different figures for sake of clarity. We notice from these two figures that the microstructure has a great effect on the flow cannot define properly their permeability. In Figure 6.6 the numerically computed K and K at each porosity level are plotted against σ/R. The effect of σ/R on K /R 2 and K /R 2 is evident. An increasing σ or also moving from a uniform array to an array showing higher degrees of disorder results in an increase of both permeability values, the transverse and the longitudinal. The results of K are in agreement with those of Chen and Papathanasiou [START_REF] Chen | Micro-scale modeling of axial flow through unidirectional disordered fiber arrays[END_REF] who also found that an increase of the microstructure's heterogeneity enhances K . This is due to the formation of large flow paths. On the other hand, our results concerning K are opposite to those of Chen and Papathanasiou [START_REF] Chen | The transverse permeability of disordered fiber arrays: a statistical correlation in terms of the mean nearest interfiber spacing[END_REF] except for φ > 0.7. They explain that the decrease in K in the range 0.45 < φ < 0.7 comes from the presence of narrow gaps which reduce the permeability and that the increase of K in the range φ > 0.7 comes from the formation of flow paths whose size is comparable to the size of fiber aggregates. In our case, we explain this increase by the simple fact that large gaps enhance K as they enhance K .

Conclusion

A method to generate unidirectional disordered fiber arrays, with a specific predefined porosity and inter-fiber spacing, in an FE context was introduced. This method is based on an advancing front method to generate the fibers' centers, on level-set functions to immerse the fibers in the FE mesh and on anisotropic mesh adaptation to properly describe the fibers' interfaces. Furthermore, permeability computations of the generated unidirectional disordered fiber arrays were performed. The flow through the fibers, which were considered as rigid discs, is governed by the Stokes equations that, when averaged, lead to Darcy's law. Hence, the permeability of a fiber array can be computed. The RVE's size, or number of fibers, was determined by studying its influence on the transverse and longitudinal permeabilities. Afterwards, we have shown that the fiber arrays are isotropic and that the transverse permeability is independent of the flow direction. Furthermore, the influence of the degree of disorder of fiber arrays on transverse and longitudinal permeabilities was studied and the results have shown that both permeabilities increase on all porosity levels as the disorder increases.

Fig. 3 .

 3 Fig. 3.1. (a) η(α) defined on the whole computational domain using a mixture law and (b) a the anisotropic remeshing at the solid-fluid interface

Fig. 3 . 2 .

 32 Fig. 3.2. The MINI-element P + 1 -P1

Fig. 4 . 1 .

 41 Fig. 4.1. Three equally spaced discs of radius R and their equilateral disposition

Fig. 4 . 2 .

 42 Fig. 4.2. Fiber distributions generated by our method in a 10 × 10 square: (a) 6130 fibers with φ = 0.51 and σ = 0.06R,(b) 6111 fibers with φ = 0.51 and σ = 0.12R, (c) 3678 fibers with φ = 0.7 and σ = 0.03R and (d) 3634 fibers with φ = 0.7 and σ = 0.4R

end for r 2 .Fig. 4 . 3 .

 243 Fig. 4.3. The L(r) functions of the fiber distributions described in Figure 4.2

Fig. 6 . 1 .Fig. 6 . 2 .

 6162 Fig. 6.1. Elementary volumes extracted from a 10 × 10 microstructure in order to determine the size of the RVE

Fig. 6 . 3 . 9 Figure 6 . 6 . 3 .Figure 6 . 4

 63966364 Fig. 6.3. The averaged permeability components at different values of σ: (a) φ = 0.51, (b) φ = 0.7 and (c) φ = 0.9

  Figure 6.4 shows two representative flow paths computed on RVEs at φ = 0.7 and Figure 6.5 illustrates the velocity computed for the boundary conditions of type 3 (Eq.(4.3)) on the same previous RVEs. The difference between the two microstructures lies in the choice of σ (σ = 0.03R in Figures 6.4(a) and 6.5(a) and σ = 0.4R in Figures 6.4(b) and 6.5(b)).We notice from these two figures that the microstructure has a great effect on the flow

Fig. 6 . 4 .Fig. 6 . 5 .

 6465 Fig. 6.4. Flow paths for boundary conditions 2 across a unidirectional disordered fiber array: (a) φ = 0.7 , σ = 0.03R and (b) φ = 0.7 , σ = 0.4R

Fig. 6 . 6 .

 66 Fig. 6.6. The effect of averaged standard deviation of the inter-fiber spacing on the normalized averaged transverse (a) and longitudinal (b) permeability