Pierre Matri
email: pmatri@fi.upm.es

Philip Carns
email: carns@mcs.anl.gov

Robert Ross
email: ross@mcs.anl.gov

Alexandru Costan
email: alexandru.costan@inria.fr

María S Pérez
email: mperez@fi.upm.es

Gabriel Antoniu
email: gabriel.antoniu@inria.fr

SLoG: Large-Scale Logging Middleware for HPC and Big Data Convergence

Cloud developers traditionally rely on purposespecific services to provide the storage model they need for an application. In contrast, HPC developers have a much more limited choice, typically restricted to a centralized parallel file system for persistent storage. Unfortunately, these systems often offer low performance when subject to highly concurrent, conflicting I/O patterns. This makes difficult the implementation of inherently concurrent data structures such as distributed shared logs. Yet, this data structure is key to applications such as computational steering, data collection from physical sensor grids, or discrete event generators. In this paper we tackle this issue. We present SLoG, shared log middleware providing a shared log abstraction over a parallel file system, designed to circumvent the aforementioned limitations. We evaluate SLoG's design on up to 100,000 cores of the Theta supercomputer: the results show high append velocity at scale while also providing substantial benefits for other persistent backend storage systems.

I. INTRODUCTION

Traditionally, cloud system designers have used purposespecific storage services for their data storage needs. These include key-value stores [START_REF] Decandia | Dynamo: Amazon's highly available key-value store[END_REF]- [START_REF] Escriva | HyperDex: A distributed, searchable key-value store[END_REF], wide-column databases [START_REF] Khetrapal | HBase and Hypertable for large scale distributed storage systems[END_REF]- [START_REF] Lakshman | Cassandra: A decentralized structured storage system[END_REF], or streaming message brokers [START_REF] Kreps | Kafka: A distributed messaging system for log processing[END_REF], [START_REF] Hunt | ZooKeeper: Wait-free coordination for internet-scale systems[END_REF]. On the contrary, high-performance computing (HPC) platforms rest on on a much more constrained set of storage primitives, typically limited to parallel file systems [START_REF] Schwan | Lustre: Building a file system for 1000-node clusters[END_REF], [START_REF] Carns | PVFS: A parallel file system for linux clusters[END_REF] or transient burst buffers [START_REF] Liu | On the role of burst buffers in leadership-class storage systems[END_REF], [START_REF] Wang | Burstmem: A high-performance burst buffer system for scientific applications[END_REF]. The low availability of local storage on the compute nodes of most supercomputers today and the lack of administrative access give little opportunity for users to deploy the storage system they need.

As the boundaries between HPC and big data analytics (BDA) continue to blur [START_REF] Reed | Exascale computing and Big Data[END_REF], new challenges arise. A critical objective set in this convergence context is to foster application portability across platforms [START_REF] Fox | Big Data, simulations and HPC convergence[END_REF]. Let us consider an application running in a cloud context, which uses a specialized storage service such as a distributed shared log [START_REF] Balakrishnan | CORFU: A shared log design for flash clusters[END_REF]. Porting this application to HPC (e.g., to leverage specific hardware capabilities) is challenging. Indeed, deploying cloud-oriented shared log services on HPC is often not possible because of the unique specificities of these platforms. While one could consider a shared log abstraction over the available parallel file system, providing the illusion of a storage paradigm atop another is extremely difficult considering conflicting sets of constraints and APIs between the two models [START_REF] Greenberg | MDHIM: A parallel key/value framework for HPC[END_REF], [START_REF] Kalmegh | Graph database design challenges using HPC platforms[END_REF]. Shared log storage is indeed one of these storage models that are both unavailable and difficult to implement on HPC platforms using the available storage primitives. Yet, in scientific applications, distributed logs could play many roles, for example for in situ visualization of large data streams, collection of telemetry events for computational steering, or data aggregation from arrays of physical sensors. A shared log enables multiple processes to append data at the end of a single byte stream. Unfortunately, in such a case, the write contention at the tail of the log is among the worst-case scenarios for parallel file systems, yielding problematically low append performance [START_REF] Tsujita | Multithreaded two-phase I/O: Improving collective MPI-IO performance on a Lustre file system[END_REF]. Thus, efficient streaming log middleware is needed that can leverage the primitives available on HPC platforms for persistence.

We present the SLoG1 shared log middleware. It provides high append velocity over distributed file systems by leveraging I/O parallelism and proxying, circumventing the inherent limitations of HPC file systems. It features pluggable backends that enable it to leverage other storage models such as object stores or to transparently forward the requests to a shared log storage system when available (e.g., on cloud platforms). SLoG abstracts this complexity away from the developer, fostering application portability between platforms. We evaluated SLoG's performance at scale on a leadershipclass supercomputer, using up to 100,000 cores. We measured append velocities peaking at 174 million appends per second, far beyond the capabilities of any shared log storage implementation on HPC platforms. For these reasons, we envision that SLoG could fuel convergence between HPC and big data. Our contributions are as follows.

• We leverage several use cases for distributed logs, (Section II), and present the design of SLoG and the diverse techniques supporting its scalability (Section III).

• We evaluate SLoG on up to 100,000 cores of the Theta supercomputer [START_REF] Harms | Theta: Rapid installation and acceptance of an XC40 KNL system[END_REF], exhibiting high append velocities on HPC platforms using both file system and object storage backends (Section IV). We review the related work (Section V) and conclude with a discussion of future research directions and perspectives (Section VI).

II. THE CASE FOR SHARED LOG MIDDLEWARE

Distributed logging is a simple storage model that naturally serves a wide variety of applications. On HPC systems it allows the collection of metadata, telemetry, or monitoring events generated by simulation that can be used either live for computational steering or offline for analytics or visualization. In the case of non-bulk-synchronous applications such as discrete event simulators, shared logging is a natural model for storing the generated events that can arrive in no predefined order. Shared logs also prove critical for real-world BDA use cases, such as large-scale monitoring systems [START_REF] Newman | MonALISA: A distributed monitoring service architecture[END_REF].

In telemetry event management, the events generated by an application carry data related to its state, its performance, or other application-specific information. Analyzing these events while the application is running can be useful, for instance to verify that the application behaves as expected or to run in situ visualization [START_REF] Dorier | Damaris/viz: A nonintrusive, adaptable and user-friendly in situ visualization framework[END_REF]. This information can also be used for in-situ analytics [START_REF] Ziegeler | In situ analysis as a parallel I/O problem[END_REF] or computational steering [START_REF] Yi | In-situ visualization and computational steering for large-scale simulation of turbulent flows in complex geometries[END_REF]- [START_REF] Butnaru | A parallel and distributed surrogate model implementation for computational steering[END_REF]. The latter can in turn influence the simulation in real time, for example increasing its resolution on a particular zone or rerunning part of the simulation with different parameters. We illustrate this in Figure 1. In this paper we focus on the application of telemetry event management to computational steering.

The case for a shared log on HPC systems. In order to deliver its real-time promise, computational steering requires collecting live telemetry data from all the nodes running the simulation. In parallel, one or multiple processes retrieve, analyze, and act on the gathered events in chronological order. This chronological order is critical for the accuracy of computational steering and is the key motivation for a distributed shared log. Moreover, a distributed shared log between the simulation and the steering helps mitigate bursting generation rates and decouple these two components.

