Nicolò Cesa-Bianchi

Tommaso R Cesari

Claire Monteleoni

Cooperative Online Learning: Keeping your Neighbors Updated

We study an asynchronous online learning setting with a network of agents. At each time step, some of the agents are activated, requested to make a prediction, and pay the corresponding loss. The loss function is then revealed to these agents and also to their neighbors in the network. Our results characterize how much knowing the network structure affects the regret as a function of the model of agent activations. When activations are stochastic, the optimal regret (up to constant factors) is shown to be of order √ αT , where T is the horizon and α is the independence number of the network. We prove that the upper bound is achieved even when agents have no information about the network structure. When activations are adversarial the situation changes dramatically: if agents ignore the network structure, a Ω(T ) lower bound on the regret can be proven, showing that learning is impossible. However, when agents can choose to ignore some of their neighbors based on the knowledge of the network structure, we prove a O( √ χT ) sublinear regret bound, where χ ≥ α is the clique-covering number of the network.

Introduction

Distributed asynchronous online learning settings with communication constraints arise naturally in several applications. For example, large-scale learning systems are often geographically distributed, and in domains such as finance or online advertising, each agent must serve high volumes of prediction requests. If agents keep updating their local models in an online fashion, then bandwidth and computational constraints may preclude a central processor from having access to all the observations from all sessions, and synchronizing all local models at the same time. An example in a different domain is mobile sensor networks cooperating towards a common goal, such as environmental monitoring. Sensor readings provide instantaneous, fullinformation feedback and energy-saving constraints favor short-range communication over long-range. Online learning algorithms distributed over spatial locations have already been proposed for problems in the field of climate informatics by McQuade and Monteleoni [START_REF] Mcquade | Global climate model tracking using geospatial neighborhoods[END_REF][START_REF] Mcquade | Large-Scale Machine Learning in the Earth Sciences[END_REF], and have shown empirical performance advantages compared to their global (i.e., non-spatially distributed) online learning counterparts.

Motivated by these real-life applications, we introduce and analyze an online learning setting in which a network of agents solves a common online convex optimization problem by sharing feedback with their network neighbors. Agents do not have to be synchronized. At each time step, only some of the agents are requested to make a prediction and pay the corresponding loss: we call these agents "active". Because the feedback (i.e., the current loss function) received by the active agents is communicated to their neighbors, both active agents and their neighbors can use the feedback to update their local models. The lack of global synchronization implies that agents who are not requested to make a prediction get "free feedback" whenever someone is active in their neighborhood. Since in online convex optimization the sequence of loss functions is fully arbitrary, it is not clear whether this free feedback can improve the system's performance. In this paper, we characterize under which conditions and to what extent such improvements are possible.

Our goal is to control the network regret, which we define by summing the average instantaneous regret of the active agents at each time step. In order to build some intuition on this problem, consider the following two extreme cases where, for the sake of simplicity, we assume exactly one agent is active at each time step. If no communication is possible among the agents, then each agent v learns in isolation over the subset T v of time steps when they are active. Assuming each agent runs a standard online learning algorithm with regret bounded by O( √ T ) -such as Online Mirror Descent (OMD)-the network regret is at most of order v √ T v ≤ √ N T , where T = v T v and N is the number of agents. Next, consider a fully connected graph, where agents share their feedback with the rest of the network. Each local instance of OMD now sees the same loss sequence as the other instances, so the sequence of predictions is the same, no matter which agents are chosen to be active. The network regret is then bounded by O( √ T ), as in the single-instance case. Our goal is to understand the regret when the communication network corresponds to an arbitrary graph G.

We consider two natural activation mechanisms for the agents: stochastic and adversarial. In the stochastic setting, at each time step t each agent v is independently active with probability q v , where q v is a fixed and unknown number in [0, 1]. Under this assumption, we show that when each agent runs OMD, the network regret is O( √ αT ), where α ≤ N is the independence number of the communication graph. Note that this bound smoothly interpolates the two extreme cases of no communication (α = N ) and full communication (α = 1). From this viewpoint, α can be viewed as the number of "effective instances" that are implicitly maintained by the system. It is not hard to prove that this upper bound cannot be improved upon: fix a network G and a maximal independent set in G of size α. Define q v = 1/α if v belongs to the independent set and 0 otherwise. Then no two nodes that can ever become active are adjacent in G, and we reduced the problem to that of learning with α non-commmunicating agents over T /α time steps. Since there are instances of the standard online convex optimization problem on which any agent strategy has regret Ω( √ T ), we obtain that the network regret must be at least of order α T /α = √ αT . Note that this lower bound also applies to algorithms that have complete preliminary knowledge of the graph structure, and can choose to ignore or process any feedback coming from their neighbors. In contrast, the OMD instances used to prove the upper bound are fully oblivious both to the graph structure and to the source of their feedback (i.e., whether their agent is active as opposed to being the neighbor of an active agent).

In the adversarial activation setting, nodes are activated according to some unknown deterministic schedule. Surprisingly, under the same assumption of obliviousness about the feedback source which we used to prove the O( √ αT ) upper bound for stochastic activations, we show that on certain network topologies a deterministic schedule of activations can force a linear regret on any algorithm, thus making learning impossible. On the other hand, if agents are free to use feedback only from a subset of their neighbors chosen with knowledge of the graph structure, then the network regret of OMD is O( √ χT ), where χ is the clique-covering number of the communication graph. Hence, unlike the stochastic case, where the knowledge of the graph is not required to achieve optimality, in the adversarial case the ability of choosing the feedback source based on the graph structure is both a necessary and sufficient condition for sublinear regret.

The extension of the OMD analysis to a multiagent setting with communication (Theorem 2), and the lower bound for the adversarial activation setting (Theorem 5) are the main technical novelties of the paper.

Related Works

The study of cooperative nonstochastic online learning on networks was pioneered by Awerbuch and Kleinberg [START_REF] Awerbuch | Competitive collaborative learning[END_REF], who investigated a bandit setting in which the communication graph is a clique, users are clustered so that the loss function at time t may differ across clusters, and some users may be non-cooperative. More recently, a similar line of work was pursued by Cesa-Bianchi et al. [START_REF] Cesa-Bianchi | Delay and cooperation in nonstochastic bandits[END_REF], where they derive graph-dependent regret bounds for nonstochastic bandits on arbitrary networks when the loss function is the same for all nodes and the feedbacks are broadcast to the network with a delay corresponding to the shortest path distance on the graph. Although their regret bounds -like ours-are expressed in terms of the network independence number, this happens for reasons that are very different from ours, and by means of a different analysis.

In their setting all agents are simultaneously active at each time step, and sharing the feedback serves the purpose of reducing the variance of the importance-weighted loss estimates. A node with many neighbors observes the current loss function evaluated at all the points corresponding to actions played by the neighbors. Hence, in that context cooperation serves to bring the bandit feedback closer to a full information setting.

In contrast, we study a full information setting in which agents get free and meaningful feedback only when they are not requested to predict. 1 Therefore, in our setting cooperation corresponds to faster learning (through the free feedback that is provided over time) within the full information model, as opposed to [START_REF] Cesa-Bianchi | Delay and cooperation in nonstochastic bandits[END_REF] where cooperation increases feedback within a single time-step. An even more recent work considering bandit networks studies a stochastic bandit model with simultaneous activation and constraints on the amount of communication between neighbors [START_REF] Martínez-Rubio | Decentralized cooperative stochastic multi-armed bandits[END_REF]. Their regret bounds scale with the spectral gap of the communication network. The work of Sahu and Kar [START_REF] Kumar | Dist-Hedge: A partial information setting based distributed non-stochastic sequence prediction algorithm[END_REF] investigates a different partial information model of prediction with expert advice where each agent is paired with an expert, and agents see only the loss of their own expert. The communication model includes delays, and the regret bound depends on a quantity related to the mixing time of a certain random walk on the network. Zhao et al. [START_REF] Zhao | Decentralized online learning: Take benefits from others' data without sharing your own to track global trend[END_REF] study a decentralized online learning setting in which losses are characterized by two components, one adversarial and another stochastic. They show upper bounds on the regret in terms of a constant representing the magnitude of the adversarial component and another constant measuring the randomness of the stochastic part.

The idea of varying the amount of feedback available to learning agents has also appeared in single-agent settings. In the sleeping experts model [START_REF] Freund | Using and combining predictors that specialize[END_REF], different subsets of actions are available to the learner at different time steps. In our multi-agent setting, instead, actions are always available while the agents are occasionally sleeping. An algorithmic reduction between the two settings seems unlikely to exist because actions and agents play completely different roles in the learning process. In the learning with feedback graphs model [START_REF] Noga Alon | Online learning with feedback graphs: Beyond bandits[END_REF][START_REF] Mannor | From bandits to experts: On the value of side-observations[END_REF], each selection of an action reveals to the learner the loss of the actions that are adjacent to it in a given graph. In our model, each time an active agent plays an action, the loss vector is revealed to the agents that are adjacent to the active learner. There is again a similarity between actions and agents in the two settings, but to the best of our knowledge there is no algorithmic reduction from multi-agent problems to single-agent problems. Yet, it should not come as a surprise that some general graph-theoretic tools -like Lemma 2-are used in the analysis of both single-agent and multi-agent models.

A very active area of research involves distributed extensions of online convex optimization, in which the global loss function is defined as a sum of local convex functions, each associated with an agent. Agents are run over the local optimization problem corresponding to their local functions and communicate with their neighborhood to find a point in the decision set approximating the loss of the best global action. This problem has been studied in various settings: distributed convex optimization -see, e.g., [START_REF] John C Duchi | Dual averaging for distributed optimization: Convergence analysis and network scaling[END_REF][START_REF] Scaman | Optimal algorithms for non-smooth distributed optimization in networks[END_REF] and references therein, distributed online convex optimization [START_REF] Hosseini | Online distributed optimization via dual averaging[END_REF], and a dynamic regret extension of distributed online convex optimization [START_REF] Shahrampour | Distributed online optimization in dynamic environments using mirror descent[END_REF]. Unlike our work, these papers consider distributed extensions of OMD (and Nesterov dual averaging) based on generalizations of the consensus problems. The resulting performance bounds scale inversely in the spectral gap of the communication network.

Preliminaries and definitions

Let G = (V, E) be a communication network, i.e., an undirected graph over a set V of N agents. Without loss of generality, assume V = {1, . . . , N }. For any agent v ∈ V , we denote by N v the set of nodes containing the agent v and the neighborhood w ∈ V | (v, w) ∈ E . The independence number α G is the cardinality of the biggest independent set of G, i.e., the cardinality of the biggest subset of agents, no two of which are neighbors.

We study the following cooperative online convex optimization protocol: initially, hidden from the agents, the environment picks a sequence of subsets S 1 , S 2 , . . . ⊆ V of active agents and a sequence of differentiable Algorithm 1 Online Mirror Descent Parameters: σ t -strongly convex regularizers g t : X → R for t ∈ {1, 2, . . .} Initialization:

θ 1 = 0 ∈ R d 1: for t ∈ {1, 2, . . .} do 2: choose w t = ∇g * t (θ t ) 3: observe ∇ℓ t (w t ) ∈ R d 4:
update θ t+1 = θ t -∇ℓ t (w t ) convex real loss functions ℓ 1 , ℓ 2 , . . . defined on a convex decision set X ⊂ R d . Then, for each time step t ∈ {1, 2, . . .}, 1. each agent v ∈ S t predicts with x t (v) ∈ X , 2. each agent v ∈ v∈St N v receives ℓ t as feedback, 3. the system incurs the loss 1 |St| v∈St ℓ t x t (v) (defined as 0 when S t ≡ ∅). We assume each agent v runs an instance of the same online algorithm. Each instance learns a local model generating predictions x t (v). This local model is updated whenever a feedback ℓ t is received. We call paid feedback the feedback ℓ t received by v when v ∈ S t (i.e., the agent is active) and free feedback the feedback ℓ t received by v when v ∈ v∈St N v \ {S t } (i.e., the agent is not active but in the neighborhood of some active agent). The goal is to minimize the network regret as a function of the unknown number T of time steps,

