Henner Jakob
email: henner.jakob@inria.fr

Nicolas Loriant
email: nicolas.loriant@inria.fr

Charles Consel
email: charles.consel@inria.fr

An Aspect-Oriented Approach to Securing Distributed Systems

Keywords: D.2.11 [Software Engineering]: Software Architectures-Domain-specific architectures Design, Languages Distributed systems, Security, Pervasive computing systems, Aspect-oriented programming

The increasing size and complexity of distributed systems create a need to raise the level of abstraction for their development. This need becomes critical for pervasive computing where non-functional properties, such as security, must be guaranteed. Architecture description languages (ADLs) propose a promising approach to coping with the size and complexity of pervasive computing systems. A system is defined by a high-level description that may be used to produce a programming framework. However, non-functional properties are not specifically addressed by existing ADL works. To address this issue aspect-oriented programming is a wellproven technique to properly modularize non-functional concerns that can be dealt with by weaving dedicated code into a program.

In this paper, we present DiaAspect, an aspect-oriented language for an ADL. ADLs are a key to our approach because they expose features enabling an accurate coordination of aspects. We demonstrate the expressiveness of Dia-Aspect with two examples of security policies in pervasive computing. We also show how, combining the knowledge of the architecture description with aspect code, improves aspect weaving in the implementation code.

INTRODUCTION

Computing platforms are rapidly shifting from desktops to environments populated with networked devices. These distributed systems consist of an ever larger range of entities that offer computing power, extensive functionalities and network capabilities. The advent of these systems provides services to users at anytime and everywhere, opening up a host of application areas, ranging from home automation to assisted living.

Non-functional properties like security typically impact every aspect of a system. To cope with these concerns, one must integrate them at an early stage of the system development, as is proposed by an ADL. An ADL allows to precisely and formally define a software architecture in terms of the elements of the system, their behavior, their externally visible properties and the relationships among them [START_REF] Bass | Software Architecture in Practice[END_REF]. However, non-functional properties have received limited attention in the ADL community in that, only specific concerns have been addressed. This situation leads architects to decorate an architecture description with non-functional concerns, obfuscating the entire design. Some middlewares provide support for non-functional concerns. For example, in the Enterprise Java Bean (EJB) component model [START_REF]Entreprise Java beans specification[END_REF], EJB containers provide support for security through encryption and authentication. Such component models are either dependent on a specific middleware or provide little development support, if any.

Aspect-oriented programming (AOP) [START_REF] Kiczales | Aspect-oriented programming[END_REF] is a software engineering approach to combining non-functional concerns with software design. AOP provides techniques and tools to systematically represent, modularize and compose concerns that are crosscutting an entire system. While the integration of AOP in programming languages has been heavily studied, its usage at the architecture level is still to be explored.

This paper.

This paper presents an approach to defining non-functional concerns alongside the description of a pervasive computing system. These concerns are automatically mapped into the system implementation, using an aspect-oriented approach. Our contributions are as follows.

• We propose to extend our lightweight ADL, named DiaSpec [START_REF] Jouve | A SIP-based programming framework for advanced telephony applications[END_REF], with annotations to express non-functional concerns, as a pervasive system architecture is being defined. A dedicated programming framework is generated from an architecture description. Additionally, non-functional annotations are mapped into aspectoriented declarations. These declarations are used by an aspect-oriented engine to inject code into the generated framework to enforce declared non-functional concerns.

• To enable aspects to be declared at the architecture level, we have developed an aspect-oriented language for ADLs, named DiaAspect. Our language provides support to manipulate the common set of architecture abstractions found in existing ADLs. This support takes the form of a join point model on top of which architectural crosscutting concerns are to be expressed. Because our model covers features common to most ADLs, we believe it is widely applicable.

• We describe how aspects at the architecture level allow powerful aspect weaving in existing middleware. This situation is due to the fact that aspect weaving is not performed on an arbitrary program, but on an automatically generated framework, whose structure is known. As a result, the aspect weaver may select the optimal join point projection depending on the architecture specification and the aspect code.

• We validate our approach by securing the communication between components and enforcing access control. First, we show how to express the distribution of certificates to components of a Remote Method Invocation (RMI) distributed system to secure the communication with the secure socket layer (SSL) protocol. Second, we illustrate how to easily enforce access control lists on the interfaces of the components of a system.

Outline.

The rest of this paper is organized as follows. Section 2 presents notions and concepts of aspect-oriented software development (AOSD) and ADL. It also introduces DiaSpec, our lightweight ADL on which we illustrate our approach. Section 3 introduces two security examples used throughout the paper to illustrate our approach. In Section 4, we describe DiaAspect, our dedicated aspect-oriented language for the DiaSpec ADL. In Section 5, we illustrate our language on the two examples introduced. In Section 6, we present our implementation. We discuss related works and conclude in Section 7 and Section 8.

BACKGROUND

Let us now introduce a few notions that are relevant to the rest of the paper. We present some ADL concepts that provide a basis to identify what is required to express aspectoriented programming for ADLs. These requirements are used to define an aspect-oriented language dedicated to ADLs, and more specifically to DiaSpec. We describe how DiaSpec descriptions are used by a programming framework generator, named DiaGen, and the distributed system-based runtime support for DiaSpec descriptions, named DiaEnv. Finally, we introduce requirements for AOSD at the architecture level.

Architecture Description Languages

"A software architecture for a system is the structure or structures of the system, which consist of elements, their externally visible properties, and the relationships among them" [START_REF] Bass | Software Architecture in Practice[END_REF]. Because a software system may not be completely represented as a single structure, it is usually modeled with multiple views. ADLs allow architects to represent and specify systems in both human and machine readable languages. A precise semantics increases the coordination of multiple participants while providing tools to perform verifications on various properties and generate code and test cases.

The vast majority of ADLs use a component-and-connector idiom to represent the software system structure. Basic connectors, such as synchronous Remote Procedure Calls (RPC) and asynchronous messages, may be used to develop new connectors. In addition, ADLs may sometimes feature code generation tools and runtime support for the deployment of distributed systems. The mandatory feature, to be supported by a middleware, is a service discovery mechanism, allowing the dynamic binding of components.

DiaSpec

We now give a brief overview of our lightweight ADL to provide some background for our aspect-oriented language, DiaAspect. Notice that DiaSpec and DiaAspect are general enough to be introduced on top of most ADLs.

DiaSpec is an ADL that is part of a generative approach to supporting the implementation of pervasive computing systems. As illustrated in Figure 1, it features the Dia-Gen compiler that takes a DiaSpec specification and produces a typed programming framework. This framework closely guides the implementation of the software components. DiaGen performs various verifications of both the architecture specification and the component implementations. It partially generates a software layer, DiaEnv, to support runtime execution. In addition to performing runtime consistency checks, DiaEnv abstracts over the communication layer, allowing to transparently deploy a system implementation over multiple middlewares. 1DiaSpec models software architectures in terms of components and connectors. As illustrated in Figure 2, it is a dedicated ADL, where component interactions are included in component declarations. It differs from most ADLs where component interfaces are defined with ports. Nevertheless, this difference has no implication on our work. DiaSpec features three built-in connectors: commands, events and sessions.

• A command connector amounts to a synchronous RPC, allowing a one-to-one interaction. • An event connector corresponds to the well-known publish/subscribe paradigm, in which a publisher sends events to receivers registered for the corresponding type.

• A session connector defines an offer/answer session negotiation protocol [START_REF] Rosenberg | An offer/answer model with session description protocol (SDP). RFC 3264[END_REF], followed by a media session.

In addition to component inheritance, a component declaration includes attributes that enable component instances to be distinguished. The example in Figure 2 specifies a simple air conditioning system in which four components must be implemented: TemperatureSensor, Fan, FireManager and AirConditioning. The first two extend Device and thus inherit the attribute location. This attribute is used to refer to specific instances of components.

DiaGen and DiaEnv provide developers with a service discovery mechanism that filters component instances according to their type and their attribute values. A dedicated implementation of this mechanism is generated in conformance with the DiaSpec description. It corresponds to a typed version of existing service discovery mechanisms.

Aspect-Oriented Programming

The decomposition of software into small, meaningful, manageable and comprehensible parts has been a core idiom of software engineering for decades. A proper separation of concerns promotes reusability, traceability, adaptation, and comprehensibility. But, the relevance of concerns may vary with respect to roles: architects, developers and administrators. It may also vary depending on the stage of the software life cycle. Moreover, concerns may be constrained by the implementation language being used. For example, the object-oriented paradigm drives the decomposition of data structures into classes.

Despite the numerous software engineering approaches to system decomposition (e.g. libraries, modules, components, etc.), achieving a proper decomposition, where every concern is correctly modularized, is not possible in practice [START_REF] Tarr | N degrees of separation: multi-dimensional separation of concerns[END_REF]. Some concerns are then spread out and mixed; these are said to be crosscutting the decomposition of the system.

AOSD is a software engineering approach that focuses on the identification and representation of crosscutting concerns, and their modularization in separate units, as well as their automated composition into a complete system. AOP [START_REF] Kiczales | Aspect-oriented programming[END_REF] is a language independent paradigm, where aspects encapsulate crosscutting concerns. An aspect associates an advice, the actual code of the concern, with pointcuts that refer to the regions in the base program where the advice is to be applied. An aspect weaver then realizes the coordination of the aspects with the base program, either statically, e.g. through static code inlining, or dynamically, e.g. using the host language's reflection mechanism.

While a significant body of work has focused on the proper decomposition of system architectures [START_REF] Araujo | Early aspects: The current landscape[END_REF], it does not explicitly focus on crosscutting concerns. Properties such as security and QoS are inherently crosscutting. Moreover, these properties must be explicitly specified at design time to reason about them.

ILLUSTRATING EXAMPLES

Security is a key challenge in the development of pervasive computing systems. Even if most middlewares provide support for the implementation of security policies (e.g. encryption and authentication) developers must properly design systems to guarantee a correct implementation of these policies. Security typically impacts every aspect of a system. As a result, it must be dealt with at an early stage of the application design.

This section introduces two examples of non-functional security concerns that crosscut a system architecture. We revisit these examples in Section 5 to illustrate the expressiveness and effectiveness of our approach.

Managing Certificates for SSL communication