Why are shared logs challenging? Parallel file systems are well known to be generally bad at handling concurrent, conflicting writes originating from different clients [START_REF] Tsujita | Multithreaded two-phase I/O: Improving collective MPI-IO performance on a Lustre file system[END_REF], [START_REF] Mclay | A user-friendly approach for tuning parallel file operations[END_REF]- [START_REF] Tessier | TAPIOCA: An I/O library for optimized topology-aware data aggregation on large-scale supercomputers[END_REF]. Unfortunately, distributed logs are by design a highly concurrent data structure in which multiple processes essentially compete to write to the tail of a unique object. Because distributed logging is not available on HPC platforms today, few or no applications are actually built around this data structure. Yet, we advocate that distributed logging could become a relevant storage model for HPC applications in the next years, as strong as it is today for BDA.

Hence, we posit that efficient distributed shared log middleware would foster the development of applications taking advantage of computational steering at scale. It needs at the same time to evade the limited write concurrency permitted by these storage systems, while abstracting this complexity away from the developer. Moreover, when other storage models are available (e.g., object storage [START_REF] Matri | Could blobs fuel storage-based convergence between HPC and big data?[END_REF]), it should take advantage of their different capabilities to further enable performance without requiring any application modification.

Forward cursor ⇒ (a) (b) (c) (d)
Fig. 1: Telemetry event collection. Simulation ranks push events to a log (a). Such events are used for online tasks (e.g., in situ visualization, computational steering) (b), or are stored for further offline processing (c). Computational steering in turn controls the simulation (d).

III. SLOG DESIGN AND IMPLEMENTATION

In this section we describe the key principles that drive the design and implementation of SLoG, large-scale shared log middleware that tackles the aforementioned challenges. SLoG adapts best practices from CORFU [START_REF] Balakrishnan | CORFU: A shared log design for flash clusters[END_REF] or Kafka [START_REF] Kreps | Kafka: A distributed messaging system for log processing[END_REF] to solve the limitations of HPC parallel file systems.

The two key design principles behind SLoG are stream partitioning (Section III-A) and write proxying (Section III-B). The former reduces write contention by distributing events to multiple partitions. The latter eliminates concurrency, channeling appends through intermediate processes each mapped to a partition. Readers are organized in groups, obtaining events from a dedicated read proxy (Section III-C). Figure 2 shows this architecture. SLoG also acts as a pass-through interface on platforms where shared log implementations are available, fueling cross-platform application portability.

A. Stream partitioning

In order to reduce write concurrency, a log is split into multiple partitions. This concept is similar to that of Kafka. When a writer generates an event, it is directed to one and only one partition. Each partition hence receives only a subset of all append operations. A partition is mapped to a single, dedicated file on the parallel file system, providing persistence as well as fault tolerance. The possible write throughput for a partition is hence determined by the maximum throughput the underlying file system can support for a single file. Increasing the number of partitions for a log increases the total throughput supported by the log by increasing write parallelism. The number of partitions for each log is set by the user at the time the log is created, and that number cannot be changed at a later time.

Multiple strategies are possible for determining to which partition each write is directed. These include uniform distribution (e.g., round robin) or a static mapping of writers to specific partitions. In practice, uniform distribution is the default on most systems [START_REF] Kreps | Kafka: A distributed messaging system for log processing[END_REF]. We adopted the same default for SLoG, thus ensuring that with p partitions, each will receive n/p events out of a total of n.

w 1 w 2 w 3 wn ... p 2 p 1 p 2 d 1 d 2 r A 2 r A 1 r A 3 r B 1 r B 2

B. Write proxying

Unfortunately, enabling clients to append directly to the files that back each partition of the system would yield problematically low performance, because of the inherent low level of concurrency permitted by parallel file systems. Indeed, while these systems excel at dealing with intensive loads of isolated writes, they are known to poorly handle conflicting operations. Yet, efficient conflict management is essential for fast file appends. The POSIX specification defines an O_APPEND mode for this purpose, but it is inherently more difficult to implement in a distributed environment considering the strong consistency semantics of POSIX I/O, and hence it has not been the focus of parallel file system optimization.

Append proxying in SLoG circumvents this issue. This principle is inherited from best practices for HPC I/O optimization. With each partition, SLoG associates a dedicated process named append proxy. Through this process all append operations are channeled from the writers to a specific partition. This eliminates the costly file system access concurrency at the expense of a small increase in append latency.

The proxy process receives append operations from the writers. It maintains internally an exclusive cursor on the file that hosts the partition it is associated with, as well as an offset counter. For each append operation, it generates a write to the file at the next free offset, along with the current timestamp, and increments the counter by the size of the event. For optimization, incoming requests can be batched and written in bulk to the file system. Also, the proxy can be used to transparently apply further data transformations (e.g., data compression).

C. Read proxying

SLoG distributes clients in reader groups. For each reader group, a dedicated read proxy process distributes the events to the readers of this group while ensuring a total delivery order between events within a group.

Reader groups serve the principal purpose of increasing the read parallelism permitted by the system. A number of partitions are associated with each reader group. The events from these partitions are distributed to all the readers the group contains. This distribution is performed observing the order of the timestamp associated with each event in the log.

Such an approach permits great flexibility in terms of grouping strategies. For example, read parallelism can be increased by mapping distinct groups to disjoint sets of log partitions. With this configuration, a partial order of the messages is guaranteed inside a reader group. It is illustrated by Figure 3a. Mapping different groups to a common set of partitions enables multiple readers to access the same events in parallel. This is the case when multiple applications access the same log simultaneously (e.g., computational steering and in situ visualization). It is illustrated by Figure 3b. Combinations of these strategies are possible in order to increase parallelism while using multiple applications.

Read proxies serve a purpose symmetric to that of write proxies. With each reader group a single read proxy process. Its role is twofold: (1) coordinate reads from the file system for all the partitions the group is mapped to and (2) distribute the events to the readers of the group. It does so guaranteeing partial delivery ordering of the events.

A read proxy holds a cursor to all files that hold each of the partitions the group is associated with. Read requests from the clients are served in timestamp order, advancing all the file cursors in parallel. As an optimization to significantly reduce read latency, the proxy reads from the file system in chunks that are cached in its memory.

If the write proxy is configured to apply data transformations such as compression, the read proxy is responsible for applying the reverse operation, hence rendering these transformations completely transparent to the client.

IV. EVALUATING SLOG AT SCALE

In this section we discuss early experiments proving the relevance of the design of SLoG. We focus on append throughput, which is harder to scale than read throughput.

A. Leveraging a high-end experimental platform

All experiments carried out in this paper are based on the Theta supercomputer [START_REF] Harms | Theta: Rapid installation and acceptance of an XC40 KNL system[END_REF], hosted at the Argonne Leadership Computing Facility (ALCF) of Argonne National Laboratory. Theta is a 9.63 petaflops system, ranked 16 in the TOP500 as of June 2017. It is composed of 3,624 nodes, each containing a 64-core Knights Landing Xeon Phi CPU with 16 GB of high-bandwidth in-package memory (MCDRAM), 192 GB of DDR4 RAM, and a 128 GB SSD. Networking leverages a Cray Aries interconnect with Dragonfly configuration.