R T = T t=1 1 |S t | v∈St ℓ t x t (v) -inf x∈X T t=1 ℓ t (x) (1)
Note that only the losses of active agents contribute to the network regret.

Online Mirror Descent

We now review the standard Online Mirror Descent algorithm (OMD) -see Algorithm 1-and its analysis. Let f : X → R be a convex function. We say that f * :

R d → R is the convex conjugate of f if f * (x) = sup w∈X x • w -f (w) .
We say that f is σ-strongly convex on X with respect to a norm • if there exists σ ≥ 0 such that, for all u, w ∈ X we have

f (u) ≥ f (w) + ∇f (w) • (u -w) + σ 2 u -w 2 .
The following well-known result can be found in [17, Lemma 2.19 and subsequent paragraph].

Lemma 1. Let f : X → R be a strongly convex function on X . Then the convex conjugate f * is everywhere differentiable on R d .

The following result -see, e.g., [13, bound (6) in Corollary 1 with F set to zero]-shows an upper bound on the regret of OMD. Theorem 1. Let g : X → R be a differentiable function σ-strongly convex with respect to • . Then the regret of OMD run with

g t = √ t η g, for η > 0, satisfies T t=1 ℓ t x t -inf x∈X T t=1 ℓ t (x) ≤ D η √ T + η 2σ T t=1 1 √ t ∇ℓ t 2 *
where D = sup g and

• * is the dual norm of • . If sup ∇ℓ t * ≤ L, then choosing η = √ 2σD/L gives R T ≤ L 2DT /σ.
A popular instance of OMD is the standard online gradient descent algorithm, corresponding to choosing X equal to a closed Euclidean ball centered at the origin, and setting g = 1 2 • 2 for all t, where • is the Euclidean norm. Another instance is the Hedge algorithm for prediction with expert advice, corresponding to choosing X equal to the probability simplex, and setting g(p) = i p i ln p i .

Stochastic Activations

In this section we analyze the performance of OMD when the sets S t of active agents are chosen stochastically. As discussed in the introduction, in this setting we do not require any ad-hoc interface between each OMD instance and the rest of the network. In particular, we make the following assumption.

Assumption 1 (Oblivious network interface). An online algorithm A is run with an oblivious network interface if for each agent v it holds that:

1. v runs an instance A v of A 2. A v uses the same initialization and learning rate as the other instances 3. A v makes predictions and updates while being oblivious to whether v

∈ S t or v ∈ u∈St N u \ S t
This assumption implies that each instance is oblivious to both the network topology and the location of the agent in the network. Moreover, instances make an update whenever they have the opportunity to do so, (i.e., whenever they or some of their neighbors are active). The purpose of this assumption is to show that communication might help OMD even without any network-specific tuning. In concrete applications, one might use ad-hoc OMD variants that rely on the knowledge of the task at hand, and decrease the regret even further. However, the lower bound proven in Section 6 shows that the regret cannot be decreased significantly even when agents have full knowledge of the graph.

We start by considering a slightly simplified stochastic activation setting, where only a single agent is activated at each time step (i.e., |S t | = 1 for all t). The more general stochastic case is analyzed at the end of this section.

We assume that the active agents v 1 , v 2 , . . . are drawn i.i.d. from an unknown fixed distribution q on V . The main result of this section is an upper bound on the regret of the network when all agents run the basic OMD (Algorithm 1) with an oblivious network interface. We show that in this case the network achieves the same regret guarantee as the single-agent OMD (Theorem 1) multiplied by the square root of independence number of the communication network.

Before proving the main result, we state a combinatorial lemma that allows to upper bound the sum of a ratio of probabilities over the vertices of an undirected graph with the independence number of the graph [START_REF] Griggs | Lower bounds on the independence number in terms of the degrees[END_REF][START_REF] Mannor | From bandits to experts: On the value of side-observations[END_REF]. The proof is included for completeness.

Lemma 2. Let G = (V, E) be an undirected graph with independence number α G and q any probability distribution on V such that

Q v = w∈Nv q v > 0 for all v ∈ V . Then v∈V q v Q v ≤ α G Proof. Initialize V 1 = V , fix w 1 ∈ arg min w∈V1 Q w , and denote V 2 = V \ N w1 . For k ≥ 2 fix w k ∈ arg min w∈V k Q w and shrink V k+1 = V k \ N w k until V k+1 = ∅. Being G undirected, we have w k / ∈ k-1
s=1 N ws , therefore the number m of times that an action can be picked this way is upper bounded by α G . Denoting

N ′ w k = V k ∩ N w k , this implies v∈V q v Q v = m k=1 v∈N ′ w k q v Q v ≤ m k=1 v∈N ′ w k q v Q w k ≤ m k=1 v∈Nw k q v Q w k = m ≤ α G
The following holds for any differentiable function g : X → R, σ-strongly convex with respect to some norm • .

Theorem 2. Consider a network G = (V, E) of N agents and assume S t = {v t } for each t, where v t is drawn i.i.d. from some fixed and unknown distribution on V . If all agents run OMD with an oblivious network interface and using g t = √ t η g, for η > 0, then the network regret satisfies

E[R T ] ≤ D η + ηL 2 2σ α G T
where D ≥ sup g, L ≥ sup ∇ℓ t * , and

• * is the dual norm of • . In particular, choosing η = √ 2σD/L gives E[R T ] ≤ L 2Dα G T /σ.
Proof. Fix x ∈ X , any sequence of realizations v 1 , . . . , v T , and any v in the support V ′ ⊂ V of the activation distribution q. Note that the OMD instance run by v, makes an update at time t only when v ∈ N vt . Hence, letting r t (v) = ℓ t x t (v)ℓ t (x) and applying Theorem 1,

T t=1 r t (v)I{v ∈ N vt } ≤ D η T v + ηL 2 2σ T t=1 I{v ∈ N vt } t s=1 I{v ∈ N vs } ≤ D η + ηL 2 2σ T v ( 2
)
where

T v = T t=1
I{v ∈ N vt }, the addends after the first inequality are intended to be null when the denominator is zero, and we used

Tv s=1 s -1 /2 ≤ 2 √ T v . Note that r t (v) is independent of v t , as it only depends on the subset of v s , s ∈ {1, . . . , t -1}, such that v ∈ N vs . Denote by Q v the probability P(v ∈ N vt ) = w∈Nv q(w) > 0. Let F t-1 be the σ-algebra generated by {v 1 , . . . , v t-1 }. Since Q v is independent of t, P v ∈ N vt | F t-1 = Q v .
Therefore, taking expectation with respect to v 1 , . . . , v T on both sides of (2), using E[T v ] = Q v T , and applying Jensen's inequality, yields

E T t=1 r t (v)Q v ≤ D η + ηL 2 2σ Q v T (3) Now, letting R T (x) = T t=1 r t (v t ), we have that E R T (x) is equal to E v∈V ′ T t=1 r t (v)I{v t = v} = E v∈V ′ T t=1 r t (v)E I{v t = v} | F t-1 = v∈V ′ q v E T t=1 r t (v)
Dividing both sides of (3) by Q v > 0, we can write

E R T (x) ≤ D η + ηL 2 2σ v∈V ′ q v T Q v ≤ D η + ηL 2 2σ T v∈V ′ q v Q v ≤ D η + ηL 2 2σ √ αT
where in the last two inequalities we applied Jensen's inequality and Lemma 2. Observing that E[R T ] = sup x∈X E R T (x) and recalling that x was chosen arbitrarily in X concludes the proof.

Note that the proof of the previous result gives a tighter upper bound on the network regret in terms of the independence number α ′ ≤ α of the subgraph induced by the support V ′ of q.

Next, we consider the setting in which we allow the activation of more than one agent per time step. At the beginning of the process, the environment draws an i.i.d. sequence of Bernoulli random variables X 1 (v), X 2 (v), . . . with some unknown fixed parameter q v ∈ [0, 1] for each agent v ∈ V . The active set at time t is then defined as S t = {v ∈ V | X t (v) = 1}. Note that, unlike the previous setting, now v∈V q v = 1 in general.

We state an upper bound on the regret that the network incurs if all agents run OMD with an oblivious network interface (for a proof, see Appendix A). Our upper bound is expressed in terms of a constant depending on the probabilities of activating each agent and such that Q ≤ 1.6(α G + 1). The result holds for any differentiable function g : X → R, σ-strongly convex with respect to some norm • . Theorem 3. Consider a network G = (V, E) of N agents. Assume that, at each time step t each agent v is independently activated with probability q v ∈ [0, 1]. If all agents run OMD with an oblivious network interface and using g t = √ t η g, for η > 0, the network regret satisfies

E[R T ] ≤ D η + ηL 2 2σ QT
for some nonnegative Q ≤ 1.6(α G + 1), D ≥ sup g, and L ≥ sup ∇ℓ t * , where

• * is the dual norm of • . In particular, choosing η = √ 2σD/L gives E[R T ] ≤ L 2DQT /σ.
In order to compare the previous upper bound to Theorem 2, consider the case q v = q for all v ∈ V . Without loss of generality, assume q > 0 (the regret is zero when q vanishes). Then

Q = Q(q) = 1 N v∈V 1 -(1 -q) N 1 -(1 -q) |Nv |
(for a proof, see Theorem 7 in Appendix A and proceed as in the proof of Lemma 3). A direct computation of the sign of the first derivative of the addends q → 1-(1-q) N 1-(1-q) |Nv | shows that these functions are decreasing in q, hence 1 = lim

q→1 -Q(q) ≤ Q ≤ lim q→0 + Q(q) = v∈V 1 |Nv | ≤ α G
where the last inequality follows by Lemma 2. Note that the lower bound Q ≥ 1 is attained if the probabilities of picking agents at each time step are all 1. In this case all agents are activated at each time step, the graph structure over the set of agents becomes irrelevant and the model reduces to a single-agent problem. The inequality Q(q) ≤ α G is not a coincidence due to the constant q. Indeed, one can prove that this is always the case, up to a small constant factor (for a proof., see Lemma 4 in Appendix A).

The previous results shows that paying the average price of multiple activations is never worse (up to a small constant factor) than drawing a single agent per time step, and it can be significantly better. A similar argument shows a tighter bound Q ≤ max{3, α G } when the activation probabilities satisfy v∈V q v = 1, which allows to recover the upper bound on the network regret proven in Theorem 2. This is consistent with the intuition that -in expectation-picking a single agent at random according to a distribution q = (q 1 , . . . , q N ) is the same as picking each v independently with probability q v . Similarly to the case |S t | = 1, the previous result gives a tighter upper bound on the network regret in terms of the independence number α ′ ≤ α of the subgraph induced by the subset V ′ of V containing all agents v with q v > 0. Note that the setting discussed in this section smoothly interpolates between the single-agent setting (q v = 1 for all v), cooperative learning with one agent stochastically activated at each time step ( v q v = 1), and beyond ( v q v < 1), where a non trivial fraction of the total number rounds is skipped.

Lower Bound for Stochastic Activations

In this section we show that, for any communication network G with stochastic agent activations, the best possible regret rate is of order Ω √ α G T . This holds even when agents are not restricted to use an oblivious network interface. The idea is that if the distribution from which active agents are drawn is supported on an independent set of cardinality α G , then the problem reduces to that of an edgeless graph with α G agents.

