domain: spi info

text:
Yosra Ben Saied
email: yosrabensaied@yahoo.fr

Alexis Olivereau

A Lightweight Threat Detection System for Industrial Wireless Sensor Networks

Keywords: Energy efficiency, hybrid approach, intrusion detection, resource constraints, wireless sensor network

Wireless Sensor Networks (WSNs) led the way to new forms of communications, which extend today the Internet paradigm to unforeseen boundaries such as eHealth, intelligent building or smart grid, to name a few. The legacy industry, however, is slower to adopt this technology, mainly for security reasons. Self-managed security systems allowing a quicker detection and better resilience to attacks, may counterbalance this reluctance. In this article, we propose a hybrid threat detection system that involves both centralized decision and local, per-cluster, work repartition and that is designed to run on top of industrial wireless sensor networks. Compared with the literature, we prove that this system is suitable for architectures mainly composed of constrained and sleeping devices, for which it achieves a fair level of autonomous security without prohibitively drawing out energy resources.

Introduction

Wireless Sensor Networks (WSNs) changed the Internet communication paradigm by introducing unattended devices in a hitherto human-centric model. WSNs filled the gap between the physical world and the Internet by offering the ability to monitor in real time a wide range of physical values through connected devices.

Initially built around very simple entities, technologies and architectures, WSNs benefitted from a strong interest from the international research community, which designed in a few years a wide variety of new components and systems. Dedicated radio technologies and routing protocols were specified, as well as novel, complex architectures. These architectures evolved from a basic, unidirectional topology made of multiple sensor nodes pushing data to a server through a sink node, to more advanced systems featuring bi-directionality. These new systems paved the way towards the emerging Machine to Machine (M2M) and Internet of Things (IoT) frameworks.

This gain of maturity of WSN technologies accelerates their adoption in the industry, and this adoption is all the quicker as WSNs answer to classical needs of industrial scenarios: monitoring of physical values, asset supervision and facilities surveillance are all key requirements in these scenarios, for which dedicated sensor nodes are available. Introducing actuators along with sensor nodes allows for more complex operations such as asset control and better production chain automation. Finally, the addition of more complex entities, such as indoor positioning nodes and mobile devices, to an existing sensor topology can give rise to new services, also profitable for the industry. For example, better worker safety can be achieved if the production chain is able to detect a worker's presence in an unauthorized area, or to determine that a certain user is not wearing adequate protection suit.

However, even though cost effective devices and energy-efficient technologies and protocols are available, the underlying security question impedes the use of WSNs in the most critical industrial scenarios. The inherent vulnerability of WSN nodes, which is due to their exposed location and their use of wireless communications, is such that a WSN has to mimic all security features from the legacy Internet, while also adding specific use cases and taking into account the strong shortcomings of the WSN nodes. Yet, strong security is often contradictory with limited resources. Furthermore, the unattended nature of WSNs and their autonomous decision taking upon a context change put them at risk of triggering harmful behavior, if misled by an attacker.

We propose to leverage on this WSN autonomy to provide it with the ability to detect new threats and react to them in the most appropriate manner. Contrary to existing work, the proposed threat detection system is lightweight enough to be run on resource-constrained sensor nodes. It greatly improves the resilience of the industrial WSN, without bringing in excessive energy con-sumption. Our solution is based on a partly centralized architecture and specifies new roles for WSN entities, in accordance with their status and capabilities.

The remainder of this paper is organized as follows. Section 2 quickly describes challenges related to the development of a threat detection framework for WSNs and highlights the limitations of existing approaches. Section 3 describes the proposed threat detection solution and its implementation on physical sensor devices. We assess its performance when compared to main state of the art solutions in Section 4. Section 5 concludes the paper.

Problem Statement and Related Work

The need for a threat detection system arises from the need to protect WSNs against attacks. Obviously, we assume that WSNs feature dedicated authentication and access control routines, key establishment functions. The security here is provided by the use of cryptographic primitives, which protects the WSN from external attacks. However, these security functions fail to protect the network against an authenticated attacker from inside the network: a sensor node owning legitimate cryptographic keys can easily launch an internal attack inside the group. It may be reluctant to make its resources available to other nodes as part of a cooperative act. Therefore, such a node may prefer to behave at times in a selfish manner and refuse to cooperate in packet forwarding in order to maximize its energy savings. On the other hand, a legitimate sensor node can also act maliciously by dropping packets, delaying the transmission or sending packets through a different route than planned. Such an internal attacker can only be detected through behavioral analysis mechanisms. The latter track the system behavior and interactions between nodes to detect threat attempts and/or occurrences. Once a security anomaly is detected, a reaction mechanism is launched to take security and service repair measures. While highly relevant to the autonomous topologies that wireless sensor networking involves, threat detection raises issues with respect to its adaptability to these domains. It challenges indeed the constrained nodes' limited energy resources by involving complex processing and high resource consumption.

Most of existing researches on threat detection systems propose to mimic models proposed for Mobile Ad hoc NETworks (MANETs) in order to counter internal attacks inside WSNs [START_REF] Da Silva | Decentralized intrusion detection in wireless sensor networks[END_REF][START_REF] Krontiris | Cooperative intrusion detection in wireless sensor networks[END_REF][START_REF] Krontiris | Towards intrusion detection in wireless sensor networks[END_REF]. These threat detection models rely on monitoring agents deployed on each sensor node. The goal of these local agents is to track the traffic within their radio range and detect misbehaving nodes causing routing disruptions. Authors in [START_REF] Roman | Applying intrusion detection systems to wireless sensor networks[END_REF] define a threat detection system based on monitoring agents. They distinguish between local agents (able to process only the packet they actively forward) and global agents. Global agents act as watchdogs by monitoring nearby traffic, especially when they determine how a forwarding node has behaved (e.g. difference, delay and loss between an incoming packet at a neighbor). In addition to regularly monitoring global agents, the authors introduce the concept of spontaneous watchdog behavior: a node determines from a sensed packet that it can be in position to act as a watchdog for this packet. This node will first identify the number n of its neighbors that could also act as watchdogs for the same packet, and turns on its watchdog behavior for that packet with a one-in-n probability. However, local passive monitoring has been shown to draw as much power from the sensor node transceiver as data reception [START_REF] Wander | Energy analysis of public key cryptography for wireless sensor networks[END_REF]. Ioannis et al. in [START_REF] Ioannis | Towards intrusion detection in wireless sensor networks[END_REF] also apply the watchdog behavior, with a higher emphasis on countering malicious monitoring nodes. To that aim, cooperation through majority is achieved, relying on an encrypted flag. That means that a node will be classified as malicious only if a majority of its neighbors flag it as such. The problem with the use of watchdog behavior is that watchdog nodes will drain their resources quickly.

Butun et al. in [START_REF] Butun | A survey of intrusion detection systems in wireless sensor networks[END_REF] review existing intrusion detection systems in WSNs and highlight the fact that IDSs designed for MANETs cannot be applied to WSNs directly. Authors consider that resource constraints characterizing the sensor nodes should be taken into account for the design of an IDS adapted to WSNs. Also, the fact that most of sensor nodes are most of the time in sleeping mode makes the operation of this security system more complex, because synchronous node actions might prove difficult, if not impossible, to achieve.

Other researches recognize that designed IDSs need to spend the least amount of energy as possible to spare enough energy for the crucial operations of the WSN. For this reason, Huang and Lee propose in [START_REF] Huang | A cooperative intrusion detection system for ad hoc networks[END_REF] an energyefficient monitoring system that fit the energy constraints of sensor nodes. They select a single node to perform monitoring at a given time within a given cluster. This node is designated through a fair election process. Authors claim that they obtain much higher energy efficiency at an equivalent security level compared with systems where all nodes perform monitoring at the same time. They also provide a detailed set of rules that allow detecting classical attacks expected to occur within a wireless sensor network. Authors in [START_REF] Bao | Hierarchical trust management for wireless sensor networks and its applications to trust based routing and intrusion detection[END_REF][START_REF] Shaikh | Group-based trust management scheme for clustered wireless sen-sor networks[END_REF][START_REF] Zhang | A trust management architecture for hierarchical wireless sensor networks[END_REF] propose hierarchical trust management systems for WSNs to detect selfish and malicious nodes. In these approaches, regular sensor nodes do not participate in the global decision making process. Only Cluster Heads (CHs) are responsible for the global decision making process and the response. The main reason for this is to reduce the energy consumption. They wanted to conserve the energy of the majority of sensor nodes, by simply assigning them as subordinates under CHs.

These cluster based detection schemes consider that a node is periodically elected to be the monitoring agent within each cluster. However, the election process could be heavy for constrained nodes -as demonstrated in what follows. In the same context, Butun et al. in [START_REF] Butun | A survey of intrusion detection systems in wireless sensor networks[END_REF] underline that clustering based IDSs may consume considerable amount of the WSN's energy through exchanged messages between nodes to periodically elect new monitoring agents.

