Paolo Ballarini
email: paolo.ballarini@centralesupelec.fr

Benoît Barbot
email: benoit.barbot@u-pec.fr

PERFORMANCE MODELLING OF ACCESS CONTROL MECHANISMS FOR LOCAL AND VEHICULAR WIRELESS NETWORKS

Carrier sense multiple access collision avoidance (CSMA/CA) is the basic scheme upon which access to the shared medium is regulated in many wireless networks. With CSMA/CA a station willing to start a transmission has first to find the channel free for a given duration otherwise it will go into backoff, i.e. refraining for transmitting for a randomly chosen delay. Performance analysis of a wireless network employing CSMA/CA regulation is not an easy task: except for simple network configuration analytical solution of key performance indicators (KPI) cannot be obtained hence one has to resort to formal modelling tools. In this paper we present a performance modelling study targeting different kind of CSMA/CA based wireless networks, namely: the IEEE 802.11 Wireless Local Area Networks (WLANs) and the 802.11p Vehicular Ad Hoc Networks (VANETs), which extends 802.11 with priorities over packets. The modelling framework we introduce allows for considering: i) an arbitrarily large number of stations, ii) different traffic conditions (saturated/nonsaturated), iii) different hypothesis concerning the shared channel (ideal/non-ideal). We apply statistical model checking to assess KPIs of different network configurations.

Introduction

Communication protocols regulate the behaviour of communicating nodes within a concurrent environment. The Open System Interconnection (OSI) model [11] defines a layered architecture for network protocols. The Medium Access Control (MAC) layer, part of the data-link layer, determines which node is allowed to access the underlying physical-layer (i.e. the medium) at any given moment in time. A MAC scheme is mainly concerned with reducing the possibility of collisions (i.e. simultaneous transmissions over a shared channel) from taking place. The basic mechanism used for reducing the likelihood of collisions, usually referred to as Carrier Sense Multiple Access (CSMA), is that, before starting a transmission, any node should sense the medium clear for a given period. IEEE 802.11 WLAN. Wireless local area networks (WLANs) are wireless networks for which either the communication is managed by a centralised Access Point (AP) or, in the case of ad-hoc, nodes communicate in a peer-to-peer fashion through a distributed coordination function. The IEEE 802.11 [START_REF]IEEE. Ieee wirless lan medium access control (mac) and physical layer(phy) specification[END_REF] is a family of standards which specifies a number of MAC schemes and the Physical (PHY) layer for WLANs. The primary MAC scheme of the standard is called Distributed Coordination Function (DCF). It describes a de-centralised mechanism which allows network stations to coordinate for the use of a (shared) medium in an attempt to avoid collision. The DCF is a variant of the CSMA/CA MAC scheme developed for collision avoidance over a shared medium using a randomised backoff procedure. Two variants of the DCF have been defined in the standard: the Basic Access (BA), which uses a single acknowledgement to confirm the successful reception of a data packet, and the Request-To-Send/Clear-To-Send (RTS/CTS), which employs a double-handshaking scheme so to reduce costly collisions on large data packets. IEEE 802.11p VANET. In the realm of Intelligent Transportation Systems (ITS) the main concern is one of improving the effectiveness as well as the safety of future transportation systems. ITS entail hybrid communication scenarios where both Inter-Vehicle-Communication (IVC), based on ad-hoc connections between moving vehicles, and Roadside-Vehicle-Communication (RVC), concerned with the exchanging of information between moving vehicles and fixed roadside infrastructures, co-exist. In order to cope with specific needs of such hybrid scenarios, an adaptation of the IEEE 802.11 MAC layer, named 802.11p Wireless Access in Vehicular Environments (WAVE) standard, has been introduced for Vehicular Ad Hoc Networks (VANETs). The 802.11p MAC is based on an adaptation of the CSMA/CA scheme, called Enhanced Distributed Channel Access (ECDA) protocol, to the case of a network whereby traffic with different level of priority, called Access Categories (AC), circulates between wireless nodes.

Our contribution. In this paper we present a formal modelling study for the analysis of performances of wireless networks using RTS/CTS based MAC protocols specifically i) WLANs (i.e. 802.11) and ii) VANETs (i.e. 802.11p).

To this aim we develop formal models in terms of Generalized Semi-Markov Process (GSMP) expressed through an high-level stochastic Petri net formalism, namely Stochastic Symmetric Nets (SSN) [START_REF] Chiola | Stochastic well-formed colored nets and symmetric modeling applications[END_REF], which we analyse through specific performance indicators defined through the HASL properties specification language [START_REF] Ballarini | Hasl: A new approach for performance evaluation and model checking from concepts to experimentation[END_REF]. The models we developed are configurable and allow for taking into account different scenario including the network dimension, different incoming traffic conditions, the possible presence of faults affecting the channel during the transmission of some packet. For both the 802.11 and 802.11p scenarios we developed and asses KPI expressed in terms of temporal logic specifications.

Related work

Performance analysis of the IEEE 802.11 family of Wi-Fi protocols has been the subject of several research studies [START_REF] Bianchi | Performance analysis of the IEEE 802.11 distributed coordination function[END_REF][START_REF] Frederico | Specification and verification of the IEEE 802.11 medium access control and an analysis of its applicability to real-time systems[END_REF][START_REF] Lyakhov | Hybrid rts/cts mechanism in wi-fi ad hoc networks with correlated channel failures[END_REF][START_REF] Kwiatkowska | Probabilistic model checking of the IEEE 802.11 wireless local area network protocol[END_REF] each of which considers specific modelling assumptions concerning e.g. the traffic model that is considered (i.e. saturated traffic, packets arrival following a Poisson law, etc), the presence/absence of errors on the channel and in case of errors how errors are modelled. In its pivotal work Bianchi [START_REF] Bianchi | Performance analysis of the IEEE 802.11 distributed coordination function[END_REF], introduced a simple analytical model that, under given constraints (i.e. finite number of stations and ideal channel condition) allow for computing the throughput of the IEEE 802.11 Distribution Coordination Function (DCF) for both the Basic Access (BA) and the RTS/CTS versions of the DCF. Taking from the two-dimensional discrete-time Markov chain (DTMC) model of Bianchi, many extensions have been considered. In [START_REF] Alshanyour | Performance of IEEE 802.11 RTS/CTS with finite buffer and load in imperfect channels: Modeling and analysis[END_REF] a 4D DTMC, inspired by Bianchi's 2D model, has been introduced for considering the case of imperfect channels, that is, taking into account that transmission over a wireless medium is affected by errors. More specifically in [START_REF] Alshanyour | Performance of IEEE 802.11 RTS/CTS with finite buffer and load in imperfect channels: Modeling and analysis[END_REF] the imperfect nature of channels is modelled through a constant bit error probability P b (i.e. the probability that an error occurs during the transmission of a single bit of data). In [START_REF] Kwiatkowska | Probabilistic model checking of the IEEE 802.11 wireless local area network protocol[END_REF] the BA version of the 802.11 MAC is analysed through probabilistic model checking but taking into account only simple modeling assumption (2-nodes network dimension, ideal channel, no specific traffic model).

RTS/CTS carrier-sensing protocol

CSMA/CA MAC schemes are based on the simple idea that a station willing to transmit data packets has first to gain access to the channel through a sensing phase depending on which the station may either start transmitting, if the channel has been sensed free along the entire sensing phase, or, if the channel has been sensed busy, refrain from doing so for a randomly chosen duration (backoff). Collisions may take place whenever at least two stations ends the sensing phase at the same time, however by employing a randomised backoff delay, the probability of collisions decreases with the number of successive collisions. There are two versions of CSMA/CA scheme, the basic access (BA), which we do not consider in this paper, and the Request-to-Send/Clear-to-Send (RTS/CTS).

802.11 RTS/CTS: with the RTS/CTS DCF the transmission of data between a sender and a receiver is regulated by means of a double handshaking scheme which employs three small-sized, control-packets to regulate the transmission of (larger) data packets (Figure 1). The sender senses the channel for a randomly chosen duration before sending a RTS control packet to the receiver, on reception of which, the receiver, replies with a CTS control packet to the sender (first handshaking). The actual transmission of DATA packets starts as soon as the sender has received the CTS. Finally when the DATA transmission is over, the receiver acknowledges the sender with an ACK control packet (second handshaking). The RTS-CTS-DATA-ACK timing sequence is illustrated in Figure 1, which also points out the latency for a successful transmission of a DATA packet over 1-hop. To decrease the probability of collisions, the carrier-sensing time t CS , EDCA is designed so that higher-priority traffic is more likely (than lower-priority) to be be granted access to the shared medium and successfully being transmitted. In practice this is obtained by associating each AC with a contention window (CW) whose size is inversely proportional to the corresponding priority level. Therefore, for example, the CW for AC_BK (lowest priority) has to be larger than that of AC_VO (highest priority). Table 2 depicts the CW characterisation for the ACs of ECDA, where aCWmin and aCWmax are parameters set by the 802.11p protocol. One can notice that AC_VO is the dominant AC as its CW is not only the minimal one but also it does not overlap to the CW of any other AC.

Background: stochastic symmetric Petri nets

To carry out the performance modelling study of 802.11 protocols we used a coloured stochastic Petri net formalism, namely the Stochastic Symmetric Net (SSN) [START_REF] Chiola | Stochastic well-formed colored nets and symmetric modeling applications[END_REF], to model networks using 802.11/802.11p MAC, and we applied the Hybrid Automata Specification Language (HASL) [START_REF] Ballarini | Hasl: A new approach for performance evaluation and model checking from concepts to experimentation[END_REF] formalism (through the COSMOS statistical model checker [START_REF] Ballarini | COSMOS: A statistical model checker for the hybrid automata stochastic logic[END_REF]) to assess relevant performance indicators against the SSN models. For the sake of space we only give a very succinct overview of the SSN and HASL formalisms referring the reader to the literature for a detailed treatment.

Stochastic Symmetric Nets. Like any Petri net formalism an SSN model is a bi-partite graph consisting of place nodes (circles) and transition nodes (bars). Places may contain (countably many) tokens and are connected to transitions through arcs which are labelled with arc functions. An SSN model describes a continuous-time, discrete-state stochastic process whose states correspond with the possible markings of the SSN (a marking gives the content of each place of the SSN). The peculiarity of SSN models is that tokens may be associated with information, hence they may have different colours, instead of being all indistinguishable as with non-coloured PN formalisms. Therefore places and transitions of an SSN are associated with a color domain (CD) built from elementary types called color classes ({C1, ..., Cn}) with cd(p), resp. cd(t), denoting the CD of place p, resp. transition t (see Table 4 and Table 3 for the CDs of 802.11p SSN model). SSN color classes are finite, non empty and disjoint sets, they may be ordered (in this case a successor function is defined on the class, inducing a circular order among the elements in the class), and may be partitioned into (static) subclasses. An SSN transition may be immediate (drawn as a thin filled in bar) or timed which means its firing delay is sampled (the moment it gets enabled) from a probability distribution that may be exponential (thick empty bar) or general (thick filled in bar). A coloured transition may be associated with a guard, i.e. a boolean-valued expression built on top of colored variables by means of the following basic predicates: x = y, x ∈ subclass, d(x) = d(y) where x and y are variables of the transition with same type, and d(x) denotes the static subclass x belongs to.A valid transition binding is an assignment of values to its variables, satisfying the predicate expressed by the guard. A pair (transition,binding) is called transition instance. Each arc connecting a place p and a transition t is labeled with an expression denoting a function arcf : cd (t) → Bag(cd (p)) where Bag(A) is the set of all possible multisets that may be built on set A. The valuation of arcf given a legal binding of t gives the multiset of colored tokens to be withdrawn from (in case of input arc) or to be added to (in case of output arc) the place connected to that arc upon firing of such transition instance. The arc expressions in SSNs are built upon a limited set of primitive functions whose domains must be color classes. a packet ready to be transmitted is invariably present) or b) its incoming traffic is given by an Poisson process with parameter λ. 4) Perfect/Imperfect channel: the wireless medium is supposed to behave either a) as a perfect channel (i.e. transmitted data are never affected by errors) or b) as an imperfect channel with some error probability (details below). 5) Vulnerable period: the radio device of each station is supposed to exhibit a certain delay for switching between transmission/reception mode. The duration of such period is called the vulnerable period and is one source of traffic collisions within the network.

