Antoine Cheron

Johann Bourcier

Olivier Barais

Antoine Michel

Comparison Matrices of Semantic RESTful APIs Technologies

Keywords: Hateoas, Semantic REST, Comparison, Linked Data, Web

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

Introduction

Today, RESTful APIs [START_REF] Fielding | Architectural Styles and the Design of Network-based Software Architectures[END_REF] have become the de-facto standard for building web applications. The main reason behind this popularity lies in the appropriate trade-off between the facility to build such applications and the benefits provided by this approach in such an opened large-scale distributed system: evolutivity, scalability and loose-coupling. However, 95% of APIs are not RESTful [START_REF] Rodríguez | Rest apis: A large-scale analysis of compliance with principles and best practices[END_REF] as they claim.

Until today, no single standard has emerged to design truly RESTful APIs. Consequently, software architects are facing the challenge of selecting the right technologies for the design and implementation of these systems. Typically, a software architect has to select the right interface description language (IDL), interchange format and framework to ease the development of such APIs.

In addition, a new trend has recently emerged to create RESTful APIs that carry their own semantics, they are called Semantic RESTful APIs [START_REF] Salvadori | A maturity model for semantic restful web apis[END_REF]. It is a vision that proposes to make fully REST-compliant APIs compatible with the Semantic Web [START_REF] Berners-Lee | The semantic web[END_REF] and Linked Data [START_REF] Berners-Lee | Linked data principles[END_REF]. From our experience at FABERNOVEL, we found that building such APIs does not require much more effort than truly RESTful systems, whereas it offers great benefits, such as loose-coupling, automated API mash-ups [START_REF] Benslimane | Services mashups: The new generation of web applications[END_REF], machine-interpretability and very powerful querying.

However, the design of semantic RESTful APIs considerably increases the complexity for the architect to choose the appropriate technology. Indeed, the specific criteria and properties to be taken into account are not explicit when choosing an IDL, an interchange format and a framework. The industrial needs are growing for proper tools to support trade-off decisions of the architect; a tool that would help him/her to understand the consequences of a design decision, i.e. the characteristics and limitations of each approach.

In this paper, we propose to fill this gap by providing three decision matrices that help architects to choose the technologies that will best fit their needs. The main contributions of this paper are:

three comparison matrices of interchange formats, interface description languages and frameworks that help choosing the appropriate set of technologies to build Semantic RESTful APIs; key features that are missing from state-of-the-art technologies to assist and make the creation of Semantic RESTful APIs more beneficial.

Using these comparison matrices, we illustrate their usage on an industrial case and draw the outline of a research road-map to ease the adoption of Semantic RESTful APIs in the industry.

The remainder of this paper is organized as follows. Section 2 provides the required background on Semantic REST APIs and the reference maturity model to choose the functionality level of an API along, with its limitation. The two following sections describe our comparison matrices and an illustration that highlights the benefits of our proposition. Finally, section 5 discusses the role of the existing frameworks to build Semantic REST APIs.

Background

This section describes the main concepts related to the design and implementation of Semantic RESTful APIs and the process of selecting an API functionality.

Semantic RESTful services Combining REST with Semantic Web and Linked Data is a promising path since it enables the description of APIs that can change without breaking client applications. These APIs advertise their available state transitions, therefore enabling automatic composition to create high level services [START_REF] Alarcon | Rest web service description for graph-based service discovery[END_REF]. Smart software agents can then automatically discover the suite of operations to realize complex workflows and even make APIs compatible with voice assistants. This is achieved by semantically enriching the data and operations of REST systems with Semantic Web ontology technologies and by linking resources to other resources.

Selecting an API functionality level

Today, Web systems offer a wide range of functionalities. For example, they may offer multiple media types or a single one, comply with the HTTP protocol or use it as a transfer protocol, or even semantically describe their resources. This diversity can make the process of comparing and selecting the minimum set of features to be implemented very time-consuming. Maturity models have been proposed as a solution to this problem [START_REF] Paulk | Capability maturity model, version 1.1[END_REF][START_REF] Salvadori | A maturity model for semantic restful web apis[END_REF].

In companies, architects use them to decide features which must be supported by their APIs. In general, a maturity model is a scale that represents the compliance of a technology with a given architecture. To reach a level, a technology must meet each constraint of the targeted level and the previous levels. Currently, the de-facto standard in the industry is the Richardson Maturity Model [START_REF] Fowler | Richardson maturity model[END_REF], which targets building REST APIs. However, we recommend using the WS3 maturity model [START_REF] Salvadori | A maturity model for semantic restful web apis[END_REF] as it combines the models proposed by Richardson, and SoHA [START_REF] Soha | [END_REF], and extends them with semantic and documentation constraints. The WS3 maturity model In [START_REF] Salvadori | A maturity model for semantic restful web apis[END_REF], authors describe the WS3 maturity model for classifying Semantic REST Web APIs. It classifies APIs along three independent dimensions: design, profile and semantic, as shown in Fig. 1. The design dimension represents the different modeling strategies adopted for designing the technical access to a Web API through four levels: (i) RPC, (ii) resources have dedicated URI and the API is stateless, (iii) operations on a resource are mapped to HTTP verbs in compliance with the protocol and (iv) the smallest data unit that can be handled by operations is the resource. The profile dimension reflects the quality of documentation that can be interpreted by software agents through two levels. The first level: interaction profile, requires the description of all available HTTP operations and how to trigger them. The second level: the domain profile, requires the description of domain specific details such as the order of operation execution, pre-and postconditions, business constraints, etc.

The semantic dimension represents the use of semantic technologies through two levels. To reach the Semantic Description level, an API must semantically describe properties and operations of resources. The next level: Linked Data, is reached when the API semantically describes relationships between resources. Usage In their paper [START_REF] Salvadori | A maturity model for semantic restful web apis[END_REF], Salvadori et al. propose to rate systems along each dimension independently, with a score going from 0 to the number of levels in the dimension. For example, a non-documented API with no semantic support that reaches level 3 of the Richardson Maturity Model will be rated D3-S0-P03 . As another example, a system that supports HATEOAS and provides a swagger-like documentation along with the data is rated D3-S0-P24 .

Discussion on the WS3 maturity level

At FABERNOVEL, we experienced two limitations to the applicability of the maturity model to a wider audience. These limitations are related to the Atomic Resources level and the granularity of the WS3 levels.

According to WS3, the Atomic Resources constraint requires that the resource is the smallest data unit handled by operations. Respecting this constraint may introduce negative properties in the API. Let us consider an API handling insurance contracts which offers read and update operations on the postal address, email address and insurance manager. Two solutions can be considered to respect the Atomic Resources constraint. The first solution is to create one resource, where every properties can be modified at once, which increases the risk of concurrent modification. With this solution, the API would have two operations. The second solution is to create one resource for each concept: contract, email address, postal address and the manager. The API would have eight operations. This solution increases dramatically the number of operations which complexifies the documentation and maintainability. Another solution would be to create one resource with four operations: (i) read, (ii) update email, (iii) update postal address and (iv) update manager. This solution lowers the concurrency risk while maintaining a reasonable complexity and offering meaningful operation names. Unfortunately this solution breaks the Atomic Resources constraint. We therefore argue that respecting this last constraint may not always lead to better API quality.

The second limitation relates to the granularity of the maturity levels. Indeed, each level implies more than one feature. This granularity allows for a coarsegrained categorization of systems. However, to precisely differentiate systems based on the features they implement, a deeper study is needed. Given two systems that reach P1, which means they describe all available HTTP operations and how to trigger them, one might also describe its authentication process and errors while not the other one. And yet they reach the same maturity level. We therefore argue for a finer grain categorization of APIs.

Comparison Matrices

We propose three detailed matrices which address the limits of WS3 identified in the previous section. The proposed matrices enable the comparison of technologies along a set of precise criteria to highlight their differences. These matrices extend the WS3 levels by adding new criteria which are used in practice (see section 3.1) and not linked to any WS3 levels.

Insights from developers and architects

We interviewed 14 developers and architects from FABERNOVEL and clients on their experience with Semantic REST technologies. Raw results and the analysis are available online5 . Our key findings are:

-Selecting the technology: 10 respondents have already built Semantic REST APIs: 30% spent more than two weeks selecting the technologies; 80% reported that the most difficult task was to understand the feature provided by each technology. -Interchange Formats: 6 out of 7 did not find a technology providing all required features (most often the missing features were the description of HTTP operations with their data model (3/8) and the Linked Data (2/8)). -Interface description languages: All respondents said that none of them provide all required features (60% said they lack the ability to describe links to other resources and business constraints; and 20% of them would like to model the resources as finite state machines (FSM)). -Frameworks: 6 out of 7 reported that no framework offered the required feature. The missing features are related to the auto-documentation of the API, the automatic generation of link and a mechanism to model resources as FSM. -Technology score: The median value of the score is 2/5. These results emphasize the difficulties in selecting technologies associated to Semantic REST APIs. They also highlight that these technologies are not yet mature and give a rough idea of the missing features.

Comparison Matrices Design Method