3 Solution Description

Overview

Considering the inadequacies of existing approaches described in the previous section, we propose an efficient threat detection system for WSNs that aims at meeting two key requirements. Considering the energy constraints of sensor nodes will be the first key requirement for the design of our system. Processing within sensor nodes will be minimized in order to increase their lifetime. So that, most of required operations are delegated to the server side and the charge on the sensor node will be kept light, which leads to extend the network lifetime. Offloading the charge from the constrained nodes to the server is not the only advantage of this approach. Having to send its observations to the central server to take decisions, a malicious node would not be in position to propagate false assessments about specific victims. With a central entity responsible for decision making, it becomes a common profit for all nodes to provide reliable evidences since false ones can globally affect decision making at the central entity, and could eventually be detrimental to the attacker itself.

The second important key requirement is that a sensor node is often in a sleeping mode, so that considering that it is able to perform synchronous actions with other nodes for threat detection as proposed in previous works would be difficult to achieve. This problem will be carefully taken into account for the design of our solution. Sensor nodes will wake up only to collect data and perform autonomous security-related actions, then revert back to sleep mode.

Network Model

The wireless sensor network is supposed to be divided into zones. Each zone contains one or more clusters and each cluster contains one or more sensor nodes. Zones and clusters have different criticalities, different security levels and different security policies. It is also assumed that sensor nodes within the same cluster can communicate with each other. In addition, we assume that the awake time is negligible compared to the sleep time. Finally, we exclude any form of synchronization between nodes.

In a WSN, there are two types of nodes: sensing devices and gateways (or servers). The sensing devices are simple nodes equipped with radio interfaces only. The gateways are equipped with an Ethernet port to communicate with remote application servers and deliver collected data. The Gateway has more computing capacities and unlimited energy resources, as compared with a constrained sensor node. Since these nodes are Ethernet-connected, assuming that they are also connected to a power supply seems a reasonable hypothesis.

Components and Roles

The security system we present in this paper is made up of the following elements:

• Threat Detection Client: the TD client is in charge of identifying threats and sending reports to the TD server.

• Threat Detection Server : the TD server receives the registration requests from sensor nodes. It chooses which sensor(s) will be in monitoring mode for each cluster by taking into account status parameters such as batteries level and available resources. The TD server updates the global network database. It receives the alarms from TD clients and transfers them to the SA Server.

• Security Adaptation Server : the SA Server receives the threats from the TD server. It decides then which is the best policy to apply accordingly and stores the new policy for sensor nodes in the security policy mailbox.

• Security Adaptation Client: After each boot, the SA Client sends a message to the SA Server to check if there are any policies to be delivered and applies those it receives in return.

• Security Service: Various security services are supposed to run on the considered sensor nodes, but two of them are especially relevant. A key management service is assumed to be able to establish security contexts between two nodes. The security level offered by this service translates into various parameters such as algorithms, key lengths or security association life-times. Another security service is the network access control enforcement, which maintains a secure connection between a sensor node and its parent node in the WSN topology. This secure connection is established upon node's entry into the WSN, and remains active as long as the node is part of the WSN. Both of these security services can raise alarms to the local threat detection client even though the node is not actively monitoring its neighborhood, for example if a peer node cannot be authenticated (by other sensor in case of key management, by parent node in case of network access control enforcement).

• Security Policy mailbox : a module that stores the generated policies. Then, it delivers them at nodes request. The use of this module is required since the nodes are not synchronized. It also reduces the overall bandwidth consumption.

• Global Network Database: contains a global view of the network and the threats detected in the past.

It is populated every time a sensor registers to the TD Server and is updated every time a new policy is applied by the SA Server.

Operation

This section describes the sequences of actions performed by the client-side and server-side entities. These sequences of actions are categorized according to the sensor mode. This mode can either be bootstrapping, normal mode, or monitoring mode.

Bootstrapping is the mode of a node when it joins the network. The newly joining node is first to send a registration request to the TD server, informing this latter about its potential monitoring abilities. Upon reception of this registration request, the TD server registers the node and sends to the node an acknowledgement (ACK), as well as a configuration message specifying whether the node should remain in normal mode or switch to monitoring mode for a specific amount of time. The decision by the TD server is based on its knowledge of current and, in some cases, foreseen contexts of the candidate monitoring nodes. This contextual information includes data relevant to the nodes resources (e.g. battery level), location, and capabilities (e.g. number of observable neighbors). With this information, the TD server is able to identify the best node in the cluster for acting as a monitoring entity, and to configure it with this role for a certain period of time. Once the monitoring delay expires, the TD server proceeds again to the identification and designation of cluster monitoring node(s).

A node switches to Monitoring Mode when ordered to do so by the TD server, either immediately after its bootstrapping, or upon reception of a configuration order, after waking up. The sequence of actions performed when in monitoring node is: 1) Once the TD Client detects a threat, it sends an alarm to the TD server that includes information about the threat. This information contains at least the IP address of the attacker, the IP address(es) of the target(s) and the type of attack;

2) Upon receiving the alarm, the TD Server reports it to the SA server, optionally after having aggregated multiple alarm messages and/or having assessed the quality of the evaluator. Next, it stores the new policies in the security policy mailbox, in order to have them be delivered to the respective SA Clients that will have to enforce them;

3) In monitoring mode, the TD client on the sensor regularly polls the security policy mailbox by sending a dedicated inquiry message to the SA server;

4) The SA server sends the requested policy if it exists. Otherwise, it replies with a message telling the node that it is not to enforce a new policy. Along with security policies, the monitoring node might be instructed to revert back to normal mode if another monitoring node is being designated within the cluster;

5) If a new policy is received, the SA Client enforces it by configuring the security services in accordance with the received policy and acknowledges it through an ACK message;

6) The SA Server receives the ACK and updates the global network database accordingly.

The monitoring mode sequence of actions is represented in Figure 1.

If a malicious node refuses to be placed as a monitoring node, the server could detect that no reports are received from its side and then it could be penalized or excluded from the network. If this node accepts the monitoring process during a period of time but lets other colluding nodes execute attacks while sending positive observations on their behalf. The server can compare its false reports with feedbacks received from other monitoring nodes selected for the same cluster in different time intervals. Hence, it will conclude its malicious behavior as a monitoring node and take the appropriate punishment decision against it.

The Normal Mode is the default mode for a bootstrapped sensor node that has not been designated as a monitoring node. In this mode, sensor nodes alternate between active and sleeping states. Upon leaving sleeping state, the node interrogates the SA server about an eventual new policy to enforce. Meanwhile, the node in normal mode may also be instructed to switch into monitoring mode, if required from the evaluation of different nodes energy levels within the considered cluster. The node then performs the task(s) for which it has left the sleeping state. An alarm may be raised by a node in Normal Mode only if one of the run tasks detects a threat. This task would then notify the TD Client through an API call.

The sequence of actions corresponding to normal mode is represented in Figure 2.

Implementation Environment

The hybrid threat detection solution presented in this paper was implemented using embedded C for sensor nodes, with Atmel studio 6 for AVR/ARM and C++ for the server side. We used Dresden Elektronik sensor devices [START_REF] Elektronik | [END_REF] that have the following specification:

• The gateways are equipped with ARM processor, Ethernet port, power over Ethernet, USB/Serial and IEEE 802.15.4 2.4GHz.

• The sensor devices are equipped with AVR Processor, battery, USB/Serial interface and IEEE 802.15.4 2.4GHz.

Both the gateway and the sensor node use 6LOwPAN protocols. The sensor node turns to the Monitor Mode 4 Performance Evaluation

Performance Evaluation Criteria

In order to assess the performance of the proposed hybrid threat detection system, it is worth detailing first all criteria that should be taken into account for the evaluation of a threat detection system. Indeed, a threat detection system can be characterized with:

1) The proportion of false positives (benevolent behaviors wrongly identified as attacks) and that of false negatives (malicious behaviors not identified as such), as compared with the successful attacks identifications. These ratios and the corresponding performance metrics assume the existence of a list of detectable attacks {A 1 , ..., A n }. Accordingly, an attack A ′ not detectable -and thus not detected -by the system would not be qualified as a false negative.

2) The exhaustivity of the detectable attacks list constitutes a second evaluation metric. Ideally, this list must be built in accordance with a risk analysis carried out on the monitored system, so that the most frequent and/or critical attacks be detectable by the threat detection system.

3) The reactivity of the threat detection system that is, the time needed to identify an attacker as malicious.

Finally, constrained environments characterized with nodes of low capabilities in terms of computing power and energy capacity exhibit additional requirements, to which a fourth evaluation criterion corresponds.