Channel error model. To model the possibility that transmissions undergo errors due to the channel we adopted the burst-noise binary channel model [START_REF] Gilbert | Capacity of a burst-noise channel[END_REF]. A burst noise channels may be subsumed by a two-states Markov chain where one state (G as in good) represents the absence of noise while the other (B as in burst errors) represents the presence of an error spike which affects the channel impeding a transmission to correctly take place. We equipped our Petri nets models with an implementation of such burst-noise binary channel model. 3) and color domains (Table 4) that are used to characterise the type of the various elements (places and tranisitions). The model consists essentially of three parts: the 802.11p core module, representing the prioritised carrier-sensing and RTS/CTS handshaking mechanism, the backoff module, representing the behaviour of a station in case of a collision and the medium module representing the state of the wireless channel. We present these 3 modules separately bearing in mind that they are part of the same SSN model. Observe that the immediate transition InternalCollision (which is also consuming tokens from Vulnearble) deals with the case of competing packets sent, with different level of priorities, by a common station sa: in this case only the packet with the highest priority stays in Vulnerable, while packets with lower priority are moved to GoingToBackoff (i.e. highest priority packet wins the internal competition by "overtaking" lower priority packets). From place Sending a token sa, sb, p, pt moves to PacketSent if the packet type is either rts, cts or data, after a delay corresponding with the kind of packets, and then, with no delay, it is further moved to WaitForResponse through immediate transition BeginWaitingForResponse which while removing a token from Medium (hence decreasing the occupation of the wireless channel) and storing the sa, sb, p, pt in place SentPacket for later retrieval, also updates the packet type of token sa, sb, p, pt to sa, sb, p, pt++ , indicating that station sa is now ready to wait for a reply message (of type pt++, i.e. cts in reply to rts or ack in replay to data) from sb. Place WaitForResponse is initialised with marking St, s1, pr1, rts indicating that each station of St is, initially, ready to get engaged into replying to an rts request (s1 and pr1 being irrelevant at this stage). As the medium gets occupied but not garbled all tokens in WaitForResponse are moved (with no delay) to Receiving and from there, as soon as the channel gets free (and assuming it hasn't been garbled before) they move on to ReadingPacket. Thus, at the end of a transmission phase, place ReadingPacket contains all tokens representing all possible combinations of responses that all stations may be engaged in. On completion of the transmission corresponding to token sb, sa, p, pt being added to SentPacket the corresponding token sa, sb, p 2 , pt , if present, is consumed from ReadingPacket GoingToBackoff (if they represent a sender that was expecting either a cts or an ack response to a previosuly sent packet rts or data). 5). The marking of BackoffMappings associates with each priority level pr i and re-transmission attempt tx j the corresponding CW expressed as the union of color subclasses ∪ k k=1 bs k (where k is a function of pr i and tx j ). Therefore, for example, for the first re-transmission attempt (i.e. tx = tx 1 ) of a priority p 1 packet the backoff value b is chosen randomly as b ∈ bs 1 = {b1 , . . . , b4 } (which corresponds with the narrowest CW i.e. the AC_V0 category) as the initial marking BM 0 contains tokens < pr1, tx1, bs1 >. Similarly for the second re-transmission attempt (of a priority p 1 packet) b is chosen randomly as b ∈ bs 1 ∪ bs 2 = {b1 , . . . , b8 } as marking BM 0 contains tokens < pr1, tx2, bs1 + bs2 > and so on. Once a sa, sb, p, b token is added to place BackoffCounter (i.e. the backoff delay b is selected) the actual backoff (sensing) period begins: b is decremented every aSlot time units (deterministic transition DecrementBackoff ) however, as soon as the medium gets occupied, the decrement of b is suspended by moving token sa, sb, p, b from BackoffCounter to PausingBackoff (immediate transition PauseBackoff ). Token sa, sb, p, b is moved back to BackoffCounter only after the channel got freed and stayed free for AIFS i time units (where i ∈ {1, 2, 3, 4} denotes the priority level p of the backoff-ing packet). The backoff for a token sa, sb, p, b ends as soon as the backoff counter has reached 1 (i.e. b = b1): at this point token sa, sb, p, b1 is removed from BackoffCounter (immediate transition ExitBackoff ) and sa, sb, p, rts is pushed back to Sense which represent the re-start of the transmission procedure for station sa to sent a priority p packet to sb.

802.11p core SSN module

802.11p backoff SSN module

802.11p medium SSN module Figure 3 (bottom right) shows the part of the SSN model representing the shared channel. It consists of 2 uncoloured places: Medium, whose marking corresponds with the number of transmitting stations, and Garbled , which contains a token whenever a collision has taken place or an error spike has occurred on the channel. Transition GettingGarbled sets the state of the medium to garbled (adding a token in place Garbled ) as soon as Medium contains at least 2 tokens (i.e. at least two stations are transmitting). Transition EndingGarbled ends the garbled state of the medium by removing a token from Garbled as soon as the Medium is emptied. Place SentPacket (domain P ) stores the information relative to an ongoing transmission in terms of a token sa, sb, p, pt where sa is the sender station, sb the receiver ,p the priority level p and tp the packet type). The medium SSN is also equipped with a subnet for modelling the presence of error spikes affecting an ongoing transmission. To study the effect of a noisy channel on the network performances it suffices that place NotErrorSpike initially contains a token in which case transition EnterinErrorSpike reproduce the occurence of errors on the channel (following an Exponential distribution with configurable rate, 0.01 in the picture), by adding a token in Medium (hence triggering a collision if a transmission is going on). The end of an error-spike is modelled by transition ExitingErrorSpike (also Exponentially distributed with a configurable rate).

analysis of 802.11/802p models

We analysed the performances of both the 802.11 and 802.11p network models 2 by means of the COSMOS statistical model checker [START_REF] Ballarini | COSMOS: A statistical model checker for the hybrid automata stochastic logic[END_REF][START_REF]COSMOS home page[END_REF]. We considered a number of key performance indicators (KPI) including 1) the throughput (THR) of a network station (i.e. the number of successfully transmitted packets per time unit); 2) the busyTimeRatio (BTR, i.e. the ratio between the time the channel is occupied by some transmitting station and the total operation time of the network). We encoded such KPIs in terms of HASL specifications for SSN models [START_REF] Elvio | Simulation-based verification of hybrid automata stochastic logic formulas for stochastic symmetric nets[END_REF] and used them to run a number of experiments to assess the impact that 1) a faulty channel, 2) the incoming traffic, 3) the network dimension have on the KPIs. In the remainder we present an excerpt of the results obtained. Figure 4 shows an example of an HASL specification (a linear hybrid automaton, LHA) we used to assess the THR for priority p traffic on the 802.11p model. The LHA in Figure 4 has an initial (l 1 ) and a final (l 2 ) location and uses a clock (t) plus a counter of successfully ). On processing a trajectory T hroughput[p] is incremented each time an instance of the CorrectButAck transition (Figure 2) with the corresponding priority p level occur (which coincides with the termination of the transmission of a priority p packet). The LHA ends measure at time T and the value of THR is obtained as T hroughput[p]/T .

Results

Impact of network dimension and incoming traffic. 2 SSN models and HASL properties used to run such experiments on COSMOS are available at https://sites.google.com/ site/pballarini/models.

3 the BTR of a network using 802.11 RTS/CTS MAC has a nominal maximum upper bound (NomMaxBTR) corresponding with an ideal (hypothetical) situation where a continuos flow of DATA packets are transmitted without inter-arrival delay and in absence of collisions. In such situation the BTR is given by the ratio of busy-time over total-time for successfully transmitting, where the latter is given the sum of delays of sequence DIFS-RTS-SIFS-CTS-SIFS-DATA-SIFS-ACK sequence, that is NomMaxBTR = T Xtime/(SENSEtime + TXtime), where TXtime = RT S + CT S + DAT A + ACK, is the total transmission time and SENSEtime = DIFS+3SIFS is the total carrier-sensing time. From parameters in Table 1, we have that NomMaxBTR = 80/88 = 0.909, meaning that the optimal channel utilisation for a WiFi based on RTS/CTS 802.11 MAC cannot trespass ∼ 90%. (h) traffic rate λ = 1000 802.11p

Figure 5: Impact of network dimension and noisy-channel on throughput and BTR for 802.11 and 802.11p networks higher the traffic the larger the difference between high priority and low priority throughput. Figure 5(f) shows the BTR as well as the ratio between the time the channel is idle (over the total observation time) and the ratio between the time the channel is garbled (collision).

Impact of faulty channel. Figure 5(c), 5(d), 5(g) and 5(h) report on assessing the effect that a faulty channel has on the network performances. The THR and BTR are measured in function of the period of an error spike (determined by the rate of transition EnteringErrorSpike). Results show that the throughput increases with the error-spike period however (for the 802.11 model under saturated regime Figure 5(c)) the THR is unaffected by the network dimension (identical plot for different number of stations) indicating that the performances gradient induced by the network dimension in case of ideal channel (Figure 5(a)) fades away in presence of a faulty channel. Conversely under a non-saturated regime (802.11 model Figure 5(d)) the error-spike period affects the throughput differently depending on the traffic arrival rates. Finally Figure 5(g) and 5(h) show the effect of the error-spike on the prioritised throughput in a 802.11p network showing that the effect of a faulty channel, in terms of the gradient between higher and lower priority throughput is more evident under a high traffic regime (Figure 5(h)) than under a low traffic regime (Figure 5(g)).

Conclusion

We presented a performance modelling study of two versions of MAC protocol for wireless networks: the 802.11 MAC and its prioritised extension 802.11p devoted to VANETS. We developed our models using a high-level stochastic Petri nets formalism which allowed us to encode the complexity of the 802.11 prioritised scheme in a model of reasonable size. The models we presented are highly configurable and allow for the analysis of the performance of networks in different respect (traffic conditions, ideal or imperfect channel, network dimension). We analysed the performance on the network models by means of statistical model checking based on the HASL specification language. Future work include the extension of the models to N -hops topologies, which would allow us to take into account the effect of routing on the performances of a given network.

Figure

  Figure 1: 1-hop carrier-sensing and packets-transmission timing in RTS/CTS 802.11

Figure 2 Figure 2 :

 22 Figure 2: SSN model of the 802.11p prioritised RTS/CTS scheme