Authentication, confidentiality and integrity are key objectives for the security of pervasive computing systems. Their effective implementation depends on numerous factors, such as the network type, the level of trust between users, etc. Let us consider the example of the RPC. It is a wellknown mechanism to perform remote computations and to exchange messages in distributed systems. RMI is the Java application programming interface that performs the objectoriented equivalent of an RPC. While RMI offers a simple programming interface, it provides no security guarantee: RMI is built on top of the Java remote method protocol, which exchanges serialized Java objects in clear.

In this situation, good practices to achieve a secure design require the use of a secure channel. SSL is a wide spread protocol used to secure transmission in distributed systems. It performs authentication and encryption. Operationally, SSL requires participants to store certificates of trusted entities. Fortunately, the RMI API provides support for RMI over This approach is suitable for the client/server model, where every client has to hold the server certificate and the server to eventually have all client certificates. However, its application to multi agent systems, where distributed entities (dis)appear and (un)register dynamically, requires additional code to manage certificates.

In the example depicted in Figure 3, three components AirConditioning, Fan and TemperatureSensor are declared. The component Fan provides the interface setSpeed required by the component AirConditioning to regulate the temperature. Upon instantiation, instances of a component must register to the Registry and pass their certificates. In this example, a new instance of TemperatureSensor registers to the server. The Registry must then verify the certificate before the service can be successfully registered.

Whenever an instance of AirConditioning wants to invoke setSpeed, it must first obtain an instance of a Fan from the service discovery service. The Registry then returns proxies on instances of Fans with their certificates. AirConditioning may then communicate with any of these Fans via SSL.

As illustrated in Figure 3, the management of certificates for SSL communications typically crosscuts multiple aspects of a distributed system, e.g. registration and discovery, in both the component code and the built-in services (DiaEnv).

Enforcing Access Control Lists

The use of SSL connections combined with signed certificates allows the middleware to enforce authentication and encryption of communications in a distributed system. Nevertheless, the security provided by SSL communication is coarse grained. In a pervasive computing system, every entity must verify each access at a fine-grained level. The fact that an entity is authenticated does not mean that it has access to all resources.