4) The energy needed to locally (at constrained nodes side) operate the threat detection system. This energy cost must be limited, so that the constant use of this system will not consume a high amount of nodes energy re-sources prohibitively fast.

The proposed solution for threat detection does not focus on mitigating false positives and false negatives. From this viewpoint, the local detection algorithm is not improved as compared to the state of the art solutions. Likewise, the weighting and aggregation of alerts at the server side is recommended, but a novel method for doing so is not part of the current proposal. Hence, our proposed system will not be evaluated according to the first evaluation criteria. Our proposed system does also not extend the list of detectable threats, as defined by the second evaluation criteria. The completeness of this list is seen as orthogonal to the current solution. We acknowledge, however, the high relevance of this problem, especially when put aside the (often neglected) nodes constraints in terms of memory capacity. Indeed, an exhaustive list of detectable attacks would require a corresponding memory space for storing attacks signatures, which would not be affordable to highly-constrained nodes.

While criteria (a) and (b) are thus not considered for assessing the performance of our proposed threat detection system, the importance we put in the reduction of energy consumption leads us to favoring criteria (d) as the most important one for doing so. Meanwhile, we will use the criteria (c) both to guarantee a proper behavior of the proposed threat detection system even though it is instantiated on highly-constrained nodes and to provide a second metric to compare our solution to those proposed in the literature.

Threat Detection Simulator

We evaluated the performance of our proposed threat detection system by means of a purposely designed discrete time simulator. This section first reviews the generic design decisions that were applied for conceiving this tool. It then considers the models that were used to implement on the simulator both the proposed system and the solution proposed in [START_REF] Huang | A cooperative intrusion detection system for ad hoc networks[END_REF], which is one of the most studied threat detection system for WSNs, to which we intend to compare our work.

The proposed simulation program Csimu (a lightweight discrete time C-based network simulator) starts with initializing a WSN-like topology made of clusters and cluster heads. It also initializes the threat detection system, in particular by assigning per-cluster initial monitoring nodes. We conceived our simulator as a program built around a main loop, with each loop execution corresponding to an elapsed (configurable) time step. A discrete time approach is therefore adopted in this tool. The main loop is made of the following operations:

• Storage of logged values including: overall energy spent, overall energy spent for all threat detection related operations, nodes compromise status.

• Per-node individual processing loop. This second loop is the most important operation. It goes through each node part of the generated topology and performs relevant nodes actions. These consist in:

• Storing per-node values relative to energy spent and energy spent for threat detection operations.

• Updating node energy level in accordance with node status during the last time step. For example, a node that spent the last time step in 'Listen' mode will have its battery level decreased with an energy E = ∆t T imestep × P Listen , with ∆t T imestep being a configurable parameter and P Listen the power for the considered node type in listening mode.

• Carrying out node tasks, e.g. uploading a measured value to a remote sink node and updating node energy level accordingly.

• Carrying out attacks, in case the node is compromised. Two types of attacks are represented: A DoS is an example of Denial of Service attack and A Comp represents a compromise attack. Each attack involves the following steps: determination of whether an attack is to occur, choice of the to-be-attacked victim, energy consequences on the attacking node (corresponding to the number of messages to send/receive and the listening durations involved by the attack), energy consequences on the victim, other outcomes on the victim (the A DoS attack has a chance to make the victim compromised), attack detection. The attack detection involves all nodes that are in position to detect the launched attacks, taking in considerations their respective locations (only neighbor nodes are susceptible to discover the attack) and modes (the probability of detecting the attack is much higher for the nodes that are in monitoring mode). In turn, when a node detects an attack, it may either be able to pinpoint the attacker, or may just notice that an attacker is likely to be present in the cluster. In accordance with the detection results and the implemented threat detection system, various actions are to be taken that are detailed in the next section.

• Threat detection system update (detailed in the next subsection) that includes the exclusion of the identified compromised nodes and, in random (nonscripted) mode, a slight probability for each node to become compromised.

Evaluated Threat Detection Systems

"All Monitor" and "Cluster Monitoring Node Election" Systems

As explained in the related work section, the "all monitor" approach considers that each sensor device acts as an independent monitoring node. Continuing in listening mode, each node tracks its neighbor's communication in order to identify possible attacks. The "Cluster monitoring node election" approach as proposed in [START_REF] Huang | A cooperative intrusion detection system for ad hoc networks[END_REF] relies on monitoring by an elected node, with election process occurring periodically. In accordance with the protocol description in [START_REF] Huang | A cooperative intrusion detection system for ad hoc networks[END_REF], we simulated the election process as a recurring procedure where:

• Each node has its energy level decreased by a decrement corresponding to: (1) a random number generation -we assumed a hash-based random number generation; (2) the computation of a hash on the computed random number and the node identifier; (3) the sending of a message containing said identifier and hash; (4) the waiting for receiving all messages from neighbor (same-cluster) nodes; (5) the actual reception of all of these (n-1) messages (n being the number of nodes within one cluster); ( 6) the verification of all received messages through the computation of (n-1) hashes; [START_REF] Ioannis | Towards intrusion detection in wireless sensor networks[END_REF] the operation of the selection function; [START_REF] Krontiris | Cooperative intrusion detection in wireless sensor networks[END_REF] the sending of a result message.

• As a result of the election, a new node becomes responsible of cluster monitoring. That is, the former monitoring node reverts back to the 'normal' behavior (sleeping mode, with periodic wake up) and the newly elected monitoring node switches to the 'monitoring' behavior (listen mode) for the next period.

Our Hybrid Solution

We implemented our proposed approach on the Csimu simulator by adding within the server processing operation a phase where the threat detection server periodically designates a new monitoring node within all clusters. To that aim, the server maintains a view of the network nodes statuses, which makes it able to designate for each cluster the node with the highest energy level as the new monitoring node. The process happens thus as follows:

• Upon waking up, each node sends to the threat detection server a message prompting for orders and informing the server about its current battery level.

The server answers accordingly. The energy costs corresponding to message sending and response waiting, receiving and processing are taken into account and decremented from the nodes energy levels.

• Periodically, the server triggers a new designation of a cluster monitoring responsible node. This does not incur any immediate action. However, the subsequent request for orders of the former monitoring node and the newly designated monitoring node will make these learn, respectively, their orders to switch into 'normal' and 'monitoring' modes.

• The protocol that was implemented in the simulator introduced a slight variation as compared to the approach described above: each cluster is initialized with no active monitoring node. It is assumed that nodes in normal mode can pinpoint -even with much higher probabilities of false positives and false negatives -malicious behaviors and consequently warns the threat detection server, which will in turn mark the cluster as suspicious and assign a monitoring node within its population. The metrics for identifying clusters (proposed approach) / nodes (both approaches) as suspicious / compromised are explained in next section.

Simulation Parameters: Assumptions and Devices Characteristics

This section reviews the simulation parameters that are used in the simulation environment to rate nodes behaviors, on one hand, and to assess energy consumption, on the other hand.

Assessment of Attack Probability

Without implementing any actual attack detection scheme, we based our simulated threat detection function on a probabilistic model of attack detection, where each behavior (A 1 , . . . , A n attacks, as well as the 'A 0 ' benevolence) can lead to a detection of any behavior within the (A 0 , A 1 , . . . , A n ) population, with different probabilities. This led us to define an attack recognition matrix M AR as:

M AR =      p(D A0 |A 0 ) p(D A1 |A 0 ) • • • p(D An |A 0 ) p(D A0 |A 1 ) p(D A1 |A 1 ) • • • p(D An |A 1 ) . . . . . . . . . . . . p(D A0 |A n ) p(D A1 |A n ) • • • p(D An |A n )     
where p(D X |Y) corresponds to the probability to detect the event X knowing that event Y occurred.

In fact, two such matrices were defined corresponding to the two modes in which a node can be, respectively normal mode (non-monitoring node involved in usual node operations only) and monitoring node. The probabilities were tuned accordingly, to reflect the fact that a node in monitoring mode is much more accurate in its identifications, whereas the probabilities of false positives (p(D Ai | A 0 ) > 0 for some i >0) and false negatives (p(D A0 | A i ) > 0 for some i >0)are much higher for a node in normal mode.

Based on this matrix, we were able to compute the probability of actual attack (whatever the A j , j∈ [1; n]) occurrence upon the detection D Ai of a given event A i .

p(Attack|D

Ai ) = p(D Ai ∩ Attack) p(D Ai ) = p (D Ai ∩ (A 1 ∪ • • • ∪ A n )) p(D Ai ) = n k=1 (p (D Ai A k ) .p(A k )) n k=0 (p (D Ai A k ) .p(A k )) = n k=1 (M AR [k, i] .p(A k )) n k=0 (M AR [k, i] .p(A k ))
In this computation, the elements M AR [k, i] are the coefficients of the MR matrix defined above. The p(A 0 ), p(A 1 ), . . . , p(A n ) probabilities of occurrence of events A 0 , A 1 , . . . , A n are however other parameters that are required, and for which only approximate values can be used. In the simulator we used values derived from the probability of node compromise, and, for a compromised node, the probability of launching a typei attack A i . In actual conditions, where these probabilities of attacks cannot be precisely known, the system should use cognitive behavior where the optimal values of p(A k ) would be statistically approached from the history of past detected events.