Figure 3 :

 3 Figure 3: SSN model of the 802.11p backoff mechanism (top left) and the medium (bottom right).

  With the 802.11p DCF the selection of the backoff duration, i.e. the selection of the number of time slots a packet has to wait before starting a new transmission attempt, depends on the priority class of the packet as well as on the backoff round, i.e. the number of times a packet has unsuccessfully went through the carrier-sensing phase without winning the contention. Therefore a model of 802.11p backoff procedure must be equipped with necessary means to take into account the priority level of each packet hence the size of the contention window for each priority level. Figure 3 (top left) depicts the SSN subnet representing the randomised selection of the backoff delay for a packet of a given priority level p. A backoff begins when a token sa, sb, p is added to place GoingToBackoff (color domain BP ) indicating that something went wrong (e.g. a collision or a lack of handshaking) during the transmission of a priority p packet from station sa towards sb. The first step then is the updating of the backoff round counter which boils down to adding of a token sa, sb, p, tx to place TxAttemptsCounter (whose color domain is SRPCount and where variable tx ∈ TxCount stores the backoff counter) through either transition GoToBackoff (enabled if the failed transmission attempt is the first one, hence initialising tx to color tx 1 ) or through transition ReturnToBackoff (enabled by any further failed transmission attempt and increasing tx by 1). Once a token sa, sb, p enters place ChoosingBackoff having, in the process, also added (or updated) a sa, sb, p, tx token to place TxAttemptsCounter , either one between the mutually exclusive transitions ChooseBackoff or DropPacket transition is enabled. Transition ChooseBackoff realises the probabilistic selection of the backoff counter value b, corresponding to the priority level p of the backoffing packet, by adding token sa, sb, p, b to place BackoffCounter . Such a selection is achieved through randomly choosing (through the test arc between transition ChooseBackoff and place BackoffMappings) a token p, tx, b among those specified by the invariant marking of BackoffMappings (Table

Figure 4 :

 4 Figure 4: The LHA for measuring the throughput terminated transmissions for a priority p packet (T hroughput[p]). On processing a trajectory T hroughput[p] is incremented each time an instance of the CorrectButAck transition (Figure2) with the corresponding priority p level occur (which coincides with the termination of the transmission of a priority p packet). The LHA ends measure at time T and the value of THR is obtained as T hroughput[p]/T .

Figure 5 (

 5 a) and 5(b) depict the estimated THR and BTR for the 802.11 model computed (for different network dimensions) in function of the traffic arrival rate λ (i.e. under a non-saturated regime) and assuming the channel is not affected by errors (ideal channel). Both the THR and the BTR exhibit an asymptotic profile (the faster the inter-arrival rate the higher the throughput, resp. the BTR) and both are upper bounded by the nominal maximum throughput under saturated regime (which is ∼ 640 packets per second), resp. the saturated nominal maximum BTR 3 Figure5(e) and 5(f), refer to the same kind of experiment for the 802.11p model (for a network with N = 2 nodes). Plots in Figure5(e) highlight the effect of prioritised management of data traffic:

  the shortest number of available slots in the CW, wins the contention. The loser, instead, goes into backoff and increases the dimension of the CW by a power of 2 (i.e. the new contention window size CW is given by CW = (CW m + 1) • 2 bc -1 where bc, the BackoffCounter, increases with the number of consecutive unsuccessful transmissions). The minimum and maximum size of contention window, respectively CW m and CW M , are set by the standard and depends on the PHY layer. For the Frequency Hopping Single Spectrum (FHSS) PHY layer they are equal to CW min = 16 (initial value of bc = 0) and CW max = 1024 (maximum value of bc max = 6). Apart from t CS , two fixed-length time intervals are also relevant in the RTS/CTS DCF, namely: Distributed InterFrame Space (DIFS), the Short InterFrame Space (SIFS), where SIF S < DIF S and which are defined by the PHY layer in the adopted networking stack (see Table1). It should be noted that a collision can take place not only when two contending stations (randomly) pick the same carrier-sensing time (i.e. t CS ), but also, as pointed out by Heindl it et al.[START_REF] Heindl | Performance modeling of ieee 802.11 wireless lans with stochastic petri nets[END_REF], because of the existence of the so-called vulnerable period, which accounts for three factors: the time of radio waves propagation through the medium (aAirP ropagationT ime), the time a station takes for accessing the medium (aCCAT ime) and the time a station takes for switching from receiving to transmitting mode(aRxT xT urnaroundT ime). As a consequence the duration of aSlot t is set by the IEEE standard to a value larger than the vulnerable period (i.e. the value aSlot t depends on the PHY-layer dependent parameters aCCAT ime and aRxT xT urnaroundT ime plus the negligible aAirP ropagationT ime).

	1: 1-hop carrier-sensing and packets-transmission timing in RTS/CTS 802.11
	also named Contention Window (CW), which is discretely slotted in a finite number of slots of fixed duration aSlot t
	(aSlot t being a parameter of the standard whose value depends on the underlying PHY layer), is randomly chosen. A
	contention takes place when at least two stations are performing carrier-sensing at the same time. The one that has
	(randomly) chosen 802.11p RTS/CTS with priority classes: The 802.11p MAC is based on an extension of the CSMA/CA RTS/CTS
	scheme, called Enhanced Distributed Channel Access (ECDA), which supports 4 priority levels, called Access
	Categories (AC), for data traffic, namely: Background (AC_BK), i.e. the lowest priority, Best effort (AC_BE), Video
	(AC_VI) and Voice (AC_VO), i.e. the highest priority.

  Typically an arc expression is a linear combination of function tuples (denoted f 1, . . . , f n ), and each element of a tuple is either a projection function, denoted by a variable in the transition color domain (e.g. sa and sb in the tuple sa, sb, p, pt appearing as the labelling in several arcs of the SSN of Figure2), a successor function, denoted x ++ were x is a variable whose type is an ordered class; a constant function, denoted C i .All returning all elements of (sub)class C i ; a complement function denoted C i .All -x where x is a variable of type C i . The dynamics of an SSN model is defined in terms of transition instances enabling and firing: a transition instance is enabled if the marking of all of its input places is compatible with it (i.e. if the marking enables the transition instance). Upon firing a transition instance modifies the state of the SSN by removing (resp. adding) tokens from its input (resp. output)

places. For example, w.r.t. the SSN in Figure

2

, assuming colour classes St = {st1 , st2 } and Pr = {pr1 }, the initial marking of place Idle (colour domain St ×Pr ) being St, Pr = { st1 , pr1 , st2 , pr1 } enables two instances of timed transition PacketArrival (given its guard sa = sb), namely st1 , st2 , pr1 , rts st2 , st1 , pr1 , rts . The firing of the first instance consumes st1 , pr1 from Idle and adds st1 , st2 , pr1 , rts to Sense.

4 Modelling of 802.11/802.11p networks

To analyse the performance of two versions of the protocol we developed two SSN models, one modelling a network whose stations use the 802.11 MAC, the other where stations that use the 802.11p MAC. The models we developed are based on the following assumptions: 1) Clicque network topology: the network consists of N stations arranged in a clique (i.e. every station can overhear every other station). 2) Traffic direction: each station has a unique target station to which it addresses its incoming traffic. 3) Incoming traffic: each station is either under a) a saturated regime (i.e.

  4.1 SSN models of 802.11/802.11p networks We developed two SSN models one for 802.11 networks the other for 802.11p's. For the sake of space we only present the 802.11p version of the model, which is somehow a generalisation of the 802.11's. The 802.11p SSN model relies on the definition of a number of color classes (Table

Table 1 :

 1 A Parameters of the SSN model of the 802.11 and 802.11p MAC Timing parameters for the SSN models of RTS/CTS

	name		meaning	value (×10µs)
	aslot	the time unit of the backoff procedure	2
	nStations		# of stations	2
	packetSizeInAslot	# of time slot for sending a DATA packet	variable
	DIFS	the length of the DCF interframe space	5
	SIFS	the length of a short interframe space	1
	vuln	delay to switch radio from RX-to-TX mode	2
	RTS	the time it takes to send a RTS packet	16
	CTS	the time it takes to send a CTS packet	11
	ACK	the time it takes to send an aknowledgement packet	11
	timeout	delay a station waits for handshaking packet	5
	CWmin	min. size of the contention window (in aslot)	15
	backoffMax	# of TX attempts before dropping a packet	6
	deadline	# of time slots for entering the livelock (end) state	5000
		AC	CWmin	CWmax
	Background (AC_BK)	aCWmin	aCWmax
	Best effort (AC_BE)	aCWmin	aCWmax
	Video (AC_VI)	[(aCWmin+1)/2]-1	aCWmin
	Voice (AC_VO)	[(aCWmin+1)/4]-1	[(aCWmin+1)/2]-1

Table 2 :

 2 Contention windows boundaries for ECDA access categories. Priority level (of a packet) P r = p{1, . . . , 4} NO T xCount max. num. of re-transmissions before dropping a packet T xCount = tx{1, . . . , 20} YES BackoffStage backoff counter domain BackoffStage = bs 1 ∪ . . . ∪ bs 9 YES bs 1 backoff counter domain stage 1 bs 1 = b{1, . . . , 4} YES bs 2 backoff counter domain stage 2 bs 2 = b{5, . . . , 8} YES bs 3 backoff counter domain stage 3 bs 3 = b{9, . . . , 16} YES bs 4 backoff counter domain stage 4 bs 4 = b{17, . . . , 32} YES bs 5 backoff counter domain stage 5 bs 5 = b{33, . . . , 64}

	name	description	definition	ordered
	P T	Packet Type	P T = {rts, cts, data, ack}	YES
	St	Station ID	St = s{1, . . . , N }	NO
	P r
				YES
	bs 6	backoff counter domain stage 6	bs 6 = b{65, . . . , 128}	YES
	bs 7	backoff counter domain stage 7	bs 7 = b{129, . . . , 256}	YES
	bs 8	backoff counter domain stage 8	bs 8 = b{257, . . . , 512}	YES
	bs 9	backoff counter domain stage 9	bs 9 = b{513, . . . , 1024}	YES

Table 3 :

 3 Color classes for the SSN model of the 802.11p protocol. St × St × P r × P T SRB Sender-Receiver-BackoffStage SRB = St × St × BackoffStage BP M ap Mapping-Priority-Backoff BP M ap = P r × T xCount × BackoffStage SRP count Sender-Recevier-Priority-Count SRP count = St × St × P r × TxCount

	name	description	definition
	SR	Sender-Receiver	SR = St × St
	SP	Station-Priority	SP = St × P r
	BP	Backoff-Priority	BP = St × St × P r
	P	Packet sending	P =

Table 4 :

 4 Color domains for the SSN model of the 802.11p protocol.

	marking	definition
		<pr1,tx1,bs1>+<pr1,tx2,bs1+bs2>+<pr1,tx3,bs1+bs2>+<pr1,tx4,bs1+bs2>
		+<pr2,tx1,bs1>+<pr2,tx2,bs1+bs2>+<pr2,tx3,bs1+bs2+bs3> +
		+<pr2,tx4,bs1+bs2+bs3>+<pr2,tx5,bs1+bs2+bs3>+<pr3,tx1,bs1+bs2> +
		<pr3,tx2,bs1+bs2+bs3>+<pr3,tx3,bs1+bs2+bs3+bs4>
		+<pr3,tx4,bs1+bs2+bs3+bs4+bs5> + <pr3,tx5,bs1+bs2+bs3+bs4+bs5+bs6>
		+<pr3,tx6,bs1+bs2+bs3+bs4+bs5+bs6+bs7> +
		<pr3,tx7,bs1+bs2+bs3+bs4+bs5+bs6+bs7+bs8>
		+<pr3,tx8,bs1+bs2+bs3+bs4+bs5+bs6+bs7+bs8+bs9> +
	BM 0	<pr3,tx9,bs1+bs2+bs3+bs4+bs5+bs6+bs7+bs8+bs9> +
		<pr3,tx10,bs1+bs2+bs3+bs4+bs5+bs6+bs7+bs8+bs9>
		+<pr4,tx1,bs1+bs2+bs3>+<pr4,tx2,bs1+bs2+bs3+bs4>
		+<pr4,tx3,bs1+bs2+bs3+bs4+bs5>+<pr4,tx4,bs1+bs2+bs3+bs4+bs5+bs6>
		+<pr4,tx5,bs1+bs2+bs3+bs4+bs5+bs6+bs7>
		+<pr4,tx6,bs1+bs2+bs3+bs4+bs5+bs6+bs7+bs8>
		+<pr4,tx7,bs1+bs2+bs3+bs4+bs5+bs6+bs7+bs8+bs9> +
		<pr4,tx8,bs1+bs2+bs3+bs4+bs5+bs6+bs7+bs8+bs9> +
		<pr4,tx9,bs1+bs2+bs3+bs4+bs5+bs6+bs7+bs8+bs9>

Table 5 :

 5 Invariant marking of place BackoffMapping.


A modelling framework for finite strain magneto-viscoelasticity

Keywords: Magneto-active soft materials, Large deformation, Multiplicative viscoelasticity, Continuum thermodynamics Prepared using sagej.cls

The main objective of this contribution is the theoretical modelling of magnetisable materials that exhibit viscoelastic properties. The current state of the art in the mathematical modelling of nonlinear magneto-mechanics in deformable media can be easily integrated within the unified framework of continuum thermodynamics that is crucial in setting the convenient forms for the constitutive laws and evolution equations. Due to the soft nature of the materials we have in mind, the finite strain range is a priori adopted and, in this first approach, only isotropic materials are considered. We adopt the nowadays well-accepted multiplicative split of the deformation gradient which, moreover, gives rise to an intermediate configuration. Herein, the novelty resides in the fact that the magnetic field vectors are transported onto the above mentioned intermediate configuration and, therefore, can be implicitly decomposed. The proposed formulation is based on the magnetic induction as main independent variable for the magnetic part of the problem. An alternative formulation based on the magnetic field as main independent variable can easily be deduced, but this latter will not be considered in this paper for the sake of clarity. A very simple model example that agrees with the laws of thermodynamics is proposed for the purpose of demonstration to study some phenomena qualitatively.

Introduction

Magneto-active polymers (MAPs) belong to the very important class of emerging new materials that deform under magnetic stimuli. They are mostly composites of soft polymer matrix impregnated with magnetically permeable particles, typically iron particles in micro-or nano-meter size. In general, the response to magnetic fields can be devided into two categories based on the property of the matrix part of the material: they can give large and prompt deformation, or they can change their mechanical properties with moderate straining. The two features have received considerable attention in recent years due to their potential applications including, for instance, sensors, actuators, and medicine, see for example [START_REF] Ginder | Magnetostrictive phenomena in magnetorheological elastomers[END_REF][START_REF] Jolly | A model of the behavior of magnetorheological materials[END_REF][START_REF] Varga | Magnetic field sensitive functional elastomers with tunable elastic modulus[END_REF][START_REF] Zrínyi | Deformation of ferrogels induced by nonuniform magnetic fields[END_REF][START_REF] Ren | Multi-functional soft-bodies jellyfish-like swimming[END_REF][START_REF] Hu | Small-scale soft-bodied robot with multimodal locomotion[END_REF][START_REF] Lu | A bioinspired multilegged soft millirobot that functions in both dry and wet conditions[END_REF] among many others.

In parallel, the mathematical modelling of the coupling of electromagnetic fields in deformable materials has also been an area of active research, see for example [START_REF] Brigadnov | Mathematical modeling of magneto-sensitive elastomers[END_REF][START_REF] Bustamante | Universal relations in isotropic nonlinear magnetoelasticity[END_REF][START_REF] Dorfmann | Magnetoelastic modelling of elastomers[END_REF][START_REF] Dorfmann | Nonlinear magnetoelastic deformations[END_REF][START_REF] Dorfmann | Nonlinear magnetoelastic deformations of elastomers[END_REF][START_REF] Ericksen | A modified theory of magnetic effects in elastic materials[END_REF][START_REF] Steigmann | Equilibrium theory for magnetic elastomers and magnetoelastic membranes[END_REF][START_REF] Vu | Nonlinear electro-and magneto-elastostatics: Material and spatial settings[END_REF][START_REF] Nedjar | A theory of finite strain magneto-poromechanics[END_REF] . Fully coupled nonlinear field theories have been developed with constitutive formulations based on augmented free energy functions.

In particular, it has been shown that any one of the magnetic induction, magnetic field, or magnetization vectors can be used as independent variable, the other two being obtained through constitutive relations, see [START_REF] Ogden | Magnetostatics: from basic principles to nonlinear interactions in deformable media[END_REF] for details about this topic. The relevant equations are based on the pioneering works of [START_REF] Brown | Magnetoelastic Interactions[END_REF][START_REF] Jackson | Classical Electrodynamics[END_REF][START_REF] Pao | Electromechanic forces in deformable continua[END_REF] among others.

This work aims at developing a continuum field theory for the behaviour of magnetically sensitive materials that experience viscoelasticity, see for example [START_REF] Amin | Nonlinear dependence of viscosity in modeling the rate-dependent response of natural and high damping rubbers in compression and shear: Experimental identification and numerical verification[END_REF][START_REF] Liao | Temperature and strain dependent large tensile deformation and tensile failure behavior of transparent polyurethane at intermediate strain rates[END_REF] . Indeed, the response to magnetic stimuli can in general not be instantaneous for all MAPs as that has been established experimentally, e.g. see the recent works in [START_REF] Pelteret | Application of metaheuristic algorithms to the identification of nonlinear magneto-viscoelastic constitutive parameters[END_REF][START_REF] Yunus | Thermal stability and rheological properties of epoxidized natural rubber-based magnetorheological elastomer[END_REF] . The finite strain range is a priori assumed due to the compliant nature of such materials that, in general, respond with large deformations. Herein, we adopt the multiplicative decomposition of the deformation gradient into an elastically relaxing part and a viscous part as nowadays widely used in the literature, see for example [START_REF] Govindjee | A presentation and comparison of two large deformation viscoelastic models[END_REF][START_REF] Hong | Modeling viscoelastic dielectrics[END_REF][START_REF] Nedjar | Framework for finite strain viscoelastic-plasticity based on multiplicative decompositions: Part I: continuum formulations[END_REF][START_REF] Nedjar | Framework for finite strain viscoelastic-plasticity based on multiplicative decompositions: Part II: computational aspects[END_REF][START_REF] Reese | A theory of finite viscoelasticity and numerical aspects[END_REF][START_REF] Sidoroff | Un modèle viscoélastique nonlinéaire avec configuration intermédiaire[END_REF] . This gives rise to an intermediate configuration on which the magnetic field vectors can be transported, and implicitly decomposed. Very few models exist in the literature that couple the magnetic field and viscoelasticity. Notice that recently in [START_REF] Saxena | A theory of finite deformation magneto-viscoelasticity[END_REF] , the magnetoviscoelasticity has also been formulated for incompressible materials based on an a priori decomposition of the magnetic induction vector, and all the continuum thermodynamic developments therein follow by starting from that choice. On the other hand, in [START_REF] Han | Coupled magnetic field and viscoelasticity of ferrogels[END_REF] , the authors did not take into account the magnetization neither for the instantaneous (equilibrium) nor for the elastically relaxed (non-equilibrium) parts of the response, reducing the magneto-viscoelasticy to its simplest form. However, the field theory developed in this paper fully couples the magnetic fields with finite viscoelasticity for compressible as well as for incompressible magnetizable materials in the spirit of an analogous development made in [START_REF] Nedjar | A finite strain modeling of electro-viscoelastic materials[END_REF] for finite strain electro-viscoelastic soft materials.

The coupling with magnetostatics is integrated within the framework of continuum thermodynamics for the correct setting of the whole set of constitutive laws together with the characterization of the dissipative phenomena through evolution equations. Herein, only isotropic MAPs are considered. In the formulation developed in this work, the magnetic induction vector is chosen to be the main independent variable for the magnetic part of the problem and, consequently, we show that the magnetic field becomes then implicitly decomposed additively into an equilibrium and a non-equilibrium parts. A simple model example is detailed for demonstrative purposes to highlight some phenomena qualitatively with parametric studies.

An outline of the remainder of this paper is as follows. In Section 2, we recall the governing equations of mechanical balance together with the specialized versions of Maxwell's equations in magnetostatics. Both of the equivalent spatial and material descriptions are considered. Then, in Section 3, the magnetoviscoelastic coupling is developed in detail and studied in Section 4 with parametric studies. Finally, conclusions and perspectives are drawn in Section 5.

Notations: Throughout the paper, bold face characters refer to second-or fourth-order tensorial quantities. In particular, 1 denotes the second-order identity tensor with components δ ij (δ ij being the Kronecker delta), and I is the fourth-order unit tensor of components I ijkl = 1 2 (δ ik δ jl + δ il δ jk ). The notation ( ) T is used for the transpose operator and the double dot symbol ':' is used for double tensor contraction, i.e. for any second-order tensors A and B, A : B = tr[AB T ] = A ij B ij where, unless specified, summation on repeated indices is always assumed. The notation ⊗ stands for the tensorial product. In components, one has (A ⊗ B) ijkl = A ij B kl , and for any two vectors Í and Î, (Í ⊗ Î) ij = U i V j . Herein, the double-striked characters will exclusively be used for vector fields related to the magnetic part of the problem, e.g. , , , Å . . . . Finally, the dot notation will always designate the material time derivative, i.e. ( ˙ ) ≡ d( )/dt.

Basic equations

When undeformed, unstressed, and in the absence of magnetic fields, the magnetically sensitive body occupies the reference configuration B 0 with boundary ∂B 0 . We identify a material particle by its position vector in the reference configuration, X ∈ B 0 , and trace its motion by its current position at time t in Prepared using sagej.cls the spatial configuration B t , x(X, t) ∈ B t . The deformation gradient is as usual defined as F = ∇ X x, where ∇ X ( ) is the material gradient operator with respect to the reference coordinates X. The Jacobian of the transformation is given by the determinant J = det F , with the standard convention J > 0.

Magnetostatics equations

The materials we consider being electrically non-conducting, the fields in magnetostatics are governed by the following specializations of Maxwell's equations in the absence of distributed currents:

curl = 0, div = 0, in B t , (1)
where and are respectively the magnetic field and magnetic induction vectors, both with respect to the spatial configuration. Inside the material, they are related by the standard relation

= µ 0 ( + Ñ), (2)
where the field vector Ñ is the spatial magnetization. The constant µ 0 is the magnetic permeability of vacuum. Note that in vacuo: = µ 0 . In (1), curl( ) and div( ) are respectively the rotational and divergence operators with respect to x. They are denoted by Curl( ) and Div( ) in the material configuration with respect to X. At the boundary of the body, in the absence of surface currents, the magnetic field and the magnetic induction vectors must satisfy the continuity conditions:

n × = 0, .n = 0, on ∂B t , (3)
where ( ) = ( ) outside -( ) inside is the jump across the boundary ∂B t of unit normal n. Notice that the continuity conditions (3) are also valid across a discontinuity within the body, see also [START_REF] Nedjar | A couple BEM-FEM method for finite strain magneto-elastic boundary-value problems[END_REF] .

Pull-back to the reference configuration gives the Lagrangian counterparts of the above laws (1):

Curl À = 0, Div = 0, in B 0 , (4)
for the Lagrangian-Eulerian interconnections

À = F T , = JF -1 , Å = F T Ñ. (5)
With these relations, the material version of the standard relation ( 2) is given by

= µ 0 JC -1 (À + Å) , (6)
Prepared using sagej.cls

where C = F T F is the right Cauchy-Green tensor which, otherwise, is a strain measure for the deformable solid. The counterpart of the continuity conditions (3), this time for À and , are

N × À = 0, .N = 0, on ∂B 0 , (7)
where N is the unit normal to the reference boundary ∂B 0 that corresponds to n through Nanson's formula nda = JF -T N dA, where da and dA are area elements on ∂B t and ∂B 0 , respectively.

Mechanical balance equations

In quasi-statics, the spatial form of the mechanical equilibrium equation together with the boundary condition can take the form div σ

+ m + ρf = 0 in B t , σn = t on ∂B t , (8)
where σ is the Cauchy stress tensor, ρ the current mass density, f the mechanical body force per unit mass, and m is the magnetic body force per unit current volume. The traction force t comprises both the mechanical traction per unit area of the boundary ∂B t , and an additional mechanical traction due to the magnetic effects. Notice that due to the magnetization, the spatial stress tensor σ is in general non-symmetric. Nevertheless, we adopt herein the nowadays well known structure where the stress is augmented with terms arising from the magnetic body force, e.g. see [START_REF] Ogden | Magnetostatics: from basic principles to nonlinear interactions in deformable media[END_REF] for details. The balance equation ( 8) is equivalently written as div σ

+ ρf = 0 in B t , σ n = 0 on ∂B t , (9)
in terms of the, this time symmetric, total Cauchy stress tensor σ defined as, see for example [START_REF] Dorfmann | Magnetoelastic modelling of elastomers[END_REF][START_REF] Dorfmann | Nonlinear magnetoelastic deformations of elastomers[END_REF] :

σ = σ + µ -1 0 ⊗ - 1 2 . 1 + Ñ. 1 -Ñ ⊗ , (10)
where the second-to-fourth terms on the right-hand side of (10) come from the magnetic body force, here given by m = (∇ x ) T Ñ, see for example [START_REF] Brigadnov | Mathematical modeling of magneto-sensitive elastomers[END_REF][START_REF] Pao | Electromechanic forces in deformable continua[END_REF] . Here and henceforth, the notation ∇ x ( ) refers to the spatial gradient operator with respect to the coordinates x. The boundary condition (9) 2 is here equivalently written in the form of a jump condition.

Pull-back of the balance equation (9) to the reference configuration gives the following useful Lagrangian form Div P + ρ 0 f = 0 in B 0 ,

P N = 0 on ∂B 0 , (11)
in terms of the total first Piola-Kirchhoff stress tensor P = J σF -T and the reference density ρ 0 = Jρ.

In particular for later use, the expression of P deduced from ( 10) is given by

P = P + µ -1 0 J -1 F ⊗ - 1 2 C : ⊗ F -T + Å. F -T -F -T Å⊗ , (12)
where P = JσF -T and, among others, use has been made of the relations (5) 2 and (5) 3 .

Continuum thermodynamics and constitutive laws

The above governing equations need to be supplemented with constitutive relations. For the magnetic part of the problem, we have at our disposal the three magnetic field vectors , À, Å, with the connection (6).

Any one of these can be used as the independent magnetic variable in the formulation of a constitutive law together with the characterization of the viscoelastic dissipation. In this work, we consider a formulation based on the use of the magnetic induction vector that we detail herein under the assumption of isothermal conditions. The balance of energy is written in the local form as, see [START_REF] Brigadnov | Mathematical modeling of magneto-sensitive elastomers[END_REF][START_REF] Dorfmann | Magnetoelastic modelling of elastomers[END_REF] for example,

ρ ė = div (σ T v) + ρf .v + w m , (13)
where e is the specific internal energy, per unit mass, v is the spatial velocity, and w m is the magnetic power given by w m = m .v -Ñ. ˙ . Using the mechanical balance (8), we get the reduced form

ρ ė = σ : ∇ x v -Ñ. ˙ , (14)
where the spatial velocity gradient ∇ x v is related to the deformation gradient through the well known

kinematic formula ∇ x v = Ḟ F -1 .
The second law of thermodynamics postulates the positiveness of the entropy production. Denoting by s the specific entropy per unit mass, we write ṡ ≥ 0 and, for isothermal processes, the dissipation D is given by

D = ρT 0 ṡ ≥ 0, (15)
where T 0 is the fixed absolute temperature. Now defining the specific free energy function ψ as ψ = e -T 0 s, using this latter together with (14) into (15) gives the classical form,

D = σF -T : Ḟ -Ñ. ˙ -ρ ψ ≥ 0. ( 16
)
Prepared using sagej.cls

Notice that in magneto-elasticity with a free energy function based on the use of the Eulerian magnetic induction as ψ(F , ), this leads to the well-known constitutive equations,

σ = ρ ∂ ψ ∂F F T , Ñ = -ρ ∂ ψ ∂ ,
for the Cauchy stress tensor and the magnetization vector. However, as the Lagrangian form is better suited for the following thermodynamic developments, inequality ( 16) must be pull-back to the reference configuration. Thus, we may regard the specific free energy function ψ, equivalently, as a function of F and , and we introduce the notation ψ for this purpose as

ψ(F , . . .) ≡ ψ(F , J -1 F . . .), (17)
where the dots mean further arguments that will eventually be introduced later on for the purpose of viscoelastic modelling. Thereby, using the relations σF -T = J -1 P , ρ = J -1 ρ 0 and

Ñ. ˙ = J -1 Å. ˙ -J -1 Å. F -T : Ḟ + J -1 F -T Å⊗ : Ḟ ,
where use has been made of (5) 2 and (5) 3 together with the well known kinematic relation J = JF -T : Ḟ , the inequality (16) becomes,

D = P + Å. F -T -F -T Å ⊗ : Ḟ -Å. ˙ -ρ 0 ψ ≥ 0. (18)
This latter can further be simplified by introducing the augmented free energy function Ω per unit reference volume defined as, see for example [START_REF] Dorfmann | Nonlinear magnetoelastic deformations[END_REF] ,

Ω(F , . . .) = ρ 0 ψ + 1 2 µ -1 0 J -1 C : ⊗ , (19)
which, replaced into (18), transforms this latter into the more convenient form

D = P : Ḟ + À. ˙ -Ω ≥ 0, (20)
where use has been made of the definition (12) for the total first Piola-Kirchhoff stress tensor P and of the relation (6) for the Lagrangian magnetic induction vector . Among others, observe that, in this form, the conjugate character between À and replaces the one between Å and in the former form (18). For similarities, see also [START_REF] Nedjar | A theory of finite strain magneto-poromechanics[END_REF] for developments in magneto-poromechanics.

Prepared using sagej.cls

The equal sign in (20) takes place only for non-dissipative processes, i.e. magneto-elasticity. In this case the thermodynamic equilibrium state is fully determined by the deformation gradient and the magnetic induction. For a viscoelastic material, however, the free energy in a non-equilibrium state differs from that in an equilibrium state. To distinguish between them, we introduce the equilibrium free energy ψ eq (F , ) so that from (20) in the case of an equal sign we get by thermodynamic arguments, see [START_REF] Bd | Thermodynamics with internal variables[END_REF] ,

P eq = ∂Ω eq ∂F , À eq = ∂Ω eq ∂ , (21)
where P eq and À eq are the total first Piola-Kirchhoff stress tensor and the Lagrangian magnetic field vector in an equilibrium state, and where Ω eq is the augmented volumetric free energy function at thermodynamic equilibrium defined as, see (19),

Ω eq (F , ) = ρ 0 ψ eq (F , ) + 1 2 µ -1 0 J -1 C : ⊗ . ( 22
)
The simple structures of the constitutive formulas ( 21) are exactly those derived for magneto-elastic modelling, see for example [START_REF] Ogden | Magnetostatics: from basic principles to nonlinear interactions in deformable media[END_REF][START_REF] Bustamante | On variational formulations in nonlinear magnetoelastostatics[END_REF] .

Kinematic assumptions and constitutive laws

The viscoelastic response is here captured by means of the nowadays well-accepted multiplicative split of the deformation gradient into an elastically relaxing part F e and a viscous part F v as illustrated in Fig. 1, for instance, see [START_REF] Reese | A theory of finite viscoelasticity and numerical aspects[END_REF][START_REF] Sidoroff | Un modèle viscoélastique nonlinéaire avec configuration intermédiaire[END_REF][START_REF] Lubliner | A model of rubber viscoelasticity[END_REF] among others:

F = F e F v . (23)
Now motivated by the decomposition of the energy in the purely mechanical theory, the nonequilibrium free energy that we denote here by ψ neq is the difference between the total free energy ψ and that of an equilbrium state ψ eq ; ψψ eq ≡ ψ neq . Herein, this latter is associated with the elastically relaxing part of the deformation gradient F e . It can moreover depend on the magnetic induction which, otherwise, is assumed to equilibrate much faster than deformation:

ψ neq ≡ ψ neq (F e = F F v -1 , = J v -1 F v ), (24)
where is the push-forward with F v of the Lagrangian magnetic induction vector to the intermediate configuration, see the illustration on the sketch of Fig. 1. It is equivalently the pull-back with F e of the spatial magnetic induction vector . Here J v = det F v and J e = det F e are the Jacobians of the viscous and the elastically relaxing parts of the deformation gradient. The precise forms of the above constitutive functions ψ eq and ψ neq are dictated by invariance and material symmetry that will be discussed below.

Prepared using sagej.cls , À, Å , , Ñ À = F v -T À = F e T Å = F v -T Å = F e T Ñ = J v -1 F v = J e F e -1 F F e F v X X x B 0 B t 1
With the additive structure of ψ, we choose for convenience to re-cast the augmented free energy Ω, Eq. ( 19), into the following form

Ω = Ω eq + ρ 0 ψ neq , (25)
where Ω eq is the augmented form already defined in (22). The constitutive equations are given by the derivatives of Ω as,

P = ∂Ω ∂F = P eq + ρ 0 ∂ψ neq ∂F e F v -T P neq , À = ∂Ω ∂ = À eq + J v -1 F v T ρ 0 ∂ψ neq ∂ À neq , (26)
Prepared using sagej.cls where P eq and À eq have been defined in [START_REF] Amin | Nonlinear dependence of viscosity in modeling the rate-dependent response of natural and high damping rubbers in compression and shear: Experimental identification and numerical verification[END_REF]. The remainder of the inequality (20) that governs the evolution of the internal variable, here F v , becomes

D = ρ 0 F e T ∂ψ neq ∂F e - ∂ψ neq ∂ ⊗ + . ∂ψ neq ∂ 1 : l v ≥ 0, (27)
where

l v = Ḟ v F v -1
is the left-rate of the viscous part of the deformation gradient, which is in general non-symmetric.

As expected, the total stress tensor is additively split into a part that is in thermodynamic equilibrium, herein written in the augmented form P eq , and an over-stress P neq that is thermodynamically not in equilibrium. For convenience, we adopt the same notation for the magnetic field vector À; with subscript notation ( ) eq for the term derived from the part Ω eq of the augmented free energy, and subscript notation ( ) neq for the term derived from the part ρ 0 ψ neq . This notation will henceforth be used for all the magnetic fields. For instance, one has then

= F -T À = F -T ∂Ω eq ∂ eq + J v -1 F e -T ρ 0 ∂ψ neq ∂ neq , (28)
for the spatial magnetic field vector, and

À = F v -T À ≡ F e T = F v -T ∂Ω eq ∂ À eq + J v -1 ρ 0 ∂ψ neq ∂ À neq , (29)
when transported to the intermediate configuration.

Now for objectivity reasons, Ω eq depends on the deformation gradient F only through the right Cauchy-Green tensor C, and for symmetry reasons, it depends on the magnetic induction vector only through the tensor product ⊗ . Likewise, the part ψ neq depends on F e only through the elastically relaxing right Cauchy-Green tensor C e = F e T F e ≡ F v -T CF v -1 , and on only through the tensor

product ⊗ ≡ J v -2 F v ⊗ F v T , both defined in the intermediate configuration, i.e. Ω ≡ Ω eq (C, ⊗ ) + ρ 0 ψ neq (C e , ⊗ ). ( 30
)
Prepared using sagej.cls

The above stress constitutive relation (26) 1 is then equivalently written as

S ≡ F -1 P = 2 ∂Ω eq ∂C S eq + F v -1 2ρ 0 ∂ψ neq ∂C e S neq F v -T S neq , (31)
for the total second Piola-Kirchhoff-type, where S neq = F v S neq F v T is the push-forward of its nonequilibrium part S neq to the intermediate configuration. The reduced dissipation (27) can alternatively be written as

D = C e S neq -J v À neq ⊗ + J v À neq . 1 : sym l v ≥ 0, (32)
where sym[ ] is the symmetric part operator of a second-order tensor. In fact, we will show later on that the above term between parentheses in ( 32) is indeed symmetric.

Characterization of the viscoelastic evolution

A local evolution equation must now be characterized such that the above inequality (32) is satisfied for any admissible process. One way to do so is to rewrite this latter in the equivalent spatial form by exploiting the kinematic relations

b e = F e F e T ≡ F C v -1 F T , with C v = F v T F v , (33)
where b e is the elastically relaxing left Cauchy-Green tensor and C v is the viscous right Cauchy-Green tensor. One obtains the following well-known relation:

sym l v = F e -1 - 1 2 £ v b e F e -T , (34)
where £ v b e = F Ċv -1 F T is the Lie derivative of b e , see for example [START_REF] Holzapfel | Nonlinear Solid Mechanics[END_REF][START_REF] Marsden | Mathematical Foundations of Elasticity[END_REF] for details. Then, replacing (34) into (32) gives the following spatial version:

D = τ neq -J neq ⊗ + J neq . 1 b e -1 : - 1 2 £ v b e ≥ 0, (35)
where use has been made of the push-forward transformations τ neq = F e S neq F e T , J e = F e and neq = F e -T À neq , together with the kinematic relation J = J e J v .

Prepared using sagej.cls Among many possibilities, the following general form can be adopted for the evolution equation:

- 1 2
£ v b e = M : τ neq -J neq ⊗ + J neq . 1 b e -1 , (36)
where the fourth-order mobility tensor M is positive definite so that the inequality (35), or equivalently (32), is satisfied automatically. This form is indeed often used in practice, see for example [START_REF] Hong | Modeling viscoelastic dielectrics[END_REF][START_REF] Boyce | On the kinematics of finite strain plasticity[END_REF][START_REF] Nedjar | An anisotropic viscoelastic fibre-matrix model at finite strains: Continuum formulation and computational aspects[END_REF] for similarities. In general, the mobility tensor may depend on various state variables. A possible physical meaning of the evolution law (36) is that one type of polymer chains relaxes in time and dissipates energy. This can be pictured by a dashpot in series with a spring. In the current state, the stress acting in the spring drives the viscous deformation in the dashpot, here the viscous strain tensor C v related to b e through (33) 1 . The key observation in this work is that the stress acting in the spring is herein coupled with stresses emanating from the magnetic effects on the intermediate configuration if the magnetic induction is accounted for in the expression of ψ neq , see (24).

A modelling example

Being isotropic functions of their respective arguments, each part of the free energy function ψ depends at most on a collection of six irreducible invariants. For the part related to the thermodynamic equilibrium ψ eq , and hence to Ω eq by ( 22), it depends on

I 1 = C : 1, I 2 = 1 2 I 2 1 -C : C , I 3 = det C = J 2 , I 4 = . , I 5 = C : ⊗ , and I 6 = C 2 : ⊗ , (37)
where, for the double tensor contraction, we recall that we use the convention stated in the Notations's paragraph of the introductory Section 1, for example

I 5 = C : ⊗ = C KL B K B L = F iK F iL B K B L .
Likewise, the non-equilibrium part ψ neq depends at most on

I 1 = C e : 1, I 2 = 1 2 I 2
see for example [START_REF] Spencer | Constitutive theory for strongly anisotropic solids[END_REF] for details on the topic of invariance. The first three invariants in each set (37) and (38) are classical in isotropic viscoelasticity. The other ones arise from the magneto-mechanics coupling in the thermodynamically equilibrated as well as non-equilibrated states for the most general case.

Prepared using sagej.cls

The following general forms for the parts Ω eq and ψ neq can be adopted as a modelling example:

Ω eq = ρ 0 ψ ′ eq (I 1 , I 2 , I 3 ) + c 1 I 4 + c 2 I 5 + c 3 I 6 magnetic coupling + 1 2 µ -1 0 J -1 C : ⊗ augmentation , ρ 0 ψ neq = ρ 0 ψ ′ neq (I 1 , I 2 , I 3 ) + c 4 I 4 + c 5 I 5 + c 6 I 6 magnetic coupling , (39)
where c 1 , c 2 , . . . , c 6 are material parameters relative to the magnetic coupling. The first volumetric terms ρ 0 ψ ′ eq and ρ 0 ψ ′ neq are related to the purely mechanical part. Furthermore, with the form (39) 2 , one can easily check that the term between parentheses in (32) is symmetric.

Applying the above forms (39) to Eq. ( 31), we obtain the following relation for the total stress tensor of the Kirchhoff-type τ = F SF T ≡ J σ:

τ = τ ′ eq + τ ′ neq + 2c 2 J 2 ⊗ + 2c 3 J 2 ⊗ b + b ⊗ + 2c 5 J e 2 ⊗ + 2c 6 J e 2 ⊗ b e + b e ⊗ + µ -1 0 J ⊗ - 1 2 . 1 , (40)
where we have introduced for convenience the notations

τ ′ eq = 2ρ 0 F ∂ψ ′ eq ∂C F T , and τ ′ neq = 2ρ 0 F e ∂ψ ′ neq ∂C e F e T , (41)
for the purely mechanical parts. Here use has been made of the transformations (5) 2 and J e = F e for the spatial induction vector , and b = F F T is the well known (spatial) left Cauchy-Green tensor. The equilibrium part of the total stress in (41) agrees with the most general form established in [START_REF] Dorfmann | Magnetoelastic modelling of elastomers[END_REF] for the case of compressible magnetoelasticity. One can also observe that a magnetostatic contribution appears in the inelastic stress through the material parameters c 5 and c 6 .

For the spatial magnetic field vector, we have from ( 39), (26) 2 and (28):

= 2c 1 J b -1 + 2c 2 J + 2c 3 J b = -Ñ eq + J v -1 J e 2c 4 b e -1 + 2c 5 + 2c 6 b e = -Ñ neq ≡ neq + µ -1 0 . ( 42
)
Prepared using sagej.cls

Comparing the latter expressions (42) with the corresponding standard relation ( 2), one concludes that c 1 , c 2 , . . . , c 6 are in fact the magnetization parameters, among which, c 4 , c 5 and c 6 are related to the intermediate configutation.

Now for the purpose of demonstration, a simple choice for the mobility tensor in ( 36) could be M = 1/η I where η > 0 is a viscosity parameter of dimension [N m -2 s] that can be assumed constant in a first approach. The evolution equation ( 36) would then become

-1 2 £ v b e = 1 η τ ′ neq b e -1 -2c 4 J e 2 b e -1 ⊗ b e -1 + 2c 6 J e 2 ⊗ +2 c 4 I 4 + c 5 I 5 + c 6 I 6 b e -1 , (43)
after having replaced the expressions of τ neq and neq from Eqs. ( 40) and ( 42), respectively, and simplifying terms. Among others, notice that both sides of ( 43) are symmetric. However, it should be stressed that in general the viscosity parameter η is not constant, e.g. see for example [START_REF] Amin | Nonlinear dependence of viscosity in modeling the rate-dependent response of natural and high damping rubbers in compression and shear: Experimental identification and numerical verification[END_REF][START_REF] Zhou | A micro-macro constitutive model for finite-deformation viscoelasticity of elastomers with nonlinear viscosity[END_REF] .

Remark 1. If the energy dissipation due to the magnetization is negligible, then by setting c 4 = c 5 = c 6 = 0, the model drastically simplifies as follows:

• The part of the the magnetic field vector relative to the intermediate configuration would vanish as

À neq = neq = 0.
• The over-stress would reduce to its purely mechanical part as τ neq = τ ′ neq for the Kirchoff-type stress tensor.

• The reduced dissipation (35) would take its classical form encountered in finite isotropic viscoelasticity.

• If, furthermore c 1 = c 2 = c 3 = 0, then the magnetization vector field vanishes as Å = 0. In this case the present modelling reduces to that of non-magnetizable magneto-viscoelasticity as developed in [START_REF] Han | Coupled magnetic field and viscoelasticity of ferrogels[END_REF] .

Homogeneous deformation and parametric studies

In this section, we will construct a simple model to study some phenomena qualitatively. To this end, it remains now to specify the purely mechanical free energy functions ψ ′ eq and ψ ′ neq so that the augmented free energy function Ω will be completely defined, see (30) and (39). Here for the purpose of demonstration, we choose a compressible Neo-Hookean type modelling for both functions as

ρ 0 ψ ′ eq = 3 8 κ eq J 4 3 + 2J -2 3 -3 + 1 2 µ eq J -2 3 I 1 -3 , ρ 0 ψ ′ neq = 3 8 κ neq J e 4 3 + 2J e -2 3 -3 + 1 2 µ neq J e -2 3 I 1 -3 , (44)
Prepared using sagej.cls

where the first terms are related to the volumetric response with κ eq and κ neq the equilibrium and nonequilibrium bulk moduli, respectively, while the second terms are related to the isochoric part of the response with µ eq and µ neq as the equilibrium and the non-equilibrium shear moduli, respectively, e.g. see also [START_REF] Nedjar | On constitutive models of finite elasticity with possible zero apparent Poisson's ratio[END_REF][START_REF] Nedjar | A finite element implementation of the isotropic exponentiated henckylogarithmic model and simulation of the eversion of elastic tubes[END_REF] for more details about the properties of the hyperelastic version of the model (44). With these forms, the purely mechanical parts of the Kirchhoff-type stress tensor are given by, see Eq. ( 40),

τ ′ eq = 1 2 κ eq J 4 3 -J -2 3 1 + µ eq dev[J -2 3 b], τ ′ neq = 1 2 κ neq J e 4 3 -J e -2 3 1 + µ neq dev[J e -2 3 b e ], (45)
where dev[ ] denotes the spatial deviatoric operator.

As a simple model problem, we study the response of a cylindrical sample to applied magnetic induction along its axial direction e 3 , i.e. (t) = B 3 (t) e 3 . For simplicity, we assume the deformation to be homogeneous. Let λ 1 , λ 2 and λ 3 be the principal stretches along the three directions e 1 , e 2 and e 3 , as the loading direction does not change in the history, the principal directions of the viscoelastic stretches coincide with those of the above total ones, and we will denote them as λ v 1 , λ v 2 and λ v 3 Furthermore, here for symmetry reasons, we also have

λ 2 = λ 1 and λ v 2 = λ v 1 .
In components, the constitutive relation (40) for the total Kirchhoff-type stress tensor τ = J σ specializes to

J σ 11 = 1 2 κ eq λ 8 3
1 λ

4 3 3 -λ -4 3 1 λ -2 3 3 + 1 3 µ eq λ 2 3
1 λ

-2 3 3 -λ -4 3 1 λ 4 3 3 + 1 2 κ neq λ e 8 3
1 λ e

4 3 3 -λ e -4 3 1 λ e -2 3 3 + 1 3 µ neq λ e 2 3
1 λ e -2 3

3

λ e -4 3

1

λ e 4 3 3 - 1 2µ 0 B 2 3 λ -2 1 λ 3 , J σ 33 = 1 2 κ eq λ 8 3
1 λ

4 3 3 -λ -4 3 1 λ -2 3 3 + 2 3 µ eq λ -4 3 1 λ 4 3 3 -λ 2 3
1 λ

-2 3 3 + 1 2 κ neq λ e 8 3
1 λ e 1 λ e -2 3

3

+ 1 2µ 0 B 2 3 λ -2 1 λ 3 + 2c 2 B 2 3 λ 2 3 + 4c 3 B 2 3 λ 4 3 + 2c 5 B 2 3 λ e 4 1 λ e 2 3 λ -4 1 + 4c 6 B 2 3 λ e 4 1 λ e 4 3 λ -4 1 , (46)
with identical transverse stresses σ 22 = σ 11 , and where λ e 1 = λ 1 λ v -1

1

≡ λ e 2 and λ e 3 = λ 3 λ v -1 3 are the elastically relaxing stretches.

Prepared using sagej.cls

Likewise, for the viscous deformation, the evolution equation ( 43) reduces in components to

λ e 4 1 log λ v 1 = 1 2 κ neq η λ e 8 3
1 λ e 4 3

3 -λ e -4 3

1 λ e -2 3 3 + 1 3 µ neq η λ e 2 3
1 λ e -2 3

3

λ e -4 3

1 1 λ e -2 3 3

λ e 4 3 3 + 2 c 4 η B 2 3 λ e 4 1 λ -4 1 + 2 c 5 η B 2 3 λ e 4 1 λ e 2 3 λ -4 1 + 2 c 6 η B 2 3 λ e 4 1 λ e 4 3 λ -4 1 , λ e 4 3 log λ v 3 = 1 2 κ neq η λ e
+ 2 c 5 η B 2 3 λ e 4 1 λ e 2 3 λ -4 1 + 4 c 6 η B 2 3 λ e 4 1 λ e 4 3 λ -4 1 . (47)
For the following analyses, we introduce the ratio β = κ neq /µ neq between the non-equilibrium bulk and shear moduli. Herein, the particular case of quasi-incompressibility in the viscous part of the response can be obtained with a relatively high value of β.

Remark 2. The (strong) hypothesis of homogeneous deformation is set here only for demonstrative purposes to focus on the intrinsic material response. At this stage, the applied magnetic induction is assumed given and the study reduces to that of a material point. Of course, in more realistic applications, the surrounding space together with the continuity conditions, i.e. ( 3) and (9) 2 , must be accounted for in a more complete structural analysis, see for example [START_REF] Nedjar | A couple BEM-FEM method for finite strain magneto-elastic boundary-value problems[END_REF] .

Magnetic induction with no deformation

In this first example, we consider that the above sample is held fixed with λ 1 = λ 2 = λ 3 = 1 during the whole loading history and we study the effect of an applied induction B 3 (t) on the purely viscoelastic part of the response, see the sketch of Fig. 2. To simplify the mathematical expressions, we introduce the following dimensionless magnetization parameters and variables:

α i = c i µ 0 i=1,...,6 , B = B 3 √ µ 0 µ neq , H = H 3 µ 0 µ neq , ( 48
)
where H 3 is the e 3 -component of the corresponding magnetic field. From (47), one can observe that the viscoelastic stretches evolve with the characteristic time η/µ neq . This latter is herein used to normalize the time t by introducing the dimensionless time t = t µ neq /η. With the constraints λ 1 = 1 and λ 3 = 1, Prepared using sagej.cls the system (47) becomes

λ v -4 1 d d t (log λ v 1 ) = β 2 λ v -8 3 1 λ v -4 3 3 -λ v 4 3 1 λ v 2 3 3 + 1 3 λ v -2 3 1 λ v 2 3 3 -λ v 4 3 1 λ v -4 3 3 + 2α 4 B 2 λ v -4 1 + 2α 5 B 2 λ v -4 1 λ v -2 3 + 2α 6 B 2 λ v -4 1 λ v -4 3 , λ v -4 3 d d t (log λ v 3 ) = β 2 λ v -8 3 1 λ v -4 3 3 -λ v 4 3 1 λ v 2 3 3 + 2 3 λ v 4 3 1 λ v -4 3 3 -λ v -2 3 1 λ v 2 3 3 + 2α 5 B 2 λ v -4 1 λ v -2 3 + 4α 6 B 2 λ v -4 1 λ v -4 3 , (49)
which is solved for λ v 1 ( t) and λ v 3 ( t), and by mere function evaluation using these latter into (46), one can also deduce the corresponding dimensionless total stress components σ 11 /µ neq and σ 33 /µ neq . Some representative solutions to the differential-algebraic system (49) are plotted in Fig. 3 for the case where the dimensionless induction is applied quickly at a relatively high rate dB/d t = 100 until it reaches the value B = 1, and then is left fixed at this value. At dimensionless time t = 0, the sample is fully relaxed with the initial conditions being λ i = λ v i = 1, i = 1, 2, 3. Here as the system (49) is independent on the equilibrium magnetization parameters, we choose to set them to zero as α 1 = α 2 = α 3 = 0, so that the dimensionless magnetic field (48) 3 is linked to the dimensionless induction by the relation

H = 2α 4 λ v 4 1 + 2α 5 λ v 4 1 λ v 2 3 + 2α 6 λ v 4 1 λ v 4 3 + 1 B. (50)
We study the influence of the non-equilibrium magnetization parameters α 4 , α 5 and α 6 (or equivalently of c 4 , c 5 and c 6 ) by choosing to activate only one of them each time. Each case is computed with three different ranges of compressibility in the mechanical viscoelastic part through the ratio β. We choose the Prepared using sagej.cls values β = 2/3, that corresponds to a zero Poisson's ratio in the limiting case of linearized kinematics, β = 4.67 that corresponds to an approximate 0.4 Poisson's ratio, also in the limiting case of linearized kinematic, and β = 500 that corresponds to a viscoelastic quasi-incompressibility within the whole strain range.

One can see from Figs. 3(a)(c)(e) that, in all cases, the viscoelastic compressibility strongly influences the decay to equilibrium of the magnetic field. For the augmented Cauchy stress, Fig. 3(b) shows an unexpected evolution with no decay. This is explained by the fact that the stress is independent of the parameter α 4 , or c 4 , see Eq. (46) 2 . On the other side, for the cases where α 5 and/or α 6 are activated, it is seen from Figs. 3(d)(f) that the stress decays to an equilibrium value as expected. This decay is particularly sharp for α 5 = 1 in the limiting case of incompressibility, i.e. the curve with β = 500 in Fig. 3(d), and likewise for the corresponding magnetic field's decay in Fig. 3(c).

Responses to cyclic loads

In this second example, the dynamic responses of the above cylindrical sample are studied. The sample is this time unconstrained with free lateral and axial stretches and we still stipulate the homogeneous assumption for the purpose of demonstration and qualitative prediction, see Fig. 4. Initially, the sample is fully relaxed at t = 0, with the initial conditions λ i = λ v i = 1, i = 1, 2, 3. In this case four equations constitute the differential-algebraic system to be solved for the stretches λ 1 , λ v 1 , λ 3 and λ v 3 : the two issued from the mechanical balance, σ 11 = 0 and σ 33 = 0 in the absence of mechanical body forces, and the two evolution equations (47). The expressions of the formers are given in (46). Prepared using sagej.cls value at stabilization seems to be the same. As a last remark, it can be seen from Figs. 5 (b) and (d) that, with α 5 , the mean stretching λ v 3 is higher than the mean total one λ 3 , while with α 6 , the mean λ v 3 is lower that λ 3 . We remember that sample experiences contraction in the e 3 -direction under the applied magnetic load, hence the stretches λ 3 and λ v 3 are lower than 1.

Conclusion and perspectives

In this paper, we have presented a general theory to model the nonlinear magneto-viscoelasticity in the finite strain range. Within the framework of continuum thermodynamics, the nowadays well know fields related to the magnetic coupling have been embedded in a sound way for a concise characterization of the whole set of constitutive equations.

The viscoelastic kinematics adopted herein is based on the multiplicative decomposition of the deformation gradient into an elastically relaxing part and a viscous part. The formulation developed in this paper is based on the magnetic induction vector as main independent variable for the magnetic part of the coupled problem and, to matters as concrete as possible, a model problem has been proposed for demonstrative purposes so as to conduct preliminary parametric studies to qualitatively show some basic phenomena revealed by the formulation.

We believe that further work has to be accomplished to optimize the present theory toward more realistic modelling of smart magneto-active materials. Experimental investigations will certainly give better knowledge of the material properties. Among others, the magnetization parameters must be quantified. Experimental tests could for example show which ones are most relevant for the modelling or, at least, whether or not the non-equilibrium ones are pertinent. On the other hand, as most of the MAPs are manufactured as anisotropic composites, the present modelling framework must be extended by taking into account the directional dependencies of the magneto-viscoelastic responses.

Figure 1 .

 1 Figure 1. Illustration of the multiplicative decomposition F = F e F v and transport of the magnetic field vectors on the intermediate configuration.

1 -

 1 C e : C e , I 3 = det C e = J e 2 , I 4 = . , I 5 = C e : ⊗ , and I 6 = C e 2 : ⊗ ,

Figure 2 .

 2 Figure 2. Fixed cylindrical sample and magnetic induction loading.

1 Figure 3 .

 13 Figure 3. (a) (c) (e) Dimensionless magnetic field evolution H( t) and corresponding (b) (d) (f) dimensionless total stress σ33( t)/µneq. In each case, only one non equilibrium magnetization parameter is activated: α4 for (a-b), α5 for (c-d), and α6 for (e-f). Three different ratios β = 2/3, 4.67, and 500 are used for each calculation. Prepared using sagej.cls

Figure 4 .

 4 Figure 4. Free cylindrical sample and sinusoidal magnetic induction loading.

1 Figure 5 .

 15 Figure 5. Total and viscoelastic stretches λ3( t) and λ v 3 ( t) in response to applied sinusoidal magnetic induction with frequency ω = 1 in (a) (c) (e), and ω = 10 in (b) (d) (f). Each time with the activation of different magnetization parameters, as indicated.

Prepared using sagej.cls

Here again, we simplify the mathematical expressions with an adimensionalization procedure as follows. While conserving the definitions (48) 1 for the dimensionless magnetization parameters α i , i = 1 . . . 6, the dimensionless variables are this time defined as

where µ = µ eq + µ neq is the instantaneous shear modulus of the material and, to simplify further, we also consider the same ratio between the couples of bulk and shear moduli of the equilibrium and nonequilibrium parts, i.e.

Now defining ζ = µ eq /µ ∈ [0, 1] as the ratio between the equilibrium and instantaneous shear moduli, we deduce the following relations for the material parameters that will also be used when adimensionalizing the equilibrium equations σ 11 = 0 and σ 33 = 0 by µ:

Last, as for the example in Section 4.1, the evolution equations ( 47) are here again normalized with the adimensional time t = t µ neq /η.

We write the applied dimensionless magnetic induction into a sinusoidal function B = B 0 sin ω t, with ω and B 0 being the dimensionless frequency and amplitude, respectively. A set of representative results are presented in Fig. 5. In all the solutions plotted, the amplitude is set to be B 0 = 2, we consider a quasi-incompressible material with β = 500, and the material parameter ζ is taken to be 0.5 so that the equilibrium and non-equilibrium contributions are comparable. Hence, we study the influence of the magnetization parameters α i , i = 1 . . . 6, together with the driving induction frequency ω.

For illustrative purposes, we choose here again to activate only one non-equilibrium magnetization parameter each time; α 5 in Figs. 5 (a-b), α 6 in Figs. 5 (c-d), and a combination of this latter with the activation of an equilibrium parameter, for instance, here with α 3 in Figs. 5 (e-f). In the absence of further investigations, we prefer to not consider the non-equilibrium parameter α 4 due to its unexpected influence seen in the first example, Fig. 3 (b). For the sake of comparison, two frequencies are used; a moderate one with ω = 1, and a relatively high one with ω = 10.

The first observation is that, for the moderate frequency, the stretches stabilize in 1-2 cycles, while for the relatively high frequency, the stabilization takes place after 10-12 cycles. In each case, the amplitude of the viscoelastic stretch λ v 3 is markedly smaller than that of λ 3 . Moreover, for each set of parameters, the amplitude of λ v 3 is markedly smaller for the high frequency than for the moderate one, but the mean Prepared using sagej.cls


Vincent Fréour
email: vincent.freour@music.yamaha.com

Hideyuki Masuda

Satoshi Usa

Eiji Tominaga

Yutaka Tohgi

Bruno Cochelin

Christophe Vergez

CONTINUATION OF A PHYSICAL MODEL OF BRASS INSTRU-MENT FOR TRUMPET COMPARISONS

Keywords: brass instrument, physical modeling, continuation, Asymptotic Numerical Method

The system formed by a trumpet player and his/her instrument can be seen as a non-linear dynamical system, and modeled by a biomechanical model of the lips, a non-linear flow equation, and an acoustic impedance of the instrument represented by a superposition of complex modes. Different numerical tools can then be used to study these models in order to understand the effects of the parameters on the behavior of the system. In particular, the modal parameters that describe the impedance are strongly dependent on the geometry and design of the instrument, and the understanding of their influence is therefore of primary interest for a musical instrument maker. An original method of continuation (path-following method) based on the combination of the Harmonic Balance Method (HBM) and the Asymptotic Numerical Method (ANM) is applied to a physical model of the {player-trumpet} system. It consists in studying the evolution of the system where one parameter of the model varies, and enables to calculate bifurcation diagrams for different configurations of the model. Especially, this allows to compare different instruments on the basis of quantitative descriptors computed from the continuation outputs.

Introduction

The system formed by the couple {player -trumpet} falls into the class of non-linear dynamical systems likely to be studied using different numerical tools such as numerical continuation methods. As a matter of fact, continuation of a {player -trumpet} model can be of much interest for precise comparisons of different instruments. Continuation of periodic solutions using the ANM (Asymptotic Numerical Method) has already been applied to wind instruments [START_REF] Karkar | Oscillation threshold of a clarinet model: a numerical continuation approach[END_REF][START_REF] Velut | Contrôle par le musicien des régimes d'oscillation des instruments de la famille des cuivres : modélisation et mesures acoustiques, analyse du système dynamique[END_REF]. Nevertheless the potential of such method for instrument comparisons is still to be determined. In this study we illustrate the interest of this approach for the categorization of Bb trumpets in the space of some performance descriptors obtained from continuation by the ANM method combined to the Harmonic Balance Method (HBM).

Physical model of the player-trumpet system

The model considered is based on one-dimensional lip model, coupled to the resonator impedance described by a series of complex modes similar to what proposed in [START_REF] Silva | MoReeSC: A Framework for the Simulation and Analysis of Sound Production in Reed and Brass Instruments[END_REF]. The coupling between the mechanical oscillator and the acoustic resonator is achieved by a Bernoulli flow equation, considering turbulent mixing in the mouthpiece with no pressure recovery [START_REF] Elliot | Regeneration in brass instruments[END_REF]. The mechanical and acoustic equations are given in system 1, where y is the vertical lip position (y 0 is the lip position at rest), ω l , Q l , µ l and b the mechanical lip parameters, s k and C k with k ∈ [1, N ] the modal parameters of the N resonances of the acoustic impedance of the instrument, Z c the characteristic impedance, u the volume flow, p the downstream pressure at the input of the instrument (in the mouthpiece), and p 0 the upstream (mouth) static pressure.

ÿ(t) + ω l Q l ẏ(t) + ω l 2 (y(t) -y 0 ) = 1 µ l (p 0 -p(t)) ṗk (t) = Z c C k u(t) + s k p k (t), ∀k ∈ [1, N ] (1)
with

p(t) = 2 N k=1 (p k (t)) and u = 2|p 0 -p| ρ b • sign(p 0 -p) • θ(y)
, where θ(y) = |y|+y 2 , b is the lip width and ρ is the air density.

The case of a negative opening of the lips is managed by introducing the Heaviside function θ(y). The modal parameters of the N modes of the impedance are extracted from measured impedances, using the high resolution method ESPRIT [START_REF] Roy | Esprit: Estimation of signal parameters via rotational invariance techniques[END_REF] and the ESTER criteria for signal enumeration [START_REF] Badeau | Esprit: A new perturbation analysis for signal enumeration in rotational invariance techniques[END_REF].

Continuation by ANM

The system 1 can be recast under the canonical form Ẋ = f (X, λ), where X is a vector of unknown, and λ one parameter of the problem (e.g. the dimensionless mouth pressure for a brass instrument model).

Adopting the notation U = (X, λ), and looking for periodic solutions represented by their Fourier series, the resulting algebraic systems writes as R(U) = 0.

The aim of a continuation method is to find the solutions of R(U) = 0. These solutions are the union of solution branches and form a bifurcation diagram. These branches may cross at bifurcation points.

Denoting U 0 = (u 0 , λ 0 ) a solution of R(U) = 0, the Asymptotic Numerical Method (ANM) is based on the search of solutions expressed as series expansions with respect to a path parameter a.

U(a) = U 0 + aU 1 + a 2 U 2 + a 3 U 3 + ... + a p U p , (2)
where p si the truncation order of the series. For details about the fundamentals of ANM, the reader is invited to refer to [START_REF] Cochelin | Méthode asymptotique numérique[END_REF][START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF].

The specificity of the approach proposed in this paper is based on combining the ANM method with the Harmonic Balance Method (HBM) for the search of periodic solutions of the system. The HBM allows to approximate the unknowns by truncated Fourier series. The new unknowns of the problem are the Fourier coefficients of each element of U. For more details about the continuations of periodic solutions using ANM and HBM, the reader is invited to refer to [START_REF] Cochelin | A high-order purely frequency based harmonic balance formulation for continuation of periodic solutions[END_REF].

Bifurcation diagram and trumpet categorization

The calculation of bifurcation diagrams of system 1 is performed using the Matlab library MANLAB developed at LMA1 . The natural frequency of the lips f l = 2πω l is set to excite the fourth regime of the instrument in open fingering (no valve pressed). In a first step the stationary solution of the system is calculated in order to identify a Hopf bifurcation from which a periodic solution emerges. This periodic solution is then followed by continuation using the method described previously. A bifurcation diagram obtained for a Bb4 on a Bb trumpet is represented in Fig. 1.

The obtained diagram is characterized by an inverse bifurcation that induces an hysteresis H. A quantity D 1 associated to the dynamic range of p can also be defined from the stable part of the solution. By extracting H and D 1 on different trumpets, and locating the corresponding points in the (H, D 1 ) 2D space, the categorization represented in Fig. 1 is obtained. For each instrument, two impedance measurements are used for calculation. The categorization obtained clearly allows to differentiate instruments in the 2D space, showing the ability of the method to provide discriminating descriptors.

Conclusions

The application of the Asymptotic Numerical Method (ANM), coupled to the Harmonic Balance Method (HBM) as implemented in MANLAB library, allows for computation of periodic solutions of physical models of trumpets with respect to the driving static mouth pressure. Descriptors obtained from the calculated diagrams show interesting potential in categorizing instruments. Nevertheless, the relation of this objective categorization to perceptual metrics (obtained from human players) remains to be assessed.

Figure 1 :

 1 Figure 1: Left: bifurcation diagram of a Bb4 for a Bb trumpet: amplitude of the mouthpiece pressure p as a function of the static mouth pressure p 0 . The dotted line indicates unstable portions of the branch, while the solid line indicates the stable branch. Right: categorization of trumpets in the (H, D 1 ) space. The different colors correspond to the different trumpets. To each trumpet, two points are associated, corresponding to two impedance measurements of the instrument.

http://manlab.lma.cnrs-mrs.fr/

Acknowledgments

The authors would like to thank Louis Guillot for his precious support on the use of MANLAB.


V Fréour

H Masuda

S Usa

E Tominaga

Y Tohgi

B Cochelin

C Vergez

Continuation of a physical model of brass instrument: application to trumpet categorization

The system formed by the couple {player -trumpet} falls into the class of non-linear dynamical systems likely to be studied using different numerical tools such as numerical continuation methods. In this study we illustrate the interest of this approach for the categorization of Bb trumpets in the space of some performance descriptors obtained from continuation by the ANM method combined to the Harmonic Balance Method (HBM).

The model considered is based on one-dimensional lip model, coupled to the resonator impedance described by a series of complex modes similar to what proposed in [START_REF] Silva | MoReeSC: A Framework for the Simulation and Analysis of Sound Production in Reed and Brass Instruments[END_REF]. The coupling between the mechanical oscillator and the acoustic resonator is achieved by a Bernoulli flow equation, considering turbulent mixing in the mouthpiece with no pressure recovery [START_REF] Elliot | Regeneration in brass instruments[END_REF]. The mechanical and acoustic equations are given in system 1, where y is the vertical lip position (y 0 is the lip position at rest), ω l , Q l , µ l and b the mechanical lip parameters, s k and C k with k ∈ [1, N ] the modal parameters of the N resonances of the acoustic impedance of the instrument, Z c the characteristic impedance, u the volume flow, p the downstream pressure at the input of the instrument (in the mouthpiece), and p 0 the upstream (mouth) static pressure.

ÿ(t) + ω l Q l ẏ(t) + ω l 2 (y(t) -y 0 ) = 1 µ l (p 0 -p(t)) ṗk (t) = Z c C k u(t) + s k p k (t), ∀k ∈ [1, N ] (1)
with

p(t) = 2 N k=1 (p k (t)) and u = 2|p 0 -p| ρ b • sign(p 0 -p) • θ(y)
, where θ(y) = |y|+y 2 , b is the lip width and ρ is the air density.

The case of a negative opening of the lips is managed by introducing the Heaviside function θ(y). The modal parameters of the N modes of the impedance are extracted from measured impedances, using the high resolution method ESPRIT [START_REF] Roy | Esprit: Estimation of signal parameters via rotational invariance techniques[END_REF].

The specificity of the approach proposed in this paper is based on combining the ANM method [START_REF] Cochelin | Méthode asymptotique numérique[END_REF][START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF] with the Harmonic Balance Method (HBM) for the search of periodic solutions of the system. The HBM allows to approximate the unknowns by truncated Fourier series. The new unknowns of the problem are the Fourier coefficients of each element of U. For more details about the continuations of periodic solutions using ANM and HBM, the reader is invited to refer to [START_REF] Cochelin | A high-order purely frequency based harmonic balance formulation for continuation of periodic solutions[END_REF].

The calculation of bifurcation diagrams of system 1 is performed using the Matlab library MANLAB developed at LMA1 . The natural frequency of the lips f l = 2πω l is set to excite the fourth regime of the instrument in open fingering (no valve pressed). In a first step the stationary solution of the system is calculated in order to identify a Hopf bifurcation from which a periodic solution emerges. This periodic solution is then followed by continuation using the method described previously. A bifurcation diagram obtained for a Bb4 on a Bb trumpet is represented in Fig. 1. The obtained diagram is characterized by an inverse bifurcation that induces an hysteresis H. A quantity D 1 associated to the dynamic range of p can also be defined from the stable part of the solution. By extracting H and D 1 on different trumpets, and locating the corresponding points in the (H, D 1 ) 2D space, the categorization represented in Fig. 1 is obtained. For each instrument, two impedance measurements are used for calculation. The categorization obtained clearly allows to differentiate instruments in the 2D space, showing the ability of the method to provide discriminating descriptors.

Figure 1 :

 1 Figure 1: Left: bifurcation diagram of a Bb4 for a Bb trumpet: amplitude of the mouthpiece pressure p as a function of the static mouth pressure p 0 . The dotted line indicates unstable portions of the branch, while the solid line indicates the stable branch. Right: categorization of trumpets in the (H, D 1 ) space. The different colors correspond to the different trumpets. To each trumpet, two points are associated, corresponding to two impedance measurements of the instrument.

http://manlab.lma.cnrs-mrs.fr/

Acknowledgments

The authors would like to thank Louis Guillot for his precious support on the use of MANLAB.