For example in Figure 4, the AirConditioning may change the speed of a Fan by using the setSpeed method. This is a valid action and thus allowed by the AuthorizationEnforcer. On the contrary, the use of the on or off method should be denied. These methods should only be available for the FireManager (Figure 2 To forbid such uncontrolled access, the system must enforce access control lists upon requests on provided operations. As shown in Figure 4, the access list enforcement logic should be externalized from components. That is, on reception of requests, the receiver must query the access controller to verify that the caller is allowed to request a particular action from the callee.

Again, the implementation of the access control enforcement impacts multiple regions of the architecture; that is, every interface entry points.

THE DIAASPECT LANGUAGE

This section presents DiaAspect, an aspect-oriented language for the DiaSpec ADL. It first describes the join point model. A join point refers to a region in the DiaSpec architecture description and/or DiaSpec generated programming framework, where aspect code is injected.

Afterwards, we present our pointcut language. Pointcuts represent the occurrence of one or more join points. A pointcut is a predicate; it may or may not match the current join point. In addition, a pointcut may expose information specific to the underlying join points at runtime.

Finally, we introduce the DiaAspect language.

The Join Point Model

The DiaAspect join point model defines the set of events of interest in a system architecture description and its associated generated programming framework. AOP events are represented as messages in DiaEnv. For the sake of conciseness, only messages of interest between components of a DiaSpec architecture specification are listed here. Each listed message corresponds to two distinct join points in our model: one at the message emitter and one at the receiver. Figure 5 illustrate our model.

Components Registration.

The register, respectively unregister, message occurs when a component instance notifies a registry of the arrival of a component instance, respectively departure, in/from the system. Both the register and unregister messages have Component Discovery.

The discover message occurs when a component instance queries a registry for a component. A discovery is parameterized with the partial signature of the component to search for, that is, the component type and the attributes on which to apply a filter. On reply, the DiaSpec registry emits a dis-cover_r message containing an array of components.

Commands.

The command message matches the call or the reception of a command. The corresponding join point takes the following arguments: the method's name, the argument's types and the return type. The command_r join point matches command reply and has a single argument, the return type of the method.

Event Subscription.

The subscribe message matches the subscription of a component instance to an event queue. Subscriptions have two arguments, the signature of the publisher and the type of events the subscribers are interested in. Note that the publisher (the parameter) is not the same component as the receiver of the subscription message. Indeed, in DiaSpec, a subscriber subscribes to the DiaEnv, not directly to the publisher. The response corresponds to the subscribe_r message that has no argument.

Events Publishing and Reception.

The publish message occurs when a component publishes an event. It is parameterized with respect to the event type and received by a registry. Consequently, the registry emits an event message to each subscriber. That message has two parameters, the signature of the publisher and the event type.

Sessions.

There are four kinds of session messages exchanged during the life-cycle of a session: offer, answer, ack and terminate. Each message has one parameter: the type of the session. That is, the type of the data that has been negotiated during the session establishment process.

The Aspect Language

Figure 6 presents the partial BNF of the DiaAspect language. We chose a syntax similar to AspectJ, an aspectoriented system for Java. The benefits of reusing an existing syntax are well-known from both a user and an implementer perspective. The latter perspective is illustrated in Section 6 where AspectJ is shown to greatly simplify the implementation of our language.

Our language associates pointcuts with advice written in Java. An advice holds the implementation of crosscutting concerns and pointcuts select join points where an advice is to be executed. We first present our pointcut language.

Then, we introduce the runtime API supporting advice implementation with the information relative to join points matched at runtime.

The Pointcut Language

In the aspect-oriented paradigm, pointcuts act as join point selectors. In addition, to capture the occurrence of one or more join points, pointcuts may expose runtime information specific to these join points. The DiaAspect language proposes six kinds of pointcuts for which we distinguish two categories:

• Architectural pointcuts. These directly relate to artifacts defined by the architect in the system specification. That is, component declaration and component relationships (connectors).

• Built-in pointcuts. These correspond to built-in services provided by the DiaSpec runtime. That is, component registration, etc.

The following lists the pointcut featured in DiaAspect. We previously stated that a join point may catch the emission and/or reception of messages between DiaSpec components. Accordingly, pointcuts in DiaAspect distinguish join points on the emitter and/or the receiver side. If a pointcut is preceded by the send, respectively recv, keyword, it matches only the join point on the emitter, respectively receiver, side. If no keyword precedes, join points on both sides are matched.

register The register pointcut (Figure 6, line 27) matches the occurrence of register join points and the corresponding register_r join points, where the registered component matches the signature given as a parameter. For example, the pointcut recv register (Service (..)) matches every register join points received, where the component instance extends the Service component.

unregister The unregister pointcut (Figure 6, line 29) captures the occurrence of unregister join points and the corresponding unregister_r join points, where the registered component matches the signature given as a parameter. For example, the pointcut send unregister (* (..)) matches every unregister join point.

discover The discover pointcut (Figure 6, line 31) matches the discover and corresponding discover_r join points given a partial component signature (that is, the type of the component requested and the attributes on which the results are filtered). For example, the pointcut recv discover (WebCam (Location loc == "Room B", ..)) matches the reception of every discovery operation for components of type WebCam with a specific location attribute.

command The command pointcut (Figure 6, line 35) captures all command and their respective command_r join points for a given method signature. For example, the pointcut command (* OnOff.*(..)) catches every call of commands defined in the OnOff interface.

subscribe The subscribe pointcut (Figure 6, line 33) intercepts all subscribe calls and their respective sub-scribe_r joinpoints for a given component signature and event type. For example, the pointcut recv subscribe (* (..), Time) matches every subscription to an event of type Time, regardless of the publisher.

publish The publish pointcut (Figure 6 line 37) matches the occurrence of publish join points for a given event type. For example, the pointcut publish (Alert) matches every published Alert event.

event The event pointcut (Figure 6 line 38) catches the occurrence of event join points for a given event type and a publisher signature. For example, the pointcut send event (Sensor (..), *) matches the sending of any kind of event following the publication by a publisher extending Sensor.

session The session pointcuts (offer, accept, ack and terminate) (Figure 6, line 38) match the occurrence of session join points matching a session type. For example, the pointcut recv session:terminate(*) catches the reception of every session termination.

Pointcuts select join points that correspond to messages exchanged between component instances. Pointcuts filter join points on their types and the values of their arguments. In our language, an aspect may further filter the messages of interest by specifying additional clauses alongside pointcuts.

from and to As stated before, join points selected by pointcuts relate to messages exchanged between DiaSpec components. Hence, pointcut arguments correspond to the specific content of these messages. In addition to filtering join points according to message type and content, one can further restrict the collected join points by using the from and to clauses. Given a component signature, the from, respectively to, clause (Figure 6 line 23) restricts collected join points to those emitted, respectively received, by components matching the signature in that clause. In pointcut && from(Service-(Location *, ..)), join point matching is limited to those emitted by components extending Service and with at least one attribute of type Location.

if To allow aspect developers to further specify the region where to trigger advice, DiaAspect features an if clause similar to the one of AspectJ. The advice only executes if the given expression holds true. That expression must be a Java boolean expression that may refer to the DiaAspect runtime API and variables bound in the pointcut definition.

The Runtime API

An advice in DiaAspect is developed in Java. To support developers in writing advice, DiaAspect features a runtime API that exposes similar features to those in AspectJ.

The proceed method is similar to the one in AspectJ. It is a virtual method that has the same signature as the pointcut matched, and is used in around aspects. Its execution evaluates the join point matched by the aspect. For example, inside an aspect on the sending of a register message, the execution of the proceed method in the advice resumes the execution of the register.

As in AspectJ, DiaAspect also provides an advice-visible variable thisJoinPoint. It exposes reflective information about the join point that triggered the advice. We extended the thisJoinPoint to expose the architecture specific information about the current join point in addition to Java semantic information about it. For example, in the case of an aspect on a command call, thisJoinPoint.getCaller-Signature() returns the signature of the caller component.

Because aspects may be developed independently of the DiaSpec implementation code, developers can not benefit from the support of a generated programming framework when writing advice. Still, to allow developers to publish events and execute commands, we expose the Processor API, that acts as a front-end to the generated framework. Hence, developers may write an advice that interacts with DiaSpec components, while still benefiting from runtime coherency checks provided by the framework.

EXAMPLES OF CROSSCUTTING CON-CERNS REVISITED WITH DIAASPECT

This section revisits the two concerns discussed in Section 3 to illustrate DiaAspect and to show that it allows such concerns to be concisely modularized at the architecture level. We first return to the management of certificates to implement SSL communications. Afterwards, we revisit the enforcement of access control lists.

Managing Certificates for SSL Communications

Figure 7 presents the DiaAspect code that implements the distribution of certificates for SSL encrypted communications. The first aspect augments the behavior of the DiaEnv registry on service discovery to pass the certificates of the discovered services to the requesting component. The second aspect intercepts discoveries to store those certificates locally. It contains a single aspect that intercepts any receiving Dia-Spec command call (line 1). We assume that a component AuthorizationEnforcer performs the actual enforcement of an access rule. As a result, the pointcut must not match when the receiver is an AuthorizationEnforcer to avoid an infinite recursion (line 2). On the first invocation, the advice code (lines 2 to 10) obtains an AuthorizationEnforcer instance using service discovery through the DiaAspect runtime API. The advice calls the isAllowed method on the AuthorizationEnforcer (line 5). Depending on the result, the advice either throws a DiaGenSecurityException, notifying the caller of the failure, or proceeds with the original command call.

Enforcing Access Control Lists

IMPLEMENTATION

Weaving is the process of coordinating the aspect code with non-aspect code, i.e. the base program. In our approach, the aspect code refers to the architecture description. However, instead of weaving the architecture description, we inject aspects in the implementation code: the programming framework generated by DiaGen and DiaEnv. This approach exploits an architecture description to adapt the weaving process to a given generated implementation support.

This section first describes the structure and organization of the Java framework generated for a DiaSpec specification. Then, we present how DiaAspect aspects are translated into AspectJ aspects, that are woven into the implementation code.

The DiaSpec Generated Framework Organization

Given a software architecture declaration, the DiaGen compiler generates a typed framework on which to develop the distributed application.

For each command interface declaration, DiaGen includes a Java interface defining those commands. Similarly, each type of events, a component may subscribe to, creates a Java interface defining a listener to that event. For each component definition, DiaGen includes an abstract class whose implementation must provide a method definition for the provided commands and event listeners. In addition, for each couple of components sharing provided/required connectors, DiaGen generates a proxy interface of the providing component that exposes a limited view of its functionalities to the requiring component. Proxies only provide the methods of the connected interfaces. For example, given the command interface OnOff defined in Figure 2, DiaGen generates a Java interface (Figure 9) that must be declared by both abstract classes Service and ServiceProxy.

To ensure implementation independence from a specific middleware, DiaGen and DiaEnv introduce an abstract layer, as depicted in Figure 10. This abstraction layer makes a DiaSpec implementation and its corresponding generated framework completely reusable to different middlewares. As The Java interface generated for the OnOff command interface declared in Figure 2 of today, DiaSpec supports local, RMI, Web Services and SIP technologies, in addition to our simulation environment DiaSim [START_REF] Jouve | DiaSim: A parameterized simulator for pervasive computing applications (demo)[END_REF].

Given this structure, a received message first passes through the specific communication layer that formats it in a unified form before passing it to the core layer. The core layer unmarshals the message content to extract the sender information, the message type, (e.g. command, command return, event, etc), and its content. The extracted information is passed to the generated framework, which in turn dispatches it to the appropriate component proxy and finally to the component implementation. Similarly, command calls and event publishing follow the inverse path. One can note that this particular implementation offers multiple regions to intercept similar join points. The next section shows how to benefit from this strategy to optimize aspect weaving.

DiaAspect Aspects Weaving

The DiaAspect language describes aspects that coordinate join point regions. These join points refer to artifacts, i.e. components, connectors and built-in services, which are defined in a DiaSpec architecture description and its runtime environment (Section 4). An advice is defined purely in Java (exception made of the proceed keyword). The compilation and weaving process must connect the architecture pointcuts with the generated framework and the implementation code. To do so, DiaAspect aspects are translated into AspectJ code that is woven into the generated code and Dia-Env.

Translating DiaAspect aspects into AspectJ code amounts to projecting the DiaAspect join point model into the programming framework generated by DiaGen. DiaAspect pointcuts are translated into AspectJ method call pointcuts on the methods of the Java interfaces generated by DiaGen. For example, the DiaAspect pointcut command (* OnOff::* (..)) intercepts any command declared in the command interface OnOff and is translated into the following AspectJ pointcut, call(* spec.package.name.interfaces.OnOff.* (..)).

The translation of the DiaAspect from and to clauses depends on the pointcut type, i.e. send or recv pointcuts. In the case of send type pointcuts, for example in send command(* *::* (..)) && from(Manager()) && to(Service) the from clause is translated to limit AspectJ pointcut matching, to callers extending the abstract Java class Manager. This is done by using the following AspectJ code this(spec.package.name.components.Manager) and the following the invocation this(spec.package.name.proxies.ServiceProxy) This results in restricting the object on which a method is called to ServiceProxy classes.

In the case of a recv, the pointcut is limited to a caller of type proxy. In addition, the object on which the method is called must extend the corresponding abstract class.

Figure 11 presents the AspectJ code generated for the Dia-Aspect code from the example on enforcing access control lists. The DiaAspect aspect in Figure 8 is translated into an AspectJ aspect. That aspect is woven around the execution of the commandReceived method that is called whenever a DiaSpec command call is received by a component instance providing that command. The pointcut limits the matching to classes, extending the Service class (the base class for every component) and excluding components of type Au-thorizationEnforcer. The advice code is left unchanged.

The generated framework is composed of multiple layers. Depending on the aspect code, it is possible to modify the aspect projection (weaving) to shortcut these layers. This is done to optimize performance of the applications. For example, consider a DiaAspect aspect intercepting all receiving calls for a given command. Instead of weaving the aspect in the generated framework layer, we can inject it directly into the communication layer (i.e. RMI, Web Service, etc). When the advice code does not call the proceed method, this optimization avoids unmarshalling the call up to the application layer. Similarly, weaving can be specialized depending on the specific communication layer and of the aspect code.

RELATED WORK

This section briefly reviews approaches related to the modeling of crosscutting concerns in architecture design and in pervasive computing.

Multi-dimensional separation of concerns [START_REF] Tarr | N degrees of separation: multi-dimensional separation of concerns[END_REF] shows how the software artifacts, corresponding to different concerns (a.k.a. hyperslices), can be merged to generate a full application. This is the approach chosen by subject-oriented programming [START_REF] Harrison | Subject-oriented programming: a critique of pure objects[END_REF], where hyperslices are pieces of code (e.g. partial class hierarchies). AOP [START_REF] Kiczales | Aspect-oriented programming[END_REF] is quite similar, but it is asymmetric: it considers the structure of a base program and it provides pointcut languages to specify where another code crosscuts the base program and the corresponding pieces should be woven. DAOP-ADL [START_REF] Pinto | DAOP-ADL: an architecture description language for dynamic component and aspect-based development[END_REF] is an XML-based architecture description language that integrates aspects as first class entities of architecture description. The interconnection of aspects with components in DAOP-ADL is specified as evaluation rules on the interfaces of the architecture. Similarly, Pessemier et al. integrate aspect entities in the Fractal ADL, by extending the component membrane to support aspect weaving [START_REF] Pessemier | Components, ADL & AOP: Towards a common approach[END_REF]. PRISMA [START_REF] Pérez | Integrating aspects in software architectures: PRISMA applied to robotic tele-operated systems[END_REF] is an ADL that integrates an aspectoriented approach directly in the component-and-connector approach: aspects are modeled as components, allowing direct reuse of consistency checks and code generation tools. In comparison to these works, our aspect language DiaAspect is not limited to the artifacts of the ADL but also allows aspects to coordinate with the built-in services, like component registration and discovery provided by the generated framework.

Ren and Taylor present an extension of xADL, a XML ADL, for modeling security at the architecture level [START_REF] Ren | A secure software architecture description language[END_REF].

Their extension permits architects to control component instantiation, interface access, and data flow by annotating the components and connectors. Vigil [START_REF] Kagal | Trust-based security in pervasive computing environments[END_REF] is a pervasive computing middleware that enforces a trust-based security policy. To this end, it exposes services such as a certificate controller, a role assignment manager and a security agent. It also augments communication stubs with policy enforcement code. In comparison to this work, our approach allows to inject the necessary code to enforce such policies without modifying the system architecture. Moreover, by specifying these changes at the architecture level, our implementation is reusable with multiple middleware.

CONCLUSION AND FUTURE WORK

This paper presented our approach to expressing nonfunctional properties of system architectures using an aspectoriented approach. The specification of these properties, e.g. security and QoS, impacts every aspect of a software system and cannot be properly expressed in the traditional component-and-connector idiom.

To overcome this limitation, we proposed DiaAspect, an aspect-oriented language dedicated to distributed ADLs and their runtime support. We developed DiaAspect on top of both the component-and-connector idiom, and common runtime services provided by relevant ADLs. We showed that DiaAspect is expressive enough to implement concrete solutions on two widespread security problems: the distribution of certificates to encrypt communication and the enforcement of access control lists in distributed systems.

We presented the implementation of the DiaAspect aspect weaver. It injects aspect code into the programming framework generated for DiaSpec specifications. We implemented our weaving process by translating DiaAspect code into aspects written in AspectJ, a well-known aspect-oriented system for the Java programming language. We also demonstrated how our translation scheme allows to modify the DiaAspect pointcuts projection to optimize the woven code according to an architecture specification and the structure of generated frameworks.

In the future, we plan to extend the DiaAspect pointcut model to capture new join points; currently, our pointcut language only captures join points related to the relationships between components or to the runtime built-in services. A number of join points would be of particular interest: component declaration, component instantiation, connector declaration, etc. Also, we plan to integrate features similar to AspectJ inter-type declaration, for the attributes of DiaSpec components and their connectors. Another line of work we intend to explore is to define a structural and runtime model of aspect injection in system specification without introducing new artifacts in the ADL. This would allow us to transparently reuse static and runtime consistency checks performed by the DiaSpec tool suite.

  == nul l | | r es. si ze( ) == 0) r et ur n new Rem ot eSer vi ceI nf o[ 0] ; Rem ot eSer vi ceI nf o[ ] r es = new Li st ( ) ; f or ( i nt i = 0; i < r es. si ze( ) ; i ++) r esAr r ay[ i ] = r es. get ( i ) ; r et ur n r esAr r ay; } cat ch ( I nvocat i onExcept i on e) { e. pr i nt St ackTr ace( == nul l | | r es. si ze( ) == 0) r et ur n new Rem ot eSer vi ceI nf o[ 0] ; Rem ot eSer vi ceI nf o[ ] r es = new Li st ( ) ; f or ( i nt i = 0; i < r es. si ze( ) ; i ++) r esAr r ay[ i ] = r es. get ( i ) ; r et ur n r esAr r ay; } cat ch ( I nvocat i onExcept i on e) { e. pr i nt St ackTr ace( == nul l | | r es. si ze( ) == 0) r et ur n new Rem ot eSer vi ceI nf o[ 0] ; Rem ot eSer vi ceI nf o[ ] r es = new Li st ( ) ; f or ( i nt i = 0; i < r es. si ze( ) ; i ++) r esAr r ay[ i ] = r es. get ( i ) ; r et ur n r esAr r ay; } cat ch ( I nvocat i onExcept i on e) { e. pr i nt St ackTr ace( == nul l | | r es. si ze( ) == 0) r et ur n new Rem ot eSer vi ceI nf o[ 0] ; Rem ot eSer vi ceI nf o[ ] r es = new Li st ( ) ; f or ( i nt i = 0; i < r es. si ze( ) ; i ++) r esAr r ay[ i ] = r es. get ( i ) ; r et ur n r esAr r ay; } cat ch ( I nvocat i onExcept i on e) { e. pr i nt St ackTr ace( == nul l | | r es. si ze( ) == 0) r et ur n new Rem ot eSer vi ceI nf o[ 0] ; Rem ot eSer vi ceI nf o[ ] r es = new Li st ( ) ; f or ( i nt i = 0; i < r es. si ze( ) ; i ++) r esAr r ay[ i ] = r es. get ( i ) ; r et ur n r esAr r ay; } cat ch ( I nvocat i onExcept i on e) { e. pr i nt St ackTr ace( == nul l | | r es. si ze( ) == 0) r et ur n new Rem ot eSer vi ceI nf o[ 0] ; Rem ot eSer vi ceI nf o[ ] r es = new Li st ( ) ; f or ( i nt i = 0; i < r es. si ze( ) ; i ++) r esAr r ay[ i ] = r es.

Figure 1 :

 1 Figure 1: The DiaSpec development process.

Figure 2 :

 2 Figure 2: A DiaSpec example specification

Figure 3 :

 3 Figure 3: Example of the distribution of certificates

Figure 4 :

 4 Figure 4: Enforcement of Access Control Lists

Figure 7 :

 7 Figure 7: DiaAspect code managing certificates in a RMI with SSL distributed system (excerpt)

Figure 8

 8 Figure 8 contains the DiaAspect code to enforce access control lists on interfaces declared in a DiaSpec architecture.It contains a single aspect that intercepts any receiving Dia-Spec command call (line 1). We assume that a component AuthorizationEnforcer performs the actual enforcement of an access rule. As a result, the pointcut must not match when the receiver is an AuthorizationEnforcer to avoid an infinite recursion (line 2).

Figure 8 :

 8 Figure 8: DiaAspect code enforcing access control lists (excerpt)

Figure 6 :Figure 9 :

 69 Figure 6: DiaAspect BNF

Figure 10 :

 10 Figure 10: DiaSpec/DiaGen generated code structure.

Figure 11 :

 11 Figure 11: AspectJ code generated for the Di-aAspect aspect from the example on enforcing access control lists (excerpt)

get ( i ) ; r et ur n r esAr r ay; } cat ch ( I nvocat i onExcept i on e) { e. pr i nt St ackTr ace( ) ; r et ur n nul l ; } cat ch ( Di aGenExcept i on e) { e. pr i nt St ackTr ace( ) ; r et ur n nul l ; } } Implementation Specific communication layer pr ot ect ed i nt [ ] get Ser vi ces( Fi l t e t r y { Li st r es = cal l ( f i l t er ) ; i f ( r es == nul l | | r es. si ze( ) == 0) r et ur n new Rem ot eSer vi ceI nf o[ 0] ; Rem ot eSer vi ceI nf o[ ] r es = new Li st ( ) ; f or ( i nt i = 0; i < r es. si ze( ) ; i ++) r esAr r ay[ i ] = r es. get ( i ) ; r et ur n r esAr r ay; } cat ch ( I nvocat i onExcept i on e) { e. pr i nt St ackTr ace( ) ; r et ur n nul l ; } cat ch ( Di aGenExcept i on e) { e. pr i nt St ackTr ace( ) ; r et ur n nul l ; } } pr ot ect ed i nt [ ] get Ser vi ces( Fi l t e t r y { Li st r es = cal l ( f i l t er ) ; i f ( r es == nul l | | r es. si ze( ) == 0) r et ur n new Rem ot eSer vi ceI nf o[ 0] ; Rem ot eSer vi ceI nf o[ ] r es = new Li st ( ) ; f or ( i nt i = 0; i < r es. si ze( ) ; i ++) r esAr r ay[ i ] = r es. get ( i ) ; r et ur n r esAr r ay; } cat ch ( I nvocat i onExcept i on e) { e. pr i nt St ackTr ace( ) ; r et ur n nul l ; } cat ch ( Di aGenExcept i on e) { e. pr i nt St ackTr ace( ) ; r et ur n nul l ; } } pr ot ect ed i nt [ ] get Ser vi ces( Fi l t e t r y { Li st r es = cal l ( f i l t er ) ; i f ( r es == nul l | | r es. si ze( ) == 0) r et ur n new Rem ot eSer vi ceI nf o[ 0] ; Rem ot eSer vi ceI nf o[ ] r es = new Li st ( ) ; f or ( i nt i = 0; i < r es. si ze( ) ; i ++) r esAr r ay[ i ] = r es. get ( i ) ; r et ur n r esAr r ay; } cat ch ( I nvocat i onExcept i on e) { e. pr i nt St ackTr ace( ) ; r et ur n nul l ; } cat ch ( Di aGenExcept i on e) { e. pr i nt St ackTr ace( ) ; r et ur n nul l ; } }

	DiaGen/DiaSpec program building chain		Provided	Extensible
				Generated	User code
					Deployment program
	"uses"	"uses"	"uses"	"uses" "uses"
				"uses"

  one argument, the signature of the component arriving/leaving. Our model distinguishes the component issuing the (un)register message from the actual component that is arriving/leaving. The register_r, respectively unregister_r, message corresponds to the acknowledgment from the Di-aSpec registry to a component, after a register, respectively unregister, message. These messages have no argument.

	Component	DiaEnv	Component
	Registration Registration
	register(comp sig)
	register r
	Discovery Discovery
	discover(comp sig)
	discover r
	Commands Commands
	command(signature)
	command r(ret type)
	Publish / Subscribe Publish / Subscribe
	subscribe(comp sig,Evt t)
		subscribe r
	publish(Evt t)
		event(comp sig,Evt t)
	Sessions Sessions
	offer(SType)
	answer(SType)
	ack(SType)
	terminate(SType)
	Figure 5: DiaAspect Join Point Model