Theorem 4. There exists a convex decision set in R d such that, for each communication network G and for arbitrary (and possibly different) online learning algorithms run by the agents, E[R T ] = Ω

√ αT for some sequence (S 1 , ℓ 1 ), . . . , (S T , ℓ T ), where S t = {v t }, v t is drawn i.i.d. from some fixed distribution on V , and the expectation is taken with respect to the random draw of the v 1 , . . . , v T .

Proof. We sketch the proof for the case |S t | = 1. Let X be the probability simplex in R d . Let G = (V, E) be any communication graph and α its independence number. We consider linear losses defined on X . Let q be the uniform distribution over a maximal independent set A = {a 1 , . . . , a α } ⊂ V . Fix now any cooperative online linear optimization algorithm for this setting. Since each active agent v t belongs to A for all t ∈ {1, . . . , T } with probability 1, it suffices to analyze the updates of the algorithm for these agents. Indeed, no other agent incurs any loss at any time-step. Since A is an independent set, each agent a i makes an update at round t if and only if v t = a i . This happens with probability q(a i ) = 1/α, independently of t. Each agent a i is therefore running an independent single-agent online linear optimization problem for an average of T /α rounds. It is well-known [7, Theorem 3.2] that any algorithm for online linear optimization on the simplex with losses bounded in [0, 1] incurs Ω T /α regret over T /α rounds in the worst case. Consequently, the regret of the network satisfies R T = Ω α T /α = Ω √ αT .

An analogous lower bound can be proven for the case of multiple agent activations per time step. Indeed, define q v = 1/α for each agent v belonging to some fixed maximal independent set and q v = 0 otherwise. This again leads to α independent single-agent online linear optimization problems for an average of T /α rounds each, and an argument similar to the one in the proof of Theorem 4 gives the result.

Adversarial Activations

In this section we drop the stochasticity assumption on the agents' activations and focus on the case where active agents are picked from V by an adversary. The goal is to control the regret (1) for any individual sequence of pairs (ℓ 1 , S 1 ), (ℓ 2 , S 2 ), . . . where ℓ t is a convex loss and S t ⊆ V , without any stochastic assumptions on the mechanism generating these pairs. For the rest of this section, we focus on the special case where |S t | = 1 for all t and denote by v t the active node at time t.

We start by proving that learning with adversarial activations is impossible if we use an oblivious network interface. We prove this result in the setting of prediction with expert advice with two actions and binary losses, a special case of online convex optimization. The idea of the lower bound is that if the communication network is a star graph, the environment is able to make both actions look equally good to all peripheral agents, even if one of the two actions is actually better than the other. This is done by drawing the good action at random, then activating an agent depending on the outcome of the draw. For a small fraction of the times the good action has loss one, the central agent is activated. Since the central agent shares feedback with all peripheral agents, we can amplify this loss by a factor of N , and thus make the good action look to all peripheral agents as bad as the bad action. Theorem 5. For each N > 3 there exists a convex decision set in R 2 and a graph G with N vertices such that, whenever N agents are run on G using instances of any online learning algorithm with an oblivious network interface, then R T = Ω(T ) for some sequence (ℓ 1 , v 1 ), . . . , (ℓ T , v T ) of convex losses and active agents.

Proof. Fix N > 3 and let X be the probability simplex in R 2 . Let G = (V, E) be the star graph with central agent a 0 , and peripheral agents a 1 , . . . , a N -1 . Because our losses are linear on X , the online convex optimization problem is equivalent to prediction with expert advice with two experts (or actions), and we may denote losses using loss vectors ℓ t = ℓ t (1), ℓ t (2) where 1 and 2 index the actions. A good action J ∈ {1, 2} is drawn uniformly at random. Denote the other one (i.e., the bad one) by J B . To keep notation tidy, we define loss vectors by ℓ t = ℓ t (J), ℓ t (J B ) . Fix any ε ∈ 0, N -1 2(N -2) . The loss vectors ℓ t are drawn i.i.d. at random, according to the following joint distribution:

P ℓ t = (0, 1) = 1 2 P ℓ t = (1, 0) = 1 2 -ε + ε N -1 P ℓ t = (0, 0) = ε - ε N -1
Recall that only a single agent v t is active at any time. At each time step t, the adversary decides whether to activate the central agent a 0 or a peripheral agent, depending on the realization of ℓ t . If ℓ t (J) = 0, then a random peripheral agent is activated. Otherwise, we set

P ℓ t = (1, 0), v t = a 0 = ε N -1 and P ℓ t = (1, 0), v t = a i = 1 /2 -ε N -1 a 1 , . . . , a N -1
Note that when v t = a 0 , then all peripheral agents receive feedback ℓ t . Similarly, when a peripheral agent is active at time t, then a 0 receives feedback ℓ t . For b 1 , b 2 ∈ {0, 1}, let E(a i , b 1 , b 2 ) be the event: agent a i receives the loss vector ℓ t = (b 1 , b 2 ) as feedback. The following statements then hold for each peripheral agent a i ,

P E(a i , 0, 1) = 1 /2 N -1 P E(a i , 0, 0) = ε N -1 - ε (N -1) 2 P E(a i , 1, 0) = 1 /2 -ε N -1 + ε N -1 = 1 /2 N -1
Hence, each instance managed by a peripheral agent observes loss vectors (1, 0) and (0, 1) with the same probability proportional to 1 /2, and loss vector (0, 0) with probability proportional to ε(N -1)/(N -2). Since the network interface is oblivious, the instance cannot distinguish between paid and free feedback (which would reveal the good action), and incurs an expected loss of 1 /2 each time ℓ t ∈ (0, 1), (1, 0) . Using the fact that a peripheral agent is active when ℓ t ∈ (0, 1), (1, 0) with probability 1 /2 + 1 /2ε = 1ε, the system's expected total loss is at least 1-ε 2 T (we lower bound the loss of the central agent by zero). Since the expected loss of J is 1 /2ε + ε N -1 T , the expected regret of the system satisfies

E[R T ] ≥ 1 -ε 2 - 1 2 + ε - ε N -1 T ≥ T 8
where we picked ε = (N -1)/(N -2) and used (N -3)/(N -2) ≥ 1 /2 in the last inequality. Therefore, there exists some sequence (ℓ 1 , S 1 ), . . . , (ℓ T , S T ) such that R T ≥ T /8, concluding the proof.

We complement the above negative result by showing that when algorithms are run without the oblivious network interface, and agents are free to use feedback only from a subset of their neighbors chosen with knowledge of the graph structure, then the network regret of OMD is O( χ G T ). The quantity χ G is the clique-covering number of the communication graph G, which corresponds to the smallest cardinality of a clique cover of G (a clique cover is a partition of the vertices such that the nodes in every element of the partition form a clique in the graph). The intuition behind this result is simple: fix a clique cover and let the agents in the same clique of the cover know each other. Now, if each agent ignores all feedback coming from agents in other cliques, then the agents in the same clique make exactly the same sequence of prediction and updates. Therefore, the effective number of OMD instances that are being run is equal to χ G .

The following result holds for any differentiable function g : X → R, σ-strongly convex with respect to some norm • . Theorem 6. Consider a network G = (V, E) of N agents, a clique cover {K 1 , . . . , K M } where M = χ G , and let K(v) be the unique element of the cover which each v ∈ V belongs to. For any sequence v 1 , v 2 , . . . ∈ V of active agents, assume each agent v ∈ V runs OMD using g t = √ t η g (with η > 0) while making updates only at those time steps t such that v t ∈ K(v). Then the network regret satisfies

E[R T ] ≤ D η + ηL 2 2σ χ G T
where D ≥ sup g, L ≥ sup ∇ℓ t * , and

• * is the dual norm of • . In particular, choosing η = √ 2σD/L gives E[R T ] ≤ L 2Dχ G T /σ.
Proof. Fix any clique K c and any v ∈ K c . Let T c be the time steps such that v t ∈ K c . Since each agent v ∈ K c ignores the feedback coming from other cliques, the nodes in K c perform exactly the same updates, and therefore make exactly the same predictions. This means that, for any t ∈ T c , the predictions in the set x t (v) | v ∈ K c are all equal to the same common value denoted by x t (K c ). Fix any x ∈ X and, for any t ∈ T c , let r t (K c ) = ℓ t x t (K c )ℓ t (x). By Theorem 1 we have that

t∈Tc r t (K c ) ≤ D η + ηL 2 2σ T c .
Therefore, recalling that r t (v t ) = ℓ t x t (v t )ℓ t (x) and using Jensen's inequality,

T t=1 r t (v t ) = χ G c=1 t∈Tc r t (K c ) ≤ χ G c=1 D η + ηL 2 2σ T c ≤ D η + ηL 2 2σ χ G T
concluding the proof.

Theorems 5 and 6 show that with adversarial activations the knowledge of the graph is crucial for learning (e.g., for achieving sublinear regret). Since χ G ≥ α G , it is not clear whether the better rate √ α G T can be proven in the adversarial activation setting when agents do not use the oblivious network interface.

Conclusions

We introduced a cooperative learning setting in which agents, sitting on the nodes of a communication network, run instances of an online learning algorithm with the common goal of minimizing their regret. In order to investigate how the knowledge of the graph topology affects regret in cooperative online learning under different activation mechanisms, we introduced the notion of oblivious network interface. This prevents agents from doing any network-specific tuning or even accessing their neighborhood structure. When activations are stochastic, we showed that sharing losses among neighbors is enough to guarantee optimal regret rates even with the oblivious network interface. Surprisingly, when activations are adversarial the situation changes completely. There exist problem instances in which any algorithm that runs with the oblivious network interface suffers linear regret. In this case knowing graph structure is not only necessary to perform optimally, but even to have sublinear regret.

Other interesting variants of this settings could be studied in the future. For example, at the beginning of each round, active agents could be allowed to ask the predictions of some of their neighbors, and base their prediction upon it. In this case, we conjecture that the optimal regret rate would scale with the dominating number δ G of the graph, which is always smaller or equal to the independence number. one agent is active at each time step, now v∈V q v = 1 in general. Before the main result, we give some definitions and prove a technical combinatorial lemma that is leveraged in the analysis.

Denote by V ′ the set of all agents v ∈ V such that q v > 0. For each v ∈ V ′ , let

c v = S⊂{1,...,N }\{v} λ S,v 1 + |S| (4)
where the convex coefficients λ S,v are defined by

N w=1 q w   u∈{1,...,N }\({v}∪S) (1 -q u )  
Let also Q v be the probability

P v ∈ w∈St N w = 1 - w∈Nv 1 -q w > 0 (5) that agent v is updated at time t -note that Q v is independent of t.
Lemma 3. Let X(1), . . . , X(m) be independent Bernoulli random variables with strictly positive parameters q 1 , . . . , q m respectively. Then, for all v ∈ {1, . . . , m},

E X(v) m w=1 X(w) = q v c v
where we define X(v)/ m w=1 X(w) = 0 when X(v) = 0. Proof. Fix any v ∈ {1, . . . , m}. Let S v be the set {1, . . . , m} \ {v} and let F v be the σ-algebra generated by

X(w) | w ∈ S v . Then E X(v) m w=1 X(w) = E E X(v) m w=1 X(w) F v = q v E 1 1 + w∈Sv X(w)
Denote the last expectation by c v . Since for all x = 0, ∞ 0 e -tx dt = 1

x , Fubini's theorem yields

c v = ∞ 0 E e -t 1+ w∈Sv X(w) dt = ∞ 0 e -t w∈Sv E e -tXt(w) dt = ∞ 0 e -t w∈Sv
q w e -t + 1q w dt = 1 0 w∈Sv (q w x + 1q w )dx

