Amir Charif

Gabriel Busnot

Rania Mameesh

Tanguy Sassolas

Nicolas Ventroux

Fast Virtual Prototyping for Embedded Computing Systems Design and Exploration

Virtual Prototyping has been widely adopted as a costeffective solution for early hardware and software co-validation. However, as systems grow in complexity and scale, both the time required to get to a correct virtual prototype, and the time required to run real software on it can quickly become unmanageable. This paper introduces a feature-rich integrated virtual prototyping solution, designed to meet industrial needs not only in terms of performance, but also in terms of ease, rapidity and automation of modelling and exploration. It introduces novel methods to leverage the QEMU dynamic binary translator and the abstraction levels offered by SystemC/TLM 2.0 to provide the best possible trade-offs between accuracy and performance at all steps of the design. The solution also ships with a dynamic platform composition infrastructure that makes it possible to model and explore a myriad of architectures using a compact high-level description. Results obtained simulating a RISC-V SMP architecture running the PARSEC benchmark suite reveal that simulation speed can range from 30 MIPS in accurate simulation mode to 220 MIPS in fast functional validation mode.

INTRODUCTION

To keep up with a fast-evolving and highly competitive embedded system industry, designers are compelled to deliver complete working solutions under tight delay and budget constraints. To make this possible, hardware architectures, as well as the full software stacks that drive them should be validated and optimized as early as possible in the design process. In this context, Virtual Prototyping has been widely adopted as a cost-effective solution for early hardware and software co-validation.

The rapid adoption of virtual prototyping solutions was greatly facilitated by the emergence of the SystemC/TLM 2.0 standard [START_REF]Open SystemC Initiative (OSCI)[END_REF], which, in addition to offering interoperability and reusability of SystemC models, provides several abstraction levels to cope with varying needs in accuracy and speed. More recently, in response to an increasing demand for simulation speed, the use Dynamic Binary Translation (DBT) for CPU modelling has gained in relevance [START_REF]Open Virtual Platforms[END_REF], [START_REF] Rodman | ARM fast models-virtual platforms for embedded software development[END_REF], [START_REF] Delbergue | QBox: an industrial solution for virtual platform simulation using QEMU and SystemC TLM-2.0[END_REF], and has effectively set a new standard for simulation performance in early prototypes.

However, as systems grow in complexity and scale, it is now more vital than ever that virtual prototyping solutions be able to wisely exploit these technologies to offer proper balance between simulation representativeness and execution speed throughout the design process. Moreover, at this level of complexity, the time required to model and explore new architectures can quickly become unmanageable, especially at the earliest stages, where it is often necessary to make heavy alterations before reaching a stable prototype.

To cope with these new difficulties, virtual prototyping solutions need to meet new requirements not only in terms of performance, but also in terms of ease, rapidity and automation of system modelling and design space exploration. This paper introduces VPSim, a feature-rich integrated virtual prototyping solution that addresses the aforementioned challenges based on three major contributions:

• A new way of integrating the QEMU emulator, making its rich CPU and peripheral model portfolio available as a collection of SystemC modules. Our method can leverage the technologies used in QEMU, such as Dynamic Binary Translation (DBT) and paravirtualization to bring outstanding simulation speeds into the more deterministic and standardized SystemC domain. • An Accuracy Control framework, that can be used to define regions of interest in both the software and the simulated hardware dynamically. A system designer can therefore enable accurate simulation only on parts of the design and portions of code that are of interest, guaranteeing the best possible trade-off between performance and accuracy given specific needs. • A generic and powerful infrastructure for dynamic platform composition and design space exploration (DSE). It allows a single compiled executable to be used to model an infinity of architectures. Combined with a Python scripting front-end, it makes it possible to simulate, explore and optimize highly complex systems in a single compact script.

The remainder of this paper is organized as follows: In Section 2, we review some related tools and methods. Section 3 provides a complete high-level view of the features and capabilities of the VPSim tool. In Section 4, details on the underlying implementation challenges are presented. The performance of our tool is evaluated in Section 5, before concluding in Section 6.

RELATED WORKS

The first generation of Virtual Prorotyping (VP) solutions used functional Instruction Set Simulators (ISSs) with more or less internal low-level details such as pipeline stages. Among such solutions, GEM5 [START_REF] Binkert | The gem5 simulator[END_REF][START_REF] Menard | System Simulation with gem5 and SystemC[END_REF] is a discrete-event simulator able to dynamically switch between different abstraction levels. Detailed CPU models with full pipelining description can be simulated at 0.1 Million Instructions Per Second (MIPS), whereas instruction-accurate CPU models can reach 1 MIPS. Other MPSoC modeling environments, such as SESAM [START_REF] Ventroux | SESAM: An MPSoC Simulation Environment for Dynamic Application Processing[END_REF] or Unisim [START_REF] August | UNISIM: An Open Simulation Environment and Library for Complex Architecture Design and Collaborative Development[END_REF] used instruction-based ISS generation librairies to support the modeling of various CPUs, reaching approximately 10 MIPS. However, this accuracy and ISA flexibility come at the cost of limited simulation speed, which hampers their capacity to address complex systems embedding several CPUs and running full-fledged OSes.

To address this complexity, a recent trend is the use of Dynamic Binary Translation (DBT). DBT consists in dynamically translating instructions of the modeled ISA (guest instructions) to host ones (usually x86), whenever needed during guest code execution, yielding very high execution speeds.

Today, all high-performance simulation environments use DBT. Industrial solutions are led by solutions like Virtualizer [START_REF] Virtualizer | Synopsys[END_REF], Vista [START_REF]Mentor, A Siemens Business[END_REF], VPS [START_REF]Cadence[END_REF], VLAB [START_REF] Vlab | ASTC[END_REF], FastModels [START_REF] Rodman | ARM fast models-virtual platforms for embedded software development[END_REF] or, for instance, Simics [START_REF] Simics | Wind River[END_REF]. However, these come at a significant cost and do not offer the degree of customization that is required when designing new architectures. In addition, they do not provide fast design space exploration capabilities.

Amongst open-source virtual platforms, OVP [START_REF]Open Virtual Platforms[END_REF] uses processor models simulated through a closed-source DBT engine named OVPSim reaching hundreds of MIPS. An OVP model can be wrapped for inclusion in a SystemC model using TLM 2.0 interfaces. However, OVPSim cannot provide performance evaluation as the models are timed relying on instruction count and neglecting memory access timings. In addition, OVP claims to support parallel multiprocessor simulation but this is at the cost of deterministic execution loss.

QEMU [START_REF] Bellard | QEMU, a fast and portable dynamic translator[END_REF] is an open source DBT based emulator supporting many CPU models but does not provide performance estimation. Several works have tried to embed QEMU within a SystemC simulation environment to provide both determinism and timing evaluation [START_REF] Gligor | Using Binary Translation in Event Driven Simulation for Fast and Flexible MPSoC Simulation[END_REF]. In [START_REF] Gligor | Using Binary Translation in Event Driven Simulation for Fast and Flexible MPSoC Simulation[END_REF], the authors wrap QEMU processors within SystemC threads and investigate several QEMU instrumentation options to take into account instruction count and data access latency. Depending on whether synchronization shall occur on every data access or on a periodic basis, the simulation speed varies from 3 to 60 MIPS. However, the cumbersome annotation is to be performed for every ISA and possible accuracy settings are limited to a few predefined configurations. GreenSocs [START_REF] Delbergue | QBox: an industrial solution for virtual platform simulation using QEMU and SystemC TLM-2.0[END_REF][START_REF] Monton | Mixed sw/systemc soc emulation framework[END_REF][START_REF] Montón | Mixed simulation kernels for high performance virtual platforms[END_REF] propose a QEMU-based framework that exploits QEMU MMIO callback mechanisms to access external Sys-temC peripherals. This allows for very fast simulation for applications with few IO communications, as the execution mostly takes place in the context of QEMU. Unfortunately, this solution runs QEMU and SystemC in separate kernel threads, requiring frequent synchronization and precluding determinism.

The solution we propose integrates QEMU by executing its CPU and peripheral models in the context of Sys-temC threads, thereby preserving the predictability of Single-Threaded SystemC simulation. The accuracy of memory accesses, including instruction fetches, can be configured at a very fine granularity, and may change dynamically during simulation.

A USER-LEVEL VIEW OF VPSIM

VPSim is a key asset in charge of virtual prototyping and DSE within SESAM, an integrated EDA framework for complex electronic systems ranging from Cyber-Physical Systems to Microservers. SESAM provides a holistic environment to address HW/SW co-design, exploration and validation through Virtual Prototyping, HW prototyping and emulation while taking into account power, temperature and reliability factors. This section provides a high-level view of the major features and capabilities of VPSim as perceived by the end-user.

VPSim is a tool that was designed specifically to accelerate software/hardware co-validation at the earliest design stages of all kinds of computer architectures. It can be used to easily compose, simulate and explore new hardware architectures, but also to run, profile and debug full software stacks on the simulated platforms.