Currently, DiaEnv supports Web Services, RMI, local and SIP deployments.


Thierry Faug
email: thierry.faug@cemagref.fr

Rémi Beguin

Benoit Chanut

Mean steady granular force on a wall overflowed by free-surface gravity-driven dense flows

Keywords:

We studied free-surface gravity-driven recirculating flows of cohesionless granular materials down a rough inclined plane and overflowing a wall normal to the incoming flow and to the bottom. We performed 2D spherical particle discrete element simulations using a linear damped spring law between particles with a Coulomb failure criterion. High-frequency force fluctuations were observed. This paper focuses on the mean steady force exerted by the flow on the obstacle versus the macroscopic inertial number of the incoming flow, where the inertial number measures the ratio between a macroscopic deformation timescale and an inertial timescale. A triangular stagnant zone is formed upstream of the obstacle and sharply increases the mean force at low incoming inertial numbers. A simple hydrodynamic model based on depth-averaged momentum conservation is proposed. This analytical model predicts the numerical data fairly well and allows us to quantify the different contributions to the mean force on the wall. Beyond this model, our study provides an example of the ability of simple hydrodynamic approaches to describe the macroscopic behavior of an assembly of discrete particles, not only in terms of kinematics, but also in terms of forces.

INTRODUCTION

The importance of granular materials in geophysics and in various industrial processes has resulted in extended research on granular flows at the frontier between physics, soil mechanics and fluid mechanics. Flows around obstacles and the force those granular flows are able to exert on the obstacle are important issues when applied to storage and conveying bulk solids [1] and also in geophysical flows [2][3][4]. Granular drag on objects was approached by the pioneering work of Wieghardt [5,6], who provided an early systematic experimental and analytical study of granular flows around immersed objects. In those free-surface flows, he showed the drag force depended relatively little on velocity and he observed a pile-up in front of the immersed object as well as depression of the surface in the wake. A recent experimental study on dense granular flows around an immersed cylinder in a vertical chute [7] showed that the mean drag is independent of the mean upstream velocity. The subject also contributed experiments on two-dimensional flows in a vertical bin around various inserts [1]. The authors recorded velocity contours and observed stagnant zones in front of the inserts.

Many recent studies have focused on the drag force on small obstacles such as cylinders in the case of rapid -dilute -granular flows [8,9], including interstitial gas effects [10,11].

In the dilute regime, the force is proportional to the square of the incoming velocity as predicted from kinetic theory for granular gas. The importance of shock waves formed in front of the obstacle has been shown for this rapid regime [12,13] and has been the subject of valuable recent studies [18][19][20][21][22][23][24].

The quasic-static regime has also been investigated by Wieghardt [5,6]. A series of recent studies examined drag forces on immersed obstacles in a rotating granular bed at very low velocities [25][26][27]. Experimental measurements showed that the drag force is linearly dependent on the cylinder diameter, quadratically dependent on the depth of insertion and independent of velocity [25]. The effect of obstacle shape on drag force and its influence on jamming behavior has been analyzed [27] as well as stick-slip fluctuations in granular drag [26]. All these studies refer to horizontal motion and showed that the average pressure is scaled as the immersion depth. A recent study investigated the average drag forces for obstacles slowly plunging into and withdrawing from shallow beds of different granular materials, and showed that the average pressure is scaled as the immersion depth exhibiting 2 (August 7, 2009) a power law behaviour rather than a linear behavior [28].

Since the studies of Wieghardt [5] and Tuzun and Nedderman [1], little has been studied on the drag of dense granular flows and the effects of the stagnant zone observed in front of the obstacle. In the present paper, dense granular flows will systematically refer to assemblies of grains moving at higher velocity than the grains in the quasi-static regime mentioned just above but at lower velocity than the grains in the rapid dilute regime. This regime refers to the so-called granular liquid regime mentioned in [29]. A typically sized roughness close to the diameter of the flowing grains is a condition to obtain this type of dense regime.