= 1 0 S⊂Sv x |S| w∈S q w   u∈Sv\S (1 -q u )   dx
Now set λ S,v = w∈S q w Sv \S (1q u ) and note that S⊂Sv λ S,v = w∈Sv (q w + 1q w ) = 1. Substituting λ S,v in the last identity gives

c v = S⊂Sv λ S,v 1 0 x |S| dx = S⊂Sv λ S,v 1 + |S|
We now give an upper bound on the regret that the network incurs if all agents run OMD with an oblivious network interface. Our upper bound is expressed in terms of a constant depending on the probabilities of activating each agent and such that Q ≤ 1.6(α G + 1). The result holds for any differentiable function g : X → R, σ-strongly convex with respect to some norm • . Theorem 7. Consider a network G = (V, E) of N agents. Assume that, at each time step t each agent v is independently activated with probability q v ∈ [0, 1]. If all agents run OMD with an oblivious network interface and using g t = √ t η g, for η > 0, the network regret satisfies

E[R T ] ≤ D η + ηL 2 2σ QT where Q = v∈V ′ (q v c v )/Q v ,
E R T (x) = E T t=1 1 w∈V X t (w) v∈V ′ r t (v)X t (v) = T t=1 v∈V ′ E X t (v) w∈V X t (w) E r t (v) = v∈V ′ q v c v T t=1 E r t (v) (7)
and the last identity follows by Lemma 3. Putting identity [START_REF] Hazan | Introduction to online convex optimization[END_REF] and inequality [START_REF] Griggs | Lower bounds on the independence number in terms of the degrees[END_REF] together gives

E R T (x) ≤ v∈V ′ q v c v 1 Q v D η + ηL 2 2σ √ T ≤ v∈V ′ q v c v Q v D η + ηL 2 2σ √ T
where in the last inequality we used Jensen inequality and v∈V ′ q v c v ≤ 1. This concludes the proof.

We now prove that the inequality Q(q) ≤ α G is always true up to a small constant factor.

Lemma 4. Let G = (V, E) be an undirected graph. For all v ∈ V , choose numbers q v ∈ (0, 1] and define c v and Q v as in (4) and (5) respectively. Then

Q = v∈V q v c v Q v ≤ α G + 1 1 -e -1
Proof. Let P v = w∈Nv q w , V 1 = v ∈ V | P v ≥ 1 , and V 0 = v ∈ V | P v < 1 . We begin by splitting the sum as follows

v∈V q v c v Q v = v∈V1 q v c v Q v + v∈V0 q v c v Q v
We upper bound the two terms separately. Since the minimum min v∈V1 Q v is attained when q v = 1/|N v | for all v ∈ N v , we can lower bound, for each v ∈ V 1 ,

Q v ≥ 1 -1 - 1 |N v | |Nv | ≥ 1 -e -1
This together with v∈V q v c v ≤ 1 yields

v∈V1 q v c v Q v ≤ 1 1 -e -1
To upper bound the sum over V 0 , we first use the inequality 1x ≤ e -x that holds for all x ∈ [0, 1]. Setting x = q w gives

Q v ≥ 1 -exp - w∈Nv q w = 1 -e -Pv
For all v ∈ V 0 , we can then use the inequality 1e -x ≥ (1e -1 )x, holding for all x ∈ [0, 1]. Setting x = P v < 1 we conclude that Q v ≥ (1e -1 )P v for all v ∈ V 0 . Finally, using c v ≤ 1 we can write

v∈V0 c v q v Q v ≤ 1 1 -e -1 v∈V q v P v ≤ α G 1 -e -1
where the last inequality follows by Lemma 2. Putting everything together gives the result.

.

  D ≥ sup g, L ≥ sup ∇ℓ t * , and• * is the dual norm of • . In particular, choosing η = √ 2σD/L gives E[R T ] ≤ L 2DQT /σ.Proof. Fixing an arbitrary x ∈ X , setting r t (v) = ℓ t x t (v)ℓ t (x), and proceeding as in Theorem 2 yields, for each v ∈ V ′ ,

	E	T t=1	r t (v) ≤	D η	+	ηL 2 2σ	T Q v	(6)
	Now we write E[R T ] = sup

x∈X E R T (x) , where

Two adjacent agents that are simultaneously active exchange their feedback, but this does not bring any new information to either agent because we are in a full information setting and the loss function is the same for all nodes.

Acknowledgments

Nicolò Cesa-Bianchi and Tommaso Cesari gratefully acknowledge partial support by the GAMENET COST Action (CA 16228) and by the MIUR PRIN grant Algorithms, Games, and Digital Markets (ALGADIMAR). Nicolò Cesa-Bianchi is also partially supported by the Google Focused Award Algorithms and Learning for AI (ALL4AI).

A Stochastic Activations: Multiple Agents

In this section we present all missing results related to the stochastic activation model with multiple activations per time step. Recall that, at the beginning of the process, the environment draws an i.i.d. sequence of Bernoulli random variables X 1 (v), X 2 (v), . . . with some unknown fixed parameter q v ∈ [0, 1] for each agent v ∈ V . The active set at time t is then defined as S t = {v ∈ V | X t (v) = 1}. Note that, unlike when only


Solenne Willemin

Paul Carminati

Sylvain Jacques

Jérôme Roger

Francis Rebillat

Identification of complex oxidation/corrosion behaviours of boron nitride under high temperature

Keywords: boron nitride, oxidation, thermo-gravimetric analysis, corrosion

To characterize the oxidation/corrosion behaviour of boron nitride (BN) at high temperature and to determine kinetics, thermo-gravimetric analyses were performed on BN flat coatings under dry (80% N 2 , 20% O 2 ) and wet conditions (10% H 2 O, 70% N 2 , 20% O 2 ), at temperatures ranging from 700 to 900 °C. The surface morphology was analysed by scanning electron microscopy. The aim of the study was to understand the complex behaviour of BN at high temperature under wet air. A model was then proposed to explain the recorded mass variations, by taking into account a competition between absorption and volatilization phenomena. Finally, an estimation of the amount of adsorbed species was extracted from these data.

Introduction

Ceramics and composites based on silicon carbide (SiC) are interesting materials for structural components in high-temperature applications. Indeed, compared with metallic alloys historically used in the aeronautic industry, ceramic matrix composites (CMCs) are lighter (enabling a reduction in energy consumption) and stiffer. Their ability to resist to severe environmental conditions such as high temperatures and oxidizing/corrosive environments makes them key candidates for aeronautic industry. Boron nitride is a promising interphase material for SiC fibre-reinforced SiC matrix composites. Its layered structure allows matrix cracks deviation from mode I to mode II, protecting the fibres from premature failure. It has a better oxidation resistance than pyrocarbon used until now [START_REF] Gallet | Oxidation Resistance of BN Coatings with Different Textures[END_REF]. Furthermore, after oxidation, boron and its compounds give a B 2 O 3 glassy phase with a low melting point of 450 °C. Boron nitride then turns out to be an efficient protective layer for fibres. However, this material is known to be very sensitive to moisture because boron oxide volatilizes, even at room temperature: when poorly crystallized, adsorption of moisture may occur in significant amount [START_REF] Gallet | Oxidation Resistance of BN Coatings with Different Textures[END_REF]. The aim of this study is to understand the complex behaviour of boron nitride at high temperature under wet air.

Experimental procedure

Sample Preparation

Boron nitride coatings were prepared by chemical vapour deposition under chemical vapour deposition/chemical vapour infiltration (CVD/CVI) conditions (low pressure). In order to remove remaining partially decomposed reactive species trapped during the CVD deposition process, thermal post-treatments at high temperature were performed. The coatings were deposited on silicon carbide (SiC) rectangular substrates. For the studied coatings (Table 1), sample dimensions were approximately 5 mm × 6 mm, with an average mass of 180 mg and a CVD boron nitride thickness of about 15 µm.

Sample Characterization

Thermo-gravimetric analyses (TGA) were performed using Setsys 1600 from Setaram. The specimens were placed in dry (80 kPa N 2 , 20 kPa O 2 ) and wet conditions (10 kPa H 2 O, 70 kPa N 2 , 20 kPa O 2 ), at temperatures ranging from 700 to 900 °C. Dry air was generated by mixing N 2 and O 2 . To obtain wet conditions, a dry gas mixture is saturated with water vapour by bubbling through a wash bottle, maintained in a thermally controlled bath. The temperature of the bath is fixed to a value corresponding to the required P H 2 O(g) . In every performed test, O 2 partial pressure was maintained constant to 20 kPa.

Elsewhere, surface morphology of oxidized coatings was characterized using scanning electron microscopy (SEM) FEI 400 FEG.

The crystallographic parameters of BN were characterized using X-ray diffraction (XRD, Bruker D8 Advance apparatus). The followed parameters were d 002 (shortest distance between two planes in the 002 orientation), L c (crystallite sizes, inversely proportional to the full width half maximum of the characteristic XRD peak) and the structural inhomogeneity (linked to the asymmetry of the XRD peak at the small angle side). Indeed, the coatings were made of a mixture of more or less crystallized BN (strong asymmetry of XRD peak). The crystallographic parameters of the studied deposits, before oxidation, are presented in Table 1. The factor ''peak asymmetry'' is calculated as the ratio between the right full width half maximum and the left full width half maximum for 2Ɵ < 2Ɵ average of the XRD peak. A value of 1 then indicates a symmetrical peak, i.e., a homogeneous coating.

Results and discussion

Behaviour Under Dry Air

In the absence of moisture, boron nitride oxidation occurs according to the following reaction:

BN + O → B O + N (1)
At temperatures below 1000 °C, B 2 O 3 volatilization is negligible because of its low saturated vapour pressure. Isothermal TGA analyses, therefore, only indicate a mass gain, induced by B 2 O 3 formation on the flat coating surface. Oxidation kinetics follows a general relationship similar to the one described earlier by Deal and Grove [START_REF] Deal | [END_REF] for the thermal oxidation of silicon. As long as the oxide thickness is thin, the relation describing the scale growth can be simplified by a linear law: oxidation kinetics is limited by reactions at the interface (limiting step). For longer duration, oxidation is later controlled by oxygen diffusion through the oxide thickness, and the kinetic rate follows a parabolic law. Two different oxidation rate constants can, therefore, be determined: the linear rate constant k i and the parabolic rate constant k p .

The TGA weight evolution of BN coatings aged at temperatures of 800, 850 and 900 °C in dry conditions (Fig. 1) clearly follows the previous description. The values of rate constants extracted from experimental curves are gathered in Table 2. The weight evolution of the samples from TGA can be converted either to BN or B 2 O 3 thickness evolution (consumption or formation, respectively). In this study, the linear rate constant k i was determined from BN thickness evolution in function of time during the very first hours of the tests (linear evolution due to an interfacial regime). The parabolic rate constant k p was determined by plotting the square of B 2 O 3 thickness evolution in function of time (parabolic evolution due a regime limited by diffusion).

As expected, the oxidation kinetics is thermo-activated. Also, Arrhenius plots (ln(k) versus 1/T) allows values of overall activation energies to be determined: 126 and 63.4 kJ/mol for linear and parabolic regimes, respectively.

Behaviour Under Wet Air