Figure 1 gives a global overview of how VPSim can be used. The user specifies a platform using a platform composition front-end, which forwards a high-level platform description to a central VPSim object named the Platform Builder. The Platform Builder uses this description, and a library of registered components to construct and run a SystemC simulation. Components can be either internal hardware models, or Proxy Components for interfacing with external subsystems. VPSim loads software binaries through the ElfLoader and BlobLoader special components. Once the SystemC simulation is launched, the user can use standard debugging tools or VPSim's built-in debug and profiling utilities, to evaluate and optimize their software. VPSim collects fine-grained per-components statistics and pushes them back to the highlevel front-end, which may use them to present performance profiles to the user or to perform design space exploration.

VPSim is also capable of operating as merely one component of a wider Cyber-Physical System (CPS) simulation, in compliance with the Functional Mockup Interface (FMI) co-simulation standard [START_REF] Blochwitz | Functional mockup interface 2.0: The standard for tool independent exchange of simulation models[END_REF]. In this mode, VPSim is packaged as a Functional Mockup Unit (FMU), that can be loaded by an FMI simulation master. This makes it possible to evaluate the cyber part of a cyber-physical system (modelled in VPSim) along with its surrounding physical environment.

Component library

VPSim includes a large library of CPU models, buses and peripherals to chose from. Common controllers such as UART, I2C and SPI, PCI-Express from various vendors (Xilinx, Renesas, Cadence, etc.) are modelled in VPSim. Network, block devices and GPUs are also made available through Paravirtualization [START_REF] Russell | virtio: towards a de-facto standard for virtual I/O devices[END_REF]. All the components that are instantiated in VPSim are fully standard-compliant SystemC/TLM 2.0 modules. CPU and Virtio models are made available to the SystemC world through a new approach for integrating QEMU into SystemC, which will be described in Section 4.

To further broaden its range of supported models, VPSim provides proxy components. These components are able to load and connect SystemC/TLM 2.0 initiator and target subsystems from any third-party EDA vendor to the system described in VPSim. The external subsystem is separately compiled in a shared library that implements simple glue functions. Through this interface, VPSim can seamlessly integrate CPU models such as the ARM Fast Models [START_REF] Rodman | ARM fast models-virtual platforms for embedded software development[END_REF], OVP (Open Virtual Platforms) CPUs [START_REF]Open Virtual Platforms[END_REF] or QBox [START_REF] Delbergue | QBox: an industrial solution for virtual platform simulation using QEMU and SystemC TLM-2.0[END_REF].

Composition and Exploration

To enable truly rapid prototyping of complex architectures, it is not enough for simulations to be fast. The time required to fully describe an architecture, configure it, modify it is also of critical importance to designers and software developpers.

Most related tools employ a common scheme, wherein the platform's top is written in SystemC, either manually or generated using a Graphical User Interface, and then compiled. Later on, the configuration of the simulation is usually described in a higher level language such as Lua [START_REF] Delbergue | QBox: an industrial solution for virtual platform simulation using QEMU and SystemC TLM-2.0[END_REF].

VPSim adopts a different philosophy: Compilation shall take place only once, and the whole simulation, including the platform to be simulated, the applications to be run, and the configuration, shall be described using a dynamic front-end. The benefits of this approach are manyfold. For instance, changes in the simulated architecture do not require recompilation. This is a true game changer at early design phases, where it is very common to make adjustments in the architecture itself. Also, the ability to describe the system hierarchy along with the components' configuration in the same environment removes the need to create a custom configuration file format for every new platform.

VPSim also provides the front-end with fine-grained simulation statistics, making it possible to perform both platform specification and Design Space Exploration within the same environment. VPSim communicates with the front-end using the XML markup language, as described in Section 4. [START_REF] August | UNISIM: An Open Simulation Environment and Library for Complex Architecture Design and Collaborative Development[END_REF].

By default, VPSim ships with a Python-based composition and exploration front-end, detailed in Section 4.4. Note that unlike other tools that use high-level description languages, e.g. GEM5 [START_REF] Binkert | The gem5 simulator[END_REF], VPSim is not tighly linked to its Python interface, nor does it know about its existence. By using XML, which is a standardized data exchange format, many new interfaces can be developped for VPSim according to specific needs. In addition, it makes it much easier to couple VPSim with existing XML-capable tools.

Debug facilities

A decisive factor in the usefulness of any virtual prototyping tool is its ability to inspect and help understand the behaviour of the system under evaluation. VPSim offers debug and control capabilities out of the box. These are described in what follows.

Standard tool support.

VPSim makes it possible to debug individual CPU cores in the system using GDB. Each CPU core in the system has a gdb enable attribute which, when set to True, enables a GDB session to control and debug the code it executes.

3.3.2

The VPSim Monitor. The system designer interacts with VPSim through a command line interface named the "VPSim Monitor". At any point during the simulation, the system can be frozen using a key stroke, bringing up a command prompt. Through the Monitor, the user can inspect/alter the memory and register contents, but also reconfigure the entire simulation. For instance, it is possible to change the debug level of any component in the system to enable more or less debug information for the rest of the simulation. It is also possible to create watches on specific address ranges to monitor memory accesses, and possibly manually override the read/written values.

Dynamic Checkpointing. When inside the VPSim

Monitor, it is possible to create one or several named checkpoints. Then, at anytime during the rest of the simulation, it is possible to roll back to the named checkpoint to investigate the cause of bugs. When a checkpoint is created, the entire user process is copied using the fork() system call. The parent process then sleeps and waits for a wakup signal. Rolling back to the checkpoint simply consists in terminating the current process after waking up the correct parent. The list of checkpoints and their associated process IDs is transported with each rollback to make sure that all checkpoints are accessible from any process instance. This non-intrusive approach has several advantages. Unlike [START_REF] Montón | Checkpointing for virtual platforms and SystemC-TLM[END_REF], our method does not require implementing a checkpointing method for each individual component, nor does it introduce changes to the SystemC kernel. Since it operates on the entire process, it works even with closed-source SystemC modules. Compared to other methods that use snapshots of the entire process, such as [START_REF] Kraemer | A checkpoint/restore framework for SystemCbased virtual platforms[END_REF], our method can operate fully in-memory, without having to save images to disk. This makes rolling back to checkpoints created within one same simulation much faster than [START_REF] Kraemer | A checkpoint/restore framework for SystemCbased virtual platforms[END_REF], where several seconds are necessary to store and load the snapshots.

Profiling and Performance Evaluation

Unless otherwise specified in global simulation parameters, VPSim attempts to run the target software as fast as possible. That is, most accesses to main memory will be emulated in a fast, untimed fashion based on the Direct Memory Interface (DMI) of TLM components. This mode of operations greatly eases and accelerates the functional validation and debugging phases. However, when optimizing the target software, it would also be helpful to get detailed information about its execution, such as cache usage and network latency.

One specificity of VPSim is its ability to dynamically define regions of interest, in both the software and the hardware. These regions of interest are simulated more precisely, i.e. bus accesses, caches, and other transactions are fully modelled. This helps provide a performance profile for a specific portion of the executed code.

The user can describe the beginning and end of the software region of interest dynamically (during execution), or statically (in global parameters). The hardware components that need to be precisely simulated can also be designated by the user. One common usage is when profiling a Linux user-space application. In such cases, the Linux boot phase is of little interest and can be quickly skipped through in untimed mode. VPSim then switches to a more accurate simulation mode once the user application is entered.

At the end of the software region of interest, statistics from all the simulated components are registered and displayed to the user (e.g. cache miss rate, network latency, executed instructions, etc.). This fine-grained control over the simulation's accuracy allows the user to make the best trade-offs between simulation speed and accuracy given their specific needs, which may evolve from one design phase to another.

KEY TECHNICAL CONTRIBUTIONS

In this section, we describe the key technical contributions that we made to enable all the features presented so far.

QEMU in VPSim

QEMU is one of the leading Dynamic Binary Translation (DBT) -based emulation solutions to date [START_REF] Bellard | QEMU, a fast and portable dynamic translator[END_REF]. In addition to unmatched execution speeds, QEMU supports a wide range of virtualized and paravirtualized hardware models, including most of the existing CPU architectures, and many peripherals.

Methods for leveraging QEMU in a SystemC simulation environment have been extensively explored [START_REF] Monton | Mixed sw/systemc soc emulation framework[END_REF], [START_REF] Montón | Mixed simulation kernels for high performance virtual platforms[END_REF], [START_REF] Delbergue | QBox: an industrial solution for virtual platform simulation using QEMU and SystemC TLM-2.0[END_REF], [START_REF] Gligor | Using Binary Translation in Event Driven Simulation for Fast and Flexible MPSoC Simulation[END_REF], [START_REF] Tlmu | Edgar E. Iglesias[END_REF], [START_REF] Chiang | A QEMU and SystemC-based cycle-accurate ISS for performance estimation on SoC development[END_REF]. Some of these methods view QEMU as merely an Instruction Set Simulator (ISS), while others make it possible, in addition to CPU clusters, to model complete subsystems within QEMU, and some peripherals in SystemC [START_REF] Delbergue | QBox: an industrial solution for virtual platform simulation using QEMU and SystemC TLM-2.0[END_REF]. In the latter approach, accesses to the peripherals modelled within QEMU are not visible to the outter SystemC world.