In the dense granular regime, when the obstacle has a typical size close to the flow depth, large stagnant zones or dead zones are able to form in front of the obstruction, while grains that are non trapped inside the dead zone continue to flow around the obstacle. These stagnant zones observed in the dense regime are not accompanied by a granular jump (large discontinuity in flow depth and velocity), as described earlier for rapid flows on smooth beds [14][15][16][17]. Furthermore, free-surface curvatures and large flow-depth gradients are observed in the vicinity of the obstacle. All these conditions result in a mean force on the obstacle that cannot be expressed as (i) a simple function of the incoming depth-averaged velocity of the flow, as is usually observed for granular flows around small immersed obstacles in the dilute regime [8,9,12,13] or (ii) a simple function of a depth of insertion (equal to the flow depth when the small object is located at the base of the flow) in the quasi-static regime [25,27,28,30]. This paper focuses on 2D free surface gravity-driven flows of dense granular materials down an inclined plane and overflowing a wall normal to the bottom. These flows correspond to a flow geometry typically encountered in geophysical flows when avalanches overtop protection dams [4]. A roughly triangular dead zone, whose length greatly depends on the slope inclination, is formed upstream of the wall and largely influences the mean force on the wall. We investigated this mean force using discrete numerical simulations combined with a hydrodynamic modeling approach. The numerical simulations method and results are presented first. Then we describe a simple hydrodynamic analytical model showing that the mean force is the sum of many contributions: the weight of the granular material (motionless grains inside the dead zone and moving grains above), the incoming momentum force, the incident pressure force and the basal friction force inside the dead zone. Based on simple arguments to calibrate the free parameters, the predictions from the proposed

3 (August 7, 2009)
analytical model are successfully compared to numerical results.

DISCRETE NUMERICAL SIMULATIONS

Numerical methods

Simulation method and contact law

Numerical simulations are carried out using the molecular dynamics method, as introduced by Cundall [31] and largely used to simulate dense granular flows (see for example [32][33][34]). We defined the normal contact force as the sum of two contributions, an elastic one N e (spring) and a viscous one N v (dashpot). In the simplest case used here, the normal elastic force N e ij is assumed to depend linearly on the overlap δ ij (the displacement of the spring):

N e ij = k n δ ij ,
where k n is the normal stiffness. The dashpot contributes a normal dissipative force N v proportional to the time derivative of the overlap:

N v ij = C ij δij
, where C ij is a normal damping coefficient that may be related to a coefficient restitution e in a binary collision of cohesionless grains [35,36]. The sum of both contributions N ij is restricted to being repulsive, i.e., tensile normal forces are not allowed, as the particles are assumed to be non cohesive. Finally, the normal contact force may be expressed as follows: The simulated system is an assembly of spheres of average diameter d and density ρ P (average mass m = 1/6ρ P πd 3 ). A small polydispersity (±10% in size) is introduced to prevent crystallization. The average particle diameter is d = 1 mm and the particle density is ρ P =2450 kg m -3 corresponding to the glass material. The normal stiffness k n is calculated from the Young modulus E [37]: k n ≈ Ed. The Young modulus of the glass material is E = 69GP a, which gives k n =6.9 10 7 N m -1 but k n was reduced to k n = 10 4 N m -1 with respect to the limit of rigid grains (no influence of the stiffness above this limit), in order to decrease the calculation time.

N ij = k n δ ij + C ij δij
It has been previously shown that the dimensionless number N 1 = k n /(ρd 3 γ2 ) (where γ is the shear rate and ρ the density) has to be greater than 10 4 in order to conclude that we are in the limit where grains behave as if they are perfectly stiff [36]. N 1 can be thought as the inverse of the square of the ratio between the shear velocity u γ = γd and the sound velocity c s = k/(ρd) [36]. The condition N 1 > 10 4 gives u γ /c s < 10 -2 (condition 1). This condition 1 controls the overlap generated by flow inertia and gives the following condition on the shear rate (with k n = 10 4 N m -1 and a density ρ = φρ P where φ = 0.55 is a typical volume fraction): γ < 860 s -1 . The dense granular flows investigated in this paper occur below this limiting value of the shear rate (as well as many similar experimental granular flows [38]). We have also to consider the overlap generated by quasistatic mechanisms such as the compression of force chains. Similarly to the number N 1 , we can define the dimensionless number N 2 = k n /(P d) for spherical grains, which can be interpreted as the inverse of the square of the ratio between the characteristic velocity of grains rearrangement u P = P/ρ (where P is the pressure) and the sound velocity c s . By analogy to the condition 1, the condition N 2 > 10 4 , which gives u P /c s < 10 -2 , has to be respected (condition 2). For dense granular flows, the typical pressure level is maximum at the bottom and may be approximated by the hydrostatic pressure P ≈ ρgh. The condition 2 gives the following condition on the flow depth (with k n = 10 4 N m -1 ): h < 75d. The granular flows investigated in the paper are below this limiting value of the flow depth (thereafter, in Fig. 1b, are provided the typical values of the flow depth). Furthermore, recent discrete simulations on dense granular flows showed that the effect of variations in k n is minimal as long as k n > 2.10 5 mg/d [39]. With our values for d and m, it gives k n > 2570 5 (August 7, 2009) N m -1 , which is also compatible with the chosen value for k n .