The Theta storage system is based on Lustre. Its 10 TB parallel file system spans 170 storage servers. Theta is also connected to the GPFS file system of the Mira supercomputer colocated in the same facility using high-speed interconnects. This configuration allows for experiments on both Lustre and GPFS from the same experimental platform.

The minimal API required for computational steering is straightforward to benchmark. Indeed a single operation, append, is performed by all event generator ranks. We develop a simple MPI-based benchmark that sequentially generates events at maximum rate from a number of generator ranks. These 1 KB events are appended to the log using that operation. The computational steering process is simulated by a single dedicated rank that reads the logs concurrently, repeatedly reading events from the log. Since we focus only on storage I/O, the data is discarded immediately by this process. We advocate that this benchmark, although simple, is well suited for evaluating the maximum performance achievable by the system. The append operation causes a very high write contention in distributed logging systems, a large number of processes concurrently push data at maximum rate to a single log puts a very high pressure on the storage system.

B. Evaluated persistence backends

Besides the Lustre and GPFS file systems, we experimented using SLoG with two additional backends: Zlog [START_REF]Zlog[END_REF], which implements the CORFU [START_REF] Balakrishnan | CORFU: A shared log design for flash clusters[END_REF] log design over the Ceph file system [START_REF] Weil | Ceph: A scalable, high-performance distributed file system[END_REF]. At its core lies a network sequencer, an atomic counter initialized with the size of the file and incremented by the size of each write.

Týr [START_REF] Matri | Týr: Blob storage meets built-in transactions[END_REF], which is an object store targeted at coping with highly concurrent I/O patterns. It features a lightweight transaction protocol that natively enables atomic operations such as appends to be performed at high speed.

These two systems provide a baseline performance compared with (1) a purpose-built shared log storage system (i.e., Zlog) and (2) an object storage system built for concurrent I/O patterns (i.e., Týr).

C. Baseline append performance with a single partition

We start evaluating the baseline performance of SLoG using a single log. This demonstrates the relevance and efficiency of write proxying.

Theta is one of the few supercomputing platforms offering access to the local SSD drives of the compute nodes, we can compare Týr and Zlog on the machines. The application uses between 1 and 8,192 cores as append clients. We deploy Týr and Zlog using between 64 and 1,536 cores for the storage servers, on different machines from the append clients (we measured this ratio to be optimal for our setup). We measure the number of appends per process per core over 10 minutes. We do not use any write batching or compression. Týr and Zlog are configured with a replication factor of 3.

We plot the results on Figure 4. Overall, they show that the performance achieved with SLoG is competitive with that of Týr or Zlog when the number of concurrent writers is low, despite the overhead of write proxying and the limitations of the underlying file system. As soon as the I/O capabilities of a single write proxy are saturated, however, the write performance offered by SLoG with both Lustre and GPFS reaches a plateau, slightly decreasing as more writers are added. Týr and Zlog exhibit similar behavior, which is due to the high write contention at the tail of the log.

D. Varying the number of partitions

The performance decrease observed with a single log when the write proxy gets saturated is the main motivation driving the parallelism model of SLoG. By evenly distributing the writes across partitions, we increase the number of proxies, hence increasing the total write capacity of the log.

Since neither Zlog nor Týr offers such a parallelism model to support high logging velocity, we implement for each a SLoG backend. For Zlog, each partition is mapped to a Zlog log. For Týr, we map a partition to a blob. Since thee systems are built for concurrent I/O, the SLoG client library accesses the storage systems directly, bypassing the write proxy.

In Figure 5 we plot the results. We how increasing the number of partitions impacts the performance of the storage. We use the same setup as in the preceding section with 8,192 clients, varying the number of partitions between 1 and 64. The results show a clear increase in the throughput per client. As expected, doubling the number of partitions roughly doubles the achieved throughput on the log by evenly splitting the load across multiple files, reducing hot spots at the tail of the log and consequently increasing throughput.

We also measure the achieved read performance during this experiment. We do not note any substantial performance degradation when increasing the number of partitions, with only a 5,3% performance drop between 1 and 8,192 partitions.

E. Scaling up to 100,000 append cores

We demonstrate the scalability of the system by varying the number of client cores from 12,500 to 100,000. Based on our previous experiments, 16 client cores per proxy yield the best overall performance. We hence vary the number of partitions proportionally between 780 and 6,250. We use both Týr and Lustre as backends, using the same configuration as in the preceding section. Týr is deployed on 20,000 storage cores.

Overall, the results depicted in Figure 6 confirm the nearhorizontal scalability that can be achieved with SLoG by increasing the number of concurrent partitions, using either file-or object-based storage as the backend. We measure a throughput of up to 174 million appends per second on Lustre and up to 241 million appends per second with Týr.

V. RELATED WORK SLoG seeks to enable HPC applications to leverage distributed shared logs over the storage primitives provided by the platform. This tackles the challenges posed by large-scale computational steering. At smaller scale, however, researchers have used comparable approaches on HPC platforms.

In [START_REF] Agelastos | The lightweight distributed metric service: a scalable infrastructure for continuous monitoring of large scale computing systems and applications[END_REF], Agelastos et al. propose a lightweight distributed service able to collect live metrics on the Blue Waters system. However, this tool is specifically built for system monitoring. Knüpfer et al. propose a comparable contribution in [START_REF] Knüpfer | Score-p: A joint performance measurement run-time infrastructure for periscope, scalasca, tau, and vampir[END_REF], describing a large-scale monitoring infrastructure capable of collecting performance metrics from a large number of compute nodes. In contrast with these two contributions, SLoG tackles a more general problem with significantly higher throughput requirements.

The challenges of parallel I/O to a single file on HPC systems have been tackled by Frings et al. in [START_REF] Frings | Scalable massively parallel I/O to task-local files[END_REF]. Indeed, their design of SIONlib enables a large number of processes to concurrently collaborate on a single file. While that work is well suited for a number of use cases, its goal is not to provide a total ordering across the events that is crucial for distributed logging. Actually, while not the target of this paper, SIONlib could be seen as complementary to SLoG by enabling all partitions to share a single file.

Also relevant to the distributed logging field are the excellent publications describing state-of-the-art streaming data systems, including CORFU [START_REF] Balakrishnan | CORFU: A shared log design for flash clusters[END_REF], Kafka [START_REF] Kreps | Kafka: A distributed messaging system for log processing[END_REF], and BookKeeper [START_REF] Hunt | ZooKeeper: Wait-free coordination for internet-scale systems[END_REF]. SLoG builds on the best practices learned from these systems to cope with the hard specificities of HPC storage systems. SLoG does not intend to replace such systems but to complement them by providing an alternative for platforms on which they are not available.

VI. CONCLUSION AND PERSPECTIVES

Distributed shared logs provide a relevant storage model for a number of applications (e.g., computational steering, in situ visualization). However, they are currently not available on HPC platforms, typically offering parallel file systems for persistent storage and only limited possibilities for installing specialized storage services such as distributed shared logs.