VPSim adopts a different approach. In VPSim, all CPUs and peripherals must implement a common consistent interface and be executed in the SystemC context. This is key to guaranteeing a certain level of predictability, inspectability and compatibility with VPSim's debugging and profiling tools. Therefore, in VPSim, models that are backed by QEMU can be instantiated, configured, and controlled in the same way as native components. For instance, the VirtioNet and VirtioBlk components, which are backed by virtio-net and virtio-blk devices in QEMU, can be instantiated, mapped to any address, connected to any interrupt line in VPSim. Accesses to these peripherals, as well as the interrupts they generate, are all visible and debuggable through the VPSim Monitor, for instance.

In our approach, no IO accesses are served within QEMU. Instead, all accesses are visible to the SystemC world and completed by SystemC models. Accesses to RAM (fetches and data) can be completed in QEMU for efficiency, but only if the target Memory module has a Direct Memory Interface (DMI) as per the SystemC/TLM 2.0 standard. QEMU's CPU execution thread and IO thread are both executed in the context of SC THREADs and SC METHODs, making simulations more predictable and controllable. In the absence of external input (Network, Keyboard, etc.), simulations in VPSim are deterministic and guaranteed to be repeatable.

While transparent to the user, internally, VPSim maintains a number of QEMU instances, as shown in Figure 2. Each QEMU instance may be used to instantiate CPUs of the same model, plus any number of peripherals. These models are then associated to proxy VPSim components, which are exposed to the user like any other components.

Dynamic Accuracy Control

In the SystemC/TLM 2.0 standard, some memory-mapped modules may have a DMI (Direct Memory Interface). That is, it is possible to get a pointer to their internal memory space, and access it directly, without simulating the entire TLM transaction. QEMU in VPSim makes good use of this feature by declaring all DMI regions as RAM regions [2], thereby enabling ultra-fast inline accesses to main memory.

However, to enable the accurate simulation mode presented in Section 3.4, it is necessary to perform the full TLM transactions. Components in VPSim possess a DMI OK flag, which is active by default. When set to false, the component will issue a DMI invalidation request for its address range. Optionally, only a subset of the entire a accessible range can be invalidated. This propagates the invalidation request upstream and will force the upstream component to initiate full TLM transactions for the invalidated range. The upstream component then sets its own DMI OK flag to false, which will provoke the invalidation of its own accessible address space by its upstream initiator, and so on until the root of the device tree is reached.

VPSim can take TLM transactions only as far as necessary to get the desired level of accuracy. For instance, if the user disables the DMI OK flag only on a Cache component, VPSim only makes sure that TLM transactions get to the designated Cache component, which then completes memory accesses using DMI pointers to the final targets. VPSim's Component interface automatically manages these pointers internally and is able to use them to complete a transaction when eligible.

In VPSim, components take the port from which these invalidation requests are received into account. Therefore, a Component may force TLM accesses only on specific ports, while allowing fast DMI simulation on the rest of the device tree. This is extremely useful when only part of the memory hierarchy is of interest. For instance, forcing TLM accesses on the instruction cache component will only force blocking TLM transactions from QEMU when fetching instructions, while data accesses remain unaffected.

Platform Builder

At the core of the dynamic platform composition presented in Section 3.2, is a central VPSim object named the Platform Builder. The Platform Builder interacts with a composition front-end, such as the Python frond-end presented in Section 4.4, and has three roles:

First, it communicates to the front-end the structure of the platform description document that it expects. Currently, VPSim supports platform descriptions in XML format. Therefore, the Platform Builder automatically generates, from the list of self-registered Components, an XML Schema Document (XSD) that describes the structure of the expected XML document. This document also includes information about all the available components in VPSim and their attributes. The front-end should use this Schema to present the components to the user. A GUI front-end might display them as boxes that can be connected to each other, for instance. The Python frond-end presented in Section 4.4 uses the Schema to dynamically generate Python classes corresponding to each component.

Second, upon receiving an XML document from the frontend, it elaborates and configures the specified platform. Platform elaboration takes place in three phases:

• Make: During this phase, all the components of the platform are instantiated and initialized. The Platform Builder checks whether all the required attributes were specified, sets the specified attributes, and assigns defaults when applicable. • Connect: All connections between components that were specified in the platform description are realized.

The availability and compatibility of the connected ports are checked during this phase. • Finalize: The Platform Builder sets some Platformwide parameters as specified in the input XML document, and, more importantly, invokes a finalize callback that Components may implement to perform finalization actions after all the platform has been elaborated.

Finally, at the end of a simulation, the Platform Builder collects per-component statistics and forwards them back to the front-end in an XML document. The front-end may use these statistics to build performance profiles, or for architectural exploration.

The Python front-end

The default and preferred front-end in VPSim uses the Python language. With this front-end, VPSim strives to provide a dynamic platform composition environment that is compact, intuitive, but also expressive enough to allow detailed platform configuration and description.

To illustrate with an example, Listing 1 shows how a 4-core ARMv8 Platform capable of booting the Linux operating system is described and simulated in VPSim. A simulated system is represented by the System class. Each component is represented by a corresponding Python class, e.g. Memory, Bus, etc. The attributes of each component can be either specified in the constructor for compactness (see Listing 1, Line 9), or as regular class members (see Listing 1, Line 3). This makes it possible to set or modify some attributes programmatically. Of course, each of these components has many more attributes that can be configured. When omitted, VPSim sets them to sensible defaults. Some relationships between components can also be inferred automatically. For instance, since the GIC400 (Generic Interrupt Controller) is the only interrupt controller in the system, it is automatically connected to the interrupt lines of the UART and the architected timers of the CPUs. Calling the simulate method interprets the entire system and launches a simulation. At the end of the simulation, per-component statistics are returned in a Python dictionary object. This model allows for as much expressiveness and modelling freedom as one would get writing a Top in SystemC. In addition, we get all the dynamicity, compactness and power of the Python language. The Python front-end of VPSim was also designed with rapid design space exploration in mind. Because several System instances can live within the same Python script, and because simulation results are made available as Python objects, complete design space exploration can be performed in one same script. To accelerate DSE, the Python frond-end makes it possible to perform the simulations concurrently. To do so, it is enough to call the simulate method with a wait=False argument, to run simulations asynchronously. This fully integrated approach saves the user from having to write ad-hoc automation scripts for each new study.

EXPERIMENTAL SETUP AND RESULTS

In this section, we describe the experimental setup and the results obtained in terms of simulation performance.

Simulated platform and experimental setup

All the experiments have been conducted on a RISC-V-based simulated platform [START_REF] Waterman | The RISC-V Instruction Set Manual[END_REF]. It is an SMP platform whose number of cores ranges from 1 to 32, running Linux kernel 4.6.2 as guest. Each simulated CPU has 32KB instruction and data caches, and all cores share an L2 cache of 2 MB and 1 GB of RAM. This platform was simulated under VPSim on a i7-4770 host machine with 16 GB of RAM running Ubuntu 18.04 LTS. We have checked that the CPU frequency remains stable at 3.9GHz under any load.

When driven by the Accellera SystemC simulation kernel, VPSim is single-threaded and relies exclusively on singlecore performance. Simulation performance is never limited by RAM, as less than 200 MB are used in any simulation scenario. The simulations were executed one at a time inside Docker containers based on Ubuntu 18.04 LTS. We have verified that Docker has no significant impact on the runtime performance of VPSim.

The simulated platform runs Linux with various benchmarks from the PARSEC [START_REF] Bienia | The PARSEC Benchmark Suite: Characterization and Architectural Implications[END_REF] suite: bodytrack, ferret, swaptions, blackscholes, fluidanimate and dedup. The standard simmedium datasets were used. The execution sequence starts with Linux boot, followed by the benchmark setup, warmup, the benchmark's Region of Interest (ROI), teardown and finally, the platform's shutdown. Simulation time starts when the sc start function is called and ends when the kernel sends the poweroff signal. Platform elaboration time is negligible with respect to the duration of a simulation.

Depending on the experiment, four accuracy modes are used:

(1) DMI: this is the fastest mode where as many transactions as possible are completed using DMI. (2) L1-ROI: DMI is used everywhere except in the ROI for transactions initiated by CPU 0. In that case, TLM transactions reach the L1-caches, which then complete the requests using DMI. (3) ACCURATE-ROI: DMI is used everywhere except in the ROI, where blocking TLM is used. (4) ACCURATE: this is the most accurate mode where all memory accesses are modelled and timed using blocking TLM transactions.