In wet air, oxidation kinetics is significantly modified: as soon as small quantities of water vapour are introduced, volatilization of B 2 O 3 can no more be neglected. Meschi et al. [4] have observed, using mass spectrometry, the presence of H 3 BO 3 , HBO 2 and (HBO 2 ) 3 gaseous species (for T ranging from about 800-1200 °C and P H 2 O ranging from 0.05 to 0.14 mmHg). From a thermodynamic approach of gas/ liquid equilibria, at temperatures above around 800-900 °C, HBO 2 is considered as the major species, whereas below this temperature the major species are mostly (HBO 2 ) 3 and H 3 BO 3 . In wet air, a competition occurs between an oxide formation on the surface of the deposit, and its volatilization into H x B y O z species [5,6]. This leads to a paralinear behaviour, which is a mixture of simultaneous parabolic and linear kinetics. It can be schematized as presented in Fig. 2. This theoretical model detailed by Tedmon [7] for Cr and Fe-Cr alloys and by Opila and Hann [8] for silicon carbide oxidation, has already been considered to describe the corrosion behaviour of CVD boron nitride [9].

Figure 3 shows the weight changes measured at temperatures of 700, 750 and 800 °C in wet air (10% H 2 O) for BN coatings. As expected, recession of the boron nitride thickness and thus oxidation/corrosion kinetics increase with temperature. However, unlikely to the previously considered model [7], the curves exhibit: (1) a critical mass gain without the formation of an oxide scale and (2) two linear parts instead of only one during the stable recession step. A model explaining the differences between the observed and expected theoretical behaviours is developed below.

Model for Wet Air

As boron nitride is known to be moisture sensitive even at room temperature, a reasonable assumption explaining this behaviour would be an absorption process in the deposit under wet air. As shown in Fig. 3, and schematized in Fig. 4, a weight gain is clearly noticeable at the very beginning of oxidation while no oxide scale is formed over the surface (corresponding to stage 0 in Fig. 4). This could correspond to oxidizing/corrosive species absorption into BN coatings. Then, a continuous mass loss follows this weight gain. The mass loss shows two successive linear domains, corresponding to stages 1 and 2 in Fig. 4, respectively. The first linear regime is expected to result in the competition between the OH species progression in the boron nitride, and the volatilization of H x B y O z species. The second one indicates a faster mass loss. It can reasonably be assumed that: (1) the BN coating is saturated with H x B y O z adsorbed species and (2) the contribution of the weight gain due to absorption of species is no longer present. Thus, this latter regime corresponds to the only volatilization of gaseous species until complete consumption of the coating: the resulting degradation process could be resumed as the active oxidation of hydrated BN molecules, ''BNOH''. This model seems to properly explain the differences encountered with the theoretical paralinear model. Such a kind of behaviour has already been observed for oxidation in air of poorly crystallized carbon (phenolic resin-based carbon): similar layered crystalline structure, sensitive to adsorption and/or intercalation phenomena [10,11].

To highlight this assumed absorption phenomenon during an oxidation test, a reverse phenomenon of desorption is momentarily created through the introduction of a step under argon at the working temperature. Thus, oxidation tests were interrupted for a short period of time at constant temperature under neutral gas (argon), during the second linear regime (Fig. 5). After this modification of the ''reactive'' surrounding gas, a mass loss identified as desorption is clearly visible (once a correction of mass fluctuations is subtracted, due to the change in nature of gases). Moreover, when oxidizing gases (zone 3 in Fig. 5) is reintroduced, a mass gain due to absorption occurs again but at a lower level. The duration of the first linear regime is also reduced. This supports the hypothesis of absorption in the deposit: as BN coating becomes thinner, the volume of the sample that can be hydrated is lower.

The chemical reactions involved during the two linear domains are shown in Fig. 6. The differences between the two linear rates are mainly due to the degree of absorption: (1) progression and (2) saturation. It can be assumed that the same molecular amount of boron nitride is consumed during regime 1 and regime 2: the BN consumption in mol/h in domain 1 (k BN1 ) equals that in mol/h in domain 2 (k BN2 ). The difference of slopes between both linear regimes is proportional to the difference of mass balances between both involved overall reactions. Thus, the difference of mass balance can directly be related to the average molecular mass of adsorbed species by one mole of BN.

In a first step, from the slope of the TGA curves, a mass consumption rate, k 1 in mg/h, is determined (Table 3). Then, a molecular consumption rate, k BN1 , which corresponds to the only BN consumption, is deduced. As the mass weight rate k 2 measured by TGA is proportional to k BN1 and to M BNOH , it is then possible to extract the molecular mass of hydrated BN molecules, BNH x O z . This average quantity of absorbed species, H x O z , is thus estimated, reported to one mole of BN (It has been checked that after desorption, the coating is only made of BN; no premature removal of nitrogen occurs during an oxidation/corrosion test).

The gaseous hydroxides responsible for the volatilization of ''B 2 O 3 '' in moisture environments (H 3 BO 3 , HBO 2 and (HBO 2 ) 3 ) indicate that the absorption capability around a boron atom is expected to be important (several chemical bonds between B and O or H atoms in each specie). Further, the layered crystalline structure of this material favours the insertion of gaseous species in the structure. Consequently, the average calculated quantities of H x O z absorbed species for one mole of BN, from TGA data seem to be conceivable. A high capability of absorption could be expected around 800-900 °C in regard to the H x B y O z gaseous species responsible to the volatilization of ''B 2 O 3 '', characterized by a higher M H x O z /M B y mass ratio [5,6].

To highlight the phenomenon of hydrated species absorption in BN coatings, XPS analyses have been conducted on as-received and oxidized samples (quenched to keep their high-temperature state). The results, presented in Fig. 7, show that after an oxidation test, the amount of oxygen trapped inside the coating is larger than the initial quantity: close to the surface, the oxygen amount presents a maximum of 35 at.%, which decreases to the stable value of 27 at.% inside the deposit. Moreover, deconvolutions of B1s peaks indicate an initial B-O bonding percentage of 8 at.%, which reaches 23 at.% after oxidation. This clearly supports the presence of hydrated BN molecules inside the coating. Theses analyses will be further completed by Auger electron spectroscopy, in order to estimate the amount of oxygen across the whole thickness of the boron nitride coating.

SEM observations of BN surface morphologies after oxidation/corrosion have been carried out. As shown in Fig. 8, cracks are observed over poorly crystallized boron nitride coatings. This cracking is expected to occur during the desorption process, due to volume variation and/or to the elimination of trapped gases. From further experiments, a coating with a higher crystallization degree shows a lower weight gain due to an absorption/intercalation process and presents no cracking after testing, respectively.

Conclusions

In this study, the oxidation/corrosion behaviour of CVD boron nitride flat coatings has been investigated. In wet air, no theoretical behaviour allows to describe accurately the oxidation/corrosion behaviour of BN coatings. Starting from a paralinear behaviour, a modelling is proposed by taking into account an absorption phenomenon of H x O z species through the boron nitride thickness. The estimated quantities of absorbed species are consistent with the nature of H x B y O z hydroxides met at such level of temperatures in moist environments.

The next step in this work will be to identify the absorbed elements in the various expected reactions and to deduce the different kinds of chemical bonding involved between them deep inside BN deposits. By quenching the system during a corrosion test, it can be expected to keep trapped a major part of absorbed species. Then, the identification of these trapped species can be carried out by X-ray photoelectron spectroscopy, Auger electron spectroscopy and glow discharge optical emission spectroscopy analyses. Another undertaken possibility is an in situ identification by mass spectrometry of gaseous species removed during the desorption step.

Tables

FiguresFig. 1

 1 Figures

Fig. 2

 2 Fig. 1 Isothermal weight gain for BN in dry air at 800, 850 and 900 °C obtained from TGA experiments

Fig. 3 Fig. 4

 34 Fig. 3 Isothermal weight loss for BN in wet air (10 kPa H 2 O) at 700, 750 and 800 °C

Fig. 5 Fig. 6

 56 Fig. 5 Isothermal weight loss for BN at 800 °C in wet air (10 kPa H 2 O). During the oxidation tests (1 and 3) two interruptions under inert gas (2 and 4) were introduced, at constant temperature

Fig. 7 XPSFig. 8

 78 Fig. 7 XPS Analyses of (a) BN before wet air (10 kPa H 2 O/20 kPa O 2 /70 kPa N 2 ) exposure and (b) BN after wet air (10 kPa H 2 O/20 kPa O 2 /70 kPa N 2 ) exposure at 800 °C

Oxidation temperature (°C) k 1 (mg h -1 cm -2 ) k 2 (mg h -1 cm -2 ) k BN1 (mol BN /h) M BNOH (g/mol) M OH (g/mol BN )

  Estimation of the amount of absorbed species in the BN coating

		Sample		BN
		L c (nm)		2.1 (±0.1)
		d 002 (nm)		0.343 (±0.001)
		Peak asymmetry factor	0.6
	Tab. 1 Crystallographic parameters of the studied CVD BN coating
	Oxidation temperature (°C)	k i (µm BN /h)	k p (µm B 2 O 3	2 /h)
	800			1.8	190
	850			3.5	224
	900			6.0	350
	Tab. 2 Rate constants for oxidation of boron nitride in dry air
	700	-0.37	-0.99	-1.5 × 10 -5	66.2	41.4
	750	-0.71	-1.96	-2.9 × 10 -5	68.0	43.2
	800	-1.30	-2.89	-5.2 × 10 -5	55.2	30.3
	Tab. 3

Acknowledgements

This work has been supported by Herakles, through a PhD grant given to S. Willemin.


Cédric Batifol
email: cedric.batifol@ec-lyon.fr

Mohamed N Ichchou
email: mohamed.ichchou@ec-lyon.fr

Marie-Annick Galland
email: marie-annick.galland@ec-lyon.fr

Hybrid modal reduction for poroelastic materials

Keywords: Poroelastic materials, Dynamic substructuring, Modal synthesis Matériaux poroélastiques, Sous-structuration dynamique, Synthèse modale

Introduction