In this paper we introduce the SLoG middleware, proposing a distributed shared log abstraction over a parallel file system. To provide high write throughput, we circumvent the low write concurrency of parallel file systems using write proxying, while at the same time partitioning the log to ensure high horizontal scalability. We show on up to 100,000 cores of the Theta supercomputer that SLoG delivers its promise of high append performance with multiple persistence backends, which include object stores in addition to parallel file systems.

We foresee that such an interface could serve in the context of HPC and big data convergence, providing users with a high-level abstraction of distributed logging that can adapt to a variety of backend systems depending on the primitives available on each platform. SLoG could therefore prove to be a strong building block for distributed logging, fueling application portability between HPC and big data platforms.

Fig. 2 :

 2 Fig. 2: SLoG high-level architecture.

Fig. 3 :

 3 Fig. 3: SLoG reader grouping strategies.

Fig. 4 :

 4 Fig.4: SLoG single-partition append performance, using up to 8,192 append client cores, with 1 KB appends.

Fig. 5 :

 5 Fig.5: SLoG multi-partition append performance, using 8,192 append client cores, with 1 KB appends.

Fig. 6 :

 6 Fig.6: SLoG horizontal scalability up to 100,000 append client cores, using up to 6,250 SLoG partitions, with 1 KB appends.

SLoG stands for shared log gateway

ACKNOWLEDGMENT

This work is part of the "BigStorage: Storage-based Convergence between HPC and Cloud to handle Big Data" project (H2020-MSCA-ITN-2014-642963), funded by the European Commission within the Marie Skłodowska-Curie Actions framework. This work was supported by the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research (Contract DE-AC02-06CH11357). It is also supported by the ANR OverFlow project (ANR-15-CE25-0003). We thank Quincey Koziol (LBNL), Noah Watkins (UC Santa Cruz), Robert Jacobs (ANL), and Jay Lofstead (SNL) for their expert assistance. The experiments presented were carried out by using resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility (Contract DE-AC02-06CH11357).


Rémi Curti
email: remi.curti@ensam.eu

Philippe Lorong

Guillaume Pot

Stéphane Girardon

How to model orthotropic materials by the Discrete Element Method (DEM): random sphere packing or regular cubic arrangement?

Keywords: Discrete element method, orthotropic behavior, elements packing, elasticity

come

How to model orthotropic materials by the Discrete Element Method (DEM): random sphere packing or regular cubic arrangement?

Introduction

Modeling of continuous media using a particle approach such as DEM (Discrete Element Method) is recent. It started in the early 2000's with aggregates modeling [START_REF] Jauffrès | Simulation of the elastic properties of porous ceramics with realistic microstructure[END_REF][START_REF] Jefferson | The elastic response of a cohesive aggregate-a discrete element model with coupled particle interaction[END_REF][START_REF] Potyondy | A bonded-particle model for rock[END_REF]. The transition from continuum towards particle models requires considering the interactions between particles not as simple contacts but as "cohesive bonds" having a more complex behavior. Two main aspects of these bonds have to be taken into account: their "force model" and their fracture conditions (if the modeled media is brittle). The force model is used to compute interaction forces between two elements/particles based on their relative positions.