With an attack detection being mapped to the probability of actual attack occurrence, we were able to increase node compromise scores by a proportional coefficient. Eventually, we defined four thresholds t 1 , t 2 , t 3 , and Figure 4: Energy consumption for threat detection operation in all three evaluated systems t 4 respectively for marking a cluster as 'to be monitored' (when the cluster score exceeds t 1 ),marking a cluster as 'no longer to be monitored' (when the cluster score falls below t 2 ), marking a node as 'likely compromised' (when the node score exceeds t 3 ) and marking a node as 'likely benevolent' (when the node score falls below t 41 ). We made the cluster and node scores decremented at each execution of the main loop, to mitigate false positives and let the system revert back to normal state after a compromised node has been excluded.

From Node Characteristics to Computation Parameters

We used the node characteristics listed in [START_REF] Saied | D-hip: A distributed key exchange scheme for hip-based internet of things[END_REF] corresponding to the following parameters: power in transmit mode (54 mW), power in receive mode (61 mW), power in listen mode (60 mW) and energy consumption for one processor cycle (8.64 nJ / cycle). Likewise, sensor nodes throughput was supposed to be equal to 75 kbps, as determined in [START_REF] Saied | D-hip: A distributed key exchange scheme for hip-based internet of things[END_REF].

Random number generation and hash function required in [START_REF] Huang | A cooperative intrusion detection system for ad hoc networks[END_REF] were assumed to be based on the SHA-1 hash function, for which the required number of cycles per byte was computed, based on eBASH [START_REF] Bernstein | eBASH: ECRYPT Benchmarking of All Submitted Hashes[END_REF], to be 53.24 on the considered processor for a less than 64-byte long message (the shorter the message, the higher the cost for hashing one byte, due to the SHA-1 fixed operations that do not depend on message length).

Simulation Results

This section reviews the obtained results in terms of energy cost and efficiency, which prove that our proposed system, though slightly slower to react to malicious operations, is more energy-efficient than the other evaluated solutions.

Energy Consumption

Figure 4 presents the overall energy consumption during all processes related to threat detection for a 5-cluster, 50node topology. This energy consumption is expressed in Joules per 10 -1 s for the three evaluated systems, namely 'All monitor' (all cluster nodes stay in monitoring mode), 'Elect' (there is only one monitoring node per cluster, which is regularly refreshed through an election carried out between the cluster members [START_REF] Huang | A cooperative intrusion detection system for ad hoc networks[END_REF]) and our proposed 'Hybrid' scheme.

As was expected, the 'All monitor' approach is the most requiring in terms of energy with the 50 nodes (actually, 45 sensor nodes, since energy consumption is not measured on gateways by the simulation environment) consuming around 135 J in 50 s (a result that corresponds to the 60 mW that each node consumes in listening).

Figure 4 also highlights the high energy cost of [START_REF] Huang | A cooperative intrusion detection system for ad hoc networks[END_REF], due to its heavy election processes. Even though only one node monitors a cluster at a time, the cost of election phases (the vertical transitions between plateaus in the curve) does not allow an overall energy gain higher than 33%, as compared with the 'All monitor' approach.

In order to highlight the difference between the local election phase of [START_REF] Huang | A cooperative intrusion detection system for ad hoc networks[END_REF] and the server-assisted monitoring node designation we proposed, we simulated in Figure 4 a synchronous approach of our hybrid approach, where all cluster nodes in all clusters perform mailbox checking (and, if required, mode change) at the same moment. Nevertheless, the vertical transitions are much smaller than those of [START_REF] Huang | A cooperative intrusion detection system for ad hoc networks[END_REF], as could be expected from the lighter requirements exhibited by the mailbox checking, as compared to local election processing. Eventually, our solution performs approximately twice better than [START_REF] Huang | A cooperative intrusion detection system for ad hoc networks[END_REF] and three times better than the 'All monitor' approach, with a final threat detection energy cost of only 44 J for 45 sensing nodes in 50 seconds (around 20W for each node). It has to be noted that the exponential-like allure of the 'Hybrid' system curve corresponds to the tuned behavior of the system, wherein clusters are initially not monitored (and consume almost no energy for threat detection) and start to get monitored only when compromised nodes reveal their malicious activities and are suspected by nodes. This leads to an acceleration of the energy consumption, which quickly reaches a cruising speed, where all clusters are monitored.

Threat Detection Reactivity

The behavior of all three simulated solutions with respect to compromised nodes identification is detailed below in Figure 5, Figure 6 and Figure 7.

Figure 5, Figure 6 and Figure 7 represent the compromise/react graph for each of the three simulated models. The dark, diamond-marked curves represent the number of unidentified compromised nodes within the system. The light, square-marked peaks represent detections of compromised nodes, and correspond therefore to a decrease of the other curve. Note that all three threat detection systems were tested in the same conditions, with compromising of the same nodes being triggered at the same time (scripted, non-random, compromising mode).

From these results, it appears that the 'All Monitor' approach performs the best with respect to speed in threat identification. Cluster election and the proposed hybrid approach, on the other hand, show almost similar results requiring more regularly a few seconds to detect and react to the presence of an attacking node.

Other Benefits of the Proposed Hybrid Approach

A synchronous embodiment of the proposed hybrid approach was depicted in Figure 4, in order to reflect the lower cost induced by the threat detection management phase (mailbox checking and required actions enforcement) as compared with that of the election approach. Beside this explanatory purpose, there is no need however to mandate synchronicity in our proposed solution: the sensor nodes may very well wake up at different times; the operation of the proposed approach would still work exactly the same. We simulated this asynchronous approach by explicitly requiring each sensor node to wake up at different times. We obtained equivalent performance with respect to reactivity to attacks; on the other hand, the resulting threat detection energy consumption curve presented, as expected, a smoother aspect (Figure 8). Obviously, the support by our proposed approach of such asynchronicity (hence, of multiple sensor nodes that are not in phase regarding their sleep/wake up periods) could not be provided by the election approach, where all nodes within a common cluster must be involved in a synchronous way in the monitoring node election process. This support of sleeping node represents an important advantage of our proposed approach.

Finally, another benefit of our hybrid approach lies in the fact that interactions with a server entity can be exploited to further improve the threat detection system, by making it dynamically able to react to new threats. Indeed, the orders from the threat detection server to its clients could be extended to carry new attack signatures, thereby making the system both more evolutive and more in line with a given threat model. Simple extensions to the threat detection system, in line with the adaptive security approach, could lead to dynamically choosing which node(s) in a given cluster receive which attack signature(s).

Conclusion

This paper presents a threat detection system for industrial wireless sensor networks that could be qualified as hybrid in that it involves a semi-centralized process. The switch from normal threat detection mode to monitoring mode is triggered by the threat detection server, which bases on regular reports from nodes and updates its decision accordingly. We showed that this system, as reactive as the most studied state of the art solutions with respect to identifying threats, performs better from the viewpoint of energy consumption, saving around 50% of energy. The hybrid approach we propose, characterized by centralized management and local instantiation of threat detection, is also more flexible in terms of extension possibilities, opening up interesting development axes for the future of autonomous security in WSNs.

Figure 1 :Figure 2 :

 12 Figure 1: Threat detection and security adaptation in monitoring mode

Figure 3 :

 3 Figure 3: Overall logical architecture and state machines

Figure 5 :Figure 6 :Figure 7 :

 567 Figure 5: All monitor

We did not use that fourth threshold, though in our simulation: instead, we considered that a compromised node was to be removed from the system; and could optionally be re-introduced after a while -but considered as an entirely new node.

Acknowledgments

This work was financially supported by the EC under grant agreement FP7-ICT-258280 TWISNet project.


Soraya Rouifed

M Cottet

M De Battista

Y-F Le Lay

Florence Piola

P Rateau

A Rivière-Honegger

Landscape perceptions and social representations of Fallopia

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Mamadou Toungara

Gregory Chagnon

Christian Geindreau
email: christian.geindreau@hmg.inpg.fr

NUMERICAL ANALYSIS OF THE WALL STRESS IN ABDOMINAL AORTIC ANEURYSM: INFLUENCE OF THE MATERIAL MODEL NEAR-INCOMPRESSIBILITY

Keywords: Abdominal aortic aneurysm, anisotropic models, numerical analysis, near-incompressibility