ACCURATE-ROI mode is as accurate as ACCURATE mode during the ROI, which is likely to be what the user is focusing on in a real world application. L1-ROI mode gives accurate information about the L1-caches of a single core during the ROI. In cases where the platform and the benchmarks are symmetric in terms of executed code, the behaviour of one core with respect to caches is likely to be representative of the rest.

The number of MIPS and the real simulation time are collected at the end of each simulation. A total of 3 experiments are conducted:

Simulation results and performance

In Figure 3a, the raw performance of VPSim in fast DMI mode on a single-core simulated platform can be observed. The speed ranges from 121 to 219 MIPS and the real simulation time from 31 to 193 seconds. The Linux boot and shutdown processes take approximately 14 seconds. The difference in observed performance mainly comes from varying benchmark profiles and the way this impacts QEMU. Indeed, QEMU performance can be hampered by conditional jumps or memory accesses requiring the use of the software MMU. However, these factors are hard to predict and there is no clear correlation between, for instance, the frequency of load/store instructions of a benchmark and the speed at which it can be simulated. In Figure 3b, the effect of platform complexity on MIPS and real simulation time is reported. While there is a slight drop in MIPS between the single and dual-core platforms, the speed remains in the noise margin up to 32 cores. This was expected as the overall simulation pattern remains the same: a single core doing the work outside of the ROI and all cores sharing the work inside the ROI. However, the real simulation time increases almost linearly with the number of cores. This can be explained by the sequential simulation of all cores, which introduces extra overhead in single-threaded phases (e.g. boot, setup, etc.). By contrast, the time required to simulate the ROI remains more-or-less constant, since a fixed amount of computation is shared between all cores. As a result, ROI simulation time becomes negligible compared to the other phases, hence the linear increase in simulation time.

Figures 3c and3d show the speedups obtained when variable accuracy is applied. The reference speed for both MIPS and real simulation time is the ACCURATE mode. It can be observed that ACCURATE-ROI introduces a MIPS speedup comprised between 1.3x and 1.5x with no loss in simulation accuracy in the ROI. The only slight difference that occurs resides in the cold caches at the beginning of the ROI. However, if the ROI is big enough, this effect is negligible. Also, simulating the caches while not recording the statistics during the first steps of the ROI could mitigate this inaccuracy. The MIPS speedup obtained with L1-ROI ranges from 3.9x to 6.6x. In the context of cache performance evaluation on a multithreaded application, this mode may be sufficient. Indeed, if the application makes a comparable usage of all cores and exhibits negligible data sharing, then there is no need to simulate all caches.

The observation on real simulation time is comparable but speedup is larger. It ranges from 2.1x to 9.2x for ACCURATE-ROI and from 6.8x to 17.8x for L1-ROI. One particularly interesting thing to note is the difference between these speedups and the respective MIPS improvements. The reason resides in the timing accuracy variations. When simulating with maximum accuracy, all memory transactions are timed and participate in the update of the internal SystemC time, which passes faster than in DMI mode as a result. Consequently, all actions triggered by timed interruptions occur more often. This concerns several Linux routines such as thread scheduling, IO polling, etc. This variation in the simulated instructions is a side effect of the simulation accuracy changes. It plays in favor of variable accuracy as simulations get even shorter.

CONCLUSIONS

This paper introduced several contributions. First, we presented an alternative approach to integrating the QEMU emulator into SystemC. Our approach makes the CPU and peripheral models of QEMU available as TLM modules. One key benefit is the ability to leverage QEMU's high-performance emulation technologies such as Dynamic Binary Translation and Paravirtualization to boost simulation speed. Results show speeds approaching 220 MIPS when simulating workloads on top of Linux. Then, we presented a method for fine-tuning the simulation to get the best performance given specific accuracy requirements. Our framework is based on a clever use of the TLM 2.0 DMI specification. We experimentally demonstrated that by focusing simulation effort only on regions of interest, speedups of over 17x could be achieved. Finally, we showed how it was possible to abstract away all the complexity SystemC architecture modelling, by introducing a sophisticated dynamic platform elaboration infrastructure. By allowing designers to specify their modelling and exploration needs at a higher level, hours of development and compilation time can be saved. All of these promising building blocks are integrated into a unique virtual prototyping solution named VPSim.

Figure 1 :

 1 Figure 1: An Overview of the VPSim Platform.

Figure 2 :

 2 Figure 2: QEMU integration in VPSim.

1 from vpsim import * 2 3

 2 c l a s s ExampleSystem ( System ) :

8 9 13 rom . i s r e a d o n l y = True 14 s e l f . s y s b u s >> rom 15 16 s 26 27Listing 1 :

 8131416261 f . s y s b u s = Bus ( l a t e n c y =10* ns ) ram = Memory ( b a s e=0x10000000 , s i z e =1*GB) 10 s e l f . s y s b u s >> ram 11 12 rom = Memory ( b a s e=0x00000000 , s i z e =100*KB) e l f . s y s b u s >> GIC400 ( b a s e=0 x f f 0 0 0 0 0 0 ) 17 18 s e l f . s y s b u s >> 19 CadenceUART ( b a s e=0x f e 0 0 0 0 0 0 , i r q =0x70 ) 20 21 f o r i in range ( 4 ) : 22 cpu = Arm64 ( id=i , model= ' c o r t e x -a57 ' ) 23 cpu . r e s e t p c = rom . b a s e 24 cpu ( ' i c a c h e ' ) >> s e l f . s y s b u s 25 cpu ( ' d ca ch e ' ) >> s e l f . s y s b u s BlobLoader ( o f f s e t=rom . base , 28 f i l e= ' u-boot . b i n ' ) 29 BlobLoader ( o f f s e t=ram . base , 30 f i l e= ' example . dtb ' ) 31 BlobLoader ( o f f s e t=ram . b a s e+0x80000 , . addParam ( param="quantum" , v a l u e =1000* ns ) 37 s t a t s = s y s . s i m u l a t e ( ) Simulating a quadcore ARMv8 System running Linux in VPSim

Figure 3 :

 3 Figure 3: (a) MIPS and real simulation time for various benchmarks on a single-core platform running Linux. (b) MIPS and real simulation time for the fluidanimate benchmark with 1 to 32 cores. (c-d) Speedup with various simulation accuracy levels compared to ACCURATE mode on a quad-core platform.


MOTS-CLÉS : Matériaux hétérogènes

Keywords: modélisation morphologique, comportement quasi-fragile Heterogeneous materials, Morphological modelling, quasi-brittle mechanism 31 èmes Rencontres

Multi-scale simulation of quasi-brittle heterogeneous materials : application to concrete like materials E. Roubin 1 1 LMT-Cachan, Ecole Normale Supérieure de Cachan, Université Paris VI, CNRS, UniverSud Paris, 61 avenue du Président Wilson, 94235 Cachan Cedex, FRANCE emmanuel. roubin@ lmt. ens-cachan. fr Prix Jeunes Chercheurs « René Houpert » RÉSUMÉ. À des échelles d'observation fines (inférieures à l'échelle centimétrique) les matériaux à matrice cimentaire ne peuvent plus être considérés comme des milieux homogènes. Ils sont, dans le cadre de cette étude, considérés comme bi-phasiques, les granulats étant inclus dans une matrice cimentaire. Deux problématiques majeures qu'entraine cette thématique de modélisation sont abordées. D'un coté la modélisation morphologique et d'un autre, la modélisation mécanique prenant en compte le caractère hétérogène et quasi-fragile du matériau. Le résultat principal de l'étude montre que la simple modélisation de la morphologie permet, avec des modèles mécaniques simples, de retrouver des comportements macroscopiques complexes. Une application est faite montrant l'impact des incompatibilitées de déplacement entre la pâte de ciment et les granulats du au retrait de séchage sur la dissymétrie du comportement en traction/compression du béton.

ABSTRACT. A 3D meso-scale model for failure of heterogeneous quasi-brittle materials is presented. At such scale, concrete can be represented as an heterogeneous material with two phases, where inclusions (aggregates) are included within the matrix (mortar). Both morphological quasi brittle failure mechanism modelling are expressed here. Finally, impact of drying shrinkage incompatibilities between aggregates and cement paste on asymmetrical response in tension/compression tests is shown.

Introduction

In view of the growing complexity of macroscopic models for heterogeneous materials, such as concrete or more generally cement-based ones, the question of scale observation becomes relevant. On the one hand, it clearly appears that most of the macroscopic behaviours of such materials (creep, shrinkage, cracking. . .) take their origin at smaller scales (mesoscopic, microscopic. . .). Hence, predictive and accurate numerical simulation of long-term behaviour related to concrete structures durability shall account for these scales. On the other hand, usage of very fine imaging procedures, mostly inspired from the medical field, with concrete-like materials has increased dramatically over the last few years. Chief among them is the X-ray tomography [ROU 10] which offers a nondestructive way to get very precise data, dealing with both morphological and mechanical aspects. Numerical simulations shall be able to use this huge volume of data for both parameter identification and validation of numerical models.