In the rigid grains limit, it has been shown that the macroscale behavior does not depend on µ (except for frictionless grains µ = 0) and e (except for the extreme values e = 0 and e = 1) for plane shear flows without gravity [34]. It has also been shown that the flow behavior (bulk density and velocity profile) depends little on the local particle friction if µ > 0.5 for steady flows down a rough inclined plane (see Fig. 9 in [32]). In our simulations, we used µ = 0.5 corresponding to a typical value for glass beads. The restitution coefficient e was also shown to only slighly influence the bulk density and the velocity profile of steady dense granular flows down rough inclined planes if e < 0.9 (see figure 10 in [32]). In the dense granular regime, the macroscopic flow quantities are known to be almost insensitive to the coefficient of restitution e [29,40]. In our simulations, we chose e = 0.5. It has also been shown that the k t /k n ratio has a very small influence [32,36]. Here we used k t /k n = 1/2, which gives k t =5 10 3 N m -1 . The commercial Particle Flow Code was used here (PFC 2D version 3.0 [START_REF]Particle Flow Code in 2 Dimensions (PFC 2D version 3[END_REF]). In our simulations, the time step was kept constant at 4 10 -6 s, which guaranteed the stability of our calculations.

Flow geometry

We defined a channel of length L c = 1700d and a reservoir of length L r = 150d. The roughness of the bottom was made with grains of a mean diameter d = 1 mm with the same properties as the grains in motion. The channel was constantly fed by releasing grains into the reservoir. Grains moved out of the reservoir by an aperture of constant height H r = 35d (see inset in Fig. 1a). After a certain simulation time, the imposed height H r at the exit of the channel led to a constant mass flow rate. We investigated two types of flow in a large range of slope inclinations (16 • ≤ θ ≤ 32 • ): (i) flows with no obstacles in order to characterize our flows regarding the existing literature reporting dense granular flows down rough inclined planes and (ii) flows overflowing a vertical wall normal to the bottom and located at a distance 900d from the exit of the reservoir. The following section presents the results.

6 (August 7, 2009)

Numerical results

Control flows without an obstacle

Dense granular flows down rough inclined planes have been thoroughly studied and an overview of recent progress and remaining unanswered questions can be found in [29,40]. This section reports a summary of the main results obtained from our discrete simulations of control flows with no obstacles. We show that these results are compatible with previous reports in the literature [40]. The existence of steady and uniform flows in a large range of slope inclinations [θ min , θ max ] has been shown for dense granular flows down a rough inclined plane [32,40,42]. The function h stop (θ) has been defined as the thickness of the granular material left by a steady and uniform flow at slope inclination θ and can be measured. The angle θ max is the angle for which no grains are able to stay on the inclined plane (h stop = 0) and the angle θ min is the angle for which h stop tends towards ∞ (no flow). We determined the function h stop (θ) numerically, which is depicted in Fig. 1a. The curve can be fitted by the

function h stop (θ) = Bd[(θ max -θ)/(θ -θ min )]
where B is a constant depending on material properties [40]. We found the following values: θ min ≈ 14 • , θ max ≈ 24 • and B = 2.3. It is found to be in agreement with previous data on 2D flows [32,40,43].

A single dimensionless number called the inertial number has been defined and may be interpreted as the ratio [40] between (i) a microscopic time scale d/ P/ρ P , which represents the time it takes for a particle of density ρ P to fall in a hole of size d under the pressure P giving the typical time scale of rearrangements, and (ii) a macroscopic time scale 1/ γ linked to the mean deformation ( γ is the local shear rate). I is the square root of the Savage number [44] or the Coulomb number [45,46] previously introduced in the literature. The averaged inertial number I may be estimated from depth-averaged velocity and thickness of the flows assuming a Bagnold-like velocity profile [40]:

I = 5ūd h 3/2 √ 2g cos θ (1)
Fig. 1b shows the velocity profiles obtained at different slope inclinations and at the position x/d = 900. The Bagnold profiles (shown in dash lines in Fig. 1b) roughly fit the numerical data, except at the base of the flow as it has been previously reported [40]. Fig. 1c shows the macroscopic inertial number I, calculated from Eq. ( 1), along the distance x from the reservoir. For steady and uniform flows, we can define a unique value of I far enough from the reservoir. Above θ max , we observed a continuous increase in the inertial number I along x showing that the flows are non-uniform flows. For these gradually accelerated flows, we defined the inertial number I at two positions: x/d = 900, corresponding to the location of the wall in the second set of simulations, and x/d = 1500 (maximum position up to which we took measurements because for x/d > 1500 the flow properties were influenced by the boundary limit due to the end of the inclined plane at x/d = 1700).

The robust scaling has been shown for the variation of the basal friction coefficient µ * versus the inertial number I for plane shear flows [34]:

µ * = µ * min + bI (2)
where µ * min = tan θ min and b is a constant. Another scaling is proposed in the literature for µ * (I) for inclined planes [47,48]:

µ * = µ * min + ∆µ I 0 /I + 1 (3)
where ∆µ = tan θ max -µ * min and I 0 is a constant typically equal to 0.3 for glass beads. The depth-averaged equations, introduced by [49] and recently revisited [29,50] in the context of shallow granular flows down an inclined plane, allow us to estimate the effective friction coefficient µ * . The acceleration is balanced by the gravity parallel to the plane, the tangential stress between the fixed bottom and the flowing layer, and a pressure force related to the thickness gradient [29,49,50]. The momentum balance is reduced to the following equation in steady 2D flow conditions:

β ∂hū 2 ∂x = tan θ -µ * -k ∂h ∂x gh cos θ. (4)
The β factor, in the acceleration term, is related to the velocity profile and is defined by

1 h h 0 u 2 (y)dy = β 1 h h 0 u(y)dy 2 .
It is generally taken to be equal to 1, which corresponds to the exact value for plug flows [49]. It can be calculated provided an assumption on the shape of the velocity profile: β = 4/3 for linear velocity profiles and β = 5/4 for Bagnold velocity profiles [51]. The k factor, in the thickness gradient term, is the ratio of the normal stress σ xx to the normal stress σ yy , classically introduced for dense granular flows [49,50].

For steady dense granular flows, k can be chosen equal to 1 (isotropic fluid-like behavior)

8 (August 7, 2009)
as shown by previous studies [32]. From Eq. ( 4), we could estimate the effective friction coefficient at a given position x 0 /d:

µ * (x 0 ) ≃ tan θ - △h △x x=x 0 - β gh cos θ △(hū 2 ) △x x=x 0 (5)
We performed our calculations with

β = 1, △x = 200d, △h = h x=x 0 -100d -h x=x 0 +100d and △(hū 2 ) = (hū 2 ) x=x 0 -100d -(hū 2 ) x=x 0 +100d
. Fig. 2 shows the variation of the effective friction coefficient µ * versus the macroscopic inertial number I for different positions x/d. All the data collapse into a single curve whatever the position x/d and the slope, which shows that the existence of a unique relationship between µ * and I is still valid for accelerating flows investigated here. We also reported the laws from Eq. ( 2) with b = 0.5 and from Eq. ( 3).

These laws should be valid for steady and uniform flows when tan θ = µ * . This is the case for the Eq. ( 2) which predicts fairly well the data for θ < θ max (low I values) but fails for larger values of θ (large I values). The Eq. ( 3) was shown to be in good qualitative agreement with the data in the sense it predicts the asymptotic saturation of the friction at large values of I (µ * → µ max ). We reported the equation prediction for three values of θ max : and φmin = 0.4.

In conclusion, steady and uniform flows (θ ≤ θ max ≈ 24 • ) as well as slightly non-uniform flows (θ ≥ θ max ) were investigated and characterized with regards to recent new insights on dense granular flows. We also measured the velocity and density profiles over depth and found good agreement with previously published data on dense granular flows down rough inclined planes [29,32,40], i.e., Bagnold-like velocity profiles and a constant volume fraction. We carried out a second set of numerical simulations for which we added a vertical wall normal to the bottom at the location x/d = 900, as depicted in Fig. 3. The obstacle height H was systematically taken equal to the flow depth h 1 of the control flow: H/h 1 = 1 (for θ ≥ θ max , h 1 was defined at the position x/d = 900). The following analyzes the force on the wall resulting from these dense granular flows.

Mean force on the obstacle Fig. 4a shows an example of how the force exerted on the obstacle evolves over a 0.4s time-duration window. The force is characterized by high-frequency fluctuations with possible high amplitudes. We systematically observed a fluctuating chain forces network formed inside the dead zone being the source of the high-frequency fluctuations. Force chains in granular media are highly fluctuating physical processes that have been widely studied in the literature (see, for example, [52][53][54][55][56]). This issue is not discussed in detail in this paper. The total normal force F n (t) on the wall, at a given time t, is the sum of each force f i n (t) mobilized at the contact points between the wall and each bead i in contact with the wall:

F n (t) = i=1...n f i n (t) ( 7
)
where n is the number of beads in contact with the wall at the given time t. The number n is closely approximated by the ratio between the obstacle height and the mean particle on the wall. When the ratio f i=imax n /F n (t) is close or equal to one (point A in Fig. 4b), it means that the total force is mainly due only one bead in contact with the wall and trapped in a force chain. Inversely, when no bead in contact with the wall is trapped in a force chain, each bead has a similar contribution to the total force F n (t) and the ratio f i=imax n /F n (t) is lower (point B in Fig. 4b).

We represented the moving average over 0.02 s and 0.

= 1 β h 1 0 1 2 ρ 1 u 2 (y)dy = 1 2 ρ 1 u 2 1 h 1 (dynamic contribution) or F hydro 0 = h 1 0 ρ 1 g(h-y) cos θdy = 1 2 ρ 1 gh 2 1 cos θ (hydrostatic contribution). u 1 = h 1
0 u(y)dy is the mean velocity over flow depth h 1 of the incoming flow of density ρ 1 . The following section proposes a hydrodynamic model to describe this behavior in greater detail and quantify the different contributions to the mean force in each regime defined by the inertial number I.

HYDRODYNAMIC MODELING

Momentum conservation in a control volume

In a steady regime, the variation of momentum in a fixed control volume of fluid V is equal to the sum of the volume forces (the weight of the granular material here) and of the external forces resulting from the elements in contact with the control volume:

S ρu(u • n)dS = - S pndS + S σ • ndS + V ρgdV (8)
We use a single integral symbol for external forces ( S ) and a double integral symbol for the volume forces ( V ) because we consider variables per unit width (in N m -1 ). n is the unit vector normal to the surface S of the control volume V , ρ the density, u the velocity, g the gravity acceleration, p the pressure and σ the stress tensor. We apply the momentum conservation on the control volume V 0 defined on the schematic view given in Fig. 5. If we consider sections (S 1 ) and (S 2 ) defined in Fig. 5, Eq. ( 8) gives:

βQ 0 m (u 2 + u 1 ) = P + P 0 + P 1 + P 2 + R -F (9)
where u 1 = 1 h 1 S 1 udS and u 2 = 1 h 2 S 2 udS are the mean velocities at sections (S 1 ) and (S 2 ) defined in Fig. 5, and also shown in Fig. 3. S 1 is the section normal to the bottom and the incoming flow at the position where the obstacle does not create disturbance in the upstream incoming flow. S 2 is the section normal to the outgoing flow at the top of the obstacle and making an angle α with the bottom. Mass flow rate conservation gives

Q 0 m = ρ 1 u 1 h 1 = ρ 2 u 2 h 2
where Q 0 m is the mass flow rate, h 1 is the thickness of the incoming flow at section (S 1 ), ρ 1 its density, and h 2 is the thickness of the outgoing flow at section (S 2 ) and ρ 2 its density. Let us define all the forces involved in momentum conservation. P is the weight of the control volume: P x = ρ 0 gV 0 sin θ (component of P in the x-axis direction) and P y = -ρ 0 gV 0 cos θ (component of P in the y-axis direction). ρ 0 is the mean density of the granular material in the control volume assumed to be close to the incoming flow density:

ρ 0 ≈ ρ 1 .
The force P 0 due to the intersitial fluid (air) at the free surface of the flow are ignored: P 0 ∼ 0. The pressure force P 1 due to the incoming fluid on section (S 1 ) is parallel to the bottom: P 1x = 1 2 kρ 1 gh 1 cos θ (P 1y = 0), where k is the earth pressure coefficient classically introduced for dense granular flows [49,50]. The pressure force P 2 resulting from the outgoing fluid on section (S 2 ) is assumed to be negligible (P 2 ≈ 0). This assumption is argued by the presence of a flying and dilute granular jet (downstream the wall) which does not exert any pressure on the more dense incoming flow at section (S 2 ) (zero pressure boundary condition). The reaction of the bottom R has two components. R y is the y-axis component and R x is the x-axis component corresponding to the mean basal friction force assumed to be proportional to the normal force (y-component of the weight of the volume friction µ zm later in the paper. F is the mean force exerted on the wall by the granular flow.

Numerical simulations showed a curvature of the free-surface of the flow in the zone of influence of the obstacle. Here, for the sake of simplicity, we assume that the free-surface can be modeled by a straight line between sections (S 1 ) and (S 2 ) in Fig. 5. The disturbance created by the obstacle in the flow is characterized by a free-surface inclined at an angle α sl and the dead zone inclined at an angle α zm with the bottom (see Fig. 5). This assumption and the relation H = h 1 allow us to estimate the volume V 0 :

V 0 = 1 2 h 1 2 + δ h cos α L -h 1 δ 2 h tan α ( 10
)
where L is the upstream distance at which the flow is no longer influenced by the obstacle (i.e., the distance between the section S 1 and the obstacle: see Fig. 5 and also the simulation picture in Fig. 3) and δ h = h 2 /h 1 is the ratio between the flow depths in sections (S 1 ) and (S 2 ). The obstacle's influence zone can be determined by the following equation:

L = h 1 δ h (cos α)(tan α sl ) (11)
The projection of momentum conservation on the x-axis and y-axis gives the following equations for F n and F t :

F n = F dyn + F pressure + F weight-f riction (12)
F t = -βρ 1 u 2 1 h 1 δ u sin α -ρ 1 gV 0 cos θ + R y (13)
with:

F dyn = βρ 1 u 2 1 h 1 (1 -δ u cos α) (14)
F pressure = 1 2 kρ 1 gh 2 1 cos θ (15)
F weight-f riction = ρ 1 gV 0 (sin θµ zm cos θ)

where δ u = u 2 /u 1 represents the variation of mean flow velocity between sections (S 1 ) and (S 2 ). As mentioned above, mass flow rate conservation gives

Q 0 m = ρ 1 u 1 h 1 = ρ 2 u 2 h 2 .
Assuming that the density is unchanged (ρ 2 ≈ ρ 1 ), we have δ u ≈ 1/δ h . The following 13 (August 7, 2009) proposes simple empirical equations to close the model and determine the angles α sl , α zm and α (defined in Fig. 5), the basal friction µ zm , and the velocity ratio δ u .

Free-surface and dead zone angles

The more inclined the slope is, the higher the free-surface angle upstream of the wall. Therefore, we can assume that α sl is a simple affine function of the slope inclination θ: α sl = aθ + b, where a and b are empirical coefficients to be determined. Two asymptotic conditions may be considered for each incoming regime: (i) the uniform -dense -regime and (ii) the non-uniform -more dilute -regime. It gives different values of the parameters a and b for each regime. We note a 1 and b 1 (respectively, a 2 and b 2 ) the coefficients in the dense (respectively dilute) regime. First, let us consider the incoming dense uniform regime for θ < θ max . When θ = θ min , no steady regime is possible. The influence zone propagates increasingly upstream of the obstacle (L → ∞), which corresponds to a value of α sl close to zero. This asymptotic condition gives: 0 = a 1 θ min + b 1 , which implies α sl = a 1 (θθ min ).

When increasing the inclination angle θ, the angle α sl increases until it reaches the critical value α c sl for θ = θ max , which gives the asymptotic condition: α c sl = a 1 (θ max -θ min ). Second, we consider the incoming non-uniform and dilute regime for θ > θ max . We can consider the situation for which the influence zone tends toward zero, which implies that the angle α sl tends toward π/2 at a critical slope inclination θ c . This asymptotic condition gives: π/2 = a 2 θ c + b 2 . α sl may be expressed as: α sl = a 2 (θθ c ) + π/2. When decreasing the inclination angle θ, the angle α sl decreases until it reaches the value α c sl defined above for θ = θ max , which gives: α c sl = a 2 (θ max -θ c )+π/2. Combining all these asymptotic conditions, we obtain the following changes in α sl according to the slope angle θ:

α sl (θ < θ max ) = α c sl θ max -θ min (θ -θ min ) ( 17
)
α sl (θ > θ max ) = α c sl -π 2 θ max -θ c (θ -θ c ) + π 2 ( 18
)
The values of the critical angles θ c and α c sl will be discussed and determined later in the paper. In section (S 2 ), the grains at the free surface have a velocity parallel to the free surface and grains at the top of the obstacle (base of the launching downstream jet) have a

α = α sl + α zm 2 (19)
Considering that the length of the dead zone formed upstream of the obstacle is identical to the length of the influence zone defined according to the free-surface (length L of the control volume in Fig. 5), the angle α zm defined in Fig. 5 can be determined by the following implicit equation:

tan α zm = tan α sl 1 + h 1 H δ h cos( α sl +αzm 2 )
+ 1

(20)

Variation in velocity and depth

Incoming flows encounter a local decreased slope created by the dead zone upstream of the obstacle, which leads to a decrease in velocity. The term δ u defined in momentum conservation (Eq. ( 9)) is then smaller than one: δ u = u 2 /u 1 < 1. We do not have precise measurements of the velocity profiles at section S 2 , which would allow us to test a possible theoretical prediction of δ u based on the assumption of Bagnold velocity profiles at sections S 1 and S 2 . Future work is needed in that sense. Here we simply propose an empirical law based on the following statement: the larger the deflecting angle α, the higher the expected velocity decrease is. Thus, we simply assume that the relative velocity reduction is proportional to the angle α: △u/u 1 = (u 1 -u 2 )/u 1 = κα, where κ is a coefficient to be determined. This equation allows us to quantify δ u :

δ u = 1 -κα (21)
Let us note that the free-surface angle α is strongly correlated to the obstacle's influence zone L: the larger α, the longer L. We could have also proposed a simple relation between δ u and L. By mass flow rate conservation and assuming an unchanged density, the flow depth is increased (δ h = h 2 /h 1 > 1) and can be approximated by:

1/δ h ≈ 1 -κα (22) 15 (August 7, 2009)
The κ coefficient is a model parameter that will be discussed later. at the beginning of the dead zone and a value which is less and less determined when the obstacle is approached. For simplicity reasons, the constant value that can be given will be discussed below.

NUMERICAL DATA COMPARED TO HYDRODYNAMIC MODELING

Calibration

To compare the predictions of the proposed hydrodynamic model to the data from discrete numerical simulations, we need to validate the empirical laws proposed to close the model and to determine the different parameters introduced in the model.

A first empirical law has been proposed for the prediction of the angle α sl . Fig. 6a gives the angle α sl versus the inclination angle θ. Each group is clearly described by a linear fit and allows us to determine the coefficients a i and b i defined previously. Identifying these coefficients to Eqs. ( 17) and ( 18), we can obtain the critical angles α c sl and θ c . We have both unknown parameters and four equations, giving the following results: α c sl ≈ θ min and θ c ≈ π/2. Eqs. ( 17) and ( 18) can be then expressed:

α sl (θ < θ max ) = θ min θ max -θ min (θ -θ min ) (23)
α sl (θ > θ max ) = π 2 - θ min -π/2 θ max -π/2 (π/2 -θ) (24)
Fig. 6b shows the angle (θα zm ), which corresponds to the mean angle of the dead zone with the horizontal (see Fig. 5), versus the slope inclination θ. (θα zm ) is remarkably close to θ min whatever the slope inclination. The angle α zm can be approximated as follows:

α zm ≈ θθ min , instead of using the implicit Eq. ( 20). The prediction of Eqs. ( 23) and ( 24)

compared to the angle α sl estimated from numerical simulations are also reported in Fig. 6b.

The linear law proposed to describe α sl is in quite good agreement with the numerical data.

The inset graph in Fig. 6b shows the numerical angles compared to proposed predictions in terms of α zm and α sl as a function of I. This plot is obtained by considering the value of the inertial number that the flow would have in absence of obstacle at x/d = 900, corresponding to the position of the obstacle.

The relation θα zm = θ min gives a simple equation (instead of Eq. ( 11), which demands knowledge of δ h and α sl ) to estimate the length L of the obstacle's influence zone. Indeed, we can consider that it is identical to the length of the dead zone (as shown in Fig. 3), which gives: L ≈ H/ tan(θθ min ). Fig. 6c shows the prediction of the latter equation compared to the length of influence measured directly from numerical simulations. Predictions are satisfying regarding the simple assumptions made. L was estimated graphically using the simulation pictures similar to the typical picture depicted in Fig. 3. It is the distance between the obstacle and the section (S 1 ), and coincides quite well with the length of the dead zone represented by grains with a velocity smaller than a threshold speed v t (see black colored grains in Fig. 3). This threshold velocity was typically chosen equal to 5% of the depth-averaged velocity of the granular flow in absence of obstacle (v t ≈ 0.05ū).

A second empirical law has been proposed to predict the velocity ratio δ u . It implies a κ coefficient difficult to determine. κ is expected to vary and display different behaviors in the different flow regimes (dilute, dense and quasi-static regimes). One approximation considers the value of κ approximated from the dilute regime for which collisional interactions are dominant. If the effects of the ambiant fluid are not considered, the only source of velocity reduction in the dilute regime stems from collisions between grains and is expected to be proportional to the restitution coefficient. Therefore, we can assume that the u 2 /u 1 ratio scales as e in the dilute regime. Furthermore, assuming that the length of the dead zone tends to vanish in the dilute regime, it corresponds to a limit angle α of π/2, giving the following limit condition according to Eq. ( 21): δ u = e = 1 -κ e (π/2). The value of κ is then equal to κ e = (1 -e)/(π/2) = 0.32 with e = 0.5. Fig. 6d shows the prediction of Eq. ( 22)

with κ = κ e = 0.32 compared to the depth ratio h 2 /h 1 measured directly from numerical simulations. The prediction from the simple law proposed for the depth ratio δ h , is not perfect, but the order of magnitude given by this κ e value is satisfying regarding the crude assumption made. The assumption of unchanged density probably has an effect on results.

Density ρ 2 at the top of the obstacle (free boundary and granular jet at section S 2 ) may be smaller than density ρ 1 , which could explain why the flow depth ratio is underestimated.

Another parameter in the model is the basal friction µ zm . We measured shear and normal basal stresses directly from numerical simulations. Fig. 7a shows the ratio between basal shear and normal stresses τ xy /σ yy versus the position x/d upstream of the obstacle (example for θ = 26 • ). Outside of obstacle's influence zone (-x/d > L/d with x < 0 in Fig. 7a), the friction is close to (tan θ) as was expected, and it decreases when approaching the obstacle with high fluctuations close to the obstacle. The variation of τ xy /σ yy is nontrivial.

We interpolated all the curves by polynomial functions, named P µ (x), which allowed us to estimate an averaged friction. µ zm is estimated from the averaged value obtained at a distance 1/3L upstream of the obstacle: µ zm = P µ (x = -1/3L) with x < 0 (Fig. 7a). It corresponds to the position of the dead zone's center of gravity assumed to be triangular in shape. Fig. 7b gives µ zm versus I. The values measured are relatively constant and found to be very close to the value calculated from the minimum angle θ min , which is also reported in Fig. 7b. In spite of the complex force chains network inside the dead zone leading to high fluctuations of the local basal friction τ xy /σ yy , it is possible to consider a continuum approach and define an averaged basal friction µ zm equal to (tan θ min ). This value seems compatible with the observed geometry of the dead zone provided that the angle of the dead zone with the horizontal (θα zm ) is equal to θ min , as discussed above (see Fig. 6b) and shown again in the inset graph in Fig. 7b.

All the parameters are now determined, which allows us to compare the predictions of the hydrodynamic analytical model to the numerical data from discrete simulations in terms of the force exerted on the wall.

18 (August 7, 2009)

Fig. 8a shows the normal force on obstacle F n versus the slope inclination θ obtained from both the hydrodynamic modeling approach and the data from discrete numerical simulations. to the normal stress σ yy , was chosen equal to 1, which corresponds to isotropic material conditions. A value derived from a Mohr-Coulomb plasticity [49,57] does not provide a better prediction, similar to the results from previous studies [21,32].

Fig. 8a also reports the different contributions to the normal forces given by Eq. (12).

At high inertial numbers I, the contribution from the momentum term (F dyn ) is dominant, proving an inertial regime. At low inertial numbers I, the pressure term (F pressure ) and the difference between weight and basal friction (F weight-f riction ) become dominant contributions.

In the dense regime (0 < I < 0.3 -0.4), Fig. 8b shows that the force scales like the hydrostatic force F hydro 0 but it is four times greater than F hydro 0 . The contribution of the term (F weight-f riction ) to the total normal force F n is dominant because of the large increase in length of the dead zone, as shown in Fig. 6c.

An important result from discrete simulations, shown by the Fig. 4a, is that the tangential force F t on the obstacle is found to be close to zero whatever the slope inclination: F t ≈ 0.

From Eq. ( 13), we can estimate the reaction of the bottom in the y-axis direction (R y ).

Fig. 8c shows the reaction of the bottom R y versus the inertial number I. The reaction R y is strongly increased at low inertial numbers because -P y = ρ 1 gV 0 cos θ is very large (the influence zone upstream of the obstacle being very long). This behavior is proof of the transition toward the quasi-static regime and prevents the occurrence of a steady flow regime when approaching the angle θ min . In this quasi-static regime, the momentum force in the y-axis direction is negligible and the reaction R y is equal to the weight of the material upstream of the obstacle in the y-axis direction: βρ 1 u 2 1 h 1 δ u sin α << R y = -P y . At higher inertial numbers, there exists a critical value I c for which the reaction of the bottom R y is equal to zero (see inset in Fig. 8c), which indicates that the momentum force in the y-axis direction is exactly balanced by the weight of the material stored upstream of the obstacle in the y-axis direction. Above I c , the three contributions are balanced according to the Eq. ( 13) and R y < 0. The fact that F t is found to negligible, and the consequences it has on the prediction of the bottom reaction R y , should be further studied by initiating experimental studies.

CONCLUSION

This paper has described the results on the mean force on a vertical wall caused by freesurface dense granular flows down an inclined plane. We presented numerical data from 2D spherical particle discrete simulations performed for a wide range of slopes. The normal force on the obstacle was characterized by high-frequency fluctuations with high amplitudes due to a fluctuating force chains network inside the dead zone formed upstream of the wall. This paper has focused on the mean normal force. We considered a hydrodynamic modeling approach based on momentum conservation in a control volume. The analytical prediction was compared to numerical data and quite accurately reproduced the behavior of the force.

Although the numerical simulations display a more complex shape of the dead zone and the free-surface upstream of the obstacle (particularly when approaching the quasi-static regime), we assumed a triangular shape for the dead zone. These simple assumptions, made to close the theoretical model, provided good results. Indeed, the mean angle of the dead zone and the length of the influence zone upstream of the obstacle were captured quite well by the simple laws proposed. The predictions of the angle α sl and of the variation in depth were not perfect and will need further investigation. However, we believe that the hydrodynamic analytical model proposed in this paper with simple empirical arguments to close the model is strikingly satisfying in the sense that it is able to predict the exact order of magnitude of the mean granular force computed from discrete numerical simulations if relevant values are provided for the parameters needed in the model. All these parameters were estimated from simple assumptions with no fitting processes. We can conclude that the dynamics 20 (August 7, 2009) of the incoming flow (u 1 , h 1 and φ 1 ) and the granular material properties (θ min = 14 • , k = 1, β = 5/4 and e = 0.5 here) made it possible to estimate the mean steady granular normal force on a vertical obstacle for uniform incoming flows (θ min < θ < θ max ) but also for gradually accelerated flows (θ > θ max ). Besides numerical simulations showed that the mean tangential force exerted on the wall is close to zero, which could initiate future experimental studies. Beyond this analytical model, this study provides a new example of the ability of hydrodynamic modeling approaches to describe the mean macroscopic behavior (kinematics and depth-averaged force) of an assembly of rigid grains. It has to be kept in mind that the continuum model proposed in this paper has been developed on the basis of stationary flow conditions (recirculating flows) with an obstacle height close to the incoming flow depth.

Neither the force due to rapid granular flows impinging high walls in the presence of a upward traveling shock wave [17], nor the peak impact force due to a granular avalanche front, can be reproduced by the continuum model proposed in this paper. Further investigations will be devoted to the situation of dense granular avalanches (finite volume) overflowing small planar obstacles, in order to see if the continuum model is able to catch the force exerted by the decelerating avalanche tail (after the occurrence of the peak impact force). We distinguished grains of individual velocity v smaller (black grains) and greater (white grains) than a threshold velocity v t , which gives the shape of the so-called "dead zone". The threshold velocity was typically taken equal to 5% of the depth-averaged velocity of the granular flow in the absence of obstacle: v t ≈ 0.05ū. The vertical dashed gray line represents the transition at θ = θ max .

FIG. 5: Sketch of the control volume V 0 (hatched zone) inside which the momentum conservation is applied. h 1 and u 1 are flow depth and mean velocity of the incoming flow at section (S 1 ). The section (S 1 ) is normal to the bottom and represents the beginning of the influence's zone of the obstacle. The wall height is H and L is the length of the obstacle's influence zone, i.e., the distance between the section (S 1 ) and the foot of the obstacle. h 2 and u 2 are flow depth and mean velocity of the flow at section (S 2 ). The section (S 2 ) is normal to the main direction of the outgoing flow of velocity u 2 . We defined the angle α between u 2 and the bottom. θ is the bottom slope. Assuming a triangular shape for the dead zone, we also defined α zm as the angle of the dead zone with the bottom and α sl as the angle of the free-surface (inside V 0 ) with the bottom.

  . Two microscopic parameters are then needed: the normal stiffness k n and the damping coefficient C ij (or restitution coefficient e) to define the normal intergranular force. The tangential component T ij of the contact force is implemented in terms of a linear spring T ij = k t u ij , where k t is the tangential stiffness and u ij the displacement of the spring. T ij is restricted to absolute values smaller than µN e ij according to friction between grains described by a Coulomb condition enforced with the sole elastic part of the normal force. Here the local particle friction parameter µ is introduced. When this threshold is reached, the tangential relative motion is regarded as sliding with the sliding friction µN e ij (directed opposite the tangential relative velocity). Finally the tangential contact force may be expressed as follows: T ij = min µN e ij , k t u ij . Two miscroscopic parameters are then needed: the tangential stiffness k t and the local particle friction parameter µ to define the tangential intergranular force. 4 (August 7, 2009)

24 • , 26 •

 26 Fig.2bshows the mean volume fraction φ3D versus the inertial number I for different positions x/d. Our simulations were performed in two dimensions with spherical beads (no disks). The mean volume fraction φ3D was calculated from the mean volume fraction φ2D assuming φ3D =2 3 φ2D if we compare a sphere of diameter d included in a cube of identical size d to a disk of diameter d included in a square of size d. Again, all the data are shown to collapse into a single curve whatever the position x/d and the slope, which supports the existence of a unique relationship between φ and I. Fig.2bpresents a comparison of the results to the following law proposed in the literature[29,48]:

  diameter: n ∼ H/d. When a bead j in contact with the wall is trapped in a relatively long force chain, the force f i=j n increases substantially and largely contributes to the total force F n (t), which gives a force peak. Such a mechanism can be illustrated in Fig.4bwhich gives the change over time of the ratio of the force f i=imax n to the total normal force F n (t), where f i=imax n is the instantaneous force which corresponds to the bead exerting the maximum force 10 (August 7, 2009)

  1 s in Fig.4a. The graph indicates that the moving average tends toward a non-time-dependent hydrodynamic force, one of the main points of the present paper. Fig.4cshows the mean steady normal force F n on the wall versus the slope inclination θ. The results show that the normal force increases sharply on low slope inclinations (low value of I) when decreasing the slope inclination.The dramatic increase in force when θ is close to θ min is proof of a transition toward a quasi-static regime and the flow stopping. The dead zone formed upstream of the obstacle tends to move upward indefinitely and the mean force is greatly increased by the weight of the granular material stored upstream of the wall. The mean force is quasi-constant at intermediate slope inclinations. The last data point seems to show a slight increase in force at the highest slope inclination which was investigated (highest value of I). Here we have a transition toward an rapid -or inertial -regime for which the force is mainly controlled by the flow velocity. These regimes are also visible in the inset in Fig.4c, which shows the rescaled normal force F n /F 0 versus the inertial number I, where F 0 is a typical flow force equal to F dyn 0

3. 4 .

 4 Basal friction in the dead zoneThe quantification of the basal friction µ zm is a crucial and non-trivial point. Let us note T b the basal friction force. If we consider grains at the base of the flow at section (S 1 ), these grains are in movement and the basal friction force is known according to Coulomb sliding condition: T b = (tan θ)P y for θ min < θ < θ max (steady and uniform flows). For slightly accelerating flows at slopes inclinations larger than θ max , the basal friction force T b at section S 1 is slightly smaller than (tan θ)P y (positive flow-depth gradient in the x-axis direction). If we consider motionless grains at the base of the flow immediately upstream of the obstacle, the basal friction is locally undetermined (the non-sliding condition gives: T b < µ zm P y ). The basal friction µ zm is then expected to vary between a value equal to µ * (µ * = tan θ for steady and uniform flows and µ * < tan θ for gradually accelerating flows)

Fig. 8b showsor F 0 = F hydro 0 . 4 (

 04 Fig. 8b shows the same results in terms of the rescaled force F n /F 0 versus the inertial number I, where F 0 = F dyn 0

FIG. 1 :FIG. 3 :

 13 FIG. 1: (a) Thickness h stop normalized by d versus the slope inclination θ: open squares represent our numerical data, plain squares represent data from [43] reported also in [40], the continuous line gives the fitted function h stop (θ)/d = B[(θ max -θ)/(θθ min )] with θ min = 24 • , θ min = 14 • and B = 2.3. Inset: sketch of the simulated recirculation system. (b) Velocity profiles measured at x/d = 900 for control flows at different slope inclinations; the dotted lines give the Bagnold profiles. (c) Inertial number I versus normalized distance x/d for different slope inclinations.

24 (

 24 August 7, 2009)

FIG. 4 : 1 2 ρ 1 u 2 1 h 1 1 2 ρ 1 gh 2 1

 41111 FIG. 4: (a) Change over time for the instantaneous normal (dark gray) and tangential (light gray) components of the force per unit width exerted on the obstacle: example for θ = 28 • . Lines (continuous line for F n and dotted line for F t ) represent the time-averaged values. (b) Change over time of the ratio of the instantaneous force f i=imax n

FIG. 6 :

 6 FIG. 6: (a) Free-surface angle α sl versus the slope inclination θ: numerical data and linear fits for both dense and dilute regimes. (b), (c) and (d): predictions from the empirical laws compared to numerical data. (b) Angles θα zm and θα sl versus the slope inclination θ; inset graph: angles α zm and α sl versus the inertial number I; (c) length of the influence zone upstream of the obstacle normalized by the obstacle height L/H versus I; (d) depth ratio δ h = h 2 /h 1 versus I (with κ = κ e = 0.32).

Acknowledgments

The authors acknowledge financial support from the French National Research Agency (ANR-MONHA). The authors are grateful to Mohamed Naaim for fruitful discussions and to David Bertrand for his experience and help in conducting the numerical simulations with the commercial code PFC.