Recently, hyperelastic mechanical models were proposed to well capture the aneurismal arterial wall anisotropic and nonlinear features experimentally observed. These models were formulated assuming the material incompressibility. However in numerical analysis, a nearly incompressible approach, i.e., a mixed formulation pressure-displacement,i s usually adopted to perform finite element stress analysis of abdominal aortic aneurysm (AAA). Therefore, volume variations of the material are controlled through the volumetric energy which depends on the initial bulk modulus κ. In this paper, an analytical analysis of the influence of κ on the mechanical response of two invariant-based anisotropic models is first performed in the case of an equibiaxial tensile test. This analysis shows that for the strongly nonlinear anisotropic model, even in a restricted range of deformations, large values of κ are necessary to ensure the incompressibility condition, in order to estimate the wall stress with a reasonable precision. Finite element simulations on idealized AAA geometries are then performed. Results from these simulations show that the maximum stress in the AAA wall is underestimated in previous works, committed errors vary from 26% to 58% depending on the geometrical model complexity. In addition to affect the magnitude of the maximum stress in the aneurysm, we found that too small value of κ may also affect the location of this stress.

Introduction

In the last two decades, numerous numerical studies [START_REF] Humphrey | Intracranial and abdominal aortic aneurysms: Similarities, differences, and need for a new class of computational models[END_REF] have been performed in order to predict the wall stress distribution in abdominal aortic aneurysms (AAA), i.e., to predict the rupture of AAA. Accurate and reliable numerical stress analysis requires to take into account the geometrical features (tortuosity, asymmetry...) of AAA, the wall thickness, the presence of the intraluminal thrombus, the heterogeneity of the wall, and the use of appropriate constitutive laws for the AAA material. The present work mainly concerns this last point. Like many soft tissues, the arterial wall is assumed to be an incompressible material namely in experimental studies. Under physiological conditions, this assumption could be justified. [START_REF] Carew | Compressibility of the arterial wall[END_REF] The majority of published numerical studies in AAAs 1 have been performed using incompressible hyperelastic isotropic models fitted thanks to classical uniaxial tensile tests. [START_REF] He | The composition and mechanical properties of abdominal aortic aneurysms[END_REF][START_REF] Raghavan | Toward a biomechanical tool to evaluate rupture potential of abdominal aortic aneurysm: Identification of a finite strain constitutive model and evaluation of its applicability[END_REF][START_REF] Thubrikar | Mechanical properties of abdominal aortic aneurysm wall[END_REF][START_REF] Xiong | Measurement and analysis of ultimate mechanical properties, stress-strain curve fit, and elastic modulus formula of human abdominal aortic aneurysm and nonaneurysmal abdominal aorta[END_REF] However, biaxial tensile tests are more appropriate to characterize the anisotropic behavior of aneurismal arterial wall and for the development of suitable material models. In this way, planar biaxial tensile tests on normal and pathologic human abdominal aortic tissues were carried out in Ref. 7. The obtained results clearly show that the degeneration of the aorta leads to an increase in mechanical anisotropy and also an increase in the circumferential stiffness for AAA tissue as compared to normal aortic tissue. In order to take into account such features, hyperelastic anisotropic constitutive models were recently proposed using a Green-strain approach [START_REF] Geest | The effects of aneurysm on the biaxial mechanical behavior of human abdominal aorta[END_REF] or an invariant approach. [START_REF] Rodriguez | Mechanical stresses in abdominal aortic aneurysms: Influence of diameter, asymmetry, and material anisotropy[END_REF][START_REF] Basciano | Invariant-based anisotropic constitutive models of the healthy and aneurysmal abdominal aortic wall[END_REF] In both cases, the strain energy density W is formulated assuming that the arterial material is incompressible. These models were then implemented in a finite element code in order to compute both the distribution and the magnitude of stresses in patient-specific geometries using shell elements [START_REF] Geest | The effects of aneurysm on the biaxial mechanical behavior of human abdominal aorta[END_REF] or in idealized AAA geometries 8 using 3D elements. To perform such simulations, a nearly-incompressible approach, i.e., a mixed formulation pressuredisplacement implemented in many finite element softwares, is usually adopted. For that purpose, the strain energy density W is split into an isochoric part W and a volumetric part U : W = W + U . The volumetric function U which depends on the jacobian J and on the initial bulk modulus κ of the material is used to control the volume variation. [START_REF] Doll | On the development of volumetric strain energy functions[END_REF] In practice, the incompressibility of the material is nearly reached for large values of κ which is used as a user-specified penalty parameter.

The main objective of the present work is to emphasize the influence of the initial bulk modulus κ values on both the magnitude and the distribution of the stress in an anisotropic AAA wall. To our knowledge, these values of κ as well as the form of the volumetric function U (J) are rarely specified or discussed in previous numerical works. For that purpose, the paper is divided into two sections. In Sec. 2, an analytical analysis of compressibility effects on equibiaxial tensile tests is carried out. Two strain energies based on polynomial or exponential functions [START_REF] Rodriguez | Mechanical stresses in abdominal aortic aneurysms: Influence of diameter, asymmetry, and material anisotropy[END_REF][START_REF] Basciano | Invariant-based anisotropic constitutive models of the healthy and aneurysmal abdominal aortic wall[END_REF] are considered and a classical quadratic volumetric function is used. Through a parametric study, the influence of the initial bulk modulus on the near-incompressibility of the mechanical response of the considered constitutive models is highlighted. Then, this influence is emphasized on idealized AAA geometrical models, in the last section. For that purpose the model proposed in Ref. 8 was implemented in a finite element code and numerical simulations on the AAA models were performed. The obtained results are finally compared to those in the literature.

Analytical Analysis of Compressibility Effects on Equibiaxial

Te n s i l e Te s ts 2.1. Material models

2.1.1. General framework
Hyperelastic models are often expressed thanks to strain elongations or strain invariants. The isotropic strain invariants I 1 , I 2 and I 3 (or J) are defined by:

I 1 =tr(C),I 2 = 1 2 tr(C) 2 -tr(C 2 ) ,I 3 = J 2 =det(C), (1)
where C = F T F is the right Cauchy-Green strain tensor and F is the deformation gradient whose principal components are the elongations λ i . For fiber reinforced materials, an anisotropic approach must be used and a fiber direction vector M is defined for each fiber orientation. This leads to the introduction of two new invariants, I 4 and I 5 which are defined as:

I 4 = MCM,I 5 = MC 2 M. (2)
Other invariants can be introduced if there are many reinforced directions. Many hyperelastic materials, as arterial tissues, are often considered as incompressible. To overcome well-known numerical difficulties associated with the incompressibility condition, a nearly incompressible approach is usually adopted. This later approach leads to use the isovolume deformation gradient F = J -1/3 F.I t s principal components are thus defined by λi = J -1/3 λ i (i =1, 2, 3) and it permits to build incompressible invariants:

Ī1 = I -1/3 3 I 1 , Ī2 = I -2/3 3 I 2 , Ī4 = I -1/3 3 I 4 , Ī5 = I -2/3 3 I 5 . (3)
The strain energy density W is thus split into two parts, an isochoric one W and a volumetric one U :

W = W ( Ī1 , Ī2 , Ī4 , Ī5 )+U (J). (4)
The volume variation is then controlled by the choice of the volumetric function U (J) which depends on the initial bulk modulus κ of the material. [START_REF] Doll | On the development of volumetric strain energy functions[END_REF] The incompressibility of the material is nearly reached when the volumetric energy U is greater than the isochoric energy W , i.e., for large values of the bulk modulus κ which in this particular case serves as a user-specified penalty parameter. This last point is discussed in the following.

Strain energy densities for artery material

Within this framework, two forms of isochoric strain energy ( W ) using exponential 8 or polynomial 9 functions were proposed in order to describe the experimentally observed hyperelastic anisotropic behavior of arteries. [START_REF] Geest | The effects of aneurysm on the biaxial mechanical behavior of human abdominal aorta[END_REF] In these models, the material is considered as a matrix material reinforced with two fiber networks (Fig. 1). Thus, two vectors are introduced to describe the two distinguished fiber orientations: M =(cosθ sin θ 0) and M ′ =(cosθ ′ sin θ ′ 0), with θ = -θ ′ . As a consequence, the strain energy density (W ) may depend on the following anisotropic strain invariants defined for the two orientations, as:

Ī4 = I -1/3 3 MCM, Ī5 = I -2/3 3 MC 2 M, Ī6 = I -1/3 3 M ′ CM ′ , Ī7 = I -2/3 3 M ′ C 2 M ′ , Ī8 = I -1/3 3 (MM ′ )MCM ′ ,I 9 =( MM ′ ) 2 . (5)
In practice, Ī1 , Ī4 and Ī6 are usually sufficient to capture the typical features of arterial response. [START_REF] Rodriguez | Mechanical stresses in abdominal aortic aneurysms: Influence of diameter, asymmetry, and material anisotropy[END_REF][START_REF] Basciano | Invariant-based anisotropic constitutive models of the healthy and aneurysmal abdominal aortic wall[END_REF][START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF] In this study, models in consideration are:

• polynomial reinforcement constitutive equation (P6 model)

W = U (J)+c 1 ( Ī1 -3) + α( Ī1 -3) 2 + k 1 ( Ī4 -1) 6 + k 3 ( Ī6 -1) 6 , (6)
where c 1 , α, k 1 and k 3 are material parameters. Basciano and Kleinstreuer 9 recently proposed a similar form without the isotropic neo-hookean part. In this particular case, let us remark that there is no initial slope in the stress-strain curve (c 1 = 0), which can lead to numerical difficulties (nonconvergence).

• Rodriguez constitutive equation (R model) 8 W = U (J)+c 1 ( Ī1 -3) + k 1 2k 2 {exp[k 2 [(1 -ρ)( Ī1 -3) 2 + ρ( Ī4 -Ī0 4 ) 2 ]] -1} + k 3 2k 4 {exp[k 4 [(1 -ρ)( Ī1 -3) 2 + ρ( Ī6 -Ī0 6 ) 2 ]] -1}, (7)
where k 2 and k 4 are dimensionless material parameters, ρ is a measure of anisotropy. Ī0 4 and Ī0 6 correspond to the initial crimping of the fibers. In the first model, the anisotropic term is assumed to contribute when either Ī4 > 1o r Ī6 > 1, or both.

Similarly, the term ( Īi -Ī0 i )( i =4 , 6), does not contribute to the strain energy whether Ī4 ≤ Ī0 4 or Ī6 ≤ Ī0 6 . This last model is equivalent to the model proposed by Holzapfel et al. [START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF] for ρ =1and Ī0 i =1. A general form for the volumetric function (Eq. ( 8)) was proposed in Ref. 10.

U/κ =((α +1) -1 J (α+1) +(β -1) -1 J (1-β) )(α + β) -1 -(α +1) -1 (β -1) -1 , (8)
where α and β are dimensionless parameters that can be used to obtain the most common volumetric functions proposed in the literature, α>0a n dβ>1. In this work, the classical quadratic function [START_REF] Simo | A framework for finite strain elastoplasticity based on maximum plasticdissipation and the multiplicative decomposition: Part I. Continuum formulation[END_REF] usually implemented in finite element analysis codes is chosen for the two material models:

U (J)= κ 2 (J -1) 2 . ( 9
)
It was obtained with α =1andβ =0inEq. (8). This volumetric function violates the condition on β, but as the deformation range investigated in this study leads to very small volume changes (J ≈ 1), all the volumetric functions proposed in Ref. 10 coincide as we can observe in Fig. 2(a). Moreover, the evolution of the derivative of this ratio with the volume change is also illustrated in Fig. 2(b), since this derivative is involved in the stress calculation. Table 1 shows the values of α and β for the different volumetric functions in Fig. 2.

In order to control the material volume variation, one can choose an arbitrary higher initial bulk modulus. However, since this modulus must be considered as a material property, linking it to the material parameters seems to be coherent. For that purpose, we choose to come back to linear elasticity thanks to the wellknown relationship between the initial bulk modulus κ, the Young's modulus and the Poisson's ratio. In the case of orthotropic materials, like the arterial wall, the initial bulk modulus can be expressed in term of the Poisson's ratio ν ij , and Young's moduli E i , E j and E k as follow [START_REF] Itskov | Elastic constants and their admissible values for incompressible and slightly compressible anisotropic materials[END_REF] :

κ = E i E j E k 3[(1 -2ν ij )E i E k + E j E k -E i E j ] , (10)
where i = j = k (i, j =1 , 2, 3). For the sake of simplicity, the classical isotropic relationship often available in finite element codes is adopted in the following:

κ = E 3(1 -2ν) , (11)
Table 1. Parameters for the volumetric function (Eq. ( 8)). 10 where E and ν are the Young's modulus and the Poisson's ratio. Due to this simplification, by convention we choose the e 1 direction (Fig. 1) as reference for determining the Young's modulus in infinitesimal deformation of anisotropic materials under consideration. This choice has no consequence on the following analysis. The Poisson's ratio in (Eq. ( 11)) is a driving representative parameter of the compressibility (which has no physical relevance) and the strict incompressibility is reached for ν =0.5. After some algebraic manipulations, it comes that:

U U 1 U 6a U 6b U 6c U 6d U 6e U 6f α 1 2 0.
• P6 model E =6c 1 . (12)
• Rm o d e l

E = 2 1+2 ρk1 c1 sin 4 θ [3c 1 +4ρk 1 (cos 4 θ +sin 4 θ) -4ρk 1 cos 2 θ sin 2 θ]. (13)
The material parameters for each model are given in Table 2.T h em a t e r i a l parameters for the exponential model were identified in Ref. 8 on biaxial tensile tests data from Ref. 7. Using these parameters, a Young's modulus of 264 kPa was obtained with (Eq. ( 13)). Let us remark that the fibers orientation, θ,i nt h e s e models was considered as a phenomenological parameter, so its values have no physical relevance. [START_REF] Rodriguez | Mechanical stresses in abdominal aortic aneurysms: Influence of diameter, asymmetry, and material anisotropy[END_REF] For P6 model, to have the same Young's modulus as for R model, c 1 =4 4 .0 kPa was adopted. Thus, for a given value of the Poisson's ratio, we have the same initial bulk modulus for the two material models (Eq. ( 11)). The parametric study in the following is achieved with the initial bulk modulus, κ.

Equibiaxial loading

To study the different constitutive equations, a loading state close to that happens in arterial walls, i.e., equibiaxial loading, is chosen. The same elongation is imposed in the two plane directions (λ 1 = λ 2 ). The stress states are determined by means of the derivative of the strain energy density: W i = ∂W/∂ Īi and U J = ∂U/∂J.The Cauchy stress formulation can be written as:

               σ 1 = 2 J [( λ2 1 -λ2 3 )W 1 +(W 4 + W 6 )( λ1 cos θ) 2 ] σ 2 = 2 J [( λ2 2 -λ2 3 )W 1 +(W 4 + W 6 )( λ2 sin θ) 2 ] σ 3 = 2 3J [(3 λ2 3 -Ī1 )W 1 + Ī4 (W 4 + W 6 )] + U J . (14)
The third equation is used to determine the third elongation λ 3 by imposing a stress free condition σ 3 = 0, the plane stress assumption. For each model, the inversion of this later equation was numerically done using Matlab r for several values of the initial bulk modulus. The evolution of the third elongation λ 3 versus the first one λ 1 for the different constitutive equations is presented in Figs. 3(a) and 3(b). These figures show that even if the initial bulk modulus is very large κ ≥ 1 × 10 17 kPa, the elongation may be very far from the incompressible case, namely for the R model. Moreover for this model, it appears that whatever the value of the initial bulk modulus the material is thickening after a given value of the deformation (beyond 10%), which is not physically acceptable. This is due to the fact that the volumetric function U is not increasing enough compared to the isochoric energy W to satisfy the incompressibility condition. For the polynomial model, when κ ≥ 4.40 × 10 6 kPa, the material response is identical to the incompressible case, in the range of the considered deformation. These results are reinforced by the evolution of J versus λ 1 ,presentedinFigs.3(c) and 3(d). It clearly appears that the volume variation can become very high, even if κ is very large. Nevertheless, the curves also highlight that the incompressibility hypothesis seems to be verified until a given level of deformation depending on the constitutive equation.

c 1 (kPa) α(kPa) k 1 = k 3 (kPa) k 2 = k 4 ρ Ī0 4 = Ī0
The near-incompressibility formulation has also a large consequence on the stress state as illustrated in Figs. 3(e) and 3(f), which show the evolution of the first principal stress σ 1 versus λ 1 for the different models. For a given value of λ 1 , an important difference can be observed between the incompressible theoretical stress, σ 1 (inc), and the stress calculated for different values of the initial bulk modulus, σ 1 (κ). In order to better quantify this difference, the relative error defined as:

Er = σ 1 (inc) -σ 1 (κ) σ 1 (inc) (15)
for κ = {4.40 × 10 4 , 4.40 × 10 6 } kPa, is plotted in Fig. 4 for the two constitutive models. All these results point out that to keep a constant volume during strain hardening, a high volumetric function is needed compared to the isochoric one.