The framework presented here is to be seen in this context, and especially within a sequenced multiscale way (as opposed to integrated ones) where the macroscopic behaviour comes from a mesoscopic description of the material. In other words, we consider the macrosocopic mechanical behaviour of such material as a complex system. It is assumed that the final result cannot be predicted directly. A way of dealing with such systems is to set several basic evolution rules and simulate the system in order to obtain the structured final result. In our case, this result is the macroscopic behaviour and the basic rules are defined by the mesoscale description. This point of view can also be seen as virtual testing procedure.

At fine scale, concrete must be represented as a heterogenous materials. Therefore, two kinds of basic rules need to be established in such framework. The first kind is a purely geometrical and topological matter with the aim of explicitly representing, as faithfully as possible, the morphology of material heterogeneities. The random aspect of this representation links our virtual testing framework with Monte Carlo experiments. The second kind of basic rules is purely mechanical and defines the local behaviours of each phases of the material and their interfaces. A quasi-brittle failure mechanism is wanted, adding a strongly non-linear aspect to the system. Both kind of representation of the meso-scale are complex (random heterogeneities and non-linear behaviour) yet one have to keep in mind that those two aspects are completely independant.

The global spirit of the whole framework is to consider that with a representative morphological description, the mechanical behaviour of each phase can be kept as simple as possible and still lead to complex macroscopic behaviour. For example, it is shown in this work that with a simple brittle tension failure mechanism at the mesoscale, macroscopic compression failure can be simulated.

Explicit representation of the meso-scale heterogeneities

In this work, the mesoscopic scale characteristics are explicitly defined. Therefore a special effort had to be made concerning the morphological representation of heterogeneities. Two characteristics shall be considered concerning this modelling. First, the geometry and topology has to fit as close as possible to reality (volume fraction, area fraction, number of aggregates, size repartition...) and secondly, the random aspect of aggregates repartition within the matrix has to be explicitly represented. Two different methods of heterogeneity morphological representation are quickly overviewed. The first one represents aggregates with non-intersecting spheres of various radii randomly placed into space (see figure 1(a)). A special packing algorithm of collective rearrangement is used here in order to increase the maximum volume fraction computable in acceptable time [BEZ 02]. This method is extremly efficient when it comes to specific size repartition of aggregates. Only little work on the radii probabilistic distribution is required. Major drawback is the extremly simple spherical aspect. The idea behind the second method is to yield the phases using levet set of correlated random fields (see figure 1(b)). In opposition to white noises, the use of correlated random fields gives a spatial structure to the resulting morphologies (hereafter called excursion sets). A probabilist link between Gaussian correlated random field characteristics (mean, variance or correlation length) and both geometrical and topological quantities (volume, surface area, Euler characteristic...) has been recently made in [ADL 08], giving the possibility of controlling the excursion sets characteristics and applying it to represent each material phase and even make it evolve through time easily if needed, opening possibilities of hydration modelling problematics for instance. This method gives several advantages such has accurate control on each geometrical characteristics or representation of random shapes, but still makes it hard to yield high fraction volume disconnected morphologies. Yet, using Gaussian related random field (like the χ 2 field) and union of random fields allows to increase the fraction volume to reasonnable values.

Efforts of morphological modelling are here made within a multi-scale linear framework using a FE model with embedded discontinuities. In order to represent these heterogeneities, those morphologies are projected onto the FE mesh, thus defining a set of discontinuities within the strain field interpolation (weak discontinuities). These kinematics enhancements lead to "non-adapted" meshes in the sense of independence between heterogeneities morphology and the underlying FE mesh.

Finite Element failure modelling

The second main aspect of this work is about modelling the quasi-brittle behaviour of concrete-like materials. During the last decades, it has been shown that, in contrast to strain softening in classical continuum material models, allowing displacement discontinuities in the kinematics field leads to a more meaningfull description of the material nonlinear behaviour [SIM 93]. In the presented framework, discontinuities in the displacement field represent what we can call "meso-scale crack". In order to handle both heterogeneities and crack-opening mechanism at meso-scale, weak discontinuities are completed with a set of strong discontinuities within the framework of local enhancement [BEN 10]. The magnitude of the displacement discontinuity in each element gives the measure of the crack opening. Strong discontinuities are activated according to a certain localization criterion and are then governed by a simple relationship between traction vector and crack opening. Such elements are represented figure 2 from blue to red, regarding the magnitude of the crack. At a larger scale, the coalescence of several "local-cracks" shall represent the localization of a macroscopic crack. In other words, the local failure mechanism leads to global macroscopic failure by upscaling structural effect. This vision of modelling the macroscopic behaviour by this natural method is the cornerstone of this work, leading to several original considerations. The main one is that this structural effect has a physical sense only if the mesoscale is physically well described. A simple tension/compression test presented in the next section point out the major effect of morphological description. However, by simulating a shrinkage mechanism, we show that adding an other kind of meso-scale description such has a pre-cracked pattern leads to more realistic macroscopic behaviour.

Natural asymmetric tension/compression behaviour

Results on both tension and compression simulations are presented here. They have to be seen has two completely independant calculations yet using exactly the same mesoscopic description in terms of morphology, material characteristics and mechanical behaviour. Hence, we show that very simple local failure mechanism are enough to retrieve much more complex features at macro-scale. For the two-phases material presented here (aggregates within a cement paste) three local desciptions are applied depending on whether an element is completely in the matrix, completely in an aggregate (in grey on figure 2) or in both (in black onfigure 2). The latter is called "interface element" and handles material discontinuities with special inherent enhancement, simulating a perfect interface. All the material characteristics are given in table 1. It is assumed that aggregates keep an elastic behaviour. On the other hand, after reaching a certain maximal stress in tension, maxtrix and interface element, by the mechanism of crack-opening dissipate a certain amount of energy before total failure. Hence, the tensile strength and the fracture energy are the two parameters governing the non linear behaviour.

Figure 2. Representation of the structural effect on the macroscopic crack for tension and compression test.

A first comparison between the results of tension and compression tests can be made regarding the shape of the macroscopic crack (see figure 2). One can first notice that in both case, macroscopic cracks find their way around elastic aggregates leading to complex geometry and even "natural" branching. However, the global shape differs completly. When in tension, a single crack is well localized following a plan perpendicular to the solicitation, in compression, nearly all the elements are broken, leading to several plans, this time parallel to the solicitation. The differences observed between the macroscopic cracks reflect the asymmetric response of the material. All the upscaled mechanical characteristics of the material at the macroscopic scale are given table 2. Regarding the necessary energy to reach macroscopic failure we can note that both are bigger than the mesoscopic one, as it has been shown in experiments. Futhermore, those energies indicate that the behaviour in tension is way more brittle than in traction. Regarding the strength and the energy needed before total failure of the material, the difference between the two tests reflects the general more brittle behaviour of concrete-like material in tension than in compression. These simple tests show that, using the same mesoscopic description (mechanically and morphologically speaking), the upscaled macroscopic behaviour naturally represents a complex behaviour such as asymmetry between tension and compression response. The structural effect due to the perturbation of heterogeneities is the only reason of the asymmetric response observed. Yet several drawbacks can be noted regarding strength values. Both are rather low compared to reality, espacially the compression one. It can be improved with a more accurate morphological representation. Ideed, results presented herein are base on a monodisperse spherical method, yielding an aggregate phase of 20% fraction volume. Calculations (not presented here) have shown that increasing this volume by adding smaller balls (by following actual grrain size repartitions), or using random field excursion sets increase significantly the compression strength yet not much the tension one (increasing additionaly the strength ration).

Young modulus

Effect of homogeneous shrinkage on the macroscopic response

We saw that the outcome of the structural effect due to morphological representation brings strong asymmetric behaviour regarding tension and compression. In this section, we add to the simulation a homogeneous shrinkage mechanism by imposing a given strain, ε rd , to the matrix before the tension or compression solicitation leading to a "pre-cracked" state of the material. Such homogeneous sollicitation can represent, for exemple, the effect of drying shrinkage incompatibilities between cement paste and aggregates [LAG 11]. It can be seen as a infinitly slow drying, leading to a diffuse crack pattern, in opposition to oriented one (that is the usual effect of differential drying [DESA 08]). The macroscopic response with (resp. without) shrinkage are shown figure 3 Tension ε rd = 0 24 800 2 17 ε rd = 5 10 -5 16 000 1.2 9

Compression ε rd = 0 24 800 9 750 ε rd = 5 10 -5 24 300 8.5 720

Table 3. Upscaled material characteristic at macro-scale with drying