The design of our comparison matrices follows a 5-step sequential process: (i) search for candidate technologies, (ii) select candidate technologies, (iii) read carefully each candidate technology, (iv) elaborate fine grain criteria to characterize and differentiate technologies, (v) verify that the elaborated criteria highlighted the differences between technologies. We looped on step (iv) and (v) to avoid duplicating criteria or hiding important details.

The research of candidate technologies (step i) was done by:

1. Searching Google and Google Scholar for Semantic REST Technologies using combinations of keywords from the set: ["web", "semantic", "restful", "rest", "service", "API", "interface", "description", "documentation", "language", "modeling", "hypermedia", "document", "format", "RDF", "datainterchange", "linked data", "hateoas", "rest api", "framework"]; 2. Searching Google Scholar for tools automating tasks from services description, using keywords: "matchmakers", "service composition", "service discovery", "rest service analysis", "automated mashups", we then selected papers and technologies from their references and the papers that cite those we selected; 3. Searching the proceedings of ICWE and WS-REST.

We selected 81 papers, standards, articles and web pages (step ii) based on abstract or introduction. We selected documents that were specifications of interface description languages or models, frameworks supporting HATEOAS features, interchange formats that support RDF or HATEOAS features, comparisons between these technologies or tools leveraging them. We considered frameworks available as programming libraries that helps implementing HTTP APIs in any programming language. We opened our research to technologies from the 1990s to today and retained those that are still available today.

Then, we read the specification of each chosen technology (step iii) and elaborated classification criteria (step iv). We included those of the H Factor6 which is a measurement of the level of hypermedia support and sophistication of a mediatype. Others were carefully designed to highlight differences between technologies, based on the core design of the technologies, the features they provide and the details of the WS3 maturity model. All the material is available online 7 .

As a final step (step v), we read the specifications again to verify results and validate that the selected criteria highlighted differences and commonalities well.

Popularity criteria We defined a popularity criteria to provide a rough idea of the community support and the likelihood of the technology to last in time. It respects the following rules: 0 -Not enough to reach 1; 1 -More than 100 questions on Stack Overflow AND (2500+ NPM weekly downloads OR 100+ maven usages); 2 -More than 400 questions on Stack Overflow AND (500.000+ total downloads OR 15.000+ NPM weekly downloads OR 500+ maven usages).

Interface Description Languages

Interface Description Languages (IDLs) provide a vocabulary to document domain, functional and non-functional aspects of an API. We identified 16 candidates that are classified according to 31 criteria in Fig 2 . Among them, 4 are meta-models from conference papers [START_REF] Haupt | A model-driven approach for rest compliant service[END_REF][START_REF] Mitra | Rapido: a sketching tool for web api designers[END_REF][START_REF] Schreier | Modeling restful applications[END_REF][START_REF] Zuzak | Formal modeling of restful systems using finite-state machines[END_REF]. The 11 others are open-source projects or W3C recommendations.

In [START_REF] Mitra | Rapido: a sketching tool for web api designers[END_REF] authors present a tool to sketch CRUD or Hypermedia APIs. On the latter mode, users sketch the application using state-machines and then obtain a description in the HAL or Collection+JSON format. [START_REF] Schreier | Modeling restful applications[END_REF] models each resource type as a finite-state-machine with deterministic transitions and conditions to inform about the availability of transitions. However, they are not modeled in more details, which make them not machine-interpretable. In [START_REF] Zuzak | Formal modeling of restful systems using finite-state machines[END_REF], authors propose to model systems as non-deterministic state machines. This method thus makes software agents unable to discover the set of messages to exchange in order to make an operation available. Haupt et al. [START_REF] Haupt | A model-driven approach for rest compliant service[END_REF] propose a multi-layered model that separates the domain model from the URI model. However, resources have a fixed model, which prevent them from having one data model per state.

It is important to note than when IDLs and interchange formats are both compatible with RDF, they can be combined to form a file format usable as data-interchange format and IDL. This has great benefits to lower the overall complexity and increase the evolvability of the system. Synthesis First, the matrix highlights the fact that most technologies help with building mature systems on the design dimension and interaction profile level of the profile dimension, D3-P1 following the WS3 categories. On the other hand on the semantic dimension, we notice that 5/16 technologies support the use of RDF vocabulary, which allows to build Linked Data APIs. As a reminder, this is required to reach full Semantic REST compliance. Moreover, by supporting the use of RDF vocabulary, IDLs can be enriched to reach a higher level of maturity.

Among the technologies, four can be distinguished by the number of criteria they meet: Hydra (18), RADL [START_REF] Soha | [END_REF], OpenAPI [START_REF] Serrano | Linked rest apis: A middleware for semantic rest api integration[END_REF] and RESTdesc [START_REF] Serrano | Linked rest apis: A middleware for semantic rest api integration[END_REF]. Ope-nAPI is the only one that has no support for RDF. Thus, it helps in building systems up to D3-P2-S0 on the WS3 scale. On the other hand, Hydra, RADL and RESTdesc support the use of RDF vocabulary, which makes these technologies better suited to build systems that are mature on the semantic dimension.

Towards HATEOAS APIs From the matrix, we notice that most technologies target the documentation of the API in a single, non-splittable file. Hence, they are not suited to provide hypermedia controls at runtime.

On the other hand, only one approach, [START_REF] Schreier | Modeling restful applications[END_REF], supports the description of the conditions that determine the availability of a link, and none makes this metadata machine-interpretable. This makes software agents unable to find a way to make an operation available when it is not.

Towards better-documented APIs Only four technologies support the description of business constraints which lowers coupling and improves user experience, e.g., with the automatic generation of forms with client-side validation.

Finally, we note that most scientific publications recommend the modeling of RESTful systems with state-machines whereas open-sourced or W3C IDL authors don't consider this design method. And yet, the use of deterministic state-machines eases the determination of the available operations of a resource.

Data-interchange formats

These formats provide a data-structure, a vocabulary and a layout to represent a resource and its meta-data at runtime. When the API does not need to send meta-data, JSON and XML are the two widely used formats in the industry.

On the other side, when the system to be built have to support a hypermedia interchange format, none is considered as a standard today. We selected 11 candidate technologies, which are classified in Fig. 3 according to 24 criteria. JSON is included for comparison purposes.

Synthesis First, from this matrix, we notice that formats can be differentiated based on their compatibility with RDF. Indeed, RDF formats (Turtle, RDF XML and JSON-LD) propose very few features by default because they can be enriched with RDF vocabularies. To depict what is achievable by combining vocabularies with a RDF format, we selected two vocabularies: Hydra and SHACL, a RDF schema validation vocabulary, that we combined with JSON-LD and evaluated them. As a result, they match 12 more criteria than JSON-LD alone. From this, we infer that combining RDF formats with vocabularies allow building mature Semantic REST systems. However, this requires additional effort to find relevant vocabularies. On the other hand, non-RDF formats help building systems that can be mature on the profile dimension but not on the semantic dimension.

Furthermore, the matrix shows that no format supports the description of constraints despite the fact that it can be leveraged to reduce coupling and improve the user-experience.

Finally, it highlights that no format advertise the state of the resource even though most scientific approaches we found describe REST APIs as state-machines.

Implementation Frameworks

Implementation frameworks are software libraries that guide developers through the implementation of Web APIs. We limit the comparison to frameworks that claim to support HATEOAS. We identified six frameworks that do so. Frameworks to build Semantic Web Services are excluded because their triple-centric approach differs too much from REST.

Among the selected papers, in [START_REF] Salvadori | A framework for semantic description of restful web apis[END_REF] authors propose Hypermedia Web API Support, a Java framework based on JAX-RS 2.0 that offers annotations to semantically describe REST APIs. The end result is the description of the whole API in a JSON-LD document enriched with the Hydra vocabulary. Unfortunately, the framework is not available in Maven Central. In [START_REF] Parastatidis | The role of hypermedia in distributed system development[END_REF] Parastatidis et al. present Restfulie, a framework that uses resources, state transitions and contentnegotiation as its core building blocks. We found 4 other frameworks that support HATEOAS features. They are all classified in Fig. 4 according to 23 criteria. Synthesis Despite the fact that only one framework enforces the Atomic Resources constraint, all frameworks allow to reach the highest level of maturity on the design dimension easily. This is because supporting the Atomic Resources constraint only requires developers to use the data model of the resource as the input of write operations and as the output of read operations. We notice that only API Platform and Restfulie offer a mechanism to model relations between resources from which links are generated, instead of adding them programmatically in the response, thus increasing maintainability. Otherwise, most frameworks do not ease the semantic and domain description of APIs. To us, this is the biggest challenge framework designers should tackle. Last, as for IDLs, most frameworks creators do not provide mechanisms to describe resources as state machines, thus not taking advantage of its benefits.

Matrices usage example

In this section, we present the service of an insurance company that manages insurance contracts to illustrate how the presented comparison matrices can be used in a real world scenario. This example is a light version of projects we have carried out at FABERNOVEL for large French insurance companies.

Domain description

