Emmanuel Moulay
email: emmanuel.moulay@irccyn.ec-nantes.fr

Marc Baguelin

Meta-dynamical adaptive systems and their applications to a fractal algorithm and a biological model

Keywords: dynamical systems, adaptive systems, biological systems, fractal algorithm. MSC: 93A05, 37F50, 92D15, 92D25

In this article, one defines two models of adaptive systems: the metadynamical adaptive system using the notion of Kalman dynamical systems and the adaptive differential equations using the notion of variable dimension spaces. This concept of variable dimension spaces relates the notion of spaces to the notion of dimensions. First, a computational model of the Douady's Rabbit fractal is obtained by using the meta-dynamical adaptive system concept. Then, we focus on a defense-attack biological model described by our two formalisms.

Introduction

In the two last decades, there has been much interest in the study and formalization of complex adaptive systems (see [START_REF] Waldrop | Complexity: The Emerging Science at the Edge of Chaos[END_REF] and [START_REF] Gell-Mann | The Quark and the Jaguar: Adventures in the Simple and the Complex[END_REF]). Many different approaches have been proposed: artificial chemistries, evolutionary formalism, cellular automata. However, from a more theoretical viewpoint, few mathematical formalisms exist for adaptive systems. Though we may cite the chapter "Categorical System Theory" proposed by A.H. Louie in [START_REF] Rosen | Theoretical Biology and Complexity: Three Essays on the Natural Philosophy of Complex Systems[END_REF] who discuss the relationship between natural and formal systems, most of the studies are mainly simulation-based (see [START_REF] Gell-Mann | The Quark and the Jaguar: Adventures in the Simple and the Complex[END_REF] and [START_REF] Holland | Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control, and Artificial Intelligence[END_REF]). With the recent evolution of physics and biology, a general mathematical formalism for adaptive systems would be very useful. Actually, behind their apparent heterogeneity, most of the adaptive systems share one important feature: they are dynamical objects whose structures are sometimes modified by a top level automation-like rule. On the basis of this observation, we have built a two-level formalism that helps us to design an algorithm for fractals and a biological model of co-evolution in a bacterium-phage system.

The main goal of this article is to present a concept of adaptive systems, general enough to be used in different fields (mathematics, physics or biology) and to apply this framework to a fractal implementation and a biological system. Besides, a new concept of spaces is developed. Since the second part of the 19 th century, various definitions of dimensions, in particular those developed by Cantor and Peano, have appeared (see [START_REF] Mandelbrot | Des Monstres de Cantor et de Peano À la Géométrie Fractale de la Nature[END_REF]). These new definitions lead to the concept of a fractal model and with it, a global view of spaces (see [START_REF] Mandelbrot | Les Objets Fractals[END_REF]). In order to study adaptive systems, we will use a sort of "adaptive" space called variable dimension space whose dimension changes. This new kind of space brings the notion of space and the notion of dimension together. The paper is organized as follows. In section 2, a very general adaptive system concept called meta-dynamical adaptive system is given by using an extension of Kalman systems. Section 3 is dedicated to the Douady's Rabbit fractal implemented as a meta-dynamical adaptive system. In section 4, the concept of variable dimension space leads to a special model of adaptive system: the adaptive differential equation. Finally, we use our concepts of adaptive systems to describe a biological example of a defense-attack model in section 5.

Meta-dynamical adaptive system

It is possible to give a general definition of adaptive systems by using the concept of meta-dynamical adaptive system developed by one of the author during his Ph.D. ( [START_REF] Baguelin | Modélisation de Systèmes Complexes Par Métadynamiques: Application À la Modélisation de Populations En Coévolution[END_REF][START_REF] Baguelin | A formalism for models with a metadynamically varying structure[END_REF]). But, in order to be able to present our formal approach, we will have a look at the main ideas behind the formalization of the concept of dynamical system by Kalman in [START_REF] Kalman | Topics in Mathematical System Theory[END_REF]. The aim of Kalman's approach is to show that some very common mathematical structures plus a few axioms can provide a very general framework where the notion of dynamical systems (of all kinds) is captured. Time is modelled as an ordered subset of the reals (to cover both the continuous and discrete paradigms). The important objects are the state set (the variables characterizing the system) and the transition function. The transition function defines the trajectory in the state set starting from an initial state. Only a few axioms are required to characterize these objects and allow them to form a "dynamical system". Among the more important ones are direction of time, consistency, composition property and causality. Let us recall the fundamental definition of Kalman (see [START_REF] Kalman | Topics in Mathematical System Theory[END_REF]). A dynamical system respecting Kalman axioms is defined by the following axioms:

1. There is a given time set T , a state set X, a set of input values U , a set of acceptable input functions Ω = {ω : T → U }, a set of output values Y , and a set of output functions Γ = {γ : T → Y },

(Direction of time)

T is an ordered subset of the reals1 , 3. The input space Ω satisfies the following conditions:

(a) (Nontriviality) Ω is nonempty, (b) (Concatenation of inputs) An input segment ω ]t1,t2] is ω ∈ Ω re- stricted to ]t 1 , t 2 ] ∩ T . If ω, ω ′ ∈ Ω and t 1 < t 2 < t 3 , there is an ω ′′ ∈ Ω such that ω ′′ ]t1,t2] = ω ]t1,t2] and ω ′′ ]t2,t3] = ω ′ ]t2,t3] ,
4. There is given a state-transition function

ϕ : T × T × X × Ω → X.
whose value is the state x (t) = ϕ (t, τ, x τ , ω) ∈ X resulting at time t ∈ T from the initial state x τ = x (τ ) ∈ X at initial time τ ∈ T under the action of the input ω ∈ Ω. ϕ has the following properties:

(a) (Direction of time) ϕ is defined for all t ≥ τ but not necessarily for all t < τ .

(b) (Consistency) ϕ (t, t, x, ω) = x for all t ∈ T , all x ∈ X and all ω ∈ Ω.

(c) (Composition property) For any t 1 < t 2 < t 3 , we have

ϕ (t 3 , t 1 , x, ω) = ϕ (t 3 , t 2 , ϕ (t 2 , t 1 , x, ω) , ω)
for all x ∈ X and all ω ∈ Ω.

(d) (Causality) If ω, ω ′ ∈ Ω and ω ]τ,t] = ω ′ ]τ,t] then ϕ (t, τ, x, ω) = ϕ (t, τ, x, ω ′ ) .

5.

There is given a readout map η : T × X → Y which defines the output

y (t) = η (t, x (t)). The map ]τ, t] → Y given by σ → η (σ, ϕ (σ, τ, x, ω)), σ ∈ ]τ, t] is an output segment, that is, the restriction γ ]τ,t] of some γ ∈ Γ to ]τ, t].
A dynamical system is referred as Σ = {T, X, U, Ω, Y, Γ, ϕ, η}. Now, inspired by the definition of dynamical systems of Kalman, we propose a formalization of our meta-dynamical adaptive system. Definition 1 A meta-dynamical adaptive system M is a composite mathematical concept defined by the following axioms:

1. The dynamical level: the suitably indexed set

Σ i,j = {T 0 , X i , U, Ω, Y, Γ, ϕ i,j }
is a dynamical system respecting Kalman axioms for all (i, j) ∈ I × J where {ϕ i,j } j∈J are transition functions on state set X i :

ϕ i,j : T 0 × T 0 × X i × Ω → X i .

2.

The meta-dynamical level: let X = i∈I X i be the set of all the possible states of the system and D = {ϕ i,j } (i,j)∈I×J be the set of all possible transition functions, then there exists a meta-dynamical time T 1 and a meta-dynamical rule 2 Φ :

T 1 × X × D → X × D such that: (a) (Temporal hierarchy) T 1 ⊆ T 0 , (b) (Coherence states/transitions) If Φ(t, x 1 , ϕ i1,j1 ) = (x 2 , ϕ i2,j2 ) for t ∈ T 1 , then x 1 ∈ X i1 ⇒ x 2 ∈ X i2 .
If t ∈ T 1 and Φ(t, x, ϕ) = (x, ϕ), Φ is said to be mute at (t, x, ϕ) else (t, x, ϕ) is a commutation point.

3.

Evolution rule between dynamical and meta-dynamical levels: let x t = ϕ(t, t 1 , x 1 , ω) with ω an input function, there exists a "meta"transition function

Ψ : T 1 × T 1 × X × D × Ω → X × D such that: (a) (Dynamical phase) If Φ is mute on (t,
Ψ(t 2 , t 1 , x 1 , ϕ, ω) = (x t2 , ϕ),
ii. else t 2 ∈ T 1 and

Ψ(t 2 , t 1 , x 1 , ϕ, ω) = Φ(x t2 , ϕ).
2 It is possible to consider input values at the metadynamical level, then Φ is defined as follows Φ : The meta-dynamical rule in point 2) can operate at instants for which the system is defined and not necessarily at all of them (see point 2.(a)); actually a higher level is usually slower. Moreover, as the meta-dynamical rule can change the state and the transition function, we have to consider that both match well: the resulting state has to belong to the state set on which the new transition function operates (point 2.(b)). In point 3. we describe how dynamics and the meta-dynamical rule combine together to make the system change with time.