Finite element modelling of poroelastic materials is still a complicated and time-consuming task. Indeed, the classical biphasic Biot model involves 6 degrees of freedom (dofs) for the (u s , U f ) formulation [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid[END_REF] or 4 dofs for the (u s , p) formulation [START_REF] Atalla | A mixed displacement-pressure formulation for poroelastic materials[END_REF]. Moreover, materials' parameters are complex-valued and frequency-dependent. This prevents from using such formulations in optimisation procedures. Reduced models are needed, and some have already been proposed. Recently, Jaouen [START_REF] Jaouen | A simplified numerical model for a plate backed by a thin foam layer in the low frequency range[END_REF] developed a plate-like model and Dazel [START_REF] Dazel | An extension of complex modes for the resolution of finite-element poroelastic problems[END_REF] extended the complex modes theory to the poroelastic case. Approximations still remain in these models and the accuracy for a wide range of materials and configurations is not obvious.

In this Note, a hybrid clamped-free (C-F) component mode synthesis (CMS) method is developed. Coupled dynamic modes are calculated for simultaneously fixed fluid dofs and free solid dofs at the junction boundary. Using such boundary conditions allows to shift closer both fluid and solid phases resonances. This assumption is widely used in more classical fluid/structure interaction problems, where the fluid (mainly air) is softer than the solid part [START_REF] Ohayon | Reduced symmetric models for modal analysis of internal structural-acoustic and hydroelastic-shloshing problems[END_REF]. In this contribution, eigenmodes are calculated, using the (u s , U f ) formulation, solving the frequency-dependent eigenproblem without simplifying it. Then, the projection subspace is spanned by a dynamic basis and fluid static boundary functions. Firstly, the complete derivation of the projection basis is detailed. Then, the last part of this paper presents a convergence study where results are compared to classical Craig and Bampton (C-B) and MacNeal (M-N) approaches, showing the effectiveness of the method.

The Clamped-Free poroelastic eigenproblem

The configuration studied is depicted on Fig. 1. A poroelastic domain Ω is subjected to harmonic forcing on its Neumann boundary ∂Ω N . Fixing conditions involve linear constraints on its Dirichlet boundary ∂Ω D . Finally, this sample is coupled to other layers on its junction boundary ∂Ω J . The whole boundary ∂Ω admits the direct summation of Eq. [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid[END_REF].

∂Ω = ∂Ω D ⊕ ∂Ω N ⊕ ∂Ω J (1)
The associated eigenproblem, using the solid displacement u s and the fluid displacement U f as kinematic unknowns, with hybrid (C-F) conditions on ∂Ω J , states as below (Eqs. ( 2) and ( 3)), with boundary equations ( 4)- [START_REF] Allard | Propagation of Sound in Porous Media[END_REF]. Find ω and (u s , U f ) = (0, 0) solution of

∇ • T s = -ω 2 ρ11 u s + ρ12 U f (2) ∇ • T f = -ω 2 ρ22 U f + ρ12 u s (3) u s = u fix and U f • n = u fix • n on ∂Ω D (4)
T s • n = 0 and T f • n = 0 on ∂Ω N (5)
T s • n = 0 and U f = 0 on ∂Ω J [START_REF] Allard | Propagation of Sound in Porous Media[END_REF] where n, T s , T f and E s stands for the unit normal vector, the solid stress tensor, the fluid stress tensor and the solid Green tensor, respectively. u fix stands for the displacement vector imposed by the fixing conditions. The stress-strain relations are described in Eqs. [START_REF] Ohayon | Dynamic substructuring of damped structures using singular value decomposition[END_REF] and [START_REF] Voss | An Arnoldi method for non-linear eigenvalue problems[END_REF] as follows.

T s = λs ∇ • u s I + 2μ s E s + λsf ∇ • U f I (7)
T f = λf ∇ • U f I + λsf ∇ • u s I ( 8
)
where I is the identity matrix. In above equations, λs ans μ s stand for the two Lamé-like coefficients of the solid phase. λf and λsf are the fluid phase and coupling Lamé-like coefficients. Finally, ρ11 , ρ22 and ρ12 are the solid, fluid and coupling relative densities. All these coefficients, except μ s , are complex-valued and frequency-dependent parameters. For the sake of clarity, the derivation of all these parameters is not recalled. The reader is referred to [START_REF] Allard | Propagation of Sound in Porous Media[END_REF] for additional informations. Let C u s (Ω) and C U f (Ω) be the admissible spaces for the solid and fluid displacements, respectively. Multiplying Eqs. ( 2) and (3) by admissible displacement variations δu s and δU f , applying the Green formula, the stress-strain relations, and considering boundary conditions lead to the variational formulation of the eigenproblem. 9) and ( 10)

Find (u s , U f , ω) ∈ C u s (Ω) × C U f (Ω) × C satisfying Eqs. (
∀(δu s , δU f ) ∈ C u s (Ω) × C U f (Ω).
λs

Ω ∇ • u s ∇ • δu s dΩ + 2μ s Ω E s : δE s dΩ + λsf Ω ∇ • U f ∇ • δu s dΩ + ω 2 ρ11 Ω u s • δu s dΩ + ω 2 ρ12 Ω U f • δu s dΩ = 0 ( 9 ) -λf Ω ∇ • U f ∇ • δU f dΩ + λsf Ω ∇ • u s ∇ • δU f dΩ + ω 2 ρ22 Ω U f • δU f dΩ + ω 2 ρ12 Ω u s • δU f dΩ = 0 ( 10
)
These formulations are discretized using the finite element method. The numerical eigenproblem is then written as in Eq. [START_REF]A hybrid method of component mode synthesis[END_REF].

λs K s int,1 + 2μ s K s int,2 -ω 2 ρ11 M s int,1 λsf K s int,c -ω 2 ρ12 M s int,c λsf K f int,c -ω 2 ρ12 M f int,c λf K f int,1 -ω 2 ρ22 M f int,1 u s U f = 0 0 (11)
where u s and U f are vectors of solid and fluid dofs. The above dimensionless matrices are discretized forms of domain integrals. As an example, the meaning of K s int,1 is detailed in Eq. ( 12).

Ω ∇ • u s ∇ • δu s dΩ ⇒ δu sT K s int,1 u s (12)
where . T is the transpose operator.

Finally, linear constraints are applied thanks to Lagrange multipliers. An incompatibles meshes linking procedure [START_REF] Ohayon | Dynamic substructuring of damped structures using singular value decomposition[END_REF] can be applied.

The non-linear Arnoldi eigensolver

The non-linear Arnoldi eigensolver is an iterative solver for sparse problems. The version used for solving the poroelastic eigenproblem is adapted from the one originally proposed by Voss [START_REF] Voss | An Arnoldi method for non-linear eigenvalue problems[END_REF]. It relies on three main steps. First, approximated eigenvectors and eigenvalues are calculated. It provides accurate first vectors forming a projection basis. Then the non-linear projected problem is solved using a dense solver. Finally, the error between calculated and exact eigenvector is estimated. If convergency is not reached, a new relevant vector is added to the search space.

Approximated problem

For the poroelastic case, the approximated eigenpairs are calculated by solving the frequency-independent eigenproblem (13), with the constant approximated poroelastic parameters expressed in Eqs. ( 14), (15).

λ s0 K s int,1 + 2μ s K s int,2 -ω 2 ρ 11 M s int,1 λ sf 0 K s int,c -ω 2 ρ 12 M s int,c λ sf 0 K f int,c -ω 2 ρ 12 M f int,c λ f 0 K f int,1 -ω 2 ρ 22 M f int,1 u s U f = 0 0 ( 13
)
ρ 12 = -φρ 0 (α ∞ -1), ρ 11 = (1 -φ)ρ s -ρ 12 , ρ 22 = φρ 0 -ρ 12 (14) λ f 0 = φP 0 , λ sf 0 = (1 -φ)P 0 , λ s 0 = (1 -φ) 2 φ P 0 + 2μ s ν 1 -2ν (15)
A few characteristic coefficients have to be introduced. ρ 0 , ρ s , φ, α ∞ , P 0 , ν are the fluid density, the density of the material of the skeleton, the porosity, the tortuosity, the static pressure, and the Poisson ratio, respectively. Problem (13) involves symmetric matrices and is solved using standard sparse eigensolvers. The calculated eigenvalues and eigenvectors are denoted τ j and B j .

Resolution of the projected problem

The eigenpair (ω j , V j ) now needs to be searched for. Denoting M(ω) the matrix involved in problem [START_REF]A hybrid method of component mode synthesis[END_REF], a first order development is written as in Eq. ( 16).

M(ω j ) ≈ M(τ j ) -θ M (τ j ) (16)
The method of successive linear problems introduced by Ruhe [START_REF] Ruhe | Rational Krylov, a practical algorithm for large sparse nonsymmetric matrix pencils[END_REF] allows to search for an eigenpair near a pole ξ j .

Obviously, the relevance of such a search depends on the chosen pole. For the poroelastic problem, one may choose ξ j = 0.9τ j . If the pole is too close to the exact eigenvalue, matrices could be ill-conditioned. This solver being not directly applicable to large sparse matrices, a projection of problem [START_REF]A hybrid method of component mode synthesis[END_REF] is to be solved. A projection basis B is spanned by appropriate vectors. The first of them is the approximated eigenvector B j . The resulting linear eigenproblem is presented in Eq. ( 17). The first approximation of the seek eigenvalue is denoted r 0 and is equal to τ j . At the lth-step, the retained value of θ is such that ξ j -(r lθ) is minimal. The (l + 1)th eigenvalue approximation, r l+1 , is then equal to r lθ . Finally, the algorithm stops when the convergence criterion (18) is reached.

B H M(r l )BY = θB H M (r l )BY, with B H B = I (17) ξ j -(r l+1 -θ) ξ j -(r l -θ) -1 < (18)
. H is the adjoint operator. is a small scalar parameter.

Error estimate and basis updating

The preceding step gives the estimated eigenpair ( ωj , Ûj ), with Ûj = BY. The error between estimated and exact eigenpair is expressed via the calculation of the residue R expressed in Eq. ( 19).

R = M( ωj ) Ûj (19)
Convergence is reached if the norm of R is sufficiently smaller than the norm of Ûj . If this is not the case, the basis B has to be expanded with the new vector b n satisfying conditions (20).

M(ξ j )b n = R, B H b n = 0, b n = 1 (20)

Construction of the hybrid (C-F) reduced poroelastic model

We consider now the response of a poroelastic sample subjected to non-zero prescribed loads and/or displacements. The admissible search space for such forced problem is denoted C for u s × C for U f . The hybrid (C-F) reduced poroelastic model relies on the direct decomposition (21).

C for u s × C for U f = 0 × C stat U f ⊕ C dyn u s ,U f (21)
where C stat U f is the space of fluid static boundary functions, and C dyn u s ,U f is the dynamic (C-F) space. For the discrete case, this dynamic (C-F) space is spanned by a dynamic basis P d calculated thanks to dynamic modes. Indeed, the dynamic eigenvectors, calculated using the non-linear Arnoldi solver, do not form a unitary matrix. Due to the frequency-dependence of the eigenproblem, eigenvectors have to be orthogonalised using the Gramm-Schmidt process. Calculated eigenvectors Ûj are first classified in ascending order in regards to eigenvalues. At iterate j , the vector Φ j is then calculated thanks to the procedure described in Eq. ( 22). The dynamic basis P d , formed by vectors Φ j , is then a modal-like basis.

D j = Ûj -P d j -1 P dH j -1 Ûj , Φ j = D j D H j D j (22)
Next, the fluid static functions matrix, P stat , is formed by static fluid functions calculated for each junction dofs. As an example, the lth static fluid function is obtained by solving the discretized linear problem expressed by Eqs. ( 23) and (24).

λ s 0 K s int,1 + 2μ s K s int,2 λ sf 0 K s int,c λ sf 0 K f int,c λ f 0 K f int,1 u s U f = 0 0 (23) U f • n = 0 on ∂Ω J \l, U f l • n = 1 (24)
where ∂Ω J \l is the junction boundary without the lth fluid junction dof U f l . Remaining Dirichlet conditions expressed in Eq. ( 4) have also to be fulfilled. Finally, the reduced model is written as in Eqs. (25).

⎡ ⎣ u s U f |∂Ω J U f |Ω\∂Ω J ⎤ ⎦ ≈ M C-F Z = ⎡ ⎣ 0 P d s I |∂Ω J 0 P stat f |Ω\∂Ω J P d f |Ω\∂Ω J ⎤ ⎦ U f |∂Ω J Z α (25)
where Z α is the vector of generalised dofs, P stat f |Ω\∂Ω J is the fluid static functions matrix restricted to fluid dofs belonging to Ω \ ∂Ω J ,

P d f |Ω\∂Ω J
is the dynamic (C-F) matrix restricted to fluid dofs belonging to Ω \ ∂Ω J , P d s is the dynamic (C-F) matrix restricted to solid dofs. The relevance of the projected model depends obviously on the number of dynamic modes retained.

Convergence study for a 1D sample

The test case involves a one-dimensional felt sample, 12 centimetres long. This poroelastic material is clamped at one side and is subjected to an harmonic loading or an imposed displacement on its other side. Reference results are obtained using a standard finite element formulation involving 962 dofs, and quadratic Lagrange shape functions. The hybrid (C-F) projection basis is compared to classical Craig and Bampton (C-B) [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF] and MacNeal (M-N) [START_REF]A hybrid method of component mode synthesis[END_REF] approaches. Errors on calculated fluid and solid displacements are estimated using formula (26).

err = | û -u| | û| ( 26
)
where û and u are reference and projected solid or fluid displacements, respectively. Results presented are calculated for a 2000 Hertz forcing frequency. Eleven dynamic modes are calculated. The projection basis is then spanned by calculated modes, their complex conjugates and relevant static functions. The number of generalised coordinates is equal to 22, 23 and 24 for (M-N), (C-F) and (C-B) projection basis, respectively. Remaining relative errors are listed on Table 1 for both phases and excitation types. The (C-F) basis proves to be the most efficient. Indeed, residual errors are small for both fluid and solid phases, which is not the case for the (C-B) and (M-N) basis.

On Fig. 2 are plotted error for all three basis, when one of the eleventh dynamic modes is not retained. Modes around the fifth one are shown to be the most relevant. Indeed, the fifth eigenfrequency is near 2000 Hertz for all Table 1 Error estimate at 2000 Hertz for an imposed displacement (a) and an imposed pressure (b). basis. For both (C-B) and (M-N) basis, solid and fluid phases have opposite behaviours around this mode. For the (C-B) case, the fluid error continue to diminish whereas the solid one increases. The free interface basis exhibits clearly two different peaks for the two different phases. The fluid displacement converges faster than the solid one due to the use of similar interface boundary conditions for both phases. In contrast, (C-F) errors behave the same way for both phases. Each calculated eigenmode is then relevant for both phases displacements reconstruction. This is not the case for the (C-B) and (M-N) basis.

Convergence study for a 2D sample

The test case involves a 2D sample, 10 centimetres wide and 1 centimetre deep. The material is subjected to a prescribed displacement on its junction boundary ∂Ω J . In the cases of the (C-F) and (M-N) basis, the constraint is corrected by using residual attachment modes. Finally, remaining boundaries are subjected to Dirichlet conditions as written in Eq. ( 27). A sketch of the studied configuration can be found on Fig. 3.

u s • n = U f n = 0 on ∂Ω D1 , u s = 0 on ∂Ω D2 , U f • n = 0 on ∂Ω D2 ( 27
)
The forcing frequency varies from 100 to 2300 Hertz. For each frequency, the relative error between direct and reduced results is calculated following Eq. ( 26

Conclusion

A reduction procedure for poroelastic finite element models has been proposed. The projection subspace is spanned by hybrid clamped-free dynamic modes and static fluid functions. The poroelastic eigenproblem is solved without simplification using an iterative process. The fluid phase is fixed to the junction boundary whereas solid dofs are not constrained. The whole process of building such a basis has been detailed. The convergence properties of the (C-F) basis have been studied on both 1D and 2D test cases. For the 1D test case, the reconstructed response is better converged, using the clamped-free basis, than with Craig and Bampton or MacNeal approaches. Each calculated hybrid mode is relevant to describe both fluid and solid dynamic behaviours. This is not the case for classical substructures' basis. In the 2D test case, the (C-F) basis clearly appears as a good compromise between accuracy and reduction ratio. Both phases errors are less than 2 percents until half the highest eigenfrequency, with almost three times less dofs than in the (C-B) basis. Next steps consist in studying convergence properties for a 3D sample, formulating coupling with others substructures and deriving a simple criterion for choosing modes to be retained.

Fig. 1 .

 1 Fig. 1. Poroelastic domain Ω and its boundary ∂Ω.

Fig. 1 .

 1 Fig. 1. Domaine poroélastique Ω et sa frontière ∂Ω.

Tableau 1

 1 Erreur résiduelle à 2000 Hertz pour un déplacement (a) et une pression (b) imposée.

2 Fig. 2 .

 22 Fig. 2. Error evolution versus non-retained modes at frequency 2000 Hertz for imposed displacement (a) and pressure (b). ! and " stand for (C-F) solid and fluid errors. £ and stand for (C-B) solid and fluid errors. 1 and 2 stand for (M-N) solid and fluid errors.

Fig. 2 .

 2 Fig. 2. Évolution de l'erreur à 2000 Hertz en fonction des modes non-retenus. (a) : déplacement imposé. (b) pression imposée. ! et " : erreurs solide et fluide avec la base (C-F). £ et : erreurs solide et fluide avec la base (C-B). 1 et 2 : erreurs solide et fluide avec la base (M-N).

Fig. 3 .Fig. 3 .

 33 Fig. 3. Sketch of the 2D test case. The shape of the imposed normal displacement is plotted on the left.Fig. 3. Configuration étudiée pour le cas test 2D. L'allure du déplacement imposé est visible à gauche.

Fig. 4 .

 4 Fig. 4. X-displacement relative errors. ! and ": Solid and fluid errors, basis (C-F). £ and : Solid and fluid errors, basis (C-B). 1 and 2: Solid and fluid errors, basis (M-N).

Fig. 4 .

 4 Fig. 4. ! et Erreurs relatives pour le déplacement selon x. " : Erreurs solide et fluide, base (C-F). £ et : Erreurs solide et fluide, base (C-B). 1 et 2 : Erreurs solide et fluide, base (M-N).


S Del Broccolo

R Palumbo1

M Campana1

Y Dobah1
email: y.dobah@bristol.ac.uk

F Scarpa1
email: f.scarpa@bristol.ac.uk

M Ouisse

M Ichchou

Sevilla (Spain) KIRIGAMI INSPIRED NATURAL FIBRE CELLULAR STRUCTURES FOR FUTURE VIBROACOUSTICS APPLICATIONS

The Aerospace segment makes large use of cellular structures in sandwich panels, which are known to possess good stiffness while saving weight. Metallic cores are currently being re placed by dif f erent classes of composites, especially because of their poor vibration transmissi bility. The inclusion of vibroacoustics constraints in the design process gives the opportunity to obtain multifunctional structures that still provide mechanical efficiency while introducing ab sorption, tunable vibration transmissibility or damping capabilities. In this work, natural . fibre prepregs were turned into periodic panels to obtain cellular structures with improved damping properties, using the Japanese art technique of cutting and folding paper, known as Kirigami. Varying the geometric characteristics of the honeycomb, numerical simulations were carried out in order to identify the best performing candidate in terms of the above mentioned proper ties. Results show that hexagonal cell honeycombs have smaller values of transmissibility than the auxetic one.

INTRODUCTION

Periodic structures have been largely investigated for their in-plane vibroacoustics performances, starting from Brillouin initial studies on the subject [START_REF] Brillouin | Wave Propagation in Periodic Structures, Electric Filters and Crystal Lattices[END_REF]. Mead demonstrated [START_REF] Mead | Wave propagation in continuous periodic structures: research contributions from southampton 1964-1995[END_REF] how periodical geometry and/or material discontiuities attenuate the propagation of in-plane waves through the structure, as a result of the interaction between incident, reflected and transmitted waves at the discontinuity points.

In this work, the out-of-plane vibration damping of Kirigami made honeycomb cores for future vibroacoustic purposes has been investigated.

The paper is organized as follows. In section 1 the Kirigamj manufacturing process of the hon eycombs is described. In section 2, the numerical simulations carried out so far are summarized.

Details about the Finite Element model used for the simulations are given, followed by the re sults obtained from the transmissibility analysis. Finally, section 3 gives an outline of the main results of the work and of the future investigations.

KIRIGAMI MANUFACTURING

Kirigami is an ancient Japanese art of cutting and folding paper. The idea of engineering Kirigami was first explored by Nojima and Saito [START_REF] Nojima | Development of newly designed ultra-light core structures[END_REF]; since then, the technique has been widely applied to both composite and thermoplastic materials [START_REF] Neville | Design of shape morphing cellular structures and their ac tuation methods[END_REF] to form a 2D sheet material into a 3D cellular structure.

In this work, the material used to manufacture the honeycomb is natural-fibre prepreg -flax fibres with polypropylene matrix -characterized by 2x2 Twill weaving. Both hexagonal and auxetic cell honeycomb were produced. In order to obtain the final honeycomb 3D shape, the Kirigami manufacturing process consisted of three main steps:

1. Corrugating+ Thermoforming: the sheet material was thermoformed between hexag onal tooling in a hot press for 20 minutes, at a temperature of 170° C and a pressure of 10 GPa. The hexagonal and auxetic shapes were here obtained using the same mould, but folding the sheet on the hexagonal tooling in different ways.

2. Cutting: a pattern of slits was laser cut into the corrugated material. The result for the hexagonal cell honeycomb is shown in Figure 1 (picture on the left).

NUMERICAL SIMULATIONS

A finite element model was deveJoped in order to simulate a shaking table test. In the test, the honeycomb is placed on the shaking table with a mass on the top of it. The mass has the double function of shifting its characteri.stic frequencies towards lower values and assuring uniform displacements at the top. Transmissibility will be calculated from the ratio of out-of-plane acceleration measured both on the upper surface of the mass and on the shaking table.

From a numerical point of view, since in the Fourier domain the ratio of accelerations is equal to the ratio of displacements, an harmonic analysis was carried out to calculate the displacements, and so the transmissibility.

Finite Element model

ANSYS ® Mechanical APDL 15.0 was used for the Finite Element modelling and the harmonic analyses. On the other hand, MATLAB ® provided an interactive development environment to perform data analysis and visualization, numerical computation and algorithm development.

In terms of the FE modelling, the material and geometric properties of the auxetic and hexagonal cell honeycombs are shown in Figure 3. Moreover, damping was assumed to be 4%.

Transmissibility analyses

The harmonic analysis was carried out in the frequency range [10> 4000] Hz, with a step of 10

Hz. The output displacements were used to calculate the transmissibility, as the ratio of out-of plane displacements on the upper surface of the mass and on the honeycomb's base.

For both the hexagonal and the auxetic configurations, a parametric study was carried out by either varying the honeycomb height or the top mass. In paiticular, varying the height, the mass is kept constant and equal to 1.5 kg, while as the mass changes the height is 15 mm. Figure 4 shows the transmissibility curves for the hexagonal cell honeycomb in function of its height.

In Figure 5 , the graphs show the peaks of transmissibility for both hexagonal and auxetic configurations and for both parametric studies.

CONCLUSIONS

The results of the numerical parametric study show how, for both the auxetic and hexagonal configurations, the peak of transmissibility decreases as either the height of the honeycomb or the mass applied on the top of it increase. Furthermore, among the two configurations, the hexagonal shows better absorption properties than the auxetic one, as the peaks of transmissi bility have larger values for the latter. Future developments will be carrying out the shaking table test on the various configurations that were analyzed in this paper in order to compare nu merical and experimental results, and a Dynamic Mechanical Analysis characterization of the material in order to verify the value for the damping. Finally, in-plane vibroacoustic properties of the honeycomb cores will be investigated.

Figure 1 .

 1 Figure 1. Hexagonal cell honeycomb: corrugated sheet after the thermoforming and laser cut ting of the slits (left) and after folding along the slits (right).

3 .

 3 Figure 2 both the hexagonal and auxetic cell honeycombs are shown.

Figure 2 .

 2 Figure 2. Hexagonal cell honeycomb (left) and re-entrant (auxetic) cell honeycomb (right) resulting from the Kirigami manufacturing.

  Shell elements (SHELL181, element size= 1.5 mm) were used throughout the models. In h density, p [m kg-3 ) 1040 h [

Figure 3 .

 3 Figure 3. Hexagonal cell -From the left: characteristics of the material; honeycomb dimen sions; geometric characteristics of the cell. w = length in x direction, d = length in y direction. emphAuxetic cell -Geometric characteristics

Figure 4 .Figure 5 .

 45 Figure 4. Hexagonal cell honeycomb -Transmissibility curve for variable top mass values.


Kévin Billon
email: kevin.billon@femto-st.fr-morvan.ouisse@femto-st.fr-emeline.sadoulet-reboul@univ-fcomte.fr

Morvan Ouisse

Emeline Sadoulet

Manuel Collet
email: manuel.collet@ec-lyon.fr

E Sadoulet-Reboul

Numerical tools for efficient simulations of wave propagation in damped periodic structures

come

INTRODUCTION

A periodic medium is a material or a structural system that exhibits spatial periodicity. The study of periodic structures has a long history in the field of vibrations ans acoustics. Newton described the propagation of sound in air [START_REF] Brillouin | Wave propagation in periodic structures[END_REF]. He assumed that sound was propagated in air like an elastic wave would be propagated along a lattice of point masses. This topic has interested researchers over the years, and a growing activity on this field is observed on the last years, with the objective of designing structures exhibiting properties that conventional ones cannot possess. The methods currently used are most of the time based on those derived from wave propagation in crystals [2], where almost no dissipation occurs. Reaching the upper scale for structural dynamics implies that damping effects have to be included in the analyses which are performed. A very detailed review of historical origins, recent progress and future outlook of this topic has been published recently by [START_REF] Hussein | Dynamics of phononic materials and structures: Historical origins, recent progress, and future outlook[END_REF]. The reader is invited to refer to this article and the following discussion by [START_REF] Mace | Discussion of Dynamics of phononic materials and structures: Historical origins, recent progress, and future outlook[END_REF] that cover the most important aspects of this topic.

STRUCTURE PRESENTATION

The structure is presented in ref. [START_REF] Wu | Waveguiding and frequency selection of lamb waves in a plate with a periodic stubbed surface[END_REF]. The system consists in an infinite periodic bidirectional waveguide shown in (Fig. 1). It is a 1 mm thick plate with periodic cylindrical pillars made of isotropic Aluminium 6063-T83 (ν = 0.33, E = 69e9[P a] and ρ = 2700[kg/m 3 ]). Due to the periodicity, the unit cell is used and the corresponding first Brillouin zone is described in (Fig. 2).

CLASSICAL FLOQUET-BLOCH METHOD

The Floquet-Bloch approach is a method commonly used for the study of periodic structures but its use in a damping context is not easy in 2D and 3D cases.

The material constitutive law is linear, elastic and isotropic. The periodicity is defined on the borders of the domain u R = e -jkxr 1 u L and v R = e -jkyr 2 v L (Fig. 2) where u R is the displacement on the right border and u L is the displacement on the left border in x axis, v R is the displacement on the right border and v L is the displacement on the left border in y axis, k x and k y are respectively the wavenumbers in the x and y directions. Tetrahedral quadratic elements (TET10) are used to mesh the plate.

The harmonic homogeneous dynamical equilibrium of the system is driven by the following partial derivative problem

ρω 2 u + ∇σ = 0 σ = C : ε, (1)
where u ∈ R 3 is the displacement, σ is the stress tensor, C is the elastic tensor and ǫ is the strain tensor.

A parametric eigenvalue analysis is performed using the Pardiso solver [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with pardiso[END_REF], two parameters (wavenumbers) are considered, namely k x = [0 π/r 1 ] and k y = [0 π/r 2 ]. The eigen frequencies are obtained by solving the problem. The wave's dispersion curves of the undamped system are plotted on all the first Brillouin zone (Fig. 3) and on the contours (Fig. 4). The frequencies defining the bandgaps can always be found by considering only the contour of the irreducible Brillouin zone [2]. So, (Fig. 4) is sufficient to observe partial bandgaps and complete bandgap in all the directions. It appears that the structure exhibits both partial bandgaps for frequencies around 2 × 10 4 and 1 × 10 5 Hz in some specific directions and a complete bandgaps between 114 and 143 kHz.

In the approach, the wavenumbers are parameters of the eigenvalue problem, the solutions being the frequencies. Therefore, this approach is not suitable for systems with frequency-dependent characteristics.

"SHIFTED-CELL OPERATOR" METHOD

The "Shifted-Cell Operator" technique consists in a reformulation of the PDE problem by "shifting" in terms of wave number the space derivatives appearing in the mechanical behavior operator inside the cell, while imposing continuity boundary conditions on the borders of the domain. The formulation leads to the following eigenvalue problem [START_REF] Collet | Floquet-bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems[END_REF] [(K -

ω 2 M ) + λ i (L -L T ) -λ 2 i H]φ r i = 0, (2)
where λ i = jk i is the i-th eigen value, φ r i denotes the right eigenvector associated to λ i , M and K are respectively the standard symmetric definite mass and symmetric semi-definite stiffness matrices, L is a skew-symmetric matrix and H is a symmetric semi-definite positive matrix. In this formulation, all matrices can depend on ω.

Continuity boundary conditions are imposed on the borders of the domain, namely u R = u L and v R = v L .

A parametric eigenvalue analysis is performed using Pardiso solver. In this problem, the pulsation ω is fixed as real parameter. The wavenumbers λ i = jk i and the associated right eigenvectors φ r i are computed by solving the quadratic eigenvalue problem. This approach allows to introduce damping effects. The problem being damped, all the wavenumbers are complex, and the distinction between "propagative" and "evanescent" waves becomes difficult. This is why sorting criteria are proposed in the following.

First of all, all the waves are shifted to the first irreducible Brillouin zone. Then, criteria are defined as:

• ratio between real and imaginary part of each wavenumbers: real(k)/imag(k) > 1;

• ratio between real and imaginary part of the velocity of energy transport [START_REF] Vladimir | Complex group velocity and energy transport in absorbing media[END_REF] v = I/E where I is the energy flux and E the total energy density, approximated from the kinetic energy

E = 2E c : real(v)/imag(v) > 1;
• ratio between real and imaginary part of the group velocity [START_REF] Moiseyenko | Material loss influence on the complex band structure and group velocity in phononic crystals[END_REF]: real(C g )/imag(C g ) > 1.

For frequency-dependent damped systems, the estimation of the group velocity is not trivial. We propose to rewrite (Eq. 2) as

A 1 (ω)ψ r i = λ i A 2 (ω)ψ r i (3)
with

A 1 (ω) = 0 I d (K -ω 2 M ) L , (4)
A 2 (ω) = I d 0 L T H , (5)
and

ψ r i = φ r i λ i φ r i . (6)
Conversely, a left-eigenvector for the same eigenvalue satisfies

ψ l i T A 1 (ω) = λ i ψ l i T A 2 (ω). (7)
Generally, left and right eigenvectors for the same eigenvalue are not equal. In the resolution of the right eigenvalue problem, the i-th mode (i ∈ N * ) is defined by its λ i ≥ 0 and its eigenvector φ r i . For each mode i, a mode -i is associated with λ -i ≤ 0 such that λ -i = -λ i and φ r -i such that φ r -i = φ l i . So, by solving the right eigenvalue problem, the adjoint solution is found too. (Eq. 2) is differentiated and multiplied by the left eigenvector such that

∂A 1 (ω) ∂ω ψ r i + A 1 (ω) ∂ψ r i ∂ω = ∂λ i ∂ω A 2 (ω)ψ r i + λ i ∂A 2 (ω) ∂ω ψ r i + λ i A 2 (ω) ∂ψ r i ∂ω , (8)
ψ l i T ∂A 1 (ω) ∂ω ψ r i + ψ l i T A 1 (ω) ∂ψ r i ∂ω = ψ l i T ∂λ i ∂ω A 2 (ω)ψ r i + ψ l i T λ i ∂A 2 (ω) ∂ω ψ r i + ψ l i T λ i A 2 (ω) ∂ψ r i ∂ω , (9)
which gives the expression of the group slowness using

λ i = jk i ∂k i ∂ω = real   -j ψ l i T ∂A 1 (ω) ∂ω -λ i ∂A 2 (ω) ∂ω ψ r i ψ l i T A 2 (ω)ψ r i   , (10)
The group velocity is the inverse of the group slowness

C g = real ∂ω ∂k i = real jφ l i T [(L -L T ) -2λ i H]φ r i φ l i T [2ωM -∂K ∂ω -λ i ( ∂L ∂ω -∂L T ∂ω ) + λ 2 i ∂H ∂ω ]φ r i . (11)
We now focus on homogeneous cases where the frequency dependency is characterized by a Young's modulus such that E = f (ω)E 0 and a constant Poisson's ratio. Hence

K = f (ω)K 0 , H = f (ω)H 0 and L = f (ω)L 0 .
Therefore, the group velocity is expressed as

C g = real ∂ω ∂k i = real   jφ l i T [f (ω)(-L 0 + L T 0 + 2λ i H 0 )]φ r i φ l i T [ω 2 ∂f ∂ω f (ω) -2ω]M φ r i   . (12)
For frequency-independent cases, a simplification is possible [9]

C g = real ∂ω ∂k i = real j[φ l i T (L -L T )φ r i -2λ i φ l i T Hφ r i ] 2ωφ l i T M φ r i . (13)
A comparison between the results obtained using the Floquet-Block method and the "Shifted-Cell Operator" is firstly performed on a conservative structure to validate the implementation of both methods and in particular the implementation of the "Shifted-Cell Operator". Fig. 5 presents the dispersion curves obtained along the Γ-X direction when φ is equal to zero. Both methods lead to similar results and the shift-cell method is thus validated. Damping effects are then introduced in the system and results obtained using the "Shifted-Cell Operator" are presented with increasing damping on (Fig. 6). An hysteretic damping factor is used as a first step. A complex Young's modulus is obtained such E * = E(1 + jη) with η is the loss factor. Damping effects mainly affect the curves at higher wavenumbers. Moreover, the damping reduces the frequency bandgap and reduces the amplitude resonance. So, from a design point of view, a compromise has to be found.

The method has another advantage, the group velocity constitutes a pertinent indicator for the branches tracking from one computational point to another (see Fig. 7 and Fig. 8). In order to do that, an algorithm allows to compare the computed group velocity ∂ω ∂k i and the post processed group velocity ∆ω ∆k i . The finite difference method is used to calculate the post-processed group velocity directly with the dispersion relation.

CONCLUSION

Periodic structures are interesting for waves propagation as they may exhibit complete or at least partial frequency bandgaps such that the associated waves cannot propagate through the structure.

The determination of the dispersion curves is thus necessary to design specific structures for an absorption purpose. The Floquet-Bloch method is a commonly used approach well-adapted for such problems but its application in presence of damping is difficult. Another approach is here presented consisting in the "Shifted-Cell Operator" based on a reformulation of the PDE problem.

The method is first applied to undamped periodic structure and the good agreement with the results given by the classical method validate the implementation. The method is then applied to a structure presenting an hysteretic behavior : the impact of the damping on the frequency bandgaps can thus be analysed.

Figure 1 .

 1 Figure 1. Infinite plate.

Figure 2 .

 2 Figure 2. Real and reciprocal lattices.

Figure 3 .

 3 Figure 3. Dispersion surfaces on the Brillouin zone.

Figure 4 :

 4 Figure 4: Dispersion relation on the contour of the Brillouin zone. Complete and partial bandgap in the Γ -X direction respectively in the blue and green regions.

Figure 5 .

 5 Figure 5. Dispersion curves when φ = 0 • -Methods comparison.

Figure 6 .

 6 Figure 6. Dispersion curves with increasing hysteretic damping.

Figure 7 .

 7 Figure 7. Dispersion curves with branches tracking.

Figure 8 .

 8 Figure 8. Group velocity with branches tracking.

ACKNOWLEDGMENTS

This work was financed by The French National Research Agency under grant number ANR-12-JS09-008-COVIA. It has been performed in cooperation with the Labex ACTION program (ANR-11-LABX-0001-01).