To manage insurance contracts, the service holds five kinds of resources: (i) thirdparties, (ii) contracts, (iii) warranties, (iv) cases and (v) services. Third-parties, for example customers or contractors, enter into contracts with the company. These contracts include warranties from the closed list that the company offers. For example a Person A has the following warranties: (i) damage coverage with a deductible of $500 and a maximum repair amount of $30.000 and (ii) premium vehicle loan in the event of immobilization of the damaged vehicle. A contract can have several cases. When an customer of the insurance has a claim the company creates a case that holds its details and the services provided to the insured. For example, Person A has a car accident, he opens the insurance's web application and reports a claim, which leads to the company opening a case. His car has been destroyed and he is expected to attend a diner. Thus, on the app, he asks for the loan of a car that he will immediately recover.

Technological constraints

The service has to communicate with both internal and external components. Internal components are front-end applications, such as mobile or web applications, and other kernel services, such as payments. External components are contractors APIs, for example taxi or mechanics companies.

In the insurance domain, there is a huge amount of business rules that determine (i) the warranties an insured can include in a contract and (ii) the available services for a case, based on the specificity of the given case and the warranties of the contract. Writing and maintaining these rules both on the server and its clients is very costly and error prone. Thus, we decided to keep these business rules on the server-side only. This constraint leads to the use of HATEOAS.

The project constitutes the core of the companys business, it should then be built with state-of-the-art technologies such as Linked Data. This enables the automatic creation of mash-ups and the use of a HyperGraphQL8 middleware to easily query the whole IS [START_REF] Tuchinda | Building mashups by demonstration[END_REF]. Moreover, considering that the contractors providing services are very diverse and numerous, the interactions with their APIs should leverage the automatic discovery and composition offered by the use of RDF semantics.

There is also a high probability that new client systems will be built in the future, the API must document its resources, resource attributes, operations, URI templates, HTTP verbs, hypermedia controls and errors in a machineinterpretable way. Moreover, because the service applies the CQRS pattern 9 we needed the IDL to enable associating an operation to its own input and output data model. Last, to minimize the disruption for software developers, we have chosen to keep the interchange formats as close as possible to what developers already know. It has therefore to be entity-centric, based on JSON and its structure had to be as close as possible to a JSON document without meta-data.

Selection of the technologies

From these constraints, we selected the set of criteria and features that the technologies should provide. These criteria are checked in the last column of Fig 2, 3, and4. We then count the number of criteria that were provided by each existing technology. Results are presented in Fig 5,6 and 7. For each matrix, the three technologies matching the highest amount of selected criteria are highlighted in green. It is important to note that the technologies do not have to match every criteria to be selected. Most of the time, missing features can be implemented afterwards, or proposed to the maintainers of the technologies. Interface Description Languages Hydra, OpenAPI and RADL are the technologies matching the highest number of selected criteria. However, none matches all criteria. Hydra lacks the ability to describe non-functional properties and mediatypes, which can be done with other RDF vocabularies. RADL lacks the ability to semantically describe resources models, operations, errors and non-functional properties, which can also be done with other vocabularies. On the other hand, OpenAPI does not support the usage of RDF vocabulary. In this project, we have chosen to setup both Hydra and OpenAPI. OpenAPI because it has most features and it is a must-have today because of its tooling and popularity. Hydra because it can be easily completed with other vocabularies and used with JSON-LD.

Interchange Formats Mason, JSON-LD + Hydra are the two technologies matching the highest number of selected criteria. JSON-LD + Hydra + SHACL is ignored as it does not match more selected criteria than without SHACL. While JSON-LD + Hydra lacks the ability to describe non-functional properties, Mason does not allow to use RDF vocabularies. Being incompatible with RDF requires a lot more effort to compensate than finding another vocabulary. This explains why JSON-LD + Hydra was preferred over Mason in this context.

Implementation frameworks API Platform, Spring HATEOAS and Hypermedia Web API Support [START_REF] Salvadori | A framework for semantic description of restful web apis[END_REF] are the three technologies matching the highest number of criteria. The latter is immediately removed from the candidates because no public library is available. In this example, API Platform should be preferred over Spring HATEOAS because it matches five more criteria than Spring. However, developers of the companies we worked with know Java and not PHP. Moreover the Spring framework is very popular with them, which compensates the need to develop some features by hand. This is why we decided to go with Spring.

Easing the selection of the technologies We have developed an open-source web recommender system10 in which the user selects the type of technology, the required criteria and finally score each criterion to indicate its importance. In return, the web application presents the list of technologies that meet the required criteria ordered by their respective score. By offering a two-step wizard, we reduce the process of identifying relevant technologies to a few minutes.

Discussion

This section provides our perspective on why no standard solution exists to meet all our criteria, and highlights the possibility of new research initiatives.

First, there is no de-facto neither IETF or W3C standard Interchange Format for building Semantic REST APIs. The technology the more likely to become a standard is JSON-LD, for which a W3C Working Group is active. In addition, none of the existing interchange formats support all the criteria described above, making it likely that new formats will emerge. For this reason, frameworks supporting Semantic REST APIs will rely on formats that are likely to evolve, which will require additional effort and costs. This reduces the likelihood of framework editors to invest time in developing such features.

To us, the second and also the most important reason is that the well-known and widely used tools do not rely on Semantic REST APIs to provide additional and useful features. Among the possible functionalities, we envision various tools to automate API testing, REST client generation, API gateways, middleware and smarter desktop REST clients. We believe that this limits the adoption of Semantic REST APIs because the cost of building these APIs is not perceived as offering a sufficient short-term return on investment.

Related Work

In [START_REF] Verborgh | Survey of semantic description of REST APIs[END_REF], and [START_REF] Salvadori | A maturity model for semantic restful web apis[END_REF], authors justify the need to provide a semantic description of REST APIs to avoid that programmers who develop client applications have to understand in depth several APIs from several providers. Based on this motivation, they survey academic approaches to add semantic to such APIs description and technique to automatically compose restful services.

In [START_REF] Serrano | Linked rest apis: A middleware for semantic rest api integration[END_REF] authors present a framework for REST-service integration based on Linked Data models. First, API providers semantically describe their REST services. Second, API consumers express data queries with SPARQL. Then, they use a middleware developed by the authors that automatically compose API calls to respond to data queries with a RDF graph. At the first step, authors needed to find and select a RDF-compatible Interface Description Language, which is precisely the kind of use-case our approach addresses. Because they could not find an existing one fitting their needs, they designed a new one by leveraging existing technologies such as MSM [START_REF] Pedrinaci | Toward the next wave of services: Linked services for the web of data[END_REF], Hydra, RAML and OpenApi.

In [START_REF] Tuchinda | Building mashups by demonstration[END_REF], Tuchinda et al. describe a programming-by-demonstration approach to build mashups by example. Instead of requiring a user to select and customize a set of widgets, the user simply demonstrates the integration task by example. Their approach addresses the problems of extracting data from web sources, cleaning and modeling the extracted data, and integrating the data across sources. It illustrates the benefits of getting meta-data on top of services to improve the definition of mashups and decrease the coupling between information system building blocks and the complexity of developing mature Semantic REST APIs. In [START_REF] Duke | Telecommunication mashups using restful services[END_REF], Duke et al. propose an approach to reduce the complexity for describing, finding, composing and invoking semantic rest services. They mainly provide an approach where they show how they can combine services when they get semantic information.

Other research efforts were done to lower the entry barrier for developing mature Semantic REST APIs. Among them is the semi-automatic annotation of web services as done by Patil et al. in [START_REF] Patil | Meteor-s web service annotation framework[END_REF]. Their contribution could help significantly increase the number of semantically described services, in case their work is open-sourced and updated to support nowadays popular technologies.

Findings Summary

In this paper, we have presented three comparison matrices that assist architects in choosing Semantic REST APIs enabling technologies that meet their needs. Through a real example, we have illustrated how the use of these matrices simplifies the choice of these technologies. As stated in the paper, technologies should be chosen not only according to the number of criteria they meet, but also according to the specific needs of the project. To facilitate this selection, we have developed an assistant available online.

We also pointed out some interesting features missing in current technologies. The description of constraints and conditions indicating the availability of state transitions is ignored by IDLs, vocabularies, interchange formats and frameworks. On the other hand, resource modeling as FSM is not available in most frameworks. More importantly, well-known tools do not take advantage of the power of Semantic REST APIs to provide additional and useful features.

Based on these findings, we identify areas for improvement in the tools around Semantic REST APIs that we believe can increase its adoption. By leveraging the semantic description and advertising of state transitions and non-functional properties, automated testing tools can become smarter, REST client libraries can lower the coupling with the server and automate tasks such as login, and middleware can automatically create responses from the composition of several APIs.

Fig. 1 .

 1 Fig. 1. WS3 Maturity Model (from [15])

Fig. 2 .

 2 Fig. 2. Interface Description Languages Comparison Matrix

Fig. 3 .

 3 Fig. 3. Data-interchange Formats Comparison Matrix

Fig. 4 .

 4 Fig. 4. Implementation Frameworks Comparison Matrix

Fig. 5 .

 5 Fig. 5. Results for interface description languages

Fig. 6 .

 6 Fig. 6. Results for data interchange formats

Fig. 7 .

 7 Fig. 7. Results for implementation frameworks

D3-S0-P0: Atomic Resources Design, no Semantic Description, no Profile description

D3-S0-P2: Atomic Resources Design, no Semantic Description, Domain Profile