Obviously, this condition depends on both the expression of the volumetric energy and the value of the initial bulk modulus κ. In the present case, a classical quadratic expression for U (Eq. ( 9)) is used. However, the volumetric functions comparison in Fig. 2 shows that another volumetric function, in the range of deformation investigated, would lead to the same conclusion. Obtained results show that a large value of κ appears to be sufficient to satisfy the incompressibility condition for the P6 model in this range of deformation. When models become strongly nonlinear, as R model, large values of κ are sufficient to ensure the incompressibility condition in a restricted range of deformation. Consequently, numerical simulations outside this range of deformation can lead to large errors on the estimations of stresses in arterial wall. This point is illustrated in the next section in the particular case of AAA.

Compressibility Effects on Abdominal Aortic Aneurysm Wall Stress Analysis

The aim of this section is to investigate the influence of the near-incompressibility, i.e., the influence of the initial bulk modulus κ, on AAA wall stress distribution.

For that purpose, the model recently proposed in Ref. 8 (Eq. ( 7)) was implemented in a finite element code and numerical simulations on idealized AAA models were performed.

Implementation and validation of the model

The model (Eq. ( 7)) was implemented in the finite element code Comsol Multiphysics r . The near-incompressibility is formulated by using a mixed formulationpressure-displacement -, with quadratic displacement and linear pressure Lagrangian's elements. A direct solver (UMFPACK) was used for all the following simulations. As in Sec. 2, the classical volumetric energy U defined by the relation (Eq. ( 9)) is used and the initial bulk modulus is still given by (Eq. ( 11)).

In order to validate this implementation, equibiaxial tensile tests were first simulated for several values of the initial bulk modulus, from 220 to 3.96 × 10 17 kPa. Figures 3(b), 3(d) and 3(f) show the evolution of λ 3 , J and σ 1 versus λ 1 .L e tu s remark that the directions (e 1 , e 2 , e 3 ) in the biaxial test correspond to (e θ , e z , e r )in the arterial wall. We can observe a perfect agreement between analytical solutions (continuous line) and finite element results (gray filled marks) in the whole range of deformation.

Finite element analysis of AAA wall stress

Geometry

Finite element analysis was performed on idealized AAA models (Fig. 5)w h o s e mathematical form was proposed in Ref. 8. Briefly, the model is characterized by the aneurysm dilatation parameter (F r = R an /R a ), its eccentricity (F e = e/((F r -1)R a )) and its aspect ratio (F ℓ = L an /R an ), where R a is the healthy aorta radius, R an the aneurysm maximal radius, e the deviation of the aneurysm axis from the the healthy aorta axis and L an defines the aneurysm length. In the following, the aneurysm aspect ratio (F r =2.5), its total length (L = 100 mm), the wall thickness (t =1.5 mm), the healthy aorta radius (R a = 10 mm) and the dilatation parameter (F r =2 .75) are assumed to be constant. This value of the dilatation parameter corresponds to a maximum diameter of 5.5 cm of the aneurysm. The CAO models were built in SolidWorks r and then imported into the finite element analysis software. In the above model, the middle section in the plane (x, y)c o n s t i t u t e sa plane of symmetry. Therefore, only the upper part of the model was used in the finite element analysis. The full model can be reconstructed after.

Mesh

As mentioned above, a pressure-displacement (U-P) formulation, inherent to the near-incompressibility assumption, was adopted. Consequenly, P2-P1 elements (i.e., quadratic and linear interpolation for the displacement and the pressure, respectively) are used. These elements are fully integrated which results in zero hourglass energy.

Mesh-influence studies indicated that the numerical results are less dependent on the number of elements along the aneurysm axis. On the other side, the number of elements in the AAA wall thickness can strongly affect the numerical results.

In order to capture the strong stress gradient within the AAA thickness (close to the inner surface of the aneurysm in the case of the anisotropic model ( 7)), several simulations were performed by increasing the number of elements from 2 to 10 (the thickness of each element being the same). The obtained numerical results (Table 3) show that the influence of the mesh becomes negligible (less than 2 %) when the number of elements is larger than 6. In order to reduce the computational time, we found that similar results can be obtained by using only three elements but not with the same size. Precisely a structured mesh was adopted (Fig. 5)w i t ht h r e e elements across the AAA wall thickness, 36 elements in the orthoradial direction and 60 elements along the AAA axis. The thickness of the element close to the inner surface of the aneurysm is much thinner (10 times) than the one close to the outer surface (Fig. 5). Using this type of mesh, we found that σ 1max = 698 kPa.

Collagen fibers orientation in nonsymmetric AAA models

The model implementation was validated in § 3.1 for equibiaxial tensile tests. The difficulty in the case of an aneurysm, namely a nonsymmetric AAA geometrical model, is the definition of the fibers orientation in the arterial wall. For that purpose, on the aneurysm internal surface, the two tangential vectors (T 1 and T 2 ) were extruded through the arterial wall thickness to make them available at every integration point. These two vectors were then used to define a local orientation of the collagen fiber. This process is accomplished by using the "Extrusion Coupling Variables" in Comsol Multiphysics r .

Boundary conditions

The following boundary conditions are applied on the AAA model (Fig. 5):

           u z =0 onΓ t and on Γ sz u y =0 onΓ sy u x = 0 at the corner F σ • n = P n on Γ i σ • n = 0 on Γ e . ( 16
)
In the proximal face Γ t (z = L/2), a free radial displacement and a zero-axial displacement (u z = 0) were imposed, instead of embedding boundary conditions which lead to stresses concentration in this region of the AAA as in Refs. 8 and 14.

A symmetry condition is applied on both the planes Γ sz (z =0 )an dΓ sy (y =0 ) . The usual systolic pressure (P = 120 mmHg) was applied in the inner surface (Γ i ) of the aneurysm and the outer side (Γ e ) is stress free.

Wall stress in an axisymmetric AAA

In this section, the first principal stress is evaluated in idealized AAA geometries by using hyperelastic anisotropic model for the arterial wall. The way the nearincompressibility assumption affects this stress distribution and its maximum value in the aneurysm is emphasized. For that purpose, the stress evolution and the material volume variation are studied when the initial bulk modulus (κ), which controls the material compressibility, varies from its classical values to very large values.

Figure 6 shows the distribution of the first principal stress σ 1 (which is closely aligned to σ θθ ) and the jacobian J in the AAA for three different values of κ.

Correspondence between the values of κ and those of ν is given in Table 4.F o rt h e clarity, one-half of the aneurysm reconstructed from 2D-axisymmetric simulations is shown on these figures. We can observe that:

• both the distribution and the magnitude of the first principal stress σ 1 are strongly affected by the value of κ.W h e nκ =4 .40 × 10 3 kPa, Fig. 6 illustrates that the maximum first principal stress σ 1max is located in the inner surface at the junction between the aorta and the aneurysm and does not exceed 393 kPa. By increasing the value of the initial bulk modulus, the location of σ 1max moves towards the center of the bulged part of the AAA and its magnitude is almost multiplied by two for κ =4.40 × 10 6 kPa. • The distribution of the volume variation parameter J (Fig. 6) is also dependent on the initial bulk modulus. As expected, its magnitude decreases with increasing κ and tends towards the incompressibility case (J =1 )i nt h ew h o l e AAA for the largest value of κ. The maximum volume variations which are located in the bulged part of the AAA are typically the order of 3.9% and 0.04% for κ =4.40 × 10 3 kPa and κ =4.40 × 10 6 kPa, respectively.

These results point out that significant errors can be made on the peak wall stress estimation (magnitude and location) even in a simple AAA geometry when the incompressibility condition is not fulfilled. Such effects could be amplified in the case of patient-specific AAA models which present more complex geometric features. The numerical values deduced from this analysis are summarized in Table 4.Be si d et h e very small volume variation (less than 1% beyond κ =4.40×10 4 kPa), the maximum first principal stress could deviate with about 24% from the incompressibility case. This is due to the strong nonlinearity of the mechanical model (Eq. ( 7)). As Fig. 7 clearly shows, the magnitude of σ 1max in the AAA wall rapidly increases with increasing κ and tends towards a plateau, when κ is larger than κ =4.40 × 10 6 kPa. This value ensures an estimation of σ 1max with an error lower than 1%, whereas a value of κ =4.40 × 10 5 kPa may lead to errors close 4%. Finally, let us remark that the maximal first principal stress obtained in Ref. 8 for the same AAA model with a normal systolic pressure is about σ 1max = 520 kPa. This value is 26% smaller than the value at the plateau (Fig. 7). The value of κ =3 .96 × 10 17 kPa was assumed to be the incompressible case since it corresponds to the limit of the numerical precision of the finite element analysis software used. Anyway, for the anisotropic model under consideration, we can observe that beyond κ =4.40 × 10 6 kPa, σ 1max reaches a stable value. 15)) in terms of the initial bulk modulus κ.

Effects of non-symmetrical AAA geometry on the wall stress