It can also include damping, plasticity or other physical behavior. The fracture condition can be either based on bond deformation, interaction forces or both. For each particle, resulting forces are calculated according to the boundary conditions of the simulated system and to other interacting elements, called its "neighbors". Then Newton's Law of Dynamic is solved with an explicit scheme leading to element's acceleration. New position of the element for the next increment is approximated by integrating its acceleration. In the present study, the integration is done thank to a velocity Verlet algorithm [START_REF] Swope | A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters[END_REF]. This approach derived from molecular dynamics [START_REF] Levitt | Computer simulation of protein folding[END_REF], was initially developed to model granular media [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. It exhibits several advantages toward FEM in order to simulate cutting or fracture of materials. It allows indeed modeling these phenomena with high robustness and as meshless methods as a discrete method, does not require mesh update nor precrack location.

The DEM has already been applied for modeling continuous materials presenting an isotropic and homogeneous behavior [START_REF] André | Modélisation par éléments discrets des phases d'ébauche et de doucissage de la silice[END_REF][START_REF] Hentz | Discrete element modelling of concrete submitted to dynamic loading at high strain rates[END_REF]. It can deal efficiently with heterogeneity of the media and inclusions when required, both in elasticity and fracture mechanics [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF][START_REF] Leclerc | On the suitability of a Discrete Element Method to simulate cracks initiation and propagation in heterogeneous media[END_REF]. It even allows to investigate the behavior of original structure composed both of granular and continuous media [START_REF] Iliev | Discrete element modeling of free-standing wire-reinforced jammed granular columns[END_REF]. Nonetheless, it is also of high interest to deal with anisotropic materials such as fiber-reinforced materials, or natural composites such as wood. DEM is suited for the modeling of dynamic and complex loadings of these materials, for instance crash [START_REF] Espinosa | Coupling continuous damage and debris fragmentation for energy absorption prediction by cfrp structures during crushing[END_REF] or machining. Artificial or natural composites are highly orthotropic materials. Their elastic modulus in fiber direction is much higher than their moduli in transversal directions, with a ratio between these moduli that can reach one hundred.

However, the numerical modeling of these materials with particle approaches, such as the DEM, is more difficult compared to isotropic media. This comes from the difficulty to define bond behavior in order to reproduce the sought-after anisotropy. The modeling of orthotropic elastic behavior with the DEM constitutes the main objective of this paper.

In the DEM, as previously said, the constitutive laws must be implemented in the bonds between discrete elements (DE) rather than in a stiffness matrix like in FEM. In the case of orthotropic behavior, several methods are already experimented in the literature for various purposes. One of the most encountered is the use of regular arrangement, where elements are intentionally lined up along the privileged directions of the material forming a regular Cartesian grid [START_REF] Eberhard | Simulation of cutting processes using mesh-free Lagrangian particle methods[END_REF][START_REF] Iliescu | A discrete element method for the simulation of CFRP cutting[END_REF][START_REF] Pfeiffer | Simulation of green wood milling with discrete element method[END_REF] or user designed regular pattern [START_REF] Roux | Influence of muscle-tendon complex geometrical parameters on modeling passive stretch behavior with the Discrete Element Method[END_REF]. By this way, the orthotropic behavior of the material is rebuilt thank to the structure of the domain. A different strategy suggested by [START_REF] Duan | Discrete element method modeling of inherently anisotropic rocks under uniaxial compression loading[END_REF] consists in using compact domains, where elements are randomly placed and packed by a compacting script [START_REF] Lochmann | Statistical analysis of random sphere packings with variable radius distribution[END_REF]. Then, bonds are generated between elements in contact except the bonds within a certain angle range around a given direction. Thereby, in this direction, the lack of bonds reduces the Young's modulus and generates anisotropy in the domain. An intermediate way of proceeding is to model fibers and matrix by using either regular or randomly packed elements [START_REF] Le | Modeling damages and cracks growth in composite with a 3D discrete element method[END_REF][START_REF] Maheo | A promising way to model cracks in composite using Discrete Element Method[END_REF]. This last approach requires elements quite smaller than the fiber radius. As the amount of elements composing a model impacts drastically the computation time, only small domains can be studied with such an approach. In the literature, to the best knowledge of the authors, the approaches set up to obtain anisotropy are not confronted.

In this paper, the elements are spheres, and the bonds behave as Euler-Bernoulli beams [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF]. Two different arrangements are compared in order to reach an orthotropic elastic behavior: regular cubic arrangement for its popularity in the literature, and a new method, based on packed arrangement, where beams Young's Moduli depend on bonds orientation. In comparison with the methods of the literature this method avoids to simply delete bonds and lower their quantity.

The ability of these two arrangements to match the expected material elastic orthotropic properties is compared along their three main directions of orthotropy, and also under off axis loading. The Hankinson's equation [START_REF] Hankinson | Investigation of crushing strength of spruce at varying angles of grain[END_REF] is used to define the direction dependent reference value of Young's modulus. Simulations were run under quasi-static loading with GranOO workbench (1.0) [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF] with the only addition of a developed plugin to affect the bonds elastic properties. The orthotropic dynamic behavior as well as fracture and anelastic behaviors are not discussed in this study.

Material and method

DEM Models

As described in the introduction, DEM requires the use of bonds to bear the material behavior. The bonds model chosen is Euler Bernoulli cylindrical beam [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF]. The beams elastic behavior depends on four parameters: Young modulus (𝐸 𝑏 ), radius ratio (𝑟b), length ( 𝑙 𝑏 ), and Poisson ratio (𝜈 𝑏 ). The radius ratio equals to the beam cross-section area divided by the mean radius of the two interacting elements (spheres). Euler-Bernoulli beams have three modes of deflection: compression/tension, flexion and torsion. The Poisson ratio allows computing the shear modulus in torsion. The radius ratio brings scalability of the model regarding the discrete element dimensions. If the size of elements is modified, the section and the length of beams will change accordingly resulting in an equivalent behavior of the built domain. In the study conducted in this article, the radius ratio is maintained equals to unity to focus only on the beam Young's moduli impact on the global domain properties. In depth properties of those bonds are explained by [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF]. In addition, numerical damping (𝜂) is added in the beams behavior in order to avoid endless oscillations around equilibrium state during simulations. The use of Euler-Bernoulli beams allows a fine tuning of the material behavior with a small number of parameter compared to commonly used spring/damper models, at the expense of the understanding of each parameter impacts due to their interdependency. Both the stiffness in compression and the stiffness in bending of the bond depend on the same parameters: its Young's modulus and its cross-section.

Regular arrangement

In cubic regular arrangements, also called simply "regular" arrangements in the following, spherical elements are regularly spaced along edges of a cubic frame (see Fig. 1). The directions of edges correspond to the material principal directions (𝑖, ⃗ 𝑗 , 𝑘 ⃗ ). Elements have all the same shape and properties. Their two parameters are their radius and their density (ρ). The density depends on the volume fraction of the model (ratio between the volume occupied by the elements and the whole represented volume). This parameter only impacts the dynamic behavior of the domain and not the static response, so its impact is not investigated in this article. The beams binding elements together are only lined up with the three main directions of the arrangement. They present three different sets of parameters depending on their direction: (𝐸 𝑏𝑖 , 𝑟b 𝑖 , 𝜈 𝑏𝑖 , 𝑙 𝑏 ) along 𝑖, ⃗ (𝐸 𝑏𝑗 , 𝑟b 𝑗 , 𝜈 𝑏𝑗 , 𝑙 𝑏 ) along 𝑗, ⃗ ⃗ (𝐸 𝑏𝑘 , 𝑟b 𝑘 , 𝜈 𝑏𝑘 , 𝑙 𝑏 ) along 𝑘 ⃗⃗⃗ . Only the initial length 𝑙 𝑏 remains mandatorily identical from one beam to another, whatever its direction is, since the initial length is the distance between two regularly spaced interacting elements centers.

In regular arrangements, beams properties can be analytically calibrated to obtain desired structure in-axis properties, for instance Young's modulus Ei, Ej and Ek, in the material principal directions:

𝐸 𝑖 • 𝑆 𝑖 = 𝐸 𝑏𝑖 • ∑𝑆 𝑏𝑖 , 𝐸 𝑗 • 𝑆 𝑗 = 𝐸 𝑏𝑗 • ∑𝑆 𝑏𝑗 , 𝐸 𝑘 • 𝑆 𝑘 = 𝐸 𝑏𝑘 • ∑𝑆 𝑏𝑘 (1)
where 𝑆 𝑖 , 𝑆 𝑗 , and 𝑆 𝑘 are the cross-section of the domain with respectively 𝑖, ⃗ 𝑗 , 𝑘 ⃗ as a normal, and ∑𝑆 𝑏𝑖 , ∑𝑆 𝑏𝑗 , ∑𝑆 𝑏𝑘 , the area occupied by beams section in this domain cross-section.

As previously said the radius is equal to 1 and thus, as all elements have the same radius, all the bonds have the same radius and section 𝑆 𝑒 = 𝑆 𝑏𝑖 = 𝑆 𝑏𝑗 = 𝑆 𝑏𝑘 .

Then equation ( 1) becomes, by isolating Eb and writing respectively 𝑛 𝑖 , 𝑛 𝑗 , 𝑛 𝑘 the amount of elements in the planes (𝑗, ⃗ ⃗ 𝑘 ⃗ ), (𝑖, ⃗ 𝑘 ⃗ ), and (𝑖, ⃗ 𝑗 ),

𝐸 𝑏𝑖 = 𝐸 𝑖 • 𝑆 𝑖 𝑛 𝑗 × 𝑛 𝑘 × 𝑆 𝑒 , 𝐸 𝑏𝑗 = 𝐸 𝑗 • 𝑆 𝑗 𝑛 𝑖 × 𝑛 𝑘 × 𝑆 𝑒 , 𝐸 𝑏𝑘 = 𝐸 𝑘 • 𝑆 𝑘 𝑛 𝑖 × 𝑛 𝑗 × 𝑆 𝑒 . (2)

Packed arrangement

An example of random sphere packing arrangement, which will be named packed arrangement for readability, is given on Fig. 2. An important difference between regular and random sphere packing arrangement is the absence of bonds predefined orientation which is, for the most of them, not close to the material principal directions (𝑖, ⃗ 𝑗 , 𝑘 ⃗ ). It is thus not as simple as for the regular arrangement to define the Young modulus of each beam. The proposed method is to deduce beams' Young's modulus from an interpolation between Ebi, Ebj and Ebk. This interpolation depends on beam's spatial orientation, which is defined by the spherical coordinates 𝜃 and φ, respectively the polar angle and the complimentary angle of the azimuthal angle (see Fig. 3). In DEM, the orientation of the bond linking two elements is defined by its branch vector: vector starting on one element center O1 and pointing out the second element center O2. Note that 𝑝 is the projection of the branch vector 𝑂 1 𝑂 2 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ onto the plane (O1, 𝑗, ⃗ ⃗ 𝑘 ⃗ ).

Four laws to set Young's modulus 𝐸 12 of the bonds of packed domains between two given discrete elements "1" and "2" are investigated in this article:

On-Off:

𝐸 12 = 𝐸 𝑏𝑖 𝑜𝑟 𝐸 𝑏𝑗 𝑜𝑟 𝐸 𝑏𝑘 (3)
if 𝑖, ⃗ 𝑗 or 𝑘 ⃗ is the closest direction.

Hankinson's equation:

𝐸 12 = 𝐸 𝑏𝑖 •𝐸 𝑝 ⃗⃗ 𝐸 𝑏𝑖 •sin 2 𝜃+𝐸 𝑝 ⃗⃗ •cos 2 𝜃 (4)
with

𝐸 𝑝 = 𝐸 𝑏𝑗 •𝐸 𝑏𝑘 𝐸 𝑏𝑗 •sin 2 𝜑+𝐸 𝑏𝑘 •cos 2 𝜑 (5)
the modulus in 𝑝 direction.

Linear:

𝐸 12 = 𝐸 𝑏𝑖 + (𝐸 𝑏𝑗 -𝐸 𝑏𝑖 ) • 𝜃 𝜋 2 + (𝐸 𝑏𝑘 -𝐸 𝑏𝑗 ) • 𝜃•𝜑 ( 𝜋 2 ) 2 . ( 6
)
Spherical:

𝐸 12 = cos(𝜑) • 𝐸 𝑏𝑖 + 𝑐𝑜𝑠(𝜃) • sin(𝜑) • 𝐸 𝑏𝑗 + sin(𝜃) • 𝑠𝑖𝑛(𝜑) • 𝐸 𝑏𝑘 . ( 7
)
Fig. 4 shows the aspect of these four laws for 𝜃 = 0.

Hankinson's equation [START_REF] Hankinson | Investigation of crushing strength of spruce at varying angles of grain[END_REF] is a very common approximation of off-axes properties in composites or wood and already used in FEM or analytical model to predict

Young's moduli according to fiber angle [START_REF] Cramer | Predicting lumber tensile stiffness and strength with local grain angle measurements and failure analysis[END_REF]. Spherical, On-Off and linear interpolations are not physics-based laws but their very simple implementation justifies their trial. Except linear interpolation, all chosen laws offer the advantage to commute: Eb depends on the values of Ebi, Ebj and Ebk but not on the order they are taken for its calculation.

Simulation process

To obtain the mechanical properties of the structure, an unidirectional compression test has been modeled (see Fig. 5 and Fig. 6). The loading direction is denoted 𝑋 ⃗⃗⃗ and can be different from material's principal directions. Four parameters are recorded during the simulation: reaction force F on the plates according to the loading direction, strain 𝜀 𝑥 , of the domain in the 𝑋 ⃗⃗⃗ direction, kinetic energy and potential elastic energy of the system. Forces and strains are then post-treated in order to calculate the global Young modulus Ex in the 𝑋 ⃗⃗⃗ direction from equation [START_REF] André | Modélisation par éléments discrets des phases d'ébauche et de doucissage de la silice[END_REF].

𝐸 𝑥 = 𝐹 𝑆 𝑥 •𝜀 𝑥 (8)
𝑆 𝑥 is the area of the domain section which is measured prior to the simulation. 𝜀 𝑥 computation is explained in section 2.5.

Data is recorded every increment. The optimal time step of one increment is automatically computed by an inbuilt Granoo algorithm according to the beams stiffness. As computations are done in dynamics, it is also verified that when modulus is calculated, kinetic energy of the domain is negligible toward its elastic deformation energy (more than four decades of difference).

In-axis loading

For in-axis loading the simulations are conducted on discretized 10-mm cubes whose frame (𝑋 ⃗⃗⃗ , 𝑌 ⃗ , 𝑍 ) and material principal directions frame (𝑖, ⃗ 𝑗 , 𝑘 ⃗ ) correspond. They are successively loaded along 𝑋 , 𝑌 ⃗ then 𝑍 by a mobile plate compressing the domain against a fixed one, without friction, until they reach 5% deformation in the loading direction. The setup is presented in Fig. 5.

Off-axis loading

Eight supplementary regular domains are designed to present a α-rotated frame (𝑖, ⃗ 𝑗 , 𝑘 ⃗ ) from 10° to 80° by 10°-steps around 𝑍 to highlight mechanical behavior when load is off-axis. The loading direction is 𝑋 ⃗⃗⃗ . For example, Fig. 6 displays a 20°-rotated regular cubic domain under off-axis loading. On the contrary to regular arrangements, the packed arrangement remains the same. However, bond properties are recalculated toward the new material frame orientation for every interpolation law given in 2.1. The packed model is also loaded in the height new cases.

Strain measurement

In order to compute 𝜀 𝑥 it is necessary, in the DEM model, to calculate an average of the relative displacements between elements. For this average, elements interacting with the compression plates are not taken into account. This choice is justified by the specimen's geometry. Indeed, rotated domains faces are "stair" shaped (see Fig. 6 above), thus, only some elements are in contact with the two plates and strains are not homogeneous when too close to the plates.

Under the hypothesis of a homogeneous strain field in the remaining portion of the domains, the displacement 𝑈 𝑥 𝑒 of one discrete element 𝑒 in the direction 𝑋 with respect to the fixed compression plate depends linearly on the deformation 𝜀 𝑥 and can be written as:

𝑈 𝑥 𝑒 = 𝑋 0 𝑒 • 𝜀 𝑥 + 𝐶𝑠𝑡, (9)
where 𝑋 0 𝑒 is the initial position of one element along 𝑋 ⃗⃗⃗ and 𝐶𝑠𝑡 is a constant whose value depends only of the position of the frame origin.

The strain 𝜀 𝑥 , which is minimizing the ordinary squared difference [START_REF] Kenney | Mathematics of Statistics[END_REF] between the approximation 𝑈 𝑥 𝑒 and the measured displacement 𝑈 𝑥 of one element is calculated as stated in [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF]:

𝜀 𝑥 = 𝑐𝑜𝑣𝑎𝑟𝑖𝑎𝑛𝑐𝑒(𝑋,𝑈) 𝑣𝑎𝑟𝑖𝑎𝑛𝑐𝑒(𝑋) = ∑ (𝑋 0 -𝑋 0 𝑚𝑒𝑎𝑛 )•(𝑈 𝑥 -𝑈 𝑥 𝑚𝑒𝑎𝑛 ) 𝑁𝑏 𝑒 ∑ (𝑋 0 -𝑋 0 𝑚𝑒𝑎𝑛 )² 𝑁𝑏 𝑒 , (10)
with, 𝑋 0 the initial position of an element, 𝑋 0 𝑚𝑒𝑎𝑛 the average position of the whole element set, 𝑈 𝑥 𝑚𝑒𝑎𝑛 the mean displacement of the elements in the set and 𝑁𝑏 𝑒 the number of element constituting it.

Model size

Regular domains are at maximum 9261 element-large (see Table 1). This number of elements guaranties at minimum 21 (√9261

3

) elements per direction. This is the minimum amount required to ensure the global properties to remain constant with changes in the number of elements [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF]. To build rotated domains, a larger regular cube is generated, then rotated from the desired angle and every element out of a given 10 mm-sided box is removed. The number of elements fitting the box depends on the angle of rotation of the domain. For this reason, rotated domains minimum element amount is 9401 (in the 40° rotated domain). Packed domains are generated with the standard GranOO compaction method. As this method uses random-based algorithm, two domains consecutively built on the same criterions with the same inputs parameters are different. Thus 10 domains were generated to evaluate the variability of the results. Since the number of elements differs between packed domains, the mean element amount and its standard deviation are calculated and presented in Table 2 below. To obtain a good compactness, the radius of each element composing the packed domain also varies within a 10% range. The distribution is uniform and the variation is centered on the mean radius which is equal to the radius of the elements of the cubic regular arrangement.

Elements and bonds properties are gathered in Table 3.

Results

Anisotropy degrees

The results gathered using equation ( 7) and the simulations described in the previous chapter are displayed in the Table 4. The focus is drawn on the two anisotropy degrees between the three main axes defined as the ratio between the Young modulus observed in two main axes: Ei/Ej and Ei/Ek. These two quantities, named the degrees of anisotropy of the domain, are the main targets to reach when modeling an orthotropic media. If the ratios are respected but all modulus are lower or higher than expected, they can be post-corrected or the calibration process must be reworked. Both solutions do not present technical issues with current algorithms. Fig. 7 displays the ratio Ei/Ej and Ei/Ek ( mean value in the case of packed domains). Notice that the y-axis is in log scale for readability. 2.17 As anticipated, for regular domains, the exact ratios are found in the three main directions of orthotropy, that is Ei/Ej = 20 and Ei/Ek = 100. At the contrary, for packed models the ratios between output degrees of anisotropy of the structure and input degrees of anisotropy in the bond Young's moduli are lower and change depending on the interpolation law chosen. The on-off law is guarantying the highest degree of anisotropy but the output ratio is still lower than the input ratio.

For the linear, spherical and Hankinson laws the results are far from the expected ratios since the degrees of anisotropy recorded are less than 15% as large as the input in every case. It is also noticeable that the relation between the ratios Ei/Ej and Ei/Ek depends on the interpolation method. Thus, the impact of the input parameters on the output degree of anisotropy is investigated in the following.

The same simulations as described previously were run again with 5 different degrees of anisotropy. Simulation parameters used and corresponding results are shown in Table 5. Ebj and Ebk are maintained constant, respectively equal to 5 MPa and 1 MPa while Ebi varies between 5 MPa and 125 MPa in order to increase the input degrees of anisotropy of the model. As a consequence, both input ratios Ebi/Ebj and Ebi/Ebk are affected. The first ratio varies between 1 and 25 while the second varies from 5 to 125. Fig. 8a displays the evolution of the degree of anisotropy of the modelled domain Ei/Ej according to the input ratio of stiffness in its beams Ebi/Ebj and Fig. 8b displays the evolution of the degree of anisotropy of the modelled domain Ei/Ek according to the input ratio of stiffness in its beams Ebi/Ebk. Fig. 8c presents the evolution of the third degree of anisotropy Ej/Ek according to Ebi. For every law, as the input degree of anisotropy is higher, the output degree of anisotropy of the domain increases, but as it grows slower than the input.

Off-axis behavior

Domain's effective Young's modulus is measured every 10°-step and the results are plotted on Fig. 9 from Ej to Ei. Each law displays a different evolution. In the following, ΔE is defined as the amplitude between E(0°) and E(90°): ΔE = E(0°) -E(90°). The effective Young's modulus of an α-rotated domain E(α) decreases very fast in the case of regular domains. E(α) has already decreased of 93% of ΔE when α reaches 40°. Afterwards, the effective Young modulus slowly decreases until E(90°). At the contrary, using either on-off, Hankinson's, linear and spherical laws, the switch between E(0°) and E(90°) is smoother.

Results obtain with the 4 different laws are compared with results obtain with the Hankinson's equation as reference. The use of Hankinson's equation was proven efficient to represent mechanical properties evolution along fiber angle in strongly anisotropic material, including Young modulus [START_REF] Radcliffe | A theoretical evaluation of Hankinson's formula for modulus of elasticity of wood at an angle to the grain[END_REF][START_REF] Yoshihara | Prediction of the off-axis stress-strain relation of wood under compression loading[END_REF].

Considering measured E(0°) and E(90°) are corrects, theoretical 𝐸 𝑡ℎ (10°) to 𝐸 𝑡ℎ (80°) are calculated with Hankinson's equation from E(0°) and E(90°):

E th (α) = E(0°)• E(90°) E(0°)•sin 2 α+ E(90°)•cos 2 α . (11)
Theoretical and mean measured effective Young's Modulus are then compared.

The indicator used to evaluate methods efficiency toward continuum mechanics is the Normalized Root Mean Squared Error [START_REF] Iliev | Discrete element modeling of free-standing wire-reinforced jammed granular columns[END_REF]:

𝑁𝑅𝑀𝑆𝐸 = √ ∑ (𝐸(α)-𝐸 𝑡ℎ (α)) 2 80°1 0°𝑛
•(max(𝐸(α))-min(𝐸(α)))

, [START_REF] Iliev | Discrete element modeling of free-standing wire-reinforced jammed granular columns[END_REF] where n is the number of summed quantities (8 for 10 °-steps). The NRMSE allows to compute a mean error between the measured and the estimated off-axis Young's moduli normalized by the range of the measured Young's moduli. The lower the NRMSE, the more accurate the chosen law. The resulting NRMSE for the 5 studied cases are displayed in Table 6. Hankinson's law comes second. Linear and spherical law's NRMSE is higher and regular arrangement NRMSE is the highest and equals 0.579.

To allow more detailed analysis, the relative error (13) evolution toward every 10°-steps rotation angle is also calculated:

𝑅𝑒𝑙𝑎𝑡𝑖𝑣𝑒 𝑒𝑟𝑟𝑜𝑟 (𝛼) = 𝐸(𝛼)-𝐸 𝑡ℎ (𝛼)
𝐸 𝑡ℎ (𝛼) .

(

) 13
The results according to the law chosen are displayed Fig. 10.

The relative error between effective Young's modulus of the regular domain and

Hankinson's equation increases with specimen angle (α) until it reaches its maximum between 0° and 20° rotation (the measurement resolution is not small enough to give more precision in the estimation). Once this maximum is reached, the error seems to decrease linearly until a 90° rotation. As hinted by NRMSE, relative errors are lower in the case of packed domains, whatever the rotation angle and the law chosen is. However, unlike NRMSE suggested, On-Off interpolation errors are higher than Hankinson's interpolation relative errors, unless for 30°. For 20° and 40°, On-Off error is as low as Hankinson's error, but the low 30° "On-off" errors compensate 10° and larger angle than 50° error.

Indeed, for larger angles, the Young's Modulus decreased enough for the squared error to be also very small once normalized. Linear and Spherical laws errors increases with α until it reaches a maximum of between 50° (approximately 5.33%) and 60° (approximately 5.28%). Then it decreases when α tends to 90°.

Discussion

Degrees of anisotropy

When focusing on two main directions of the considered orthotropic material, degrees of anisotropy of the global domain obtained using regular cubic domains equal to the input anisotropy ratio of the beam Young's Moduli. However, for packed domains, the global domain degree of anisotropy measured is always lower than the input. Moreover the output degree of anisotropy tends to reach an asymptote. Increasing significantly the input ratio does not affect significantly the domain anisotropy anymore. The value of the asymptotes varies depending on the Young's Modulus computation law chosen. This value is also affected by the stiffness of the bond in the third main direction (see Fig. 8c). Because of the previous phenomenon, it is impossible to obtain a high degree of anisotropy if the domain is packed with the Young's moduli calculation strategy tested.

Having very large differences between beams elastic properties to obtain a given global structure degree of anisotropy leads to two consequences. On one hand it increases the computational time as the time constant between two computational increments is directly related to the highest beam stiffness of the domain. Thus, to simulate a phenomenon lasting a given time, a simulation using packed domain instead of regular domain will last longer. On the other hand it can also lead to unphysical local effects leading to easy buckling of the structure because beams with extremely different Young's Modulus interact with each other.

Off-axis effective Young's moduli

Cubic regular arrangement present the exact degrees of anisotropy in the main directions of orthotropy but the effective modulus of the structure when loaded off the main directions is inaccurate compared to packed domains. Indeed, the relative error varies between 6% and 24% depending on the angle formed by the load and the main direction of orthotropy in the case of regular domains. The effective Young's modulus is over-estimated by this law. For a small angle, the relative error between the effective modulus and the theoretical modulus calculated can be ten times smaller for packed domains. In those last types of domains, the relative error, whatever the chosen law is, between effective modulus and its theoretical value is always lower than 5.4%. The most consistent law is Hankinson's equation, while on-off law displays very good results as well. Linear and spherical laws present slightly larger relative errors and tend to under-estimate the effective Young's modulus, especially around 50° to 70° rotation toward one main direction. Thus, as soon as the model is used for complex geometry, load state, displacement or strain fields, this last strategy of model arrangement becomes way more relevant than regular arrangement-based models (in the case of low degree of anisotropy).

Conclusions

To model orthotropic media at mesoscopic or macroscopic scale, with DEM, by the way of Euler-Bernoulli beams, the two proposed approaches led to two different behaviors. Regular cubic arrangement displays perfect stiffness in its main directions of orthotropy, but its response is not as good when loaded offaxis. It suffered up to 24% relative error with Hankinson's equation. In the case of packed domains, four laws to establish beams Young's modulus were evaluated.

Unlike cubic regular arrangements, it was not possible to reach a significantly high degree of orthotropy with any of these four laws. However, in every case their off-axis behavior was very good leading only to 7% relative error in the worst tested scenario. On-off law, even if very simple, seemed to present a good compromise between reachable anisotropy which was the highest of the four laws and off axis behavior which was the second best fitting with Hankinson's equation. Thus, the applications of the two approaches are distinct. Regular domains are mandatory when modeling materials with anisotropy degree is high (starting approximately from 10); but, in the case of low anisotropy, packed domains are preferable due to their advantage in off-axis behavior. For very small degrees of anisotropy (approximately up to 2) Hankinson's law would be the privileged choice. For higher anisotropy (between 2 and 8), On-off law becomes the optimized alternative. Linear and Spherical does not bring competitive benefits towards the two previous laws.

This last approach might be improved in order to increase their reachable orthotropy degree. For instance, the use of Timoshenko beams instead of Euler-Bernoulli model, adding shear deflection to the bonds would be an interesting possibility to investigate. Exploring new laws to determine beams Young's Modulus might also be a path to follow in order to both reach larger degrees of anisotropy and better off-axis behavior. Every tested law but "On-off" was monotonic, to go further, experimenting nonmonotonic laws may perhaps be a solution to rise up the asymptotic anisotropy degree reachable. In the case of either Euler-Bernoulli or Timoshenko beams, adapting the beams section could be one more solution to increase the domains anisotropy as it would change the ratio between tensile and bending potential energy of the system. In this last case, extra attention should be focused on the interaction between beam Young's moduli and sections since both those variables impact the structure global stiffness.

FiguresFig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .

 123456789 Figures

Fig. 10 .

 10 Fig. 10. Relative Error between the effective Young's modulus ratio measure and estimated for tested laws and each rotation of the domain

Table 1 Regular

 1

	domain global characteristics
	Parameter	Value
	Amount of elements	9261
	Amount of bonds	26420
	Volume fraction	0.523 (π/6)

Table 2 Packed

 2

	domains global characteristics
	Parameter	Value
	Amount of elements (mean)	11,670.7
	Amount of elements (standard deviation) 63.75

Table 3

 3 Simulation micro parameters: Elements and bonds properties. In the case of packed domains, the element and beam radius indicated are their mean value.

	Parameter	Value
	Elements (type)	spheres
	Radius (10 -3 • 𝑚)	0.2381
	Density (𝑘𝑔 • 𝑚 -3 )	1000
	Bonds (type)	Euler-Bernoulli beams
	Ebi (MPa)	100
	Ebj (MPa)	5
	Ebk (MPa)	1
	Radius (10 -3 • 𝑚)	0.2381
	Damping factor	0.2
	Poisson ratio	0.3
	Length (10 -3 • 𝑚)	0.4762

Table 4

 4 Degrees of anisotropy measured according to the method of orthotropic generation, with Ebi = 100 MPa, Ebj = 5 MPa and Ebk = 1 MPa

	Law	Ei/Ej	Ei/Ek
	Regular arrangement	20	100
	On-off	3.56	7.03
	Hank	1.55	3.27
	Linear	2.18	2.39
	Spherical	2.02

Table 5

 5 Input degree of anisotropy applied to the beams and output degree of anisotropy of the domain measured according to the method of orthotropic generation

	Ebi/Ebj	1	5	10	15	20	25
	Ei/Ej (Regular arrangement) 1.00	5.00	10.00 15.00 20.00 25.00
	Ei/Ej (On-Off)	1.00	2.07	2.75	3.21	3.56	3.85
	Ei/Ej (Hankinson)	1.00	1.38	1.49	1.53	1.55	1.57
	Ei/Ej (Linear)	1.01	1.69	1.97	2.10	2.18	2.24
	Ei/Ej (Spherical)	1.00	1.60	1.84	1.95	2.02	2.06
	Ebi/Ebk	5	25	50	75	100	125
	Ei/Ek (Regular arrangement) 5	25	50	75	100	125
	Ei/Ek (On-Off)	2.22	4.46	5.69	6.49	7.03	7.51
	Ei/Ek (Hankinson)	1.96	2.82	3.10	3.21	3.27	3.31
	Ei/Ek (Linear)	1.47	2.05	2.25	2.34	2.39	2.42
	Ei/Ek (Spherical)	1.38	1.89	2.05	2.13	2.17	2.20

Table 6

 6 NRMSE values calculated for the different method of orthotropic generation

	Law	NRMSE calculated
	Regular arrangement	0.579
	On-off	0.164
	Hank	0.288
	Linear	0.326
	Spherical	0.390
	The lowest NRMSE is shown by the On-off interpolation method: 0.164.

Acknowledgement

The authors would like to acknowledge GranOO staff members, especially Jean-Luc Charles and Damien André for the training they provided.

Compliance with Ethical Standards:

Conflict of Interest: The authors declare that they have no conflict of interest.

Funding: These works were conducted thanks to the support of the region Bourgogne Franche-Comté and were funded by the ANR-10-EQPX-16 XYLOFOREST-XYLOMAT and the French MESRI.