T 1 × X × D × Ω 1 → X × D where Ω 1 = {ω : T 1 → U 1 }
To take up notions mainly used in social science, our dynamical rule can be seen as a heterarchical level and our meta-dynamical level as a hierarchical level. Let us recall that a heterarchy is a network of elements which share the same "horizontal" position level in a hierarchy. Each level in a hierarchical system is composed of a heterarchy which contains its constituent elements (see [START_REF] Jen | Stable or robust? What's the difference?[END_REF]).

Remark 2 In the second case of (3.a), when t 2 ∈ T 1 , we deliberately consider that (t 2 , Φ (x t2 , ϕ)) is mute. We do not consider the case where the meta-dynamics would "rebound" and have several commutations at the same time. If the system has several commutations, it is always possible to consider this set of commutation as one, with the final state of the last commutation (provided that, we know from 3.(c) that the number of commutations is finite). Axiom 3.(c) is to avoid Zeno-style system with an infinite number of commutations in a finite amount of time (e.g. with a quantity x(t) = sin 2π 1-t with a commutation each time x(t) = 0 on [1 -ǫ, 1 + ǫ]).

Remark 3

The index set J is associated with the change of dynamics in the same state space. We recognize here the framework of hybrid systems (see [START_REF] Zeigler | Theory of Modelling and Simulation: Integrating Discrete Event and Continuous Complex Dynamic Systems[END_REF]).

Remark 4 It is also important to differentiate between the continuous or discrete dynamics and the continuous or discrete meta-dynamics because the confusion is easy. The first case is well known. The second one shows the difference between a meta-dynamical time T 1 which is continuous (see section 3) and a meta-dynamical time which is discrete (see section 5).

Example 5 On Fig. 1, we can see a meta-dynamical system in action. At t 5 = t c , we pass from a 2-dimensional state set to a 3-dimensional one. The dynamics are continuous, so T 0 = [t 0 , t end ] and the meta-dynamical time set T 1 is discrete:

T 1 = {t 1 , t 2 , t 3 , . . . , t 8 }. Since the meta-dynamics is mute on t 1 , t 2 , instant from T 0 ϕ t t t t X end 1 6 7 8 1 ϕ t t t t t t X = t c point commutation 1 0 2 4 5 0 3 0 Figure 1:
Illustration of some axioms of a metadynamical system t 3 and t 4 , the evolution from t 0 to t 5 is a purely dynamical phase, ended by a commutation (Axiom 3.(a)). Evolution from t 5 to t end is also purely dynamical (the meta-dynamics is mute on t 6 , t 7 and t 8 . The junction between the two of them is made by using the concatenation rule. As there is only one commutation, the system is defined. Now, it is possible to expand some classical notions of dynamical system, as the notion of trajectory.

Definition 6 Let M be a meta-dynamical adaptive system with the meta-transition function Ψ, if for t ∈ T 1 ,

Ψ(t, t 0 , x 0 , ϕ 0 , ω) = Φ(x t , ϕ t ) with x t ∈ X t then (X t , ϕ t ) is a meta-state of M in t.
The sequence (X ti , ϕ ti ) 0≤i≤n,0≤n≤+∞ of meta-state of M such that for all t ∈ [t i , t i+1 [, (X t , ϕ t ) = (X ti , ϕ ti ) and (X ti , ϕ ti ) = X ti+1 , ϕ ti+1 is called meta-orbit of M.

If n < +∞, (X tn , ϕ tn ) is said to be an absorbent meta-state.

The sequence (t i , X ti , ϕ ti ) 0≤i≤n,0≤n≤+∞ satisfying the same conditions is called a trajectory of M.

One can extend definition 1 to stochastic systems.

Definition 7 Let us consider the dynamical rule of definition 1 with the following meta-dynamical rule

pΦ : T 1 × X × D → (X × D → [0, 1])
where pΦ is a probability distribution on X × D which represents the probability

pΦ (t, x, ϕ) • (x f , ϕ f ) that (x, ϕ) becomes (x f , ϕ f ) at t.
The stochastic evolution rule between dynamical and meta-dynamical levels is given by

pΨ : T 1 × T 1 × X × D × Ω → (X × D → [0, 1])
where pΨ is a probability distribution on X × D which represents the probability

pΨ (t 1 , t 2 , x, ϕ, ω) • (x f , ϕ f )
that (x, ϕ) becomes (x f , ϕ f ) between t 1 and t 2 with the input function ω. The properties of pΦ and pΨ are the same as Φ and Ψ given in 1. Such a system is called a stochastic meta-dynamical adaptive system.

The interest of definition 7 is its generality. It is for example well adapted to model a large number of complex adaptive systems, in particular the biological model we develop in section 5.

Algorithm for the Douady's Rabbit fractal

The goal of this example is to give an algorithm based on definition 1 allowing to describe the Douady's Rabbit Fractal. In this example, the set I and thus the family {X i } i∈I of definition 1 is not countable. The Riemann Sphere S 2 is mapped one-to-one onto the extended complex plane C ∞ = C∪ {∞} by stereographic projection. Let us recall some definitions. Let K = R or C and B be the unit ball of K.

If f : K → K is a function, x ∈ X is periodic of period n ∈ N if f n (x) = x and for all k ∈ {1, .., n -1}, f k (x) = x. For x periodic of period n, the cycle O(x) is the set O(x) = x, f (x), ..., f n-1 (x) and its cardinal is n. Moreover, if f is differentiable at x, x is stable, quasi stable or unstable if (f n ) ′ (x) < 1, (f n ) ′ (x) = 1 or (f n ) ′ (x) > 1.
If x is stable, x is attractive if there exists an interval V strictly containing x so that for all

x ∈ V f n (x ′ ) → n→+∞ x.
Now, let us recall the definition of the Julia set (see [START_REF] Alexander | A History of Complex Dynamics: From Schrüder to Fatou and Julia[END_REF]). Let R be a non-constant rational function on S 2 . The Fatou set of R is the maximal open subset of S 2 on which {R n } n∈N is equicontinuous where

R n = R • . . . • R. The Julia set J R of R
is the complement of the Fatou set on S 2 . The filled in Julia set K R of a function R is all the points which are not attracted to the super-attracting fixed point at infinity, that is

K R = {z ∈ C : R n (z) → ∞} .
This closed set includes the Julia set as its boundary, J R = ∂K R . The escape set I R of a function R is all the points that "escape" to infinity, that is

I R = {z ∈ C ∞ : R n (z) → ∞} .
If R is a polynomial of degree 2, J R is called a quadratic Julia set. The following result can be found in [START_REF] Brolin | Invariant sets under iteration of rational functions[END_REF]: if P is a polynomial of degree d ≥ 2, then J P is closed and dense within itself.

Here, we are interested in quadratic Julia sets with

P c (z) = z 2 + c (1)
where c ∈ C. For small values of c, the Julia set is distorted by varying degrees from the unit circle, in these cases the Julia set has an infinite length. For large values of c, the Julia set becomes an infinite set of totally disconnected points, often said to be dust like (in the sense of Cantor). In the quadratic case there are only these two possible -connected and disconnected-types. The Mandelbrot set is the space containing the value c for which the associated Julia set is connected. It is generated by the quadratic sequence

z n+1 = z 2 n + c z 0 = c (2)
A complex point z = a + ib ∈ C will be denoted by (a; b). The Douady's Rabbit fractal is a Mandelbrot set with c = (- , so c ≃ (-0.12256; 0.744862) (see [START_REF] Douady | Etude dynamique des polynômes complexes, (première et deuxième partie)[END_REF]). As the points get closer to the Douady's Rabbit fractal, the speed of convergence becomes slower. With the choice of c, c 2 + c 2 + c = 0 so the origin is an

3 2 + 1 2 (a + b); √ 3
attractive cycle 3 of period 3 of P c O(0) = 0, c, c 2 + c .
The boundary points move chaotically. Thus, the idea is to change the points by using the speed of convergence as an adaptive value. The principle of our algorithm is to have a set of points evolved to the fractal boundary. For this, one gives a weight to each point. This weight varies with a dynamics which "rewards" the most adapted points (the points which are less attracted by limit values) and "penalizes" the least interesting ones. When the most efficient points reach a certain threshold, they are allowed to be multiplied in their neighborhood. The least efficient points disappear when they reach a minimal threshold. This system of thresholds which changes the dynamical structure is our meta-dynamics. Dynamical level (DL): with the formalism of definition 1, we have T 0 = N. The index set I is the set of all finite sets of points of C (so i is a set of points

3 Let K = R or C, f : K → K be function, x ∈ X is periodic of period n ∈ N if f n (x) = x and for all k ∈ {1, .., n -1}, f k (x) = x. For x periodic of period n, the cycle O(x) is the set O(x) = x, f (x), ..., f n-1 (x) and its cardinal is n. Moreover, if f is differentiable at x, x is stable, quasi stable or unstable if (f n ) ′ (x) < 1, (f n ) ′ (x) = 1 or (f n ) ′ (x) > 1. If x is stable, x is attractive if there exists an interval V strictly containing x so that for all x ∈ V f n (x ′ ) → n→+∞ x.
of C and I is not countable). For all i ∈ I, X i = N card(i) is the set of the point weights of i. At each given point z ∈ C, we attach a selective value

µ n : C → R + z → min |P n c (z)| , |P n c (z) -c| , P n c (z) -c 2 -c , 1 |P n c (z)| 0, c, c 2 + c is the attractive cycle of period 3, 1 P n c (z)
is referring to the attraction to the infinity seeing that

P n c (z) → n→+∞ ∞ is equivalent to 1 P n c (z) → n→+∞ 0.
One increases the weight of the points which are close to the boundary. In order to do this, we organize a "competition" by comparing their mutual slowness of convergence. So, the weights of the points z ∈ i are given by

ω z (t + 1) = ω z (t) + q∈i,q =z δ (z, q) with δ (z, q) =    1 if µ t (z) > µ t (q) 0 if µ t (z) = µ t (q) -1 if µ t (z) < µ t (q) ω z (0) = 0
where t ∈ N. Using the notation of section 2, we have the following transition function

ϕ i t + 1, t, {ω z (t)} z∈i = {ω z (t + 1)} z∈i
for all t ∈ N. For each state space there is only one associated transition function, so the use of the index set J is unnecessary. Meta-dynamical level (ML): with the formalism of definition 1, we have T 1 = kN. When the weight of a point reach an upper threshold M > 0, the point is allowed to give birth to a new point, randomly in its neighborhood. When the weight of a point reaches a lower threshold m < 0, the point is removed. This can be modelled by the meta-dynamical rule

Φ t, {ω z (t)} z∈i(n) , ϕ i(n) = {ω z ′ (t)} z ′ ∈i(n+1) , ϕ i(n+1) where i (n + 1) = {z + ̺ǫ : ̺ is a random point in B, z ∈ i (n) and ω z (t) > M } ∪ {z ∈ i (n) : ω z (t) ≥ m}
such that ǫ > 0 given. The weight of the reproduced point is shared between itself ω z (t) and the new neighbouring point ω z+̺ǫ (t). The other weights are kept equal. Let us sum up the evolution rule:

. . . → ϕ i(n) (t) → DL ϕ i(n) (t + 1) → DL . . . → DL ϕ i(n) (t + k) → ML ϕ i(n+1) (t + k) → . . .
As the time increases, the selective value becomes more and more accurate, i.e. i (n) tends to a set of points belonging to the Julia space J Pc or to the empty set when n → +∞. This is an interesting point because the calculations are concentrated on the fractal boundary. In a lot of classical algorithms, the K Pc (see [START_REF] Peitgen | The Science of Fractal Images[END_REF]). Indeed, this meta-dynamical adaptive system produces a cloud of points which gather round the Julia space J Pc called the Douady's Rabbit fractal (see Fig. 2).

Though we have developed this algorithm in the setting of Julia sets, the same framework can be used to explore many complex frontiers, for example other fractal structures where the boundary properties are largely unknown. Indeed, this kind of algorithm is really interesting because the calculation is concentrated on a specific region.

Adaptive differential equations

First, let us recall the classical notion of Hausdorff fractal dimension (see [START_REF] Peyrière | Calculs de dimension de hausdorff[END_REF]). Let E be a subspace of a metric space M and ρ a positive number, one defines R ρ as the set of all coverings (B i , ρ i ) i of E by balls B i with diameter 0 < ρ i < ρ. For each positif number α, one denotes:

H α ρ (E) = inf i ρ α i : (B i , ρ i ) i ∈ R ρ . H α (E) = lim ρ→0 H α ρ (E) is called the α-dimensional Hausdorff measure of E and belongs to [0, +∞]. Then, let dim(E) = inf {α > 0 : H α (E) = 0} ,
it is the Hausdorff dimension of E. On the one hand, the Hausdorff dimension is defined for all metric spaces. On the other hand, in the case of a classical space (non fractal), it is identical to its topological dimension (for example the Hausdorff dimension of R n is n). In the case of a simple linear fractal, such as fractals with internal homothetia obtained by an homothetic iteration with constant ratio, the Hausdorff dimension is equal to the homothetic dimension dim h (E) given by:

dim h (E) = ln(n) ln(k) = log k (n)
where n is the number of subsets obtained during the homothetic process of reduction with ratio 1 k (see [START_REF] Besicovitch | Sets of fractional dimensions (v): On dimensional numbers of some continuous curves[END_REF]). For more information on dimension theory the reader may refer to [START_REF] Tricot | Douze définitions de la densité logarithmique[END_REF] or [START_REF] Hurewicz | Dimension Theory[END_REF].

The variable dimension space is defined as follows:

Definition To illustrate these definitions, we can consider the following example:

Example 9 Let 1 4 ≤ λ ≤ 1 2
, one defines four similarities from the family K c of the compact subset of the square c = [0, 1]

2 with value in K c by

• s 1,λ (x) = λx, • s 2,λ (x) = λ   1 2 -λ -λ -1 4 λ -1 4 1 2 -λ   x + λ 0 , • s 3,λ (x) = λ   1 2 -λ λ -1 4 -λ -1 4 1 2 -λ   x + 1 2 λ -1 4 and • s 4,λ (x) = λx + 1 -λ 0 .
Then, one defines the function

Ω λ : K c → K c , x → s 1,λ (x) ∪ s 2,λ (x) ∪ s 3,λ (x) ∪ s 4,λ (x).
Ω λ being contracting, one defines F (λ) as the fixed point4 of this function for the Hausdorff distance. In [START_REF] Mandelbrot | Les Objets Fractals[END_REF], one may find the following result: the dimension

1/2, λ-1/4 0, 0 λ, 0 1-λ, 0 1, 0 Figure 3: Van Koch snowflake generator of F (λ) is d (λ) = ln4 ln( 1 λ )
. So, (F, d) is a continuous increasing variable dimension space on Λ = 1 4 , 1 2 . A similar construction allows us to turn continuously a n-dimension space into a (n + 1)-dimension space.

An interesting application of the above concept arises when the variable dimension space is the state space of the solutions of a differential equation over a period of time. The union of the solutions then "moves" on a variable dimension space. We talk about an adaptive differential equation.

Definition 10 Let X be a Banach space of finite dimension called the possible state space. Consider

1. a subdivision {t i } i∈N of R + , 2. an application d : R + × X → N, (t, y) → d (t, y) such that d (0, y) is given, 3. an application g : R + × R dim(X) → R dim(X) ,
such that one has the system

ẋi (t) = f i (t, x i (t)), t ∈ [t i , t i+1 [ , x i (t) ∈ R d(ti,yi) , i ∈ N x i (t i ) = g (t i , y i ) ( 3
)
where y 0 = x 0 (0), y i = lim

t→t - i x i-1 (t) for i ∈ N * , ẋi (.) is the right derivative of x i (.) and {f i } i∈N : [t i , t i+1 [ × R d(ti,yi) → R d(ti,yi)
is a family of applications. Such a system is called an adaptive differential equation. A trajectory of the system (3) is a family

x (t) = {x i (t) : t ∈ [t i , t i+1 [} i∈N .
One may notice that dim(g (t i , y i )) = d (t i , y i ). This concept of adaptive differential equation is not a succession of ordinary differential equations because the initial condition of each system i depends on the system i -1. The ordinary differential equations represent a case where d is constant. The adaptive differential equation is a reductive approach of the meta-dynamical adaptive system to a system described by a differential equation whose metadynamics is discrete and where dim X < +∞.

Application to a biological model

Here, we want to model the influence of the bacterium-phage interaction on the co-evolution of the populations of bacteria and phages. In the following we will need the definition of the Hamming distance. Let n ∈ N, the Hamming distance is the function d H : {0, 1} n → N defined by

d H (s 1 , s 2 ) = nc k=0 s k i -s k j
where s 1 , s 2 ∈ {0, 1} n . d H represents the number of differing bits between the two binary strings

s 1 = s 1 1 s 2 1 . . . s n 1 and s 2 = s 1 2 s 2 2 . . . s n 2 .
The attack of a bacteria population by phages is assumed to be done by the lysis process: a phage hangs on the surface of a bacterium cell, injects its DNA in it and then forces the bacterium to yield its own replicas inside the cell. When the cell is full, it bursts, releasing a huge quantity of copies of the infecting phage. The efficiency of the attack (i.e. the probability of success of the infection), depends on the couple bacterium-phage. One of the other characteristics of bacterial and phages populations are their high variabilities. They frequently mutate, creating new populations with new properties. Such a system has two dynamics to be taken into account: the dynamics of the populations of bacteria and phages and the meta-dynamics of evolution geared by mutations and extinctions. It is without any doubt a model which is not in the scope of the classical theory of dynamical systems, but our formalism applies well to it.

Dynamical level: it is made up of a set of ordinary differential equations. Let us consider the following system {S, {B i } 0≤i≤n b , {P j , I 1,j , I 2,j , I 3,j } 0≤j≤np },

with S the concentration of nutrient, B i the different bacteria strains with 0 ≤ i ≤ n b and P j and I k,j the different phages strains with 0 ≤ j ≤ n p . To each population B i and P j is associated a binary string s b i ∈ {0, 1} n and s p i ∈ {0, 1} n . These binary strings code for the properties of attack (for the phages) or defence (for the bacteria) facing the infection (see [START_REF] Baguelin | Modélisation de Systèmes Complexes Par Métadynamiques: Application À la Modélisation de Populations En Coévolution[END_REF] for the biological discussion). The populations are then characterized by the different concentrations and the two lists of bit strings. The set L of all the pairs of lists is taken as our second index (the one called J in definition 1).

The model is described by a modified version of Mosekilde equations (see [START_REF] Mosekilde | Topics in Nonlinear Dynamics[END_REF]). This set of ordinary differential equations describes the interactions of bacterial populations B i and phage populations P j in a chemostat. B i and P j also symbolize the concentration of these populations and S the concentration of the nutrient. The process of infection of bacteria by phages is modelled by three infection stages I 1,j , I 2,j and I 3,j . One associates with the ecosystem (4) the state space X 1+n b +4np = R 1+n b +4np + . Over a period of time without appearance or disappearance of any strain of bacteria or phages, the dynamical evolution of the system (4) is modelled by the evolution of the concentrations of S, the different bacteria B i and phages P j , and three infection stages I k,j . With the formalism of definition 1, that means that T 0 = R and for a fixed l ∈ L the transition function ϕ 1+n b +4np,l is the integration of the following set of differential equations

                                   dBi dt = νSBi κ+S -B i np j=1 αω ij P j -ρB i dI1,j dt = P j n b i=1 αω ij B i -3 I1,j τ -ρI 1,j dI2,j dt = 3 τ (I 1,j -I 2,j ) -ρI 2,j dI3,j dt = 3 τ (I 2,j -I 3,j ) -ρI 3,j dPj dt = 3 βI3,j τ -P j n b i=1 αB i + np j=1 3 k=1 αI k,j -ρP j dS dt = ρ(σ -S) - n b i=1 νγSBi κ+S (5)
with 0 ≤ i ≤ n b , 0 ≤ j ≤ n p , ρ the rate of dilution (ρ = 0.0045 min -1 ), κand ν respectively the saturation term and the growth from the Monod equation (κ = 10µg.ml -1 , ν = 0.024 min -1 ), α the theoretical adsorption constant depending on phage and bacterium size (α = 10 -9 ml. min -1 ), τ a time constant (τ = 30 min), β the number of copies of phage j released during the burst of the infected bacterial cell (β = 100), σ the continuous supply of substrate (σ = 10µg.ml -1 ), γ the amount of nutrient consumed in each cellular division (γ = 0.01ng) and finally ω ij the probability of infection of B i by P j which depends on the similarity between bit string s b i (attached to bacterium population B i ) and s p j (attached to phage population P j ) as follows

ω ij = 1 - d H (s b i , s p j ) n c 2
with n c the size of the binary string. Here, we are not interested in the identification of the biological dynamical level which can be found in [START_REF] Baguelin | Modélisation de Systèmes Complexes Par Métadynamiques: Application À la Modélisation de Populations En Coévolution[END_REF]. With the given size of a binary string n c , there exists a finite number of possible dynamical systems, here differential equations. Indeed, the number of possible different populations of bacteria is equal to the number of parts of the set of binary strings of size n c . So, there are 2 2 nc possible populations of bacteria. For the same reason, one deduces that the possible number of populations of phages is the same and thus the total number of possible state spaces is 2 2 nc × 2 2 nc = 2 2 nc +1 . So, with the notation of definition 1, it means that card (J) = 2 2 nc +1 . If we take for example n c = 10, one has 2 1025 ≃ 3.6×10 308 possible state spaces. Theoretically, one can consider a 3.6 × 10 308 dimensional space to embed the system [START_REF] Besicovitch | Sets of fractional dimensions (v): On dimensional numbers of some continuous curves[END_REF]. Nevertheless, when it comes to numerical simulations, such a big system is impossible to deal with.

One thus understands the need of an adaptive system to describe the system [START_REF] Besicovitch | Sets of fractional dimensions (v): On dimensional numbers of some continuous curves[END_REF].

Meta-dynamical level (meta-dynamical adaptive system point of view): this is the main difference with the model of Mosekilde given in [15, Chapter A] which is not evolutive.

Proposition 11 Consider a small interval of time ∆t, the adaptive changes of the system (4) are given by the following mechanisms

pΨ(t + ∆t, t, δ, ϕ 2+n b +4np,l1 ) • θ 1 , . . . , θ k , ϕ 1+k δ +n b +4np,l2 = e -λ(t) λ (t) k k! (6)
and

Ψ(t + ∆t, t, δ, ϕ 2+n b +4np,l1 , m δ ) = ϕ 2-ǫ δ +n b +4np,l2 (7)
where 1. formula ( 6) is the probability that the species δ (a bacterium B or a phage P ) gives birth to k ≥ 1 mutante strains θ 1 , . . . , θ k on [t, t + ∆t[ with λ (t) = δ(t) δe p e , p e the probability that a small group of mutant species of size δ e gradually replacing the species of the parent population δ. ϕ i,j is the function defined by equations ( 5) for j = l 1 , l 2 and k δ = k if δ is a bacterium and k δ = 4k if δ is a phage, 2. formula [START_REF] Douady | Etude dynamique des polynômes complexes, (première et deuxième partie)[END_REF] gives the determinist rule of the extinction of the species δ which depends on a given threshold m δ , ǫ δ = 1 if δ is a bacterium and

ǫ δ = 4 if δ is a phage.
pΨ is the stochastic transition function of definition 7 which governs the meta-dynamical rule of the appearance. Ψ is the transition function of definition 1 which governs the meta-dynamical rule of the extinction.

Proof. Consider the probability p Bi (t, k) that the population B i gives birth to k mutante strain on [t, t + ∆t[ and the probability p e that a small group of mutante bacteria of size B e is gradually replacing the bacteria of the parent population B i . Such a reasoning gives a binomial probability for

p Bi (t, k) = C k n(t) p k e (1 -p e ) n(t)-k
with n (t) = Bi(t) Be . When n (t) is large, one may approximate the binomial probability by the Poisson probability

p Bi (t, k) ≃ e -λ(t) λ (t) k k! (8)
with λ (t) = Bi(t) Be p e . Suppose that the birth of all the populations to a Hamming distance of one (only one bit is different) is equiprobable. The birth of a mutant strain results in the change of a group of bacteria (which is a part of a parent population) of size B e , in a population group with new characteristics. If this population group already exists, the mutant population is added to it. The same mechanism governs the modelling of the phage mutation. For the extension, there exists a threshold m δ below which the extinction of the population is certain. Every population under a given threshold (different for bacteria m B and phages m P ) is removed from the system. Thus,the extinction is determinist.

We have defined the macroscopic birth of a mutant population as an event occurring on an interval [t, t + ∆t[. There exists a time set T 1 where the system [START_REF] Besicovitch | Sets of fractional dimensions (v): On dimensional numbers of some continuous curves[END_REF] commutes. This commutation depends on the state of the system (5) and on the coefficient p B et p P respectively coefficient of the Poisson law of the bacterium and the phage. In order to simplify the model and to make it computable, one may suppose that T 1 = t 0 + i∆t even if in practice, the mutations have no reasons to be periodically defined. This allows to define, with the previous notations, the meta-dynamical rule

pΦ(t 0 + i∆t, δ, ϕ 2+n b +4np,l1 ) • θ 1 , . . . , θ k , ϕ 1+k δ +n b +4np,l2 = e -λ(t) λ (t) k k! Φ(t 0 + i∆t, δ, ϕ 2+n b +4np,1 , m δ ) = ϕ 2-ǫ δ +n b +4np,2 .
Then, at each step ∆t there are four possible commutations:

• birth of a new bacterial strain: a variable B n b +1 is added, the dimension of the system (5) increases by one,

• birth of a new phagical strain: four variables are added: P np+1 , I 1,np+1 , I 2,np+1 and I 3,np+1 , the dimension of the system (5) increases by four,

• extinction of a bacterial strain: the concerned variable is removed, the dimension of the system (5) decreases by one,

• extinction of phagical strain: variables of the concerned phage are removed, the dimension of the system (5) decreases by four.

The different possible state spaces resulting from a commutation are given on figure 4.

If more than one event occurs at each moment t ∈ T 1 , one composes the possible change of state spaces (for example, the extinction of a phage combined with the birth of a bacterial strain decreases the dimension of the system (5) by 4 -1 = 3).

Meta-dynamical level (adaptive differential equation point of view): we have described the biological model as a stochastic meta-dynamical adaptive system by using definition 7. This modelling corresponds to a stochastic view of the system (4) where the space of all stochastic realizations is infinite. There exists a family of points {t i } i∈N where the system (5) commutes. We have supposed that this commutation is given by t i = t 0 + i∆t. Then, at each step ∆t the system may commute. To see a possible evolution of the system (4), we choose a stochastic realization at each point t i (only for the appearance of a strain because the extinction of a strain is determinist) g(t i , y i ) where y i = lim t→t - i x i-1 (t) with x i-1 (t) the solution of the system (5) on [t i-1 , t i [. It describes the determinist evolution of the system (4). By extension, for us a stochastic realization is a function g : R + × R dim(X) → R dim(X) defined at least on {(t i , y i )} i∈N . For a given realization g, the system (4) may be modelled by an adaptive differential equation whose equations are given by the system [START_REF] Besicovitch | Sets of fractional dimensions (v): On dimensional numbers of some continuous curves[END_REF]. This modelling belongs to the variable dimensional space R n b +4np+1 +

where n b +4n p +1 depends on g and follows the rule given on figure 4. A detailed study of this biological model with implementation can be found in [START_REF] Baguelin | Preliminary approach to the simulation of coevolution in an ecosystem[END_REF].

In this example, one sees that the framework of the adaptive differential equations corresponds to the case of a transition function defined by a set of differential equations and a meta-dynamical time reduced to a discrete subset of R.

Conclusion

We started our research with a biological system whose dynamics changes in different dimensions. Seeing that there was no mathematical framework to describe such a system, we have developed a mathematical tool following Kalman's dynamical system called meta-dynamical adaptive system which was appropriate to give a constructive algorithm for some fractals. It was also adapted to describe and analyze our biological system. However, we have found our tool too general and we decided to develop a special tool for differential equations. The new system called "adaptive differential equation" is not a succession of differential equations because the initial condition of each system depends on the previous system and gives the new dimension of the following system. This model allows to describe a system of changing differential equations, in particular the stochastic realization of a stochastic meta-dynamical adaptive system. The last tool we use is the "variable dimension space". This new kind of space links the notion of space and dimension in a changing dynamics. We think that our work will contribute to understand the huge number of complex systems where the espace of exploration is too big to be investigate with classical means.

2

  (a -b)) where a = 3 25+ √ 621 54 and b = 3 25- √ 621 54

Figure 2 :

 2 Figure 2: The Douady's Rabbit fractal

Figure 4 :

 4 Figure 4: Evolution of the state space with n b and n p respectively the number of bacterial and phage strains before the transition.

  x t , ϕ) for all t ∈ [t 1 , t 2 [∩T 1 , then Ψ is defined between t 1 and t 2 and

	i. if t 2 / ∈ T 1 , then

  is the set of acceptable "meta"-input functions.(b) (Concatenation rule) If there exists t 2 ∈]t 1 , t 3 [ such that Ψ is defined between t 1 and t 2 and between t 2 and t 3 then Ψ is defined between t 1 and t 3 and

	Ψ(t 3 , t 1 , x 1 , ϕ, ω) = Ψ(t 3 , t 2 , Ψ(t 2 , t 1 , x 1 , ϕ, ω), ω).
	(c) (Stopping rule) Ψ is defined between t 1 and t 2 if, respecting previous
	axioms, there is only a finite number of commutation points in [t 1 , t 2 ].

  [START_REF] Gell-Mann | The Quark and the Jaguar: Adventures in the Simple and the Complex[END_REF] Let M be a metric space and Λ a parameter space. One defines two maps d : Λ → [0, +∞] , λ → d(λ) and F : Λ → 2 M where F (λ) verifies dimF (λ) = d(λ) and 2 M is the family of non-empty subsets of M.

F (Λ) is a set of variable dimension spaces. If F (Λ) is a totally ordered set for the inclusion, one calls (F, d) a variable dimension space. d is the dimension function and d(F (λ)) ≤ d(M) for all λ ∈ Λ. If Λ is a topological space then F is a set valued function and we may take account of the regularity of d. The variable dimension space (F, d) is continuous if d is continuous. Moreover, if Λ is an ordered space then (F, d) is increasing (respectively decreasing) if d is increasing (respectively decreasing).

Such a general definition could include exotic sets such as fractal cantor sets, for instance, in practice, the sets used are part of R or N.

To prove this, one uses two classical results: the family of subsets of a Banach space endowed with the Hausdorff distance is also a Banach space and if the similarities are contracting in K 2 , then the function Ω λ is also contracting in K c with the Hausdorff distance[START_REF] Tricot | Courbes et Dimension Fractale[END_REF].

Acknowledgement

The authors would like to thank Jacques LeFèvre to have given a start to our study of evolutionary ecosystems and Jean-Pierre Richard at the head of the SyNeR team for his help in producing this article.


T Dbouk

F Perales

F Babik

R Mozul

A DF-IBM/NSCD coupling framework to simulate immersed particle interactions

Keywords: Immersed granular flows, Fluid/structure interaction, Direct-forcing, Immersed boundary method, Non-smooth contact dynamics

Immersed granular flows a re p resent w idely i n d ifferent d omains u nder d ifferent f orms ( at various s cales) s uch a s i n nature (rivers, muds, atmosphere, blood...), and in many industrial applications (detergents, cosmetics, etc...). Studying such flows properly requires one to represent well the physics behind their dynamics: the fluid/solid interactions (FSI), the solid/solid interactions (SSI) and the coupling mechanisms at various scales.

In this work, a new coupling framework to simulate immersed granular flows has been developed. The FSI has been modeled using a direct-forcing immersed boundary method (DF-IBM) and implemented in the parallelized "PELICANS" C++ library. In this DF-IBM, all the mathematical equations, including the direct-forcing term, are discretized, both in space and time, and solved iteratively via a finite-volume a nd p rojection methods o n E ulerian G rids. A sharp-edge i nterface, that c an b e s moothed, i s used to represent the fluid/solid transition. The modeling of the multiple SSI at the grain's scale is based on the Non-Smooth Contact Dynamics (NSCD) approach developed in the "LMGC90" open-source library. The coupling of the two softwares "PELICANS" and "LMGC90", called Xper, provides an efficient framework to simulate and study dense immersed granular flows by taking into account, both advanced contact laws between grains, and hydrodynamic interactions.

We address in this paper the effects of imposing a fluid-ring numerically (or fluid-mesh-cells) around two settling solid disks on modifying their dynamics. The DF-IBM approach implemented in Xper is validated, on a 2D flow over a stationary rigid cylinder benchmark, and on the settling of a rigid buoyant sphere in an incompressible laminar fluid at different Reynolds numbers. The numerical results are in good agreement with experimental and numerical data from the literature.

Introduction

Sands of a sea, rocks of a river, dust particles of an atmosphere, red cells of a blood, solid particles of cosmetic and detergent products are all some examples on rigid grains or particles of a high concentration immersed in a fluid. Whenever these immersed granular particles are subject to a flow in motion, interactions arise between the particles themselves, called Solid/Solid interactions (SSI) and also between the particles and the surrounding fluid, called Fluid/Solid Interactions (FSI). Studying such flows experimentally is extremely important to understand and characterize their rheology. However, in many situations, the experimental setups may be a difficult task and can be replaced or completed by numerical investigations.

Experimentally speaking, many scientists performed a lot of measurements over the years in an attempt to understand the rheology behind such flows. For example, some researchers [START_REF] Cate | Particle imaging velocimetry experiments and lattice-Boltzmann simulations on a single sphere settling under gravity[END_REF] measured the settling velocity of a rigid sphere in a viscous fluid. Others conducted rheometry measurements for the flows of highly concentrated non-Brownian suspensions of hard spheres. N. Blaise [START_REF] Blaise | The dam break problem for a hyperconcentrated suspension[END_REF] characterized the wavefront dynamics in a dam break problem for a hyperconcentrated suspension for a wide variety of particle concentrations. Recently, F. Blanc et al. [START_REF] Blanc | Tunable fall velocity of a dense ball in oscillatory cross-sheared concentrated suspensions[END_REF] measured the fall velocity of a dense ball in an oscillatory cross-sheared concentrated suspension.

Numerically speaking, many methods have been developed successfully over the years to account for the presence of rigid bodies in a fluid as we will see in the coming section. These methods modeled immersed granular flows at different scales: microscopic, mesoscopic and macroscopic. However, the choice of a method is a major challenge and depends on many factors such as the computational time needed. The microscopic modeling approach is the most costly one with an extreme computational effort but the most complete one from a modeling point of view. In fact, there are many lacks in the mesoscopic and macroscopic approaches for predicting rheological behavior of the flow such as the highly concentrated suspensions where the microstructure plays a vital role. For these reasons, an objective for us is the numerical modeling on a microscopic scale in an attempt to constitute later macroscopic rheology laws. In this paper, for the first time to our knowledge, we show with numerical examples the effect of imposing/non-imposing a numerical fluid-ring (or fluid-mesh-cells) around a particle on the dynamics of two settling particles in a vertical channel. Most studies follow same or similar numerical strategies in modeling the contact laws between the particles, but without any numerical tests, illustration or validations on the effect for example of this imposing/non-imposing of a numerical fluid-ring around the particles on their dynamics. We show in our work that it has a great effect on the dynamics at the microscopic scale, which might thus regenerate larger kinematics modifications in the macroscopic behavior of an immersed granular flow.

The present paper is organized as follows. The numerical fluid-structure interactions modeling and the validation are presented in Section 2. Section 3 presents numerical coupling framework and its associated benchmark. The effects of the forcing/not-forcing of many methods have been developed around two settling disks in contact on the dynamics are investigated in Section 4.

The fluid-structure interactions modeling

The direct-forcing-immersed boundary method

In the present study, the Direct-Forcing Immersed Boundary Method (DF-IBM) is the numerical method used to account for the presence of rigid bodies immersed in a fluid (Fig. 1). The flow field is solved on an Eulerian grid discretization that does not need to conform to the solid boundaries. The solids are discretized on a Cartesian grid and their boundaries are represented by a set of Lagrangian points.

The original Immersed Boundary Method (IBM) goes back to the early works of Peskin [START_REF] Peskin | Flow patterns around heart valves: A numerical method[END_REF][START_REF] Peskin | Numerical analysis of blood flow in the heart[END_REF][START_REF] Peskin | The fluid dynamics of heart valves: Experimental, theoretical, and computational method[END_REF] and many others who added their finger-prints like [START_REF] Mohd-Yusof | Combined immersed-boundary/b-spline methods for simulations of flow in complex geometries[END_REF][START_REF] Fadlun | Combined immersed boundary finite difference methods for three-dimensional complex flow simulations[END_REF]. The history of the IBM is not the main task of this work and thus we leave it to many authors in literature who tackled it such as in [START_REF] Mittal | Immersed boundary methods[END_REF]. Briefly speaking, the DF-IBM is defined by modeling the fluid with the classical Navier-Stokes equations and then by imposing a velocity penalization term and a velocity time integral in the momentum equation only at mesh points that may define a material point of an immersed object. The general idea is to penalize and force-directly a no-slip boundary condition (BC) on the object's surface surrounded by the fluid. This methodology permits fast and efficient numerical approximations on Cartesian Eulerian Grids and may avoid body-fitted unstructured meshes and costly remeshing procedure.

The incompressible fluid dynamics, including an immersed rigid body (Fig. 1) of volume Ω s and boundary ∂Ω s , are governed by the transient directly-forced Navier-Stokes equations ( 1)-( 2): and

(1 -η)ρ f ∂u ∂t + (1 -η)[ρ f u • ∇u -∇ • σ -f ext ] = ρ f η (u s -u) ∆t (1)
∇ • u = 0 ( 2
)
where σ is the fluid stress tensor defined as:

σ = -pI + τ τ = µ[∇u + (∇u) T ] (3)
and f ext an external applied force field vector (such as the gravitational force), u the fluid velocity vector, u s the immersed rigid body velocity, ρ f the fluid density, p the pressure, µ the fluid dynamic viscosity, and η is a sharpedged coefficient defined as the following:

η = 1 -ϵ ∀ -→ X (x, y, z) ∈ Ω s (4) η = ϵ ∀ -→ X (x, y, z) ∈ Ω f (5)
with ϵ a parameter such as ϵ ≪ 1.

The boundary conditions and the initial conditions are respectively: u = u boundary on the boundary of Ω f and u(t 0 ) = u 0 inside Ω f at the initial time t = t 0 . Note that this formulation can be viewed as a penalization method.

In this work ϵ is taken equal to 10 -12 . However, the effect of the parameter ϵ on the accuracy of this method is minor as explained in literature by [START_REF] Fadlun | Combined immersed boundary finite difference methods for three-dimensional complex flow simulations[END_REF]. The authors have established rigorous estimates of the error induced by such penalization methods and showed numerically their efficiency. Note that if η = 1, the Eq. ( 1) reduces to u = u s which is a direct forcing of the fluid velocity field on and inside the immersed rigid boundary. Note that the shape of the function η is not the issue of this work and the form of this coefficient is left for investigation in a further work.

The Navier-Stokes equations (1)-( 2) are discretized temporally using an Explicit Euler Scheme and spatially via a Marker-And-Cell scheme [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF] on staggered grids. The numerical system is solved iteratively by a Pressure Schur Complement method (PSC) inherited from the classical Discrete Projection Methods provided by Chorin [START_REF] Chorin | Numerical solution of the Navier-Stockes equations[END_REF] and Van Kan [START_REF] Van Kan | A second-order accurate pressure-correction scheme for viscous incompressible flow[END_REF]. We will describe the method briefly, however a detailed description of the whole numerical technique can be found in [START_REF] Turek | Efficient Solvers for Incompressible Flow Problems: An Algorithmic Approach in View of Computational Aspects[END_REF].

A semi-discretization of Eqs. ( 1)-( 2) in a semi-implicit manner with respect to time results in:

(1 -η n )ρ f ∂u n+1 ∂t + (1 -η n )ρ f (u n • ∇)u n+1 = (1 -η n )(-∇ p n + ∇ • τ (u n+1 ) + f n ext ) + ρ f η n (u n+1 s -u n+1 ) ∆t f (6)
∇ • u n+1 = 0 (7)
where ∆t f is the fluid time step, superscript n the quantities at a time t n , n + 1 the quantities at time t n+1 , p n+1 the fluid pressure, u n+1 the velocity approximations. Note that the fluid pressure and the velocity approximations have to be corrected later in order to satisfy the Eq. [START_REF] Mohd-Yusof | Combined immersed-boundary/b-spline methods for simulations of flow in complex geometries[END_REF].

The term ∂u n+1 /∂t is an approximate of the temporal velocity derivative which depends on order scheme choice. The projection method is described in four fractional step scheme as follows:

1. Knowing, η n , p n and u n , find an intermediate fluid velocity field ũn+1 such that:

(1 -η n )ρ f ũn+1 -u n ∆t + (1 -η n )ρ f (u n • ∇) ũn+1 = (1 -η n )(-∇ p n + ∇ • τ ( ũn+1 ) + f n ext ) + ρ f η n (u n+1 s -ũn+1 ) ∆t f . (8)
2. Find ϕ such that:

∆ϕ = ρ f ∆t f ∇ • ũn+1 . (9)
3. Correct the fluid velocity with the equation:

u n+1 = ũn+1 - ∆t f ρ f ∇ϕ. (10)
4. Correct the fluid pressure such that:

p n+1 = p n + ϕ. (11)

The drag and lift forces acting on an immersed rigid body

The total fluid force, denoted F F S I , applied on the immersed boundary of a rigid body is defined as:

F F S I =  ∂Ω s σ • nd A ( 12
)
where n denotes the unit normal vector to the surface ∂Ω s . Dividing Eq. ( 1) by (1 -η) and with no external forces, we obtain:

ρ f ∂u ∂t + [ρ f u • ∇u -∇ • σ ] = ρ f (1 -ϵ) ϵ u s -u ∆t f . (13)
In the rigid body, the first and the second terms to the left of Eq. ( 13) become negligible, which leads to the following:

 ∇ • σ = ρ f (1 -ϵ) ϵ u s -u ∆t f  ∀ -→ X (x,y,z)∈Ω s . ( 14
)
Integrating Eq. ( 14) over the domain of the body gives:

 Ω s ∇ • σ dΩ =  Ω s ρ f (1 -ϵ) ϵ u s -u ∆t f dΩ . (15)
Now, applying the Gauss-Ostrogradsky divergence theorem for the term to the left of Eq. ( 15) and using the definition of the force in Eq. ( 12) lead to the following:

 Ω s ∇ • σ dΩ =  ∂Ω s σ • nd A = F F S I =  Ω s ρ f (1 -ϵ) ϵ u s -u ∆t f dΩ . (16)
Eq. ( 16) permits to compute numerically the total fluid force applied on an immersed rigid body of volume Ω s . This method is denoted here as a Penalization Term Integration Technique (PTIT). In fact, it is based on the definition of the coefficient η as a function of ϵ in Eq. ( 5) and on previous works in literature on generalized Brinkman type Models (porous medium with variable permeability) [START_REF] Caltagirone | Sur I'interaction fluide-milieu poreux; application au calcul des efforts exerces sur un obstacle par un fluide visqueux[END_REF][START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier -Stokes equations[END_REF]. Note that in this PTIT technique, we do not impose any empirical formulations like the ones that can be used for example for writing the drag coefficient as function of Reynolds number. Some authors in the literature [START_REF] Zohdi | Embedded electromagnetically sensitive particle motion in functionalized fluids[END_REF] implied empirical formulations for the drag coefficient to handle immersed particulate flows at various flow regimes. However, this limits the solver to handle only some predefined shapes for the immersed bodies (i.e. spheres, squares, cylinders) because the drag coefficient laws for such simple shapes had been well constituted in the literature for general flows at different Reynolds numbers.

The drag and lift forces are defined, respectively, as the following (Fig. 2): The drag C D and lift C L coefficients of the forces acting on a cylinder of diameter D are defined as:

F D = (-(F F S I ) x , 0, 0); F L = (0, -(F F S I ) y , 0). (17)
C D = F D 1 2 ρ f ū2 D ; C L = F L 1 2 ρ f ū2 D ( 18
)
where ū is the characteristic velocity arising from the definition of Reynolds number Re at the scale of a cylinder:

Re = ρ f ū D µ . ( 19
)
The DF-IBM (Eqs. ( 1)-( 2)) and the drag and lift computations are implemented inside the parallelized PELICANS code as a new solver [17].

Validation on a stationary immersed rigid body

In what follows, the numerical framework and its implementation are validated on the transient incompressible 2D flow over a stationary cylinder for two Reynolds numbers, 20 and 100. The numerical results obtained are compared to the Test case 2D-1 (Re = 20 at ū = 0.2 m/s) and Test case 2D-2 (Re = 100 at ū = 1 m/s) in the benchmark results provided by [START_REF] Schäfer | Benchmark computations of laminar flow around a cylinder[END_REF]. Fig. 3 shows the geometry being used with respective dimensions: (Cx = Cy = 0.2 m; D = 0.1; H = 0.41 m and L = 2.2 m). The same boundary conditions (BC) in [START_REF] Schäfer | Benchmark computations of laminar flow around a cylinder[END_REF] were imposed: a homogeneous Dirichlet condition on the walls, a Dirichlet velocity field on the inlet (with a parabolic profile) and a Neumann BC on the outlet.

2D flow over a stationary cylinder at Re = 20

Table 1 shows the parameters for Re = 20. Fig. 4(a)-(d) shows the velocity field, the pressure field and their contours at Re = 20. After 10 s, the flow is already in a steady state. The Fig. 4(e)-(f) shows the drag and lift coefficients computed along the 10 s. The Table 2 presents the comparison between the reference values provided by [START_REF] Schäfer | Benchmark computations of laminar flow around a cylinder[END_REF] and the results.

The obtained values at D/∆ x = 10 for the pressure-difference

∆ p =  p(C x -D 2 , C y ) -p(C x + D 2 , C y )
 and drag coefficient are in good agreement with the benchmark values provided by [START_REF] Schäfer | Benchmark computations of laminar flow around a cylinder[END_REF]. However, a slight difference in the drag is observed and an over estimated value is obtained for the lift coefficient. This may be due to the different schemes (spatial or temporal), different methods being used in different solvers and to the sensitivity of the numerical results on the diameter-to-space-step (DSS) ratio D/∆ x as it is shown in the next section.

Effect of the diameter-to-space-step (DSS) ratio D/∆ x

The sensitivity of the pressure-difference, drag and lift coefficients computations to six different values of the DSS ratio D/∆ x = [5 : 30] is studied. We emphasize that, whenever the D/∆ x is increased, the time step ∆t f is reduced and fixed in the solver to ensure a numerical stability with a 0.2 Courant-Friedrichs-Lewy (CFL) coefficient fulfilling the condition: CFL < 1. We recall that the mesh is always uniform with ∆ y = ∆ x . The Figs. 5 and6 shows clearly the effect of the DSS ratio on the numerical results and the convergence is achieved at 20. The values in Table 3, as expected, show that our computed pressure-difference, drag and lift coefficients values approach the reference benchmark ones [START_REF] Schäfer | Benchmark computations of laminar flow around a cylinder[END_REF] as the DSS ratio increases.

2D flow over a stationary cylinder at Re = 100

Fig. 7(a)-(d) shows the velocity field, the pressure field and their contours at Re = 100. We can see that the symmetry in the flow is broken and vortices are generated behind the wake region of the immersed cylinder which is an expected phenomenon for flows at moderate to high Reynolds numbers. The Fig. 7(e)-(f) shows the drag and lift coefficient versus time at Re = 100. Oscillations are observed due to the vortex generations behind the wake of cylinder. For the benchmark, the maximum values for both drag and lift coefficients are requested. The values are shown in Table 4. By increasing the DSS ratio, our results become close to the benchmark results obtained by [START_REF] Schäfer | Benchmark computations of laminar flow around a cylinder[END_REF].

Immersed rigid body dynamics

The numerical framework

There are several numerical approaches in the literature to account for friction or contact mechanics modeling between the particles. For detailed analyses of such contact mechanics modeling involving friction laws, we mention the early works of [START_REF] Oden | Nonlocal and nonlinear friction laws and variational principles for contact problems in elasticity[END_REF][START_REF] Martins | Existence and uniqueness results in dynamics contact problems with nonlinear normal and friction interfaces[END_REF], and [START_REF] Jean | The non-smooth contact dynamics method[END_REF]. However, our developed FSI solver implemented inside the PELICANS software is coupled to the contact mechanics SSI solver called LMGC90 [22]. The SSI solver LMGC90 relies on the Non Smooth Contact Dynamics approach of [START_REF] Jean | The non-smooth contact dynamics method[END_REF][START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] where it solves for rigid bodies motion (Newton's second law) taking into account the multiple solid to solid interactions based on rigid friction laws. An advantage in this contact mechanics solver is that one can treat polydisperse cases, with total particle/particle contacts-solvings thanks to the NSCD contact/friction laws and the powerful numerical algorithms already implemented in this LMGC90 package. We should mention that the LMGC90 [22] software has been developed over the last 25 years and was very well validated numerically and experimentally on many dry granular media dynamics applications. As we mentioned above, the particles polydispersity in LMGC90 is advantageous for further future developments. In fact, in many previous studies (if not most), for describing the contact laws between the particles, many authors assume that the particles are spherical restricted and remain spherical even after impacts [START_REF] Zohdi | Computation of strongly coupled multifield interaction in particle-fluid systems[END_REF][START_REF] Zohdi | Simulation of coupled microscale multiphysical-fields in particulate-doped dielectrics with staggered adaptive fdtd[END_REF], which limits further possibilities of applying the method to a wide variety of applications.

The dynamical system in our work of immersed bodies in a fluid is given by the following: where F T and T T are the total forces and torques respectively acting on an immersed body, m b the mass, u s = ∂x s ∂t and ω b the linear and angular velocities respectively, g the gravitational acceleration, I b the immersed body inertia matrix, F C and T C the SSI forces and torques respectively, F F S I the FSI force (Eq. ( 16)), and T F S I = r × F F S I the FSI torque with r being the local position relative to the immersed body centroid. The local velocity, on the Eulerian grid, in the immersed body (the right term of Eq. ( 1)) is defined by

F T = m b ∂u b ∂t T T = I b ∂ω b ∂t ( 20
)
F T = F F S I + F C + m b g T T = T F S I + T C (21)
u s = u b + ω × r.
The integration of the system (20) relies on a Leap-frog three-time-level scheme. Note that the solid time step is defined by ∆t s = β∆t f where β ≥ 1 is a relaxation coefficient accounting for the weak coupling between the two FSI and SSI solvers and their convergence.

The numerical platform

Fig. 8 shows the numerical algorithm of the Xper platform. At a time t, upon knowing the different fields: the state η, the fluid velocity, the solid velocity and the pressure, the F F S I force is computed in PELICANS. Then, the new fluid velocity, the pressure fields and the new F F S I force are all computed for a next time step t + ∆t f . The sequence is repeated until the solid time t + ∆t s is reached. Now, the latest F F S I is passed to the LMGC90 solver to solve the dynamics of the immersed solid by computing a new solid velocity including the effects of any contacts with other immersed solids. The latest solid velocity is then passed again to the solver PELICANS.

One of the great advantages in Xper platform is the parallelism. This decreases effectively the computation time needed to conduct a simulation, especially when treating later thousands of immersed bodies in a fluid.

Validation on the 3D settling of a sphere

In what follows, the coupling and the associated software, Xper, are validated on the 3D simulations of a rigid sphere. The results are compared to the experiments of Ten Cate et al. [START_REF] Cate | Particle imaging velocimetry experiments and lattice-Boltzmann simulations on a single sphere settling under gravity[END_REF].

The domain is rectangular with L x = L y = 10 cm and L z = 16 cm and the diameter of the rigid sphere is D = 1.5 cm. The density of the sphere is ρ s = 1120 kg/m 3 . For the fluid, four Reynolds number Re = 1.5, 4.1, 11.6 and 32.3 are considered. The associated parameters are presented in Table 5. The time step respects the CFL criteria:

∆t f = 4 × 10 -5 s.
The sphere is launched at the coordinates of (0.5L x , 0.5L y , 0.75L z ). A value of 20 for β and a 0.5 friction coefficient are used.

The mesh is cartesian and uniform ∆ x = ∆ y = ∆ z = 0.091 cm (about 2.13 millions of cubic cells).

The results for the 3D settling velocity at different Reynolds values are shown in Fig. 9. We observe clearly that Xper captures well the experimental settling velocity magnitudes of the experiment of Ten Cate et al. [START_REF] Cate | Particle imaging velocimetry experiments and lattice-Boltzmann simulations on a single sphere settling under gravity[END_REF] t = 0 s and t = 4.5 s. These results validate the coupling and the developed Xper platform. Indeed, our results are also in very good agreement with the simulations of Wachs [START_REF] Wachs | Peligriff, a parallel dem-dlm/fd direct numerical simulation tool for 3d particulate flows[END_REF] who used a different numerical method (a DEM-DLM/FD method). However, a difference is observed between the numerical and experimental results near the lower wall. This might be explained either by experimental artifacts when using the Particles Image Velocimetry (PIV) technique near the walls, and/or might be due to the lack of numerical precision in capturing well lower velocity magnitudes near the wall due to very small scale lubrication forces. A fact is that in our work, no lubrication force is modeled when the sphere approaches the bottom wall. However, a detailed study about the effects of lubrication forces can be found in the recent works of Izard et al. [START_REF] Izard | Modelling the dynamics of a sphere approaching and bouncing on a wall in a viscous fluid[END_REF]. The latter authors introduced the lubrication force modeling into the discrete element method, which is not captured usually by the fluid solver at a very small scale. This effect of lubrication forces near the wall will be addressed in a precise study in the future using our developed solvers in Xper platform where a very fine mesh near the lower wall will be mandatory.

Application: effect of the fluid ring on the settling of two cylinders initially aligned vertically

In the coming sections, we mean by the term "fluid-ring" any fluid mesh cells that are imposed numerically around mesh cells representing a solid region. When it comes to SSI, most numerical methods in the literature impose a numerical fluid ring around the immersed solids to avoid numerical convergence difficulties. However, this fluid ring cannot represent the real cases. In that manner, the developed platform Xper is used here to compute contacts dynamics problem of two immersed rigid cylinders in an attempt to clarify this important issue and to emphasize the effect of imposing a very small numerical fluid ring on changing the dynamics between solids after contacts.

The application deals with a 2D settling of two cylinders initially aligned vertically. The domain is rectangular with L = 0.32 cm and H = 3.2 cm. The cylinders are of D = 2R = 0.08 cm with a density ρ s = 1200 kg/m 3 . The parameters of the fluid are µ = 0.004313 Pa s and ρ f = 1000 kg/m 3 .

A uniform mesh is used with (∆ x = ∆ y = D/20) and the time step ∆t f = 10 -5 s insures a C F L ≪ 1. A value of 20 for β and 0.5 friction coefficient is used. Fig. 10 shows a rigid cylinder with a fluid ring of ∆R thickness. Two cases are tested to study the effect of friction on the contacts dynamics between the two cylinders: ∆R = 0.1R where friction is absent and ∆R = 0 where friction is present.

In each case, we investigate two initial alignment off-centerline and on-centerline.

Off-centerline

The two cylinders are initially aligned vertically off-centerline at the following centers: (0.25L, 0.96H) and (0.25L, 0.96H-2D).

Fig. 11 shows snapshots of the particle motion at different time. The principal interactions between neighboring particles are described as the well-known phenomenon of drafting (t = 0 s), kissing (t = 0.5 s) and tumbling (t > 2.5 s). First, the trailing particle approaches the leading one ("drafting"). After this, the two particles get into contact at approximately t = 0.5 s ("kissing"). At this time, the gap between the particles is equal to zero without fluid ring or is equal to ∆R with fluid ring. Following the kissing regime, the particles turn to the horizontal position and fall ("tumbling") at t > 2.5 s (for more details, see [START_REF] Ghosh | Numerical simulations of particle sedimentation using the immersed boundary method[END_REF]).

Fig. 12 presents the x and y positions of the two cylinders that were initially aligned vertically off-centerline. The dynamics behaviors obtained by imposing a fluid ring and without fluid ring are small and similar.

0n-centerline

The two cylinders are now initially aligned vertically on-centerline at the following centers: (0.5L, 0.96H) and (0.5L, 0.96H-2D).

Fig. 13 shows snapshots of the particle motion at different time. In this case, only the drafting mechanism is involved. The behaviors are different since the first particle appears to touch the down wall. With friction, i.e. without fluid ring, the particles are pushed away from the down wall. With fluid ring, the particles are sucked all the way to the down wall.

Fig. 14 presents the x and y positions of the two cylinders that were initially aligned vertically on-centerline. We observe a greater effect on the dynamics behavior by imposing a fluid ring, compared to the cases when the two cylinders were initially aligned vertically off-centerline.

The interaction of the particles settling in tandem can be fundamentally altered due to the presence of the fluid ring. The Figs. 12 and 14 emphasize an important new finding. This finding implies that imposing at all even a very small fluid ring, the obtained dynamics cannot correspond well to the "real" phenomena.

Conclusion

A Direct-Forcing Immersed Boundary Method (DF-IBM) is used and implemented in the parallelized numerical framework named PELICANS. This was to take into account for the presence of stationary and moving rigid bodies that are immersed in a fluid (Fluid-Structure Interactions). The numerical implementation is validated by simulating the transient 2D flow around a stationary cylinder at two different Reynold's numbers (20 and 100). Good results are obtained in computing the drag and lift coefficients, and the pressure difference using a PTIT (Penalization Term Integration Technique) inside the immersed body. These results have been characterized by studying the effect of the diameter-to-space-step ratio (DSS) on numerical precision. We found that a DSS of 20 was fair to obtain results that are close to the benchmark data provided by Turek and Schäfer [START_REF] Schäfer | Benchmark computations of laminar flow around a cylinder[END_REF] for the 2D flow over a stationary cylinder.

For multiple contacts (Solid-Solid Interactions) between immersed rigid bodies, the PELICANS framework was coupled to the Contact-Mechanics framework named LMGC90 and the resulting package is named Xper.

The 3D settling of a sphere has been simulated at different Reynolds numbers (between 1.5 and 32.2) using Xper to test the validity of the numerical package as a whole. Satisfactory results are obtained for the settling velocity values compared to the experimental and numerical data found in the literature [START_REF] Wachs | Peligriff, a parallel dem-dlm/fd direct numerical simulation tool for 3d particulate flows[END_REF][START_REF] Cate | Particle imaging velocimetry experiments and lattice-Boltzmann simulations on a single sphere settling under gravity[END_REF]. This indicates the primary validity of Xper in simulating immersed granular flows in 3D.

Finally, the effect on the contacts dynamics is addressed by imposing numerical fluid cells around two identical settling cylindrical particles in a stationary fluid as it is usually done by many authors in literature. We found that, imposing this numerical fluid ring (of thickness of the order of a one mesh space) around the particle (as done usually by most authors using different numerical methods to avoid numerical difficulties) to neglect friction alters the particles trajectories. This is a new finding that implies that imposing at all even a very small fluid ring, the obtained dynamics cannot correspond well to the "real" phenomena because the friction due to direct solid-to-solid contacts is present. Moreover, a Coulomb's friction coefficient of 0.5 has shown to have a great effect on the particles dynamics after contact.

This can be crucial in defining the macroscopic behavior when more than two particles may enter in contact during an immersed granular flow based on a certain microstructure orientation.

Further work is required to study the effect of introducing a level-set method, i.e. a smooth edge for the Fluid-Solid Interface, on the numerical results of the DF-IBM and on the computation time performance. Further work is also required to study the effect of a range of Coulomb's friction coefficient values on the particles dynamics.

Numerical simulations involving immersed granular flows, i.e. thousands of particles, will be conducted in the near future using Xper in an attempt to constitute macroscopic rheology laws taking into account the physics at the microscopic scale.

Fig. 1 .

 1 Fig. 1. An immersed boundary.

Fig. 2 .

 2 Fig. 2. Drag and lift forces acting on an immersed cylinder Ω s .

Fig. 4 .

 4 Fig. 4. 2D flow over a stationary cylinder at Re = 20 and D/∆ x = 10.

Fig. 5 .

 5 Fig. 5. Drag, lift and pressure-difference coefficients versus the DSS ratio.

Fig. 6 .

 6 Fig. 6. Transient drag and lift coefficients at DSS = 30.

  (a) Velocity field. (b) Pressure field. (c) Velocity contour. (d) Pressure contour. (e) Drag coefficient. (f) Lift coefficient.

Fig. 7 .

 7 Fig. 7. 2D flow over a stationary cylinder at Re = 100 and D/∆ x = 10.

Fig. 8 .

 8 Fig. 8. Numerical algorithm implemented in the Xper code.

Fig. 10 .

 10 Fig. 10. Rigid cylinder with an imposed numerical fluid ring of ∆R.

Fig. 11 .

 11 Fig. 11. 2D settling of two cylinders initially aligned vertically and off-centerline at t = 0 s, t = 0.5 s, t = 1 s, t = 2.5 s, t = 4 s, t = 5 s and t = 6.5 s.

  Normalized x-positions versus time. (b) Normalized y-positions versus time.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. 2D settling of two cylinders initially aligned vertically and off-centerline.

  (a) Normalized x-positions versus time. (b) Normalized y-positions versus time.

Fig. 14 .

 14 Fig. 14. 2D settling of two cylinders initially aligned vertically and on-centerline.

Table 1 2D

 1

	Uniform mesh	D/∆ x	∆t f	µ	ρ

flow over a stationary cylinder: parameters at Re = 20. f ∆ x = ∆ y = 0.01 10 0.01 s 0.001 Pa s 1 kg/m 3

Table 2

 2 Steady values for pressure-difference, drag and lift coefficients at Re = 20.

		Reference values [18]	Present solver (D/∆ x = 10)
	C D	[5.57:5.59]	5.3273
	C L	[0.0104:0.0110]	0.02017
	∆ p	[0.1172:0.1176]	0.10570

Table 3

 3 Values for pressure-difference, drag and lift coefficients at Re = 20 and t = 10 s.

	Present solver D/∆ x	C D	C L	∆ p
	5	5.3398	0.2267	0.1017
	10	5.3273	0.02017	0.10570
	15	5.3439	0.1611	0.10864
	20	5.4456	0.0142	0.1140
	25	5.4460	0.0164	0.1141
	30	5.4752	0.0140	0.1156
	Reference values [18]	[5.57:5.59]	[0.0104:0.0110]	[0.1172:0.1176]
	Range variation [18]	[5.50:5.81]	[0.0104:0.0329]	[0.1054:0.1230]

Table 4 Values

 4

	Present solver D/∆ x	C D	C L	∆ p
	5	2.97	0.574	2.0959
	10	2.788	0.441	2.2158
	15	2.911	0.9316	2.3128
	20	2.980	0.715	2.3438
	25	2.976	0.720	2.3537
	30	3.076	0.920	2.3710
	Reference values [18]	[3.22:3.24]	[0.99:1.010]	[2.46:2.50]
	Range variation [18]	[3.08:4.12]	[0.72:1.61]	[2.13:3.14]

for pressure-difference, drag and lift coefficients at Re = 100 and t = 50 s.

  between

	Reynolds	ρ f	µ
	1.5	970	0.373
	4.1	965	0.212
	11.6	962	0.113
	32.3	960	0.058
		Fig. 9. Sphere settling velocity versus time.

Acknowledgments

The author appreciated and thanks V. Topin, F. Dubois, J.-C. Latché and B. Piar for their useful scientific discussions.