The above results were obtained for an axisymmetric AAA model while due to the surrounded organs, especially the spine, the abdominal aortic aneurysm does not usually preserves its initial fusiform geometry. The loss of symmetry in the aneurysm is measured by the eccentricity parameter (F e ) which varies from 0 (axisymmetric model) to 1 ("fully" nonsymmetric model). 6). However, the magnitude of the first principal stress is strongly modified with both the aneurysm eccentricity and the initial bulk modulus. Figure 8 shows the modification of σ 1 pattern and the increase of σ 1max with the aneurysm eccentricity. This increase is nearly linear for F e ≤ 0.75 and nonlinear after (Fig. 9). In our best knowledge, such nonlinearity in σ 1max increase with the aneurysm eccentricity was not observed before, precisely when using hyperelastic anisotropic material model. We also reported values of σ 1max obtained by Rodriguez et al. [START_REF] Rodriguez | Mechanical stresses in abdominal aortic aneurysms: Influence of diameter, asymmetry, and material anisotropy[END_REF] in Fig. 9.

With the same geometrical and mechanical models, these authors found a linear increase of σ 1max with the aneurysm eccentricity increase. According to their results, σ 1max increases with 33% from an axisymmetric aneurysm to a fully nonsymmetric one. By contrast, our results exhibit an increase of σ 1max of about 135% when κ =4.40 × 10 6 kPa. Thus, for this value of κ and for the fully non-symmetric AAA model, in comparison to Ref. 8 the error committed (Eq. ( 15)) on the maximum stress is about 58%. When the same analysis is performed with κ =4.40 × 10 4 kPa, we found an increase of 109% for σ 1max .I nt h i sc a s e ,f o rF e ≤ 0.5, our results are the same order of magnitude of those in Ref. 8,s e eF i g .9. Therefore, it seems that their results would correspond to an initial bulk modulus of about 4.40 × 10 4 kPa. However, even with this value of κ, the discrepancy between the two values of σ 1max ,f o rF e =1 .0, is about 38%. This could be explained by the fact that in Ref. 8, σ 1max was only considered in the bulged region of the aneurysm. The material model (Eq. ( 7)) used in our finite element analysis was identified in Ref. 8 thanks to experimental data on AAA tissues from Ref. 7. However, whatever the geometrical model studied here, the obtained maximum stresses are very large compared to those reported from Ref. 7, about 120 kPa. This remark remains valid for many previous works. [START_REF] Geest | The effects of aneurysm on the biaxial mechanical behavior of human abdominal aorta[END_REF][START_REF] Rodriguez | Mechanical stresses in abdominal aortic aneurysms: Influence of diameter, asymmetry, and material anisotropy[END_REF][START_REF] Geest | The effects of anisotropy on the stress analyses of patient-specific abdominal aortic aneurysms[END_REF] Nevertheless, our results are the same order of magnitude as those obtained by Rissland et al. [START_REF] Rissland | Abdominal aortic aneurysm risk of rupture: Patient-specific fsi simulations using anisotropic model[END_REF] in patient-specific AAA models (with an anisotropic material model and with a uniform AAA wall thickness of 2 mm). Namely in the fully nonsymmetric idealized AAA model (F e =1.0), we found σ 1max = 1654 kPa while these authors reported a maximum Von Mises stress of 1555 kPa. Finally, let us remark that these peak values of the first principal stress are about 30% higher than those measured by Raghavan et al. [START_REF] Raghavan | Ex vivo biomechanical behavior of abdominal aortic aneurysm: Assessment using a new mathematical model[END_REF] in orthoradial direction of AAA tissues (σ ultimate ≈ 1200 kPa, in uniaxial tensile tests).

Concluding Remarks

Accurate and reliable numerical stress analysis of AAA requires suitable constitutive models for the arterial material. For that purpose, hyperelastic anisotropic models were recently proposed to capture the main features of the experimentally observed mechanical response of arterial tissues. In these models, arterial material are usually assumed to be incompressible and material parameters are identified in this way. In practice, a nearly incompressible approach, i.e., a mixed formulation pressure-displacement, to perform finite element simulations is usually adopted. Consequently, the strain energy is split into an isochoric part W and a volumetric part U , which depends on J and on the initial bulk modulus κ of the material. This parameter is used to control the material volume variation. [START_REF] Doll | On the development of volumetric strain energy functions[END_REF] In this work, an analytical analysis of the influence of κ on the mechanical response of two invariant-based anisotropic models was first performed in the case of an equibiaxial tensile test. In each model, a quadratic volumetric function usually implemented in many finite element softwares was considered. Analytical results clearly show that, in that case, large values of κ are sufficient to ensure the incompressibility condition, i.e., to estimate wall stress with a reasonable precision, in a restricted range of deformations only. Such results were then confirmed through a finite element stress analysis of idealized AAA geometries by using the model proposed in Ref. 8. Numerical results in the axisymmetric AAA model show that too small values of κ may lead to important error on the magnitude of the maximum first principal stress and also on its location. In the present case, it is shown that a value of the initial bulk modulus larger than 4.40 × 10 6 kPa is sufficient to estimate wall stresses with an error smaller than 1%. This value of κ is 10 to 100 times larger than values used in some previous numerical analyzes. [START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF][START_REF] Peyraut | Modélisation des tissus biologiques en hyperélasticitéa n i s o t r o p e-étude théorique et approche éléments finis[END_REF] In comparison to the maximum diameter criterion usually adopted for the aneurysm potential rupture, the maximum stress in the AAA wall is considered to be a more reliable indicator. [START_REF] Fillinger | Prediction of rupture risk in abdominal aortic aneurysm during observation: Wall stress versus diameter[END_REF] The above results illustrate that the aneurysm eccentricity could hugely contribute to this stress increase. However, the wall stress being mainly evaluated numerically, a particular care may be taken to include appropriate numerical, mechanical and geometrical parameters in order to obtain more realistic estimations. As we observed, from an axisymmetric to a fully nonsymmetric AAA models, the maximum stress could increase with about 135% (from 702 kPa to 1655 kPa). When the near-incompressibility condition is not fulfilled, the relative error on σ 1max varies from 24% to 44%. To obtain convenient numerical predictions of the aneurysm wall stress, the above underlined errors must be fixed. Moreover, material models have to be identified on average patient data.

Fig. 1 .

 1 Fig. 1. Fibers orientation in the arterial wall.

Fig. 2 .

 2 Fig. 2. (a) Evolution of the ratio of the volumetric function to the initial bulk modulus with the volume change; (b) evolution of the derivative of the ratio of the volumetric function to the initial bulk modulus with the volume change. 10

Fig. 3 .Fig. 4 .

 34 Fig.3. Evolution of λ 3 , J and σ 1 versus λ 1 (= λ 2 ) during an equibiaxial tensile test for different values of the initial bulk modulus: (a,c,e) P6 model, (b,d,f) R model, in this last case the gray filled marks correspond to results from the finite elements analysis.

Fig. 5 .

 5 Fig. 5. AAA geometrical model and boundaries designation and mesh model.

Fig. 6 .

 6 Fig. 6. 3D representation of an axisymmetric simulations results: Distribution of σ 1 and J in the AAA. The arrow indicates the location of σ 1max .

Fig. 7 .

 7 Fig.7. Evolutions of σ 1max and errors on its values Er (Eq. (15)) in terms of the initial bulk modulus κ.

Figure 8 exhibits σ 1

 1 distribution for F e = {0.0, 0.5, 1.0}. Unlike the symmetric model, σ 1max is located in the posterolateral surface, close to the junction of the bulged part and the healthy aorta. The risk of the aneurysm rupture is very high in these two locations as reported in Refs. 4 and 15. Here, only results for κ =4.40 × 10 6 kPa are presented since σ 1max location does not change beyond the classical value of the initial bulk modulus (4.40 × 10 4 kPa, Fig.

FFig. 8 .

 8 Fig. 8. Distribution of the first principal stress σ 1 in nonsymmetric AAA geometry when the initial bulk modulus κ =4.40 × 10 6 kPa.

FeFig. 9 .

 9 Fig. 9. σ 1max evolution versus the aneurysm eccentricity, results from Ref. 8 were reported.

Table 2 .

 2 Material constants: P6 and R refer to P6 model and R model.

Table 3 .

 3 Influence of the mesh discretization (regular mesh) within the AAA thickness (axisymmetric model).

	Number of elements	2	3	4	6	8	10
	σ 1max (kpa)	565	624	652	676	685	690

Table 4 .

 4 σ 1max in the AAA versus κ (or ν).

	κ(kPa)	νJ max(%)	σ 1max (kPa) Er(% )
	4.40×10 3	0.49000	3.90	393	44
	4.40×10 4	0.49900	0.70	539	24
	4.40×10 5	0.49990	0.40	675	4.0
	4.40×10 6	0.49999	0.40	698	0.4
	3.96×10 17	0.5 -1 × 10 -16	0.40	701	-