https://github.com/AntoineCheron/comparison-matrices-semantic-rest-api-techno

http://amundsen.com/hypermedia/hfactor/

https://github.com/AntoineCheron/comparison-matrices-semantic-rest-api-techno

https://www.hypergraphql.org/

https://martinfowler.com/bliki/CQRS.html

https://antoinecheron.github.io/morice/


Jérôme Jacob
email: jerome.jacob@univ-amu.fr

Orestis Malaspinas

Pierre Sagaut

A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann method-based large eddy simulation

Keywords: Large eddy simulation, Lattice Boltzmann method, circular cylinder Reynolds 3900

A new Lattice Boltzmann collision model for large eddy simulation (LES) of weakly compressible flows is proposed. This model, referred to as the Hybrid Recursive Regularised Bhatnagar-Gross-Krook (HRR-BGK) model, is based on a modification of previously existing regularised collision models defined with the BGK Lattice Boltzmann method (LBM) framework. By hybridising the computation of the velocity gradient with an adequate Finite Difference scheme when reconstructing the non-equilibrium parts of the distribution function, a hyperviscosity term is introduced in the momentum equation, whose amplitude can be explicitly tuned via a weighting parameter. A dynamic version of the HRR-BGK is also proposed, in which the control parameter is tuned at each grid point and each time step in order to recover an arbitrarily fixed total dissipation. This new collision model is assessed for both explicit and implicit LES considering the flow around a circular cylinder at Re = 3900. The dynamic HRR-BGK is observed to yield very accurate results when equipped with Vreman's subgrid model to compute the target dissipation.

Introduction

Large eddy simulation (LES) is now a mature simulation technique for high-resolution unsteady simulation of turbulent flows [START_REF] Garnier | Large-eddy simulation for compressible flows[END_REF][START_REF] Sagaut | Large eddy simulation for incompressible flows: an introduction[END_REF][START_REF] Sagaut | Multiscale and multiresolution approaches in turbulence[END_REF]. This technique, which is under development since the early 1960s, is now implemented in almost all computational fluid dynamics (CFD) tools, including commercial software. Among the key issues faced when developing LES-based simulation tools, one must cite the development of accurate models for subgrid scales and adequate numerical schemes that guarantee numerical stability whose induced dissipation does not overwhelm the effect of physical subgrid models. To this end, a huge amount of subgrid models and numerical methods have been proposed since the 1980s. This issue is far from being a trivial one, since it is now accepted that the discretisation errors and the subgrid models are implicitly coupled in a genuinely nonlinear way which escapes a full mathematical analysis. Therefore, finding an optimised pair (numerical scheme and subgrid model) is a very difficult task, since it has been observed that some partial error cancellation may occur, leading to: (i) a kind of super-convergence of LES results (an unexpected high accuracy is obtained in such a case since modelling errors and numerical errors partially balance each other) and (ii) the counter-intuitive result that in some cases an increase of the order of the numerical scheme or grid refinement may lead to a decrease of the global accuracy of the results [START_REF] Meyers | Is plane channel flow a friendly case for testing of large-eddy simulation subgrid-scale models?[END_REF][START_REF] Meyers | Optimal model parameters for multi-obective large-eddy simulations[END_REF][START_REF] Meyers | A computational error assessment of central finite-volume discretizations in large-eddy simulation using a Smagorinsky model[END_REF]. A promising way to solve this problem is to develop stabilised numerical schemes whose leading error terms mimic explicit subgrid models, leading to the implicit large eddy simulation (ILES) approach [START_REF] Grinstein | Implicit large eddy simulation: computing turbulent fluid dynamics[END_REF]. Using such a scheme, one may expect to guarantee numerical robustness and physical accuracy at the same time in a better controlled way.