We can note that the pre-cracked state has a more important effect in tension than in compression, leading to a bigger ratio between maximum stress in tension and in compression. The brittle behaviour of tension test is more sensitive to the shrinkage solicitation. It helps the localization of the one macroscopic crack leading to even more brittle behaviour. In opposition, the little effect it has on the compression response represents the higher complexity of the failure mechanism. However, the framework lacks mesoscopic physical effect modelling, leading to questionable results. For example, it as been shown that the increasing capillary tension due to shrinkage, increase the compression strength as well [SEL 09].

Conclusion

Most of problems related to concrete structures durability (such as reinforced concrete corrosion or delayed ettringite formation) are directly linked to mass transfers. Thus they are clearly affected by any mechanical degradation in the sense that any set of meso-or macro-cracks shall increase those transfers by several orders of magnitude. Here the detailed knowledge of the meso-scale cracks pattern (cracks openings and orientations) is a useful tool for such matter. For exemple, dealing with water or gaz permeability issues, a macroscopic crack in a material leads to strong anisotropy due to its orientation (figure 2) and linking it to a macroscopic permeability tensor is not an easy task. The knowledge of the crack provided by the framework presented here gives the basis of a new upscaling permeation scheme [JOU 11].

Figure 1 .

 1 Figure 1. Morphological modelling using non-intersecting spheres packing or random fields.

Figure 3 .

 3 Figure 3. Representation of the drying shrinkage mechanism on the macroscopic behaviour for tension and compression tests.

Table 1 .

 1 Material characteristics of the mesoscale

	Aggregate	80 000	-	-
	Matrix	20 000	3	5
	Interface	-	3	5

Young modulus Tensile strength Fracture energy [M P a] [M P a] [J].[m] -2 (a) Macroscopic crack in compression test (b) Macroscopic crack in tension test

Table 2 .

 2 Upscaled material characteristic at macro-scale

			Strength Fracture energy
		[M P a]	[M P a]	[J].[m] -2
	Tension	24 800	2	17
	Compression	24 800	9	750

èmes Rencontres de l'AUGC, E.N.S. Cachan, 29 au 31 mai 2013


Sylvain Losserand

Gwennou Coupier

Thomas Podgorski

Migration velocity of red blood cells in microchannels

teaching and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

Blood is a dense suspension of essentially RBCs within plasma, which is a Newtonian fluid. RBCs have the particularity to be highly deformable in physiological conditions, which leads to complex flow patterns in confined environment due to several mechanisms.

One of the most famous features, first observed by Poiseuille [1], is the presence of a cell-free plasma layer (CFL) near vessels walls [2][3][4][5][6] inducing a decrease of the apparent viscosity referred to as Fåhraeus-Lindquist effect [7] as well as a decrease of the hematocrit in small vessels compared to large ones [8,9]. Several phenomena can be invoked to explain this cell free layer, one being geometrical: the presence of a capillary wall implies that the centers of RBCs must lie at least one RBC half-thickness away from the wall. This means that, on average, there will be more RBCs near the center of the capillary than very near the wall. The second one is hydrodynamical: at low-Reynolds number the deformability of RBCs allows a symmetry breaking that leads to transverse migration from the wall to the center of the channel [10,11]. This symmetry breaking can also be obtained through the elastic deformation of the glycocalyx brush covering the endothelium [12]. Finally, RBC aggregation can have an influence on the CFL [5]. Transverse migration also plays a role in the so-called margination effect: contrary to RBCs which are concentrated around the centerline, other blood elements like white cells or platelets are marginated against the wall [13][14][15][16][17][18][19], which helps them to accomplish their function (e.g. in immune response, inflammation...). The centering of RBCs also influences the repartition of RBCs at bifurcations: due to the heterogeneity of the RBC concentration in vessels and the existence of a cell free layer, a heterogeneous distribution of the hematocrit takes place in the microcirculation [20][21][22][23][24][25][26][27].

The importance of understanding this migration is not limited to in vivo microcirculation but also has impor-tant possible applications in the design of Lab-on-Chip devices for particle separation. Geislinger et al. [28] presented a device capable to separate a dilute suspension of RBCs and platelets: the less deformable platelets migrate less than RBC and by controlling the position of the two outlets it is possible to separate the two populations. Important works are also done on a device named deterministic lateral displacement (DLD), which is made of a dense network of small pillars. The trajectory of particles inside this network is highly dependent on their properties, resulting in the possibility to separate different entities like RBCs [29], or different types of white cells [30].

In confined situations where the channel size is equal to a few cell radius, the effect of shape changes is strong and non-monotonous time variation of position can be observed. Numerical simulations explored the migration dynamics and the asymptotic position of simplified models of cells like 2D capsules [31] and 2D vesicles [32,33] that can only qualitatively be compared to experiments. To date, although results are available on the shape and the position of RBCs in very confined channels [34,35], there exist no experimental measurement and no quantification of migration dynamics of RBCs in confined Poiseuille flow (but in the pioneering work of Goldsmith [10]). The main objective of this paper is therefore to describe the migration of isolated RBCs for different confinements by using a simple general scaling law which could be used to validate numerical simulations, design future microfluidic devices or apply to some physiological situations where the hematocrit is low due to strongly heterogeneous distribution of cells at bifurcations in the microcirculation [20][21][22][23][24][25][26][27].

Transverse migration in simple shear flow Among biomimetic models, giant vesicles are closed lipid membranes exhibiting in certain conditions, dynam-ical properties that are similar to those RBCs, though they do not have cytoskeleton providing in plane shear elasticity. Their migration dynamics was theoretically described by Olla in [11]. For a vesicle with a fixed ellipsoidal shape, the drift velocity when the distance to the wall is large compared to vesicle size reads

ẏ = U γ R 3 0 y 2 , ( 1
)
where y is the distance of the center of mass to the wall, γ the shear rate of the unperturbed flow, R 0 the typical size of the vesicle and U a dimensionless parameter that depends on vesicle properties: the viscosity contrast λ between the inner and outer fluids and the reduced volume ν = 6 √ πV/S 3/2 where V is the vesicle volume and S its surface area. In the following, R 0 will be taken as the radius of the sphere of the same volume V as the considered particle: R 0 = (3V/4π) 1/3 . For RBCs of typical volume 90 µm 3 [36], this leads to R 0 = 2.8 µm.

This law has been confirmed later by 2D [37] or by 3D [38][39][40] numerical simulations and by theoretical analysis [40]. In the past years, several experiments in microgravity condition were conducted to measure this drift velocity [41,42]. The overall scaling for the drift velocity as well as the value of the prefactor U are in agreement with simulations and theory (even though Olla assumes a fixed shape), for λ = 1 [39,40] or higher values of λ [42]. Note that in Ref. [39] the y -α scaling with α = 2 seems to be valid for y/R 0 2 and that for 1 y/R 0 2, α is smaller than 2, in agreement with the experimental near-wall study of Ref. [43] where α = 1 is suggested.

Finally, similar y -2 scalings for the lift velocity were found for elastic capsules through numerical simulations [44,45].

In physiological conditions, in terms of external fluid (plasma) viscosity and shear rates, isolated RBCs are far from a tank-treading regime [46][47][48][49][50] where they would adopt a fixed orientation relative to the flow. It was recently shown, however, that the scaling proposed for vesicles is still valid and a mean drift velocity (which corresponds to a drift velocity averaged over one body rotation period) characterized by U R 3 0 = 0.36 µm 3 was measured [51]. We are aware of no direct theoretical or numerical validation of this measurement.

Transverse migration in a channel

In a channel, the shear rate of the unperturbed flow is not uniform. This sole symmetry breaking is sufficient to induce migration of lipid vesicles even in the absence of walls [40,52,53]. In a realistic channel, wall lift and shear rate gradient contribute to migration in a complex way, in particular near the centerline where cell shape or dynamic depends a lot on its position close to this zero shear rate centerline. Over a wide range of confinement degree(2 < w/R 0 < 10, where w is the half width of the channel), reduced volumes and viscosity ratios, it was shown experimentally in Ref. [54] that for lipid vesicles the migration velocity may be written as

ẏ = ξ R δ+1 0 γ(y) (y -y w ) δ , ( 2
)
where δ is close to 1 and ξ a dimensionless parameter that depends on the vesicle properties, similar to U for the drift under simple shear rate and y w is the position of the center of mass when the particle is as close as possible to the wall. For the quasispherical vesicles considered in [54], y w ∼ R 0 but in general, it may depend on particle deformability. This scaling was confirmed by 2D numerical simulations [54].

The alternative empirical law that we propose in this paper,

ẏ = ξ R δ+1 0 γ(y) y δ , (3)
is formally simpler and allows for comparison between different situations with no need to take into account the detail of the near-wall interactions. This formulation is the one generally used in the literature [11,17,55] and we shall consider it in the following. Indeed, RBCs are so deformable that y 0 can be very small and discussing its very precise value is not relevant from an experimental point of view, due to the dispersion in cell properties, initial cell orientation and the precision we can reach in its determination.