LES has been implemented within the Lattice Boltzmann method (LBM) framework [START_REF] Chen | Lattice Boltzmann method for fluid flows[END_REF][START_REF] Guo | The lattice Boltzmann method and its applications in engineering[END_REF][START_REF] Kruüger | The lattice Boltzmann method. Principles and practice[END_REF][START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF][START_REF] Succi | The lattice Boltzmann equation for fluid dynamics and beyond[END_REF] since the 1990s [START_REF] Ansumali | Kinetic theory of turbulence modeling: smallness parameter, scaling and microscopic derivation of Smagorinsky model[END_REF][START_REF] Chen | Extended Boltzmann kinetic equation for turbulent flows[END_REF][START_REF] Chen | Expanded analogy between Boltzmann kinetic theory of fluids and turbulence[END_REF][START_REF] Chen | Analysis of subgrid scale turbulence using the Boltzmann Bhatnagar-Gross-Krook kinetic equation[END_REF][START_REF] Chen | A large-eddy-based lattice Boltzmann model for turbulent flow simulation[END_REF][START_REF] Eggels | Direct and large-eddy simulation of turbulent fluid flow using the lattice-Boltzmann scheme[END_REF][START_REF] Filippova | Multiscale lattice Boltzmann schemes with turbulence modeling[END_REF][START_REF] Girimaji | Boltzmann kinetic equation for filtered fluid turbulence[END_REF][START_REF] Hou | A lattice Boltzmann subgrid model for high Reynolds number flows[END_REF][START_REF] Krafczyk | Large-eddy simulations with a multiple-relaxation-time LBE model[END_REF][START_REF] Malaspinas | Consistent subgrid scale modelling for lattice Boltzmann methods[END_REF][START_REF] Premnath | Dynamic subgrid scale modeling of turbulent flows using lattice-Boltzmann method[END_REF][START_REF] Teixeira | Incorporating turbulence models into the lattice-Boltzmann method[END_REF][START_REF] Weickert | Investigation of the LES WALE turbulence model within the lattice Boltzmann framework[END_REF], and both classical LES and ILES approaches have been proposed. The most common way to incorporate an explicit subgrid viscosity model is to modify the relaxation time in the collision model in order to recover the targeted effective total viscosity, defined as the sum of the molecular viscosity of the fluid and the subgrid viscosity. But the need to develop stabilised LBM schemes to handle high Reynolds number flows has also led to the definition of collision models with some ILES capabilities, such as Entropic methods [START_REF] Bösch | Entropic multi-relaxation models for simulation of fluid turbulence[END_REF][START_REF] Bösch | Entropic multi-relaxation models for turbulent flows[END_REF][START_REF] Ilio | Simulation of turbulent flows with the entropic multirelaxation time lattice Boltzmann method on body-fitted meshes[END_REF][START_REF] Dorschner | Transitional flows with the entropic lattice Boltzmann method[END_REF], Cascaded models [START_REF] Geier | De-aliasing and stabilization formalism of the cascaded lattice Boltzmann automaton for under-resolved high Reynolds number flow[END_REF][START_REF] Geier | Cascaded digital lattice Boltzmann automata for high Reynolds number flow[END_REF][START_REF] Geller | Turbulent jet computations based on MRT and cascaded lattice Boltzmann models[END_REF] and regularised methods [START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF][START_REF] Latt | Hydrodynamic limit of lattice Boltzmann equations[END_REF][START_REF] Latt | Lattice Boltzmann method with regularized pre-collision distribution functions[END_REF][START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization[END_REF][START_REF] Mattila | High-order regularization in lattice-Boltzmann equations[END_REF][START_REF] Montessori | Regularized lattice Bhatnagar-Gross-Krook model for two-and three-dimensional cavity flow simulations[END_REF]. While it has been observed in numerical experiments that these methods may lead to satisfactory results in turbulent flow simulations without adding an explicit subgrid model, the direct link between their build-in dissipation and the physical subgrid one is not fully understood.

This paper proposes a new improved recursive regularised LBM method based on the hybridisation of the computation of the velocity gradients when evaluating the coefficients of the reconstruction of the regularised non-equilibrium component of distribution functions, leading to the definition of an Hybrid Recursive Regularised Bhatnagar-Gross-Krook (HRR-BGK) model. In this new model, the numerical stabilisation can be explicitly tuned in order to recover a perfectly controlled amount of local dissipation. To get an accurate ILES method, the free parameter appearing in this new collision kernel can be set dynamically at each grid point and time step to a value which allow to recover exactly the same dissipation as the one provided by any classical LES subgrid model, according to an explicit formula given below. The scope of the paper is restricted to low-Mach number athermal flows.

This paper is organised as follows. The basic regularised BGK collision model used to develop the new collision scheme is reminded in Section 2. The new dynamic RR-BGK collision kernel is then discussed in Section 3. The analysis of its build-in dissipation and the bridge with ILES are presented in Section 4. Key features of the Vreman's subgrid model, which is used in this paper for explicit LES and tuning of the ILES simulations are reminded in Section 5. The approach is then assessed considering the turbulent flow at Re = 3900, which is a well-documented validation case. This test case is selected because it allows for a wall-resolving simulation without wall model in a flow configuration in which the subgrid model plays a very important role. The decoupling from the effects of a wall model is an important point here, since wall models may have a deeper impact on the solution than the subgrid model, masking the effect of the later. Here, the subgrid closure is known to have a leading effect on the development of the separated shear layers and the formation of the cylinder wake. It is worth noting that this is not the case considering plane channel flows on LBM grids. As a matter of fact, wall-resolving LES of plane channel flow would be based on a uniform grid with x + = y + = z + = 1 -2 near the wall, which is finer than classical DNS grid resolution requirement, yielding a poor interest for the assessment of the present subgrid closure approach.

Recursive regularised BGK LBM

The Boltzmann equation discretised in the velocity space reads

∂ t f i (x, t) + (ξ i • ∇)f i (x, t) = i , ( 1 )
where {ξ i } q-1 i=0 is the discrete velocity set, i is the collision operator. The most widely used collision operator in the lattice Boltzmann community is the BGK model, where is a relaxation towards a polynomial approximation of the Maxwell-Boltzmann distribution equilibrium, noted f (0

) i = - 1 τ f i -f (0) i , ( 2 )
where f

(0) i is the equilibrium distribution function f (0) i = w i ρ 1 + ξ i • u c 2 s + 1 2c 4 s H (2) i : uu , ( 3 )
with c s the speed of sound (which is a constant) and H (2) i = ξ i ξ ic 2 s I (I is the identity matrix). This model is accurate for low-Mach number, weakly compressible athermal flows. It suffers from a lack of numerical stability in high Reynolds number flows. Therefore several successful attempts have been made to make it more stable. Among these the regularised approach (see [START_REF] Latt | Hydrodynamic limit of lattice Boltzmann equations[END_REF][START_REF] Latt | Lattice Boltzmann method with regularized pre-collision distribution functions[END_REF]) is particularly appealing due to its simplicity. The collision operator of this model reads i = -

1 τ f (1) i , ( 4 )
where

f (1) i = - 1 2c 4
s H (2) i : P (1) , [START_REF] Meyers | Optimal model parameters for multi-obective large-eddy simulations[END_REF] with P (1) = q-1 i=0 H (2) i (f if (0) i ) is the deviatoric stress tensor. In Ref. [START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization[END_REF], a recursive procedure is introduced to add higher order terms to make the regularisation model even more stable, which was later used and generalised in Refs. [START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF][START_REF] Mattila | High-order regularization in lattice-Boltzmann equations[END_REF]. It is based on the relation between the equilibrium and non-equilibrium Hermite coefficients of the distribution function f i . By expanding f i in Hermite polynomials up to an order N (see [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF]) one gets

f i = w i N n=0 1 c 2n s n! H (n) i : a. ( 6
)
Then by using the standard Chapman-Enskog expansion [START_REF] Chapman | The mathematical theory of nonuniform gases[END_REF], one can decompose f i into two parts

f i = f (0) i + f (1) i , ( 7 )
where f (0) i is the equilibrium distribution function and f

(1) i is the off-equilibrium distribution function with f (0) i f (1) i . These two distribution functions can also be expanded in Hermite polynomials

f (0) i = w i N n=0 1 c 2n s n! H (n) i : a (n) 0 , ( 8 )
f

(1) i = w i N n=0 1 c 2n s n! H (n) i : a (n) 1 , ( 9 )
where a

(n) 0 and a (n)

1 are the Hermite coefficients of the equilibrium and off-equilibrium distributions respectively,

a (n) 0 = q-1 i=0 H (n) i f (0) i , ( 10
)
a (n) 1 = q-1 i=0 H (n) i f (1) i . (11)
In the athermal case, the equilibrium Hermite coefficients of order n are given by a

(n) 0 = a (n-1) 0 u, with (12)
a (0) 0 = ρ. ( 13
)
Using the Chapman-Enskog expansion as in [START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization[END_REF], one asymptotically first recovers the weakly compressible Navier-Stokes equations

∂ t ρ + ∇ • (ρu) = 0, (14)
ρ(∂ t u + u • ∇u) = -∇ • P, ( 15
)
where

ρ = i f i , ( 16
)
u = i ξ i f i , (17)
P = i (ξ i -u)(ξ i -u)f i = P (0) + P (1) = Ip -2ρνS, (18)
where S = (∇u + (∇u) T )/2 is the strain-rate tensor, p = c 2 s ρ is the pressure (this is the perfect gas law) and ν = c 2 s τ is the kinematic viscosity. Another very interesting property of the Chapman-Enskog expansion of Equation ( 2) is that the off-equilibrium coefficients are also related between themselves by a recursive relation (as are the equilibrium ones) and are given by

a (n) 1,α 1 •••α n = a (n-1) 1,α 1 •••α n-1 u α n + u α 1 • • • u α n-2 a (2) 1,α n-1 α n + perm(α n ) , (19)
where 'perm(α n )' stands for all the cyclic index permutations of indexes from α 1 to α n-1 (α n is never permuted), and where

a (2) 1,αβ = -2ρτ c 2 s S αβ , (20)
with

S αβ = (∂ α u β + ∂ β u α )/2
is the strain rate tensor. Until this point we made the hidden assumption that the discrete Hermite polynomials H (n) i have the same orthogonality relations between themselves that the continuous Hermite polynomials up to an arbitrary order n. In reality only, a limited number of Hermite polynomials possess this property: for the standard lattices D3Q15, D3Q19 and D3Q27 this is true up to order n = 2 (see [START_REF] Shan | The mathematical structure of the lattices of the lattice Boltzmann method[END_REF] for example). In [START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization[END_REF], one uses the D3Q27 lattice because of its better isotropy and therefore higher resemblance with the continuous case at the expense of a higher computational cost. To lower the amount of computations needed per grid point, we use here the D3Q19 quadrature.

The Hermite polynomials used here are a bit different than for the D3Q27 case. Here we will use only combinations of third-order Hermite polynomials which have the correct orthogonality properties. These are H (3) α,xxy + H (3) α,yzz , ( 21)

H (3) α,xzz + H (3) α,xyy , (22)
H (3) α,yyz + H (3) α,xxz , (23)
H (3) α,xxy -H (3) α,yzz , (24)
H (3) α,xzz -H (3) α,xyy , (25)
H (3) α,yyz -H (3) α,xxz . (26)
These Hermite polynomials being added to the set of H (0) , H (1) and H (2) the equilibrium distribution function is given by

f (0) i = w i ρ + ξ i • (ρu) c 2 s + 1 2c 4 s H (2) i : a (2) 0 + 1 2c 6 s (H (3) i,xxy + H (3) i,yzz )(a (3) 0,xxy + a (3) 0,yzz ) + 1 2c 6 s (H (3) i,xzz + H (3) i,xyy )(a (3) 0,xzz + a (3) 0,xyy ) + 1 2c 6 s (H (3) i,yyz + H (3) i,xxz )(a (3) 0,yyz + a (3) 0,xxz ) + 1 6c 6 s (H (3) i,xxy -H (3) i,yzz )(a (3) 0,xxy -a (3) 0,yzz ) + 1 6c 6 s (H (3) i,xzz -H (3) i,xyy )(a (3) 0,xzz -a (3) 0,xyy ) + 1 6c 6 s (H (3) i,yyz -H (3) i,xxz )(a (3) 0,yyz -a (3) 0,xxz ) (27)
and the off-equilibrium is

f (1) i = w i 1 2c 4 s H (2) i : a (2) 1 + 1 2c 6 s (H (3) i,xxy + H (3) i,yzz )(a (3) 1,xxy + a (3) 1,yzz ) + 1 2c 6 s (H (3) i,xzz + H (3) i,xyy )(a (3) 1,xzz + a (3) 1,xyy ) + 1 2c 6 s (H (3) i,yyz + H (3) i,xxz )(a (3) 1,yyz + a (3) 1,xxz ) + 1 6c 6 s (H (3) i,xxy -H (3) i,yzz )(a (3) 1,xxy -a (3) 1,yzz ) + 1 6c 6 s (H (3) i,xzz -H (3) i,xyy )(a (3) 1,xzz -a (3) 1,xyy ) + 1 6c 6 s (H (3) i,yyz -H (3) i,xxz )(a (3) 1,yyz -a (3) 1,xxz ) , (28)
where

a (2) 1,αβ = i H (2) i,αβ (f i -f (0) i ) = -2ρτ c 2 s S αβ , (29)
a (3) 1,αβγ = u α a (2) 1,βγ + u β a (2) 1,γ α + u γ a (2) 1,αβ . (30)
With these two modifications to the regularised model (see Equation 4), we have now the recursive regularised BGK Lattice-Boltzmann method (RR-BGK).

A new improved RR-BGK collision operator

In this paper, we propose a modified version of the regularisation procedure, which relies on the hybridisation of the computation of the velocity gradients with finite difference schemes in the computation of the recursive reconstruction parameters discussed in the preceding section. The deviatoric stress tensor in the above equation is replaced by a stabilising hybrid stress tensor P(1) which reads P(1) = P (1) σ

-(1 -σ )2ρτ c 2 s S FD , 0 ≤ σ ≤ 1, (31)
where σ = 1 is equivalent to the RR-BGK model.

Using the simple second-order centred finite difference approximation

g(x + x) -g(x -x) 2 x = ∂ x g(x) + 1 6 x 2 ∂ 3 x g(x), (32)
the finite-difference strain rate tensor S FD can be evaluated as

S FD αβ ∼ = 1 2 ∂ α u β + 1 6 x 2 ∂ 3 α u β + ∂ β u α + 1 6 x 2 ∂ 3 β u α , ∼ = 1 2 ∂ α u β + ∂ β u α + 1 6 x 2 (∂ 3 α u β + ∂ 3 β u α ) , ∼ = S αβ + 1 2 1 6 x 2 (∂ 3 α u β + ∂ 3 β u α ) , ( 33
)
where S αβ denotes the exact velocity gradient tensor.

The collision operator then reads

= - 1 τ f (1) i , (34)
where

f (1) i = f (1) i σ -(1 -σ ) ρτ c 2 s H (2) i : S FD , (35)
where f (1) i is given by Equation [START_REF] Bösch | Entropic multi-relaxation models for turbulent flows[END_REF]. A more efficient way to compute this last term is to evaluate it through the recursive relation using a following modified expression for the coefficients:

f (1) i = f (1) i ã(2) 1,αβ , (36)
with f

(1) i still given in Equation ( 28) and where ã(2)

1,αβ is ã(2) 1,αβ = σ i H (2) i,αβ (f i -f (0) i ) + (1 -σ )(-2ρτ c 2 s S αβ ) = σ a (2) 1,αβ + (1 -σ )(-2ρτ c 2 s S αβ ).
(37) The actual Boltzmann equation we will solve hereafter is given by

∂ t f i + ξ i • ∇f i = - 1 τ f (1) i σ -(1 -σ ) ρτ c 2 s H (2) i : S FD . ( 38
)
Performing the Chapman-Enskog expansion on this equation, one gets at the lowest order

∂ t f (0) i + ξ i • ∇f (0) i = - 1 τ f (1) i σ -(1 -σ ) ρτ c 2 s H (2) i : S FD . ( 39
)
Taking the zeroth-, first-and second-order moments of this last equation, one gets

∂ t ρ + ∇ • (ρu) = 0, ( 40
)
∂ t (ρu) + ∇ • (ρuu) = -∇p, (41)
∂ t (ρuu) + ∇ • (ρuuu) + ∇(ρu) + (∇(ρu)) T = - 1 τ P(1) . ( 42
)
Using the first two equations, the third can be rewritten as

-2ρc 2 s τ S + O(Ma 3 ) = P(1) , (43)
where Ma = u/c s 1 is the Mach number. The Mach number is shown to scale like Ma ∼ t/ x (see [START_REF] Latt | Hydrodynamic limit of lattice Boltzmann equations[END_REF] among others). Therefore using the diffusive limit for the scaling of t, t ∼ x 2 , which is the limit consistent with incompressible fluids, one can safely neglect the O(Ma 3 ) term which scales like x 3 .

Therefore the previous equation can be rewritten as

σ P (1) αβ = -2ρτ c 2 s S αβ + (1 -σ )2ρτ c 2 s S FD αβ , P (1) αβ = -2ρτ c 2 s S αβ + 1 -σ σ ρτ c 2 s x 2 6 (∂ 3 α u β + ∂ 3 β u α ). ( 44
)
With this relation, one recovers the equivalent Navier-Stokes equation

ρ ∂ t u α + u β ∂ β u α = -∂ α p + 2∂ β (μS αβ ) - (1 -σ ) x 2 6σ ∂ β (μ(∂ 3 α u β + ∂ 3 β u α )), (45)
where μ = c 2 s ρτ . After performing the standard time-space discretisation using the trapezoidal approximation of Equation [START_REF] Mattila | High-order regularization in lattice-Boltzmann equations[END_REF], one gets the following numerical scheme

fi (x + ξ i , t + 1) = f eq i (ρ, u) + 1 - 1 τ f (1) i , (46)
where fi =

f i + 1 2τ f (1) i , f (1)
i is given by Equation ( 36) and τ = τ + 1/2. The present model is discretised in space using embedded subdomains with uniform mesh with a grid spacing ratio of two between two successive subdomains (see Figure 1). For all the grid points far from the transition region, the collision-propagation algorithm is applied and S FD is computed using centred second-order finite differences schemes. As shown in Figure 1, the coarse grid ( ) overlaps the fine grid (•) making it possible to apply the collision-propagation algorithm on the transition grid points belonging to the coarse grid ( ). Because of a lack of neighbour points, the propagation step cannot be applied at the fine grid points (•) located in the transition layer. On these grid points, the macroscopic values (ρ and u) and the deviatoric stress tensor (P (1) ) are interpolated from the coarse grid neighbours located in the transition layer and S FD is computed using non-centred first-order finite difference schemes.

Associated kinetic energy dissipation, equivalent artificial viscosity and dynamic collision kernel

The extra dissipation of kinetic energy associated to the use of the new hybrid collision kernel introduced above, denoted ε σ , can be easily obtained considering the evolution equation for the kinetic energy 1 2 u • u derived from Equation [START_REF] Silvis | Physical consistency of subgrid-scale model for large-eddy simulation of incompressible turbulent flows[END_REF]. After some algebra, one obtains

ε σ = ν σ |∇ 2 u| 2 , ν σ = 1 -σ 6σ x 2 c 2 s τ , ( 47
)
where ν σ is the σ -dependent artificial hyperviscosity associated to the hybrid RR-BGK collision operator. It is worth noting that the induced dissipation originates in the leading error term of the finite difference scheme used to compute the velocity derivatives. Therefore, the use of a second-order centred scheme leads to a bi-Laplacian-based hyperviscosity, i.e. a ν σ ∇ 4 u correction in the momentum equation. Higher order artificial dissipation operators can be easily defined considering higher order centred finite difference schemes to compute S ij . As a matter of fact, using a pth-order scheme will introduce an hyperviscous operator proportional to ∇ 2p u. An increase in the order of the hyperviscosity is associated to an concentration of the induced dissipation at small resolved scales and to a possible way for smart control of spurious wiggles. The above expression for the induced dissipation ε σ offers a smart way to obtain an ILES scheme by tuning σ in order to recover the same dissipation as the one provided by a subgrid scale model, ε sgs . Considering a generic subgrid viscosity ν t (the same development holds for any RANS eddy viscosity model), one has ε sgs = ν t |∇u| 2 . Therefore, equalising the artificial and the subgrid dissipation leads to

ν t |∇u| 2 = ν σ |∇ 2 u| 2 , ( 48
)
whose solution is

ν σ = L 2 VK ν t , (49)
where L VK = |∇u|/|∇ 2 u| can be interpreted as a renormalised extended definition of the von Kármán length scale [START_REF] Menter | The scale-adaptive simulation method for unsteady turbulent flow predictions. Part 1: theory and model description[END_REF]. The associated value of σ is

σ = 1 6ν t L 2 VK x 2 c 2 s τ + 1 . ( 50
)
These expressions can be further simplified considering subgrid viscosity models that can be written as ν t = c sgs x 2 /τ sgs where τ sgs and c sgs are the subgrid time scale and the subgrid model constant, respectively. The expression for the tuning parameter is

σ = 1 6 c sgs L 2 VK c 2 s τ τ sgs + 1 . ( 51
)
It is worth noting that this last expression does not explicitly involve the mesh size x, leading to an easy implementation in multiresolution algorithms.

The last relation can be inverted to find the value of the Smagorinsky constant that yields the same dissipation as the regularised collision kernel:

c sgs = 1 σ -1 c 2 s τ τ sgs 6L 2 VK . (52)
In the case, the regularised collision kernel is used in the ILES mode, the influence of the implicit subgrid dissipation can be evaluated thanks to the subgrid viscosity parameter s defined as [START_REF] Geurts | A framework for predicting accuracy limitations in large-eddy simulation[END_REF]

s = ε σ ε σ + ε ν , ( 53
)
where ε ν = ν|∇u| 2 denotes the molecular viscosity-induced dissipation. After some algebra, one obtains

1 s = 1 + ε ν ε σ = 1 + ν ν σ L 2 VK = 1 + 6σ ν (1 -σ ) x 2 c 2 s τ L 2 VK . ( 54
)

Vreman's explicit subgrid model

The explicit subgrid model selected in the present paper for both explicit LES and ILES simulations is the one proposed by Vreman [START_REF] Vreman | An eddy-viscosity subgrid-scale model for turbulent shear flow: algebraic theory and applications[END_REF], which is observed to have very interesting self-adaptation features without relying on a test filter-based procedure or a prognostic equation for a subgrid quantity. More specifically, this model is observed to behave in a very satisfactory manner in fully developed turbulent shear flows, and also in transitional flows and near wall region without the use of dynamic procedure, test filter and semiempirical stabilisation step such as averaging or clipping. It is also computationally very efficient, since it does not involve the computation of eigenvalues. A complete analysis of its properties, including symmetry preservation, is available in [START_REF] Silvis | Physical consistency of subgrid-scale model for large-eddy simulation of incompressible turbulent flows[END_REF][START_REF] Trias | Building proper invariants for eddy-viscosity subgrid-scale models[END_REF].

The subgrid tensor is defined as

R ij = u i u j -u i u j , (55)
where the bar symbol denotes the LES filter. It is modelled as

R ij -2ν t S ij + R kk δ ij , ( 56
)
where viscosity is defined as

ν t = c B β α ij α ij , (57)
with

α ij = ∂ ūi ∂x j , β ij = 2 α ij α ij ( 58
)
and

B β = β 11 β 22 -β 2 12 + β 11 β 33 -β 2 13 + β 22 β 33 -β 2 23 . ( 59
)
The cutoff length is taken equal to the mesh size x in the present LBM-based simulations, which are carried out on grids with cubic cells. The model constant is taken equal to c = 2.5C 2 S , where C S = 0.18 is the Smagorinsky constant. In the classical LES mode, this model is implemented by modifying the relaxation time τ in the BGK collision kernel according to

τ = ν + ν t c 2 s + t 2 , ( 60
)
where ν denotes the fluid molecular viscosity.

Application to the flow around a cylinder at Re = 3900

Case description

The new Hybrid RR-BGK collision operator is assessed considering the flow around a cylinder at Re = U 0 D/ν = 3900 and Ma = U 0 /C 0 = 0.0585 with U 0 the inflow velocity, C 0 the speed of sound and D the cylinder diameter. Four different computational cases are defined, in order to assess both LES and ILES capabilities of the method (see Table 1). Case 1 and Case 2 correspond to ILES simulations with constant values of the weighting parameter σ . Case 3 is an explicit LES simulation, in which σ = 0.99, i.e. the total amount of numerical dissipation is expected to remain very small. The last test case (Case 4) corresponds to an ILES simulation with a dynamic evaluation of the parameter σ , Vreman's subgrid model being used to evaluate the targeted dissipation. The computation of σ in that case follows the procedure described in Section 3 for ρ and u for the transition in resolution fine grid points. The computational domain (see Figure 2) is extended 9.5 diameter upstream the cylinder and 49.5 diameter downstream the cylinder. A non-reflective subsonic outflow [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF] is used for the outlet and free slip conditions are imposed in the top and bottom boundaries 14.5 diameter far from the cylinder leading to a blockage ratio of 3.3%. The spanwise extent of the computational domain is taken equal to 4D and is used with periodic boundaries. This spanwise size is known to prevent spurious effects due to unphysical correlations that may be induced by periodic boundary conditions. It is shown in Breuer [48] that no differences are obtained on average quantities using π D or 2π D for the spanwise size.

The computational grid (see Figure 2) is composed of several embedded volumes of uniform cartesian mesh. The grid size is reduced from D/2.5 far from the cylinder to D/80 (finer than in Parnaudeau et al. [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF]) around the cylinder and in its wake. A layer of 15 grid points with x = D/80 is applied around the cylinder and extended 1.5 diameter downstream it (x = 2D) to ensure a good representation of the recirculation bubble. The refinement areas with x = D/40 and x = D/20 are applied on 10 grid point layers around the cylinder and extended 3.5 and 5.5 diameter downstream the cylinder (x = 4D and x = 6D). The distance of the first node to the cylinder is between 0 and 0.0177D with an averaged value of 0.00865D which is lower than in Alkishriwi et al. [START_REF] Alkishriwi | A large-eddy simulation method for low Mach number flows using preconditioning and multigrid[END_REF] and Ouvrard et al. [START_REF] Ouvrard | Classical and variational multiscale LES of the flow around a circular cylinder on unstructured grids[END_REF]. The use of cubic cells leads to 320 grid points in the spanwise direction for the minimal grid size area and a total of 10 million grid points were used in this computation. The computational time was around 47 CPU hours per vortex shedding period on 48 processors. Statistically steady state was reached after 93T * with T * = D/U 0 and all the statistics presented in Section 6.2 were computed on 84T * (≈ 17 vortex shedding periods). the Strouhal number St associated to the main frequency of aerodynamic forces exerted on the cylinder are very accurately recovered in all cases. The predicted values of the rms value of the fluctuations of the lift coefficient C l , the base pressure coefficient Cp b and the normalised recirculation bubble length L r /D exhibit more dispersion, but always remain within the range of variation of previous DNS (Direct Numerical Simulation), LES and DES (Detached Eddy Simulation) results. But it is worth noting that a too large value of L r is found in cases 1, 2 and 3 showing that the transition to turbulence in the separated shear layers present in the formation region is delayed when compared to Parnaudeau's experiments. On the overall, case 4 (ILES with dynamic tuning of the weighting parameter σ ) yields a very good prediction of all parameters, followed by case 2 (ILES with a small constant value of σ ).

Results

Key parameters related to the mean flow are reported in

A deeper insight into the results is obtained considering Figure 3, which displays the streamwise evolution of the mean longitudinal velocity along the symmetry axis in the cylinder wake. It is observed that in all cases the maximum amplitude of the reverse flow in the recirculation bubble is accurately predicted when compared to experimental data. It must also be noticed that the velocity in the far wake is also recovered in all cases, showing that both domain size and boundary conditions are adequately chosen since there is no spurious mass/momentum leaks on upper and lower boundaries and that outflow boundary conditions do not induce unphysical mass/momentum flux. As mentioned above, case 4 leads to an almost perfect agreement with experimental results. Moreover it is observed in Figure 3 that the transition in resolution located at x/D = 2 and x/D = 4 does not include spurious effects on the velocity fields which demonstrates the good implementation of the present method. The topology of the wake is now investigated looking at vertical profiles of the mean velocity at different location in the cylinder wake (see Figures 4 and5). Once again, the simulation based on the dynamic version of the new HRR-BGK collision model (case 4) is in almost perfect agreement with experimental data. An important point is that a V-shaped profile is recovered on the mean longitudinal velocity profile in agreement with experimental data, while many less accurate LES predict a U-shaped profile [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF]. Discrepancies observed in the three other cases are mainly due to the error on the prediction of the length of L r , but the global topology of both the mean recirculation bubble and near wake is recovered in all cases.

Since the size of the recirculation bubble is governed by the transition process in the separated shear layers, it is interesting to analyse the resolved Reynolds stresses. Vertical profiles of the longitudinal and vertical resolved Reynolds stresses at several locations in the wake are displayed in Figures 6 and7. Once again, the very good accuracy of the case 4 is observed. In that case, the shear layer dynamics is very well recovered, since both the maximum value and the shear layer thickness and spreading rate are accurately predicted. In other cases, the separated shear layer spreading rate is underpredicted, leading to the prediction of a too long recirculation bubble. This is further confirmed looking at the streamwise evolution of the resolved longitudinal Reynolds stress u u along the wake symmetry line that is displayed in Figure 8. Case 4 is in very good agreement with experimental data, since both the location and the amplitude of the peak located near the end of the recirculation bubble are satisfactorily predicted. In the three other cases, peaks are damped and shifted downstream, which is coherent with previous comments dealing with the damping of the shear layer dynamics in these simulations.

The effects of the dissipation on the flow structures are illustrated in Figures 9 and10, which display iso-surfaces of instantaneous Q criterion. While it is seen that the flow physics is qualitatively well predicted in all cases (laminar boundary layers along the cylinder with laminar separation and transition in the separated shear layers, large-scale roll-up of the shear layers and appearance of small scale worm-like vortices in the wake), some subtle differences can be detected. As a matter of fact, it can be seen that small-scale structures are mode developed in case 4 and that the separated shear layer transition occurs earlier in that case. The same worm-like vortices are observed in Figure 11 which display the iso-value of σ = 0.97 (Case 1 value) and σ = 0.985 (Case 2 value) for the case 4. The σ parameter tends to decrease down to 0.5 or lower values for a few grid points in areas where dissipation is needed to ensure a good representation of physical phenomena or stays at values close to 1 when no dissipation is needed.

A last quality criterion may be obtained analysing the frequency content of the wake. This is done looking at the frequency spectrum of both streamwise and normal velocity fluctuations in the wake at x/D = 3, see Figures 12 and13, respectively. Power spectra are computed using 3 sequences of 10 vortex shedding cycles with 50% of overlapping using the periodogram technique of Welch [START_REF] Welch | The use of fast fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms[END_REF] and Hanning window. It is seen that the main frequency peak is recovered in all cases, showing that the vortex shedding frequency f vs is accurately predicted. More interestingly, the existence of a secondary peak at a frequency three times larger than the primary peak, i.e. f = 3f vs in the normal velocity spectrum is accurately captured in case 4, showing the very good quality of this simulation [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF]. An inertial range with a slope is observed in all cases, whose width is about 1 decade, as in high-resolution LES and experiments presented in Ref. [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF].

Concluding remarks

A new regularised collision model for Lattice Boltzmann-based LES of weakly compressible flows, referred to as the HRR-BGK collision model, has been presented. It relies on the modification of previously existing recursive regularised BGK models, which consists of hybridising the computation of the velocity gradient with a centred finite difference evaluation when reconstructing the regularised non-equilibrium part of the distribution functions. The resulting effect is the introduction of a stabilising hyperviscosity term, whose amplitude can be explicitly tuned via a control parameter σ . The resulting model can be used in both LES and ILES modes. The new collision model has been assessed considering the flow around a cylinder at Re = 3900. The results obtained with a fixed uniform value of σ exhibit a delay in the transition of the separated shear layers, a phenomenon which is also observed when an explicit subgrid viscosity term is added. On the opposite, very satisfactory results are obtained when using the dynamic version of the HRR-BGK collision kernel equipped with the Vreman's subgrid model. It is important noting that ILES based on the dynamic HRR-BGK is not strictly equivalent to a classical BGK coupled to the Vreman's model, since the scale-byscale distributions of the total dissipation are not equivalent. While the original Vreman model is associated to a classical Laplacian-based dissipation, the present ILES method introduces an bi-Laplacian-based dissipation. The use of higher order dissipation is known to preserve large-scale and inertial-range dynamics [START_REF] Borue | Forced three-dimensional homogeneous turbulence with hyperviscosity[END_REF][START_REF] Borue | Self-similar decay of three-dimensional homogeneous turbulence with hyperviscosity[END_REF], and has been successfully used in several ILES methods for Navier-Stokes equations, e.g. Refs. [START_REF] Berland | Filter shape dependence and effective scale separation in large-eddy simulations based on relaxation filtering[END_REF][START_REF] Daude | A high-order finite-difference algorithm for direct computation of aerodynamic sound[END_REF][START_REF] Hickel | On implicit subgrid-scale modeling in wall-bounded flows[END_REF][START_REF] Hickel | An adaptive local deconvolution method for implicit LES[END_REF][START_REF] Stolz | An approximate deconvolution procedure for large-eddy simulation[END_REF][START_REF] Stolz | An approximate deconvolution model for large-eddy simulation with application to incompressible wall-bounded flows[END_REF][START_REF] Stolz | The approximate deconvolution model for large-eddy simulations of compressible flows and its application to shock-turbulent-boundary-layer interaction[END_REF], with a few existing extensions to the LBM framework [START_REF] Malaspinas | General regularized boundary condition for multi-speed lattice Boltzmann models[END_REF][START_REF] Marié | Adaptive filtering for the lattice Boltzmann method[END_REF][START_REF] Ricot | Lattice Boltzmann method with selective viscosity filter[END_REF][START_REF] Sagaut | Toward advanced subgrid models for Lattice-Boltzmann-based large-eddy simulation: theoretical formulations[END_REF].

The significant differences between the dynamic version of the present ILES method and non-dynamic ones and also the classical explicit LES are due to the fact that most key features of the flow are governed by the transition in separated shear layers. This phenomenon is very sensitive to viscous and hyperviscous damping. A fully turbulent flow is less sensitive to fine details of the dissipative mechanisms, and much smaller differences would certainly occur in such a flow, as observed in ILES results based on Navier-Stokes equations.

Figure 1 .

 1 Figure 1. Visualisation of the grid arrangement around transition in resolution.

Figure 2 .

 2 Figure 2. View of the computational domain and the grid points.

Figure 3 .

 3 Figure 3. Mean streamwise velocity in the wake centreline: Parnaudeau et al. [49] PIV 1 • and PIV 2 *, present case 1 , present case 2 , present case 3 , present case 4 .

Figure 4 .

 4 Figure 4. Mean streamwise velocity at x = 1.06D (top), x = 1.54D (middle) and x = 2.02D (bottom) in the wake of the cylinder: Parnaudeau et al. [49] PIV 1 • and PIV 2 *, present case 1 , present case 2 , present case 3 , present case 4 .

Figure 5 .

 5 Figure 5. Mean normal velocity at x = 1.06D (top), x = 1.54D (middle) and x = 2.02D (bottom) in the wake of the cylinder: Parnaudeau et al. [49] PIV 1 • and PIV 2 *, present case 1 , present case 2 , present case 3 , present case 4 .

Figure 6 .

 6 Figure 6. Vertical profiles of the streamwise resolved Reynolds stress u u at x = 1.06D (top), x = 1.54D (middle) and x = 2.02D (bottom) in the wake of the cylinder: Parnaudeau et al. [49] PIV 1 • and PIV 2 *, present case 1 , present case 2 , present case 3 , present case 4 .

Figure 7 .

 7 Figure 7. Vertical profiles of the normal resolved Reynolds stress w w at x = 1.06D (top), x = 1.54D (middle) and x = 2.02D (bottom) in the wake of the cylinder: Parnaudeau et al. [49] PIV 1 • and PIV 2 *, present case 1 , present case 2 , present case 3 , present case 4 .

Figure 8 .

 8 Figure 8. Variance of the streamwise velocity in the wake centreline: Parnaudeau et al. [49] PIV 1 • and PIV 2 *, present case 1 , present case 2 , present case 3 , present case 4 .

Figure 9 .

 9 Figure 9. Iso contours of normalised Q criterion (Q * = Q(D 2 /U 2 0 ) = 1) coloured by velocity magnitude: Case 1 (top left), Case 2 (top right), Case 3 (bottom left) and Case 4 (bottom right).

Figure 10 .

 10 Figure 10. Iso contours of normalised Q criterion (Q * = Q(D 2 /U 2 0 ) = 100) coloured by velocity magnitude: Case 1 (top left), Case 2 (top right), Case 3 (bottom left) and Case 4 (bottom right).

Figure 11 .

 11 Figure 11. Iso contours of the σ value coloured by the velocity magnitude for σ = 0.97 (left) and σ = 0.985 (right) for the Case 4.

Figure 12 .

 12 Figure 12. Power spectra density of the streamwise velocity at x = 3D: present case 1 , present case 2, present case 3 , present case 4 .

Figure 13 .

 13 Figure 13. Power spectra density of the normal velocity at x = 3D: present case 1 , present case 2 , present case 3 , present case 4 .

Table 1 .

 1 Computation parameters.

	Case	Subgrid model	σ
	Case 1	-	0.97
	Case 2	-	0.985
	Case 3	Vreman	0.99
	Case 4	Vreman	Dynamic

Table 2 .

 2 Table 2 with reference experimental and numerical results. It is observed that the mean drag coefficient C d and Overview of numerical and experiment results.

	Case	Model	C d	C l	St	L r /D	-Cp b
	Parnaudeau et al.	PIV			0.208	1.51
	[49]
	Dong et al. [52]	DNS			0.206-0.210	1-1.18	0.93-1.04
	Ma et al. [53]	DNS	0.84-1.04		0.203-0.219	1-1.59
	Alemi et al. [54]	LES-Smag	0.92-1.01	0.07-0.14	0.205-0.225
	Alkishriwi et al. [50]	LES	1.05		0.217	1.31
	Breuer [48]	LES-Smag	0.969-1.486			0.397-1.686	0.687-1.665
		LES-DynSmag 1.016-1.071			1.197-1.372	0.941-1.011
	Mani et al. [55]	LES	0.99		0.206		0.86
	Franke and Frank	LES	0.978-1.005		0.209	1.34-1.64	0.85-0.94
	[56]
	Kravchenko and	LES	1.04		0.210	1.35	0.94
	Moin [57]
	Lysenko et al. [58]	LES-TKE	0.97	0.09	0.209	1.67	0.91
		LES-Smag	1.18	0.44	0.19	0.9	0.8
	Meyer et al. [59]	LES	1.05-1.07		0.21-0.215	1.18-1.38	0.92-1.05
	Ouvrard et al. [51]	LES-Smag	0.99	0.125	0.218	1.54	0.85
		LES-Vreman	0.92	0.054	0.227	1.83	0.78
		LES-WALE	1.02	0.219	0.221	1.22	0.94
		ILES	0.92	0.052	0.225	1.85	0.77
	Parnaudeau et al.	LES			0.208	1.56
	[49]
	Abrahamsen Prsic	LES	1.0784-1.2365 0.1954-0.4490 0.1956-0.2152	1.27
	et al. [60]
	Wormon et al. [61]	LES-WALE	0.99	0.108	0.21	1.45	0.88
	Zhang et al. [62]	LES	1.001-1.098	0.125-0.345	0.21-0.22
	D'Alessandro et al.	SA DES	1.205-1.2776 0.428-0.6140	0.204-0.215	0.7172-0.85	1.077-1.289
	[63]
		SA IDDES v2 -f DES 0.9857-1.2553 0.1088-0.5719 0.205-0.214 0.7270-1.6780 0.8290-1.2570 1.0235-1.4106 0.1458-0.8283 0.205-0.222 0.5137-1.4270 0.8780-1.4688
	Case 1	ILES	0.936	0.044	0.212	2.05	0.779
	Case 2	ILES	0.973	0.077	0.210	1.835	0.828
	Case 3	LES-Vreman	0.954	0.048	0.209	2.04	0.779
	Case 4	ILES	1.047	0.165	0.212	1.425	0.925

Note: Smag stands for Smagorinsky model, DynSmag for Dynamic Smagorinsky model, TKE for Turbulent Kinetic Energy model, SA for Spalart Allmaras and IDDES for Improved Delayed DES.

Acknowledgments

This work was carried out using the ProLB solver. Professor Eric Lamballais is warmly acknowledged for providing experimental data.

Funding

This work was supported by the French project CLIMB, with the financial support of BPIFrance (Project No. P3543-24000), in the framework of the program 'Investissement d' Avenir: Calcul Intensif et Simulation Numérique'. This work was performed using HPC resources from GENCI-TGCC/CINES (Grant 2017-A0012A07679).

Disclosure statement

No potential conflict of interest was reported by the authors.