We show here that scaling law such as Eq. 3 can appropriately describe RBC migration in a confined channel and provide values for δ and ξ that are valid in a wide range of parameters.

EXPERIMENTAL SET-UP AND METHODS

The considered microfluidic channels are straight and of rectangular cross section. The width and the height of the section are respectively denoted 2w and 2h, the useful length of the channel is 2.7 cm. RBCs flow along the Ox direction and migrate laterally along Oy. The walls are located at y=0 and y = 2w in the Oy direction and at z = -h and z = h in the Oz direction (Fig. 1). Before the channel of interest, the flow is established over a long time resulting in a preliminary centering of the RBCs in the Oz direction (z = 0), depending on the width of the channel the preliminary centering in Ox is not necessarily complete. The suspension then reaches a T intersection: the control of the flow rate in the daughter branches permit to move the separation line, thus allowing to direct slightly more than the CFL of the incoming branch in Blood is supplied by the Etablissement Français du Sang (EFS) and is taken from healthy volunteers. Blood samples are first washed 3 times in a solution of phosphate-buffered saline (PBS tablet Sigma ™ : 10 m phosphate buffer, 2.7 m potassium chloride and 137 m with pH 7.4). Then samples are suspended at a low volume fraction (0.1%) so that RBC/RBC interactions become negligible in an isodense solution made of 70mL of Optiprep ™ (iodixanol aqueous solution), 130mL of ultrapure water and a PBS tablet [56]. The buffer has a density of 1.112 ± 0.001 which prevents the sedimentation of the RBCs and a viscosity of η e = 1.9 × 10 -3 Pa s at 20 • C which is a little higher than plasma viscosity (1.54 mPa s at 25 • C [57]).

Furthermore, for one experiment, density-fractionated RBCs were prepared using discontinuous gradients with Optiprep. The gradient was built up in four layers of 2 mL, containing ultra pure water and respectively 40% (1.128 g mL -1 ), 37% (1.118 g mL -1 ), 33% (1.106 g mL -1 ), 31% (1.099 g mL -1 ) of Optiprep, buffered with a PBS tablet. 2 mL of blood sample initially washed using the method explained above are layered on the top of the gradient solution following a centrifugation at 3,000 rpm for 30 minutes at room temperature. After the centrifugation, two subpopulations are extracted and suspended in a density matched solution. The first population corresponds to RBCs having a density between 1.099 and 1.106 and represent the first 5% of the total population, the second population has a density included between 1.118 and 1.128 representing the last 10% of the total population.

Images are acquired with a 32X objective mounted on an inverted microscope ( Olympus ™ IX71) using a CMOS camera (ImagingSource ™ , DMK UX174) at a rate of 30 frames/s. RBCs are not followed directly, measurements are made on windows of size L=350 µm at different positions x; about 200 RBCs are observed by windows (see Fig. 1). A custom Python software using the library OpenCV ™ [58] determines the center of mass of the projected shapes of RBCs. Using the trajectories of the RBCs obtained in each window, velocities v p (y, 0) are computed using a simple tracking program. Depending on the confinement, the probability density function of the position y of the RBCs is calculated on a part or on the entire window. To reduce the influence of the outliers the median trajectory y(x) is chosen to describe the trajectory of the RBCs. The dispersion y d is defined as the mean between the position of the 33 rd and 66 th percentile. All lengths dimension are rescaled by w and are denoted with a ∼.

RESULTS

Longitudinal velocity -Fig. 2 shows the measured velocity v p of RBCs as a function of their position ỹ (grey The maximum velocity is at ỹ = 1 as for the fluid but the experimental plots show that v p converges to 0 before ỹ = 0. A straightforward explanation is the influence of the finite size of a particle compared to the size of the channel. There is no simple solution for this complex flow due to the presence of walls. Furthermore, the deformability and the complex shape of the RBC eliminate the possibility of a simple analytical solution.

A description of v p is necessary in the continuation of this paper to translate the migration velocity (Eq. 3) into a description of the trajectory (Eq. 9). By considering RBCs as solid spheres, an analytical solution was proposed in 2D Poiseuille [59] flow but this approach cannot be easily extended to 3D flow. We approximated v p by the mean velocity of the flow without RBC averaged on a square area of size (2r) 2 and of center of mass y and z = 0.

The mean velocity is thus written as an integral:

v p (y, z) = 1 4r 2 z+r z-r y+r y-r v x (Y, Z)dY dZ, (4)
with v x the velocity of the fluid without RBCs in a rectangular channel [START_REF] Bruus | Theorical microfluidics[END_REF]:

v x (y, z) = v 0 ∞ n=1,3... a n   1 - cosh b n (y -w) cosh (b n w)   cos (b n z), (5)
where

v 0 = v max ∞ n=1,3... a n 1 - 1 cosh (bnw) , (6)
with v max = v x (w, 0) the maximal velocity of the fluid and

a n = (-1) n-1 2 n 3 , b n = nπ 2h . (7)
This yields:

v p (y, z) = v 0 h πr 2 ∞ n=1,3... a n n sin (b n r)   r - sinh (b n r) cosh b n (y -w) b n cosh (b n w)   . ( 8
)
Despite our rough assumption, the fit using Eq. 8 with v max as a free parameter and r = R 0 , gives a correct approximation of the experimental values (see Fig. 2). The model brings to light the effect of the confinement on the difference between v p and v x : an increase of the channel width induces a reduction of the difference between the two velocities. Discrepancies with the model appear near the wall and come from the fact that the model does not take into account the non negligible influence of the RBC on the flow. This influence being stronger near walls, the discrepancies should be stronger near walls as shown by Goldman [START_REF] Goldman | [END_REF]. v max obtained from the fit using Eq. ( 8) is between 0.5 and 3 mm s -1 which is in the range of typical values observed in microcirculation. A capillary number Ca = η e γe R 0 /µ s can be introduced, with γe = v max /w a characteristic shear stress and µ s = 3 µN m -1 for the RBC shear elasticity [62]. Ca is found to be between 1 and 7.

Transversal velocity -Fig. 3 shows an experimental result for w = 10.6 µm. RBCs move away from the wall and migrate towards the center. Eq. 3 describes the migration velocity and from this equation the theoretical trajectory y th (x) obeys the differential equation:

∂y th ∂x = ξ R δ+1 0 γ(y th , 0) y δ v p (y th , 0) . (9)
γ and v p are respectively obtained from Eq. 5 and Eq. 8. Both expressions are proportional to v max , therefore the trajectory does not explicitly depend on v 0 . Consequently, the potential inaccuracies in its determination through the cell velocity will not affect the determination of the lift parameters. The data are fitted using Eq. 9. A first step consists in integrating it numerically and in a second step the residual of a cost function is computed. The problem has three free parameters ξ, δ and the initial value of y denoted y 0 . The minimization of < |(ỹ -ỹth / ỹd | > is based on a global method named differential-evolution [63]. In Fig. 3 different fits of the same dataset are shown, all of them giving reasonable results. The value of y 0 found from minimization is pretty constant and appears to be independent of the two other parameters. In contrast, δ and ξ seem correlated as an increase of one induces the increase of the other one.

By fixing ỹ0 = ỹ(0), the residual can be explored in the (δ, ξ) parameter space, as shown in the insert of Fig. 6) which is obtained from the fitting of experimental data using Eq. ( 8). The dashed line represents the velocity of the fluid vx without RBC at z = 0 using Eq. ( 5) and the full line represents the fitted velocity using Eq. ( 8). FIG. 3. Trajectory of RBCs in a channel with w = 10.6 µm (red). Fitting curves (black) using different couples of parameters (δ, ξ, y0), in 1 ○ all parameters are free, giving at the end (0.14, 4.1 × 10 -3 , 0.38), 2 ○, 3 ○, 4 ○ are fits using a fixed δ giving respectively (0, 3.7 × 10 -3 , 0.38), (1, 7.7 × 10 -3 , 0.37), ( 2 Thus parameter values (δ, ξ) present in this zone give an error inferior to y d . This zone varies from confinement to confinement and no clear law relating the 2 parameters can be extracted. In more simple geometry [51] or with simpler objects [41,42,54], these two parameters were instead shown to be independent of each other. We consider here that the apparent correlation for a given confinement is the result of the specific migration function combined with the noise in the data that prevents a more accurate fit. Aiming at finding a (δ, ξ) couple of values that gives the best fit over a wide range of parameters, we thus take into account the whole set of results in different channels. For each experiment, for a fixed value of y 0 = y(0) and δ ranging from 0 to 2 with a step of 0.05, the value of ξ minimizing the residual is determined. By considering all the channels a mean < ξ > and a normalized mean square error ∆ξ = < (ξ-< ξ >) 2 >/ < ξ > as a function of δ can be computed. From the minimum of the mean square error (see Fig. 4, top left) we extract the most probable value for δ and ξ: δ s = 1.30 ± 0.05 and ξ s = 1.1 × 10 -2 ± 0.2 × 10 -2 .

The RBC trajectories ỹ(x) are shown in Fig. 4 for five different confinements. In our range of parameters, RBCs move from the wall towards the center of the channel in agreement with the results of [34] and the migration velocity decreases as RBCs approach the centerline. The fit obtained from (3) with δ = δ s , ξ = ξ s and y 0 free give good results even at high confinement where the lateral migration is mainly the result of the deformation of the RBC.

We also investigatde the influence of RBC mechanical properties by isolating the most and less dense fractions of RBCs in a density gradient (see Experimental Set up and Methods). Fig. 5 shows that less dense RBCs migrate faster than denser RBCs, the comparison of their respective ξ obtained from the fit of Eq. ( 9) when δ = δ s is fixed show a clear difference: 1.4 × 10 -2 for the less dense and 6.5 × 10 -3 for the denser ones. Denser RBCs are usually considered more aged than less dense RBCs [64]. Different interpretations can be found The red color represents the median trajectory of the less dense RBCs (points) and its fitting curve using Eq. ( 9) with a fixed δs, free ξ and ỹ0 (dashed line). The blue color denotes the trajectory and its respective fit for the denser RBCs. The black curve is the general trajectory found previously using δs and ξs and ỹ0 the mean between ỹ(0) of the two experimental curves. The insert represents the median trajectory of rigidified RBCs in red and the general trajectory in black.

in the literature about the effect aging on RBCs, either through a dehydration of RBCs that leads to an increase of hemoglobin viscosity, or to modification of membrane properties [65,66]. Netherless, there are strong indications that older cells have a higher density and lower deformability [67] (eventually leading to their elimination in the spleen [68]). Our study also indirectly suggests that denser cells are less deformable.

For comparison, the insert in Fig. 5 shows that artificially hardened RBCs using glutaraldehyde (at a concentration > 0.1% so that proteins are completely crosslinked and RBCs become rigid [69]) lose their capability to migrate and behave like rigid particles. The deformability is thus an important parameter in the migration dynamics.

DISCUSSION

An exponent different from δ = 2 may be interpreted in several ways. Except when the vesicle is close to the center, its migration velocity in an unbounded Poiseuille flow is constant (δ = 0). An exponent between 0 and 2 might be seen as a kind of average between the contributions of wall-induced lift and migration due to nonconstant shear rate. Still, the weight of these contributions is not the same depending on the position of the particle [70,71], so it is not clear whether such a scaling with a δ exponent independent from the cell radius to channel size ratio would hold.

In Ref. [17], 3D RBCs were simulated in a 2D Poiseuille flow; the cells' membrane have the same properties as that of red blood cells but the viscosity of the inner fluid is taken to be equal to that of the external fluid. In addition, the capillary number is quite large such that the dynamics of the cells resembles that of vesicles, with tank-treading like motion accompanied by some small oscillations. They find that for a channel width of 12 and 18 cell radius, a scaling law such as Eq. 3 is found with δ in the range 1.2-1.3 on average, with ξ in the range 0.04-0.07. We thus find a similar range for δ, but smaller values for ξ. This may be explained by the value of the viscosity contrast, which is not 1 in reality. Increasing the viscosity contrast generally leads to a decrease of ξ in Poiseuille flow [54] or of U in simple shear flow [11,42]. In Ref. [51], we found a factor 9 for the U values between the physiological case and a viscosity contrast close to 1. The similarity in the range for δ shows that the found values well describe the interplay between wall and shear gradient effects, in that range of confinement. To the contrary, our results are in contradiction with the model proposed in [18]: it is proposed that because of shape change with the shear rate, therefore with the y position, ξ may also depend on y [18]. This may somehow be transcripted into a different exponent δ, as long as the migration velocity is confronted to numerical or experimental data that spans on a limited range and/or are associated with noise. In Ref. [18], δ is set to 2 and the proposed modeling implies that ξ is a decreasing function of y, which implies that the velocity decreases much more strongly with y then in Eq. 3 with ξ constant and δ ∼ 1.3.

The structure of a RBC suspension is the result of two counterbalancing effects, lift and multiple collisions between cells, that lead to transverse diffusion [51]. This controls in particular the width of the cell free layer near the walls.

In Ref. [18], a good agreement with the literature is found for the gap size for a suspension simply sheared between two walls, using two approaches: direct numerical simulations of cells, or a theoretical approach solving the lift+diffusion equation. The simulations are run with elastic capsules whose mechanical properties are such that they are essentially in a tank-treading regime. Agreement with the theory is obtained by considering a fitting parameter U/f of order 0.5, where U is the dimensionless prefactor of Eq. 1 and f that of the diffusive flux written as -f R 2 0 γΦ ∂Φ ∂y . This makes sense, since this value of the parameter is quite close to that obtained for lipid vesicles in tank-treading regime: in [41] we found U to be of order 0.1 and in [42] f is shown to be of order 0.06, which makes U/f ∼ 1.7. It was also argued in [18] that the crowded environment in which these interactions take place justify to consider an effective viscosity of the carrying fluid higher than that of the plasma. Indeed, from [51], one can see that for RBCs in plasma, U/f 2 0.016/2.8 0.006, which would certainly not be sufficient to create a cell free layer. Of course, it is clear that a cell lifting in a crowded environment will not lift the same way as an isolated cell, because the neighbors may modify its dynamics and response to flow, and will screen the wall. Considering tank-treading capsules for mimicking the collective behavior of red blood cells as in [18] or [17] appears thus as relevant for obtaining accurate predictions, though the choice of the effective viscosity contrast must be debated, as in Ref. [3]. Still, it hides the subtle modification of lift mechanisms by cell-cell interaction processes.

By dimension analysis in a 2D Poiseuille flow, it is possible to estimate a typical time describing the time needed for RBC to migrate frome the wall (y = 0) to the distance corresponding to the CFL (y = d) :

τ L = wd δs+1 v max ξ s R δs+1 0 . (10)
Shear-induced diffusion also contributes to equilibration, with an equilibration time deduced from the expression of the diffusive flux below,

τ D = w 3 f v max Φ 0 R 2 0 , (11)
where Φ 0 is a typical volume fraction in cells.

Using results from Fedosov [2], for a channel of typical size w = 5 µm, with Φ 0 = 30% and a CFL of size d = 2 µm, the distance needed for all RBC to pass through the CFL is τ L v max = 200 µm and the one to obtain an equilibration is τ D v max = 150 µm. The found typical distances are of the same scale and indicate that, in the microcapillary network, except in the smalles capillaries, equilibration is not likely to happen between two successive bifurcations, which are often very close to each other. [1,26,72,73]. Non stationary dynamics must, therefore, be considered when looking at network flow. The factor 2 that we found in the ξ values for the denser (less deformable) and less dense (more deformable) cells also indicates that uneven distribution of cells depending on their age may also occur in the microvasculature.

FIG. 1 .

 1 FIG. 1. (a) Sketch of the experimental channel. The separation line of the incoming flow is indicated by the violet dashed line. (b) Trajectory of a RBC migrating at high confinement (one single window). (c) Position of different RBCs in successive windows for a smaller confinement. RBCs are chosen to represent the median trajectory (see bottom middle panel in Fig. 4).

FIG. 2 .

 2 FIG.2. Grey points: measured velocity of RBCs, rescaled by vmax from Eq. (6) which is obtained from the fitting of experimental data using Eq. (8). The dashed line represents the velocity of the fluid vx without RBC at z = 0 using Eq. (5) and the full line represents the fitted velocity using Eq. (8).

, 1 . 6 ×

 16 FIG.3. Trajectory of RBCs in a channel with w = 10.6 µm (red). Fitting curves (black) using different couples of parameters (δ, ξ, y0), in1 ○ all parameters are free, giving at the end (0.14, 4.1 × 10 -3 , 0.38),2 ○, 3 ○,4 ○ are fits using a fixed δ giving respectively (0, 3.7 × 10 -3 , 0.38), (1, 7.7 × 10 -3 , 0.37), (2, 1.6 × 10 -2 , 0.36). The surface plot in insert shows the residual |(ỹ -ỹ th )/ ỹd | in the (δ, ξ) parameter space. Dark blue color represents the minimal residual.

3 .

 3 The blue zone corresponds to < |(ỹ -ỹth / ỹd | >≤ 0.5.

FIG. 4 .FIG. 5 .

 45 FIG.4. Top left: mean square error ∆ξ as a function of δ. Other plots represent the median trajectory ỹ(x) of RBCs for 5 different confinements (red) and their respective fitting curves (black) found from Eq. (9) using fixed δs and fixed ξs with ỹ0 a fitting parameter.

ACKNOWLEDGMENTS

This work was partially supported by CNES (Centre National d'Etudes Spatiales). The authors' team is part of LabEx Tec 21 (Investissements d'Avenir -grant agreement ANR-11-LABX-0030). S. L. thanks LabEx Tec 21 for a Proof of Concept fellowship.