Guillaume Sahuquet

Ronny Gueguen

Lilian Fontalvo

Samuel Mer

Adrien Toutant

Françoise Bataille

Gilles Flamant
email: gilles.flamant@promes.cnrs.fr

Particle Flow Heterogeneity in Fluidized-Particles Multitube Solar Receiver

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.

• Molten salt is the commonly used heat transfer fluid in solar towers [1], but has limited operating temperatures (220-565 °C) that limit the efficiency of the power block to 42% [2].

• To overcome this issue and reduce costs, new particle receivers are under development worldwide as the fluidized particle-in-tube of the french National Center for Scientific Research.

• Tests at large scale of the Next-CSP H2020 project [3] with forty 3 m-long tubes (Fig. 1) under concentrated solar flux revealed the existence of heterogeneity in particle mass fluxes between the tubes.

• One of the possible explanation is a difference in the air velocity due to temperature heterogeneity and the existence of an interaction between the tubes. A small cold mock-up was constructed to study the hydrodynamics and behavior of the particle suspension flow at PROMES laboratory.

• Objective: Represent this phenomenom at ambient temperature by injecting different air velocities inside a 2-tubes mock-up.

• Material: olivine with a Sauter diameter of 61 μm (group A of the Geldart Classification)

• Fluidization in a dispenser thanks to a homogeneously distributed air flow rate.

By applying an overpressure, the fluidized bed passes into two 3.70 m-long transparent glass tubes of 4.5 cm internal diameter and is collected onto two weighting devices to measure each particle mass flux G P .

• A secondary air flow rate q ae is applied at an injection nozzle of 1.5 mm diameter at the bottom of each tube. The particle volume fraction and the fluidization regime of the mixture in the tube is directly dependent on the air velocity in the tube [4].

• Analysis method: Particle mass fluxes G P (in kg/(m 2 .s)) of each tube are studied and compared for homogeneous and heterogeneous air flow rates injected q ae (in sm 3 /h).

Conclusion

Fig. 2: a) 3D and b) schematic of the two tubes cold mock-up.

• Different air velocities in fluidized-particles multitube receiver lead to heterogeneous particle mass fluxes.

• An increase of air velocity in one tube increases the total particle mass flux that reduces the pressure in the dispenser with a linear correlation. → Existence of an interaction between the tubes, through the pressure of the dispenser.

• An increase of the total particle mass flux leads to a decrease of the pressure in the dispenser, with a linear correlation as shown in Fig. 4 for all conducted experiments.

• Experimental tests were conducted as follows:

 Same air flow rates q ae for both tubes A and B (Fig. • The other tube is impacted by this diminution of pressure and its particle mass flux decreases.

• This interaction can create a circulation stop in some cases, that has to be taken into account to prevent risks for an industrial application.

Fig. 3 :

 3 Fig.3 : Test 1 -a) Air flow rate (in sm 3 /h) and b) particle mass flux (in kg/(m².s)) for tube A (in blue) and B (in orange) for different time intervals. c) Pressure in the dispenser (in Pa) through time (in s).

Fig. 4 :

 4 Fig.4 : Decrease of the pressure in the dispenser (in Pa) as a function of the total particle mass flux (in kg/(m².s)).

  3.a) until the mixture flows (intervals 1-3). → Same particle mass fluxes G P in the tubes (Fig.3.b). → The pressure in the dispenser decreases (Fig.3.c).  Then only the air flow rate of tube A is increased (intervals 4 and 5) → G P A increases while the pressure in the dispenser decreases and thus G

Fig. 1 :

 1 Fig.1: a) Thémis solar tower (Targasonne, France) with b) the Next-CSP fluidized bed 40-tubes solar receiver.


Ali Ghaddar
email: ali.ghaddar@bitasoft.com

Dalila Tamzalit
email: dalila.tamzalit@univ-nantes.fr

Ali Assaf
email: ali.assaf@bitasoft.com

Abdalla Bitar
email: abdalla.bitar@bitasoft.com

Variability as a Service : Outsourcing Variability Management in Multi-tenant SaaS Applications

Keywords: SaaS, multi-tenant, variability

come   L'archive ouverte pluridisciplinaire

Introduction

In recent years, the tendency towards outsourcing IT resources that are not the key competencies of an enterprise has caused the appearance of new services providers type. These new services providers develop and maintain, on their infrastructures, such outsourceable IT resources while offering their access to enterprises through the web. The business model of such resources offering aims at providing a pay-as-you-go payment model for their use, where these resources can be provisioned and un-provisioned on demand. Such new outsourcing principle is gaining wide acceptance, especially in the small and medium enterprises (SME) segment. The reason for this acceptance is purely economic. Indeed, the majority of SME consider taking less risks by contracting the implementation of parts of their IT resources to a more experienced provider and find this option more cost effective. The term Cloud Computing summarises these outsourcing efforts and provisioning of services on demand. This term has emerged as the run-time platform to realise this vision [START_REF] Kabbedijk | Variability in multi-tenant environments: architectural design patterns from industry[END_REF][START_REF] Sengupta | Engineering multi-tenant software-as-aservice systems[END_REF]. More concretely, Cloud Computing is viewed as a stack containing these new types of services to be provided. That is to provide an Infrastructure as a Service (IaaS), also to provide a development Platform as a Service (PaaS), to finally provide a Software as a Service (SaaS). In the following, we will focus on this last service type.

SaaS can be defined as: "Software deployed as a hosted service and accessed over the Internet" [START_REF] Chong | Architecture strategies for catching the long tail[END_REF]. In the model of this service, the provider offers a complete application, ready to be used by prospective customers. These customers subscribe to the application and use it via the web through a simple browser. Full details on the adopted implementation technology and the application deployment server are transparent. This type of solution is offered by SaaS providers as an alternative to the traditional software applications that require installation on the client side, as well as significant investments in terms of materials and physical resources. However, from a SaaS provider perspectives, the real benefits of SaaS begins when it is possible to host multiple customers on the same application instance, without the need of a dedicated application to be deployed and separately maintained for each customer. This approach is called multi-tenant, where tenant means a customer organisational structure, grouping certain number of users [START_REF] Bezemer | Enabling multi-tenancy: An industrial experience report[END_REF]. In fact, application tenants could be distributed across different geographic zones and belong to different industry verticals, and thereby, their requirements from the application may differ, increasing the need to support these specific requirements in addition to common ones. This inevitably generates an additional and important variability [START_REF] Svahnberg | A taxonomy of variability realization techniques[END_REF] management concern, which obviously distract SaaS providers focus from their core competences.

According to our experience, variability management in multi-tenant applications generally implies: (1) modeling such variability (explicitly documenting the variable features and the possible alternatives), (2) binding the variability model and storing each tenant customisations and choice of alternatives, (3) adapting the appropriate application behavior at run-time for a tenant, according to its stored customisations. This last key concern is tied to the application's architecture. Our work is conducted in a SOA context where the targeted multi-tenant applications are designed and developed following the service oriented architecture (SOA) model. So, such application adaptation usually concern its look and business-processes (services composition an orchestration). As these variability requirements may be distracting and disturbing issues for SaaS providers, we are convinced that they will gain to outsource, if their is a possibility, such variability management concern to a specific service provider. For this purpose, we introduce in this paper the concept of Variability as a Service (VaaS) model (in the same spirit of IaaS and PaaS models in terms of outsourced resources). The VaaS model implies the necessity of a new type of service providers: the VaaS providers. A VaaS provider that follows the VaaS model we propose, will permit to its customers (i.e. multi-tenant SaaS providers) to model and to resolve their multi-tenant applications variability on his infrastructure. We will explain our approach and its context in the following. The remainder of the paper is structured as follows.

Section 2 presents the needs of modeling variability in terms of concepts for modeling and managing; it mainly identifies the limitations that the variability modeling techniques from software product line engineering encounters to realise our objectives from variability outsourcing. Section 3 presents our suggested architecture, process and meta-model for the variability outsourcing and the VaaS model. Section 4 presents a case study from the food industry that we have revisit in order to apply the new approach. In Section 5, discussions and future works are detailed. Section 6 shows related works and Section 7 concludes the paper.

Needs of Variability Modeling

When many variable features in a software system are identified, variability modeling will become an important requirement to express these features. Generally, variability modeling consists on explicitly documenting, in the software system, what does vary and how does it vary, allowing thus to improve traceability of the software variability, as well as to improve variability communication with the customers [START_REF] Pohl | Software product line engineering: foundations, principles, and techniques[END_REF]. However, since software systems are composed from different development artifacts, a variability concern may impact and cut across these different artifacts (i.e. in SOA-based systems, the same variability in the business-process may have effects on the GUI and may also requires the development and deployment of additional services). Modeling variability in each artifact separately and differently is a very bad idea because it makes us loose all the sense of variability tracking, and it makes representing dependencies between variability in these different artifacts very hard. In addition, the different way of modeling variability in each development artifact will complicates the reasoning about variability in the entire system. Therefore, in order to avoid these problems and to ensure strong expressiveness and management of variability, it is necessary to guarantee its representation in a separated, central and uniform model, while maintaining its relation with impacted artifacts.

Orthogonal Variability Models

Based on the need of tracing dependencies between variability in different development artifacts and to maintain the uniformity in variability definition and modeling, orthogonal variability models such as the OVM have been introduced in the software product line engineering (SPLE) [START_REF] Bayer | Consolidated product line variability modeling[END_REF][START_REF] Pohl | Software product line engineering: foundations, principles, and techniques[END_REF] domain. These orthogonal variability models represent variability in one central and separated model, and then link each variability information to one or many development artifacts. In [START_REF] Ghaddar | Decoupling variability management in multi-tenant saas applications[END_REF], we have shown how OVM can be used to model variability in multi-tenant SOA-based applications. One of our main goals was to show that concepts for variability modeling from the SPLE can be used in a multi-tenant SaaS context. These OVM modeling concepts are resumed in the following:

-Variation point (VP) and variant: A VP represents one or many locations in the software where a variation will occur, indicating the existence of different alternatives, each of them may result in a different software behavior. A Variant represents an alternative of a VP. -Artifact dependency: this concept relates the defined variability in the variability model to different development artifacts, indicating which artifacts are representing VPs and which ones are realising variants. -Constraints dependency: constraints in OVM can be operated between variants, between VPs and between variants and VPs: an exclude constraint specifies a mutual exclusion; e.g., if variant1 at VP1 excludes variant2 at VP2, the variant2 can not be used at VP2 if variant1 is used at VP1. A requires constraint specifies an implication; i.e. if a variant is used, another variant has to be used as well. -Internal and external variability: these two concepts separate between the variability that is only visible to the developers (internal ), and the variability that is communicated to the customers (external ).

However, according to our objectives from variability outsourcing, the OVM is not well-suited. In fact, OVM encounters certain limits that prevent its adoption as a variability modeling solution on a VaaS provider infrastructure. These limits result from the differences between the SPLE domain and the multi-tenant SaaS model. In SPLE, the software adaptation to customers requirements is realised through a customisation task that is usually done by the same organisation that has built the software. While in the multi-tenant SaaS model, since all tenants are using the same software instance, it is more practical to give these tenants the control to customise the software to their needs (interventions from the SaaS provider still always required to guide the tenants through the customisation of the application).

Enabling the Software Customisation by Tenants

The direct customisation of the software by tenants is one of our objectives since it relieves, in many cases, the SaaS providers from doing this heavy customisation tasks each time a new tenant subscribe to the application. However, such customisation capabilities implies that SaaS providers restrict, for certain tenants, the possibility to choose certain variants which, for example, cannot be selected in the application test period (tenants usually pass through an application test period to determine if the software can provide the service they expect). In addition, some variants may be highly expensive and they must be restricted for certain tenants (in use period) with limited budget. The possibility to define, depending on the tenants contexts, the required restriction conditions for variants and VPs is not the focus of OVM, as it assume that the software customisation will be always done by the software creator. In addition, another main reason to propose variability as a service is minimising the number of interventions that SaaS providers have to make for managing variability, and especially, providing concrete data values for variants. Therefore, our variability modeling solution gives SaaS providers the possibility to define variants with free values, thus tenants can provide the informations they like, such as their own application logo, title, background color etc., without having to always use those offered and predefined by the SaaS providers. Once again, this option is not supported by the OVM. Finally, since the application will be customised directly by individual tenants, it will be important to specify, in the variability model, the order in which the variation points must be bound (Some VPs may depend on each other, which is not always intuitive for tenants). In OVM, such ordering could be done by exploiting the constraints dependencies between VPs, but since constraints in OVM were not designed for this purpose, using them to order the VPs binding will be a workaround. We suggest in the next section, a new variability meta-model that covers the limitations of OVM and others SPLE variability modeling techniques, by bringing additional variability modeling concepts, allowing the SaaS providers to give their tenants certain customisation controls and capabilities in a simple, secure and flexible manner.

3 Variability as a Service : Architecture, Process and Meta-model

In this section we will explain in details the VaaS model and its core element: the VariaS component. The VariaS component provides different interfaces for SaaS providers to create and manage their applications variability models, as well as to resolve variation places in the applications. The VariaS component is deployed and maintained on a VaaS provider infrastructure. The high-level architecture and different actors of VaaS are presented in figure 1.

Architecture and Variability Outsourcing Process

As depicted in figure 1, the VaaS architecture defines certain steps that must be performed by SaaS providers as well as by tenants in the variability outsourcing process. These steps are ordered in time and they are divided under Specification and Execution steps. All mentioned steps ((step 1), (step 2)etc.) and elements ((a), (b), etc.) are shown in the figure.

Specification. First of all, the SaaS providers have to instantiate the variability meta-model (a) offered by the VaaS provider (step 1). The variability model resulting from such instantiation is stored in a variability models repository (b).

Prospective tenants of the SaaS application (having outsourced its variability management) must bind its variability model (choosing the variants that meet their needs), before they can start its use (step 2). Such variability binding has to be made through the SaaS provider application and then redirected to the VariaS component, since it is more easy and safe that tenants continue to communicate only with their SaaS providers. In addition, such binding (through an integrated variability binding tool) has to be controlled thus tenants will not be able to select non authorised variants and will not have access to certain VPs depending on their contexts and details (h). However, after binding the application variability model, tenants variants choices are saved as customisation documents (c) (step 3). Each customisation document concern one tenant, it indicates their variants chosen at the variation points. Having done so, the tenants can start using their SaaS application and expect from it to behave as they have customised it. Therefore, the application behavior must be adapted at run-time according to these customisations.

Execution. In order to adapt the multi-tenant application behavior according to its defined variability model and the tenants customisations, the SaaS developers must request the VariaS component at each variation place they have identified in the application. Each variation place must have a corresponding variation point in the application variability model that is stored on the VaaS provider side, thus the VariaS component can identify such variation (step 4). These variability resolution requests are in specific forms and use specific protocols defined by the VaaS provider. When a request is received, it will be handled by a variability resolver (d). Such resolver accesses the stored variability models as well as the tenants customisation documents in order to identify the variability action that must be performed (step 5). This action will result in a specific value which is either equivalent to the variant value (selected by the tenant) or deduced from it in case that additional computations based on the variant value must be performed. The variability resolution results will be returned to the application that sent the resolve variation request(step 6). These results must always have the same data format expected by the SaaS application at the corresponding variation place. The non respect of this format may break down the application. In this case, the additional computations to perform on a variant value would be, for example, to transform its data format to the one expected by the application. The expected data format from the resolution results of each variation point must be defined in the variability model.

Variability Meta-model

A variability meta-model is essential for the VaaS model, since it defines how SaaS providers describe their applications variable parts, and how they indicate restricted variants, those accessible to different customisations possibilities. In addition, the variability meta-model has to provide a technical and code level solution for SaaS providers, to help them resolving their application variations while still independent from these applications implementation details. Thus, the VaaS model, through such a meta-model, deals with applications based on standard technologies, such as web-services, and use common applications code artifacts such as texts, files and expressions. In section 2, we have discussed the need of an orthogonal variability model, and shown the limitations that the SPLE modeling techniques encounters for realising our variability outsourcing vision. In the following, new variability modeling concepts related to the VaaS model as well as reused concepts from OVM will be presented. Figure 2 shows our variability meta-model.

The variability modeling approach we present, aims to provide certain level of flexibility to SaaS providers, making them able to deal with complex variability modeling situations in a simple manner. Complex modeling situations such as constraining (by restricting or obligating) the choice of certain variants as well as the binding of certain VPs depending on the tenants contexts are needed. Being able to do so, the Activation Condition concept is introduced. Such concept allows SaaS providers to define conditions, which if they evaluate to true, their related Constraints will be activated. In this way, and depending on the conditions evaluation results, the same variant may be restricted for some tenants and enabled for others. A condition semantic is expressed by a SaaS provider the same as in If-Then-Else statements in programing languages. These conditions evaluate on informations such as tenant identity, status, already chosen variants and other informations may be useful to limit or extend the software usability for tenants depending on the SaaS provider commercial offer.

The variation point, variant and development artifact concepts are reused from OVM. As mentioned above, variability informations (variants and VPs) are represented in an abstract and artifacts-independent manner, while their relations with the impacted development artifacts is maintained through so called Artifact Dependency. However, for each variant (realised by one or many artifacts) its important to indicate its corresponding value in each artifact, thus such value can be treated by our service and returned to the SaaS application when needed. The same variant can have different values in different artifacts, each value must have a type, thus our variability service can determine the variability action to perform on resolution requests, as well as the appropriate GUI in which the SaaS provider can provide these values. For example, if a SaaS provider creates a value instance for a particular variant with the Web-Service type, automatically, a relevant GUI will be opened to provide the web-service description, end point address, operation to invoke, input and output massages types etc. The same, if a value instance is created with a File type, the GUI opened would be a simple file input. These is beneficial for SaaS providers, as free values can be added by tenants through dedicated and simple GUI's. Finally, the Role concept has been added and associated to a variation point in order to differentiate between VPs to bound by the tenants and others to bound by the SaaS providers. In addition, such role concept is beneficial in case one or many application resellers exist. In fact, the SaaS provider may find more efficient that application resellers make decisions on certain VPs that are related to a given segment specificities which some prospective tenants are parts of. Generally, such resellers are more experienced in these segments standards and rules. In the following, we present a case study that shows an instantiation example of these variability concepts.

Case Study: a Food Industry Application

In this section, we will revisit a food-industry application (FIA for short) that we have developed as presented in [START_REF] Ghaddar | Decoupling variability management in multi-tenant saas applications[END_REF]. FIA is a SOA-based multi-tenant application. It allows its food industry tenants looking for foods production and quality improvement to predict the future expenses and benefits from their potential recipes, before executing a real and expensive manufacturing process. FIA evaluates, by relying on an internally developed simulation service, the recipes manufacturing time and cost, as well as their quality characteristics such as nutritional value, taste and smell. FIA has also an integrated shipping service for foods, from the warehouse of the foods supplier to the tenants manufactures. Figure 3 left side shows the FIA business process.

The back-end business logic of the application is implemented using a Business Process Execution Language (BPEL) [START_REF] Jordan | Web services business process execution language version 2.0[END_REF] based solution, and the front-end for the tenants is a Web application which gathers data and hands it over to the BPEL engine for processing the recipes simulation requests. When a tenant user accesses the application, he can manage his existing recipes or decide to create a new one. In this last case, the user has to provide the foods composing the new recipe as well as their respective percentages in it. The user must also describe the recipe cooking process, by sending a pre-negotiated XML structure, or he can rely on a FIA integrated cooking process drawing tool. When finishing, the user must click on simulate button and wait for results. When the application receives the recipe details, the BPEL process invoke the foods supplier service, asking for foods prices and quality informations. These informations, as well as the recipe cooking process allows the simulation service to calculate exactly the final recipe cost and quality. However, when the simulation ends the process sends back a simulation report to the user containing the recipe simulation results. In case user validate the results, the process saves the recipe details in the database and invoke the shipping service thus the foods (composing the recipe) will be shipped from the warehouse of the foods supplier to the tenant manufacture. In order to detect possible variations in the application we have presented its business process to some selected prospective tenants. Three variations have been detected: 1) some tenants are interested by the software but they have asked for another foods supplier with higher foods quality. 2) Some other tenants asked for another foods shipper because the existing one is costly. Taking into account this shipping variation, we have decided to add an economic shipping service to the application but with a higher shipping time. The costly shipping service is very rapid, but it does not supports international shipments. On the other side, we have found that high quality foods suppliers have their own shipping services and that will excludes the need of the shippers proposed by our application. 3) Finally, some tenants want their own company-specific logos and specific titles for the application. Beside these tenant-oriented variations, a segment variation is added to indicate the different simulation models that the application supports. Currently, we are supporting the French and the American models which are the most used worldwide. However, these FIA variations was modeled in first time by relying on the OVM model, and the application customisation was always done by a member of our development team. In figure 4, the FIA variations are modeled following the VaaS approach.

Resolving FIA Variations

In this section we only focus on the business-process variability, as it requires additional computations to perform by the Varias component in order to resolve each service in the process there may be different alternatives which implement the service with different implementation logics or quality attributes. In this case, the variability occurs on selecting the most appropriate service depending on the tenants requirements. If a SaaS provider defines these different services alternatives as variants at a particular variation point, the VariaS component will be responsible of invoking the appropriate service (being able to invoke the appropriate service variant requires providing the informations needed for its invocation, such as end point address, input and output massages, operations etc.) and compensating the data mismatch at the input and output messages of different services, if exists. The services invoker (e) and the data transformers (f) elements in the VariaS component (see figure 1) are designed for this purpose.

As depicted in figure 3, the two variation points at different services invocations are resolved by sending requests to the Varias component. These requests are in a fixed form which abstract for SaaS providers the concrete services addresses. In addition, such abstraction is also beneficial for SaaS providers since they always interact with variable services in the same manner and expect the same output massages format.

Discussion and Future Work

In this paper, we presented the first stages of our work: VaaS or outsourcing Variability as a Service with its associated architecture and meta-model. We also presented the VaaS model in its specification and execution steps (see section 3.1) consider SOA applications. We see them as a set of several artifacts, possibly variable. It is important to outline two main hypothesis we voluntarily considered as implicit in this paper:

-Specification: a first implicit hypothesis of our work is that variability models of an application and its business design are defined at quite the same time.

So the application and its variability models are both defined from scratch and the variability concern is considered at the first stages of the design of the application. This is what we called early variability in [START_REF] Ghaddar | Decoupling variability management in multi-tenant saas applications[END_REF]. -Execution: a second implicit hypothesis of our work in its current stage is that the executed application has a stable design and a stable variability. This hypothesis implies that only variants need to be detected, uploaded, executed and potentially evolved. More precisely, in the execution step, we consider only source-code variable artifacts, those executed at run-time. The FIA application is an illustration example through its executable artifacts (i.e. Business-process, web-services, GUI etc.). We also consider that (ii) the application design and its variability models, essentially the View Points, are stable.

At this stage, we aim to stabilise the VaaS concept and VariaS component before exploring the following important future issues:

-Validation: we successfully developed multi-tenant applications by including and managing variability through VP and variants. This development

Conclusion

In this work, we have shown the importance and the complexity of supporting variability in SOA-based multi-tenant applications. We have also motivated the need of outsourcing the variability management to a VaaS provider. Architecture and meta-model supporting our approach have been also provided. In addition, we have revisited a multi-tenant application case study and transforming its existing OVM variability model into VaaS model. Our approach aims to decrease the variability management complexity, and to relieve the SaaS providers looking forward to adopt a multi-tenant solution, from developing an expensive variability solution beforehand.

Fig. 1 .

 1 Fig. 1. High-level VaaS Architecture

Fig. 2 .

 2 Fig. 2. Variability Meta-model of VaaS

Fig. 3 .

 3 Fig. 3. FIA Business-process and its Interaction with the Varias Component

them. In fact, as previously mentioned, we consider multi-tenant applications as SOA-based. In SOA, the application business-process is realised by composing and orchestrating available web-services achieving the desired functionality. For has been successfully tested and achieved. The next validation step is to implement the VaaS architecture and its VariaS component. -Non-code artifacts: rather than supporting only code , we aim to generalise the approach to other artifacts like requirements, design and test. -VaaS and existing applications: the objective is to enable the VaaS architecture able to support variability as well for new applications as for existing ones. We thus have to look how one can re-engineer an existing application in order to inject variability in terms of specification and management. -Evolution: applications are in essence evolutive. Managing evolution in a VaaS approach is our ultimate objective. It can be seen through three main concerns: (i) variability evolution, (ii) co-evolution of the application and its variability models and (iii) the co-evolution of the variability meta-model and its variability models.

6 Related Work

SOA Variability Management

Several authors studied the variability concern in the context of SOA systems.

In [START_REF] Chang | A variability modeling method for adaptable services in service-oriented computing[END_REF] authors identify four types of variability which may occur on SOA. In [START_REF] Sun | Modeling and managing the variability of web service-based systems[END_REF] authors present a framework and related tool suite for modeling and managing the variability of Web service-based systems. They have extended the COV-AMOF framework for the variability management of software product families. In [START_REF] Jiang | Pattern-based variability management in web service development[END_REF] authors describe an approach to handle variability in Web services, while [START_REF] Koning | Vxbpel: Supporting variability for web services in bpel[END_REF] focus on variability in business-process by proposing a VxBPEL language. However, all these approaches focus on variability management in SOA systems as a part of the software providers responsibilities which is different from our VaaS and variability management outsourcing approach.

SaaS and Multi-tenancy

In [START_REF] Jansen | Customization realization in multitenant web applications: case studies from the library sector[END_REF] authors provide a catalog for customisation techniques that can guide developers when dealing with multi-tenancy, they have identified two types of customisation: Model View Controller (MVC) customisation and system customisation. In [START_REF] Bezemer | Multi-tenant saas applications: maintenance dream or nightmare?[END_REF] the authors propose some architectural choices to make when building multi-tenant applications, while in [START_REF] Bezemer | Enabling multi-tenancy: An industrial experience report[END_REF] authors discuss their experiences with re-engineering an existing industrial single-tenant application into a multitenant one. In [START_REF] Mietzner | Variability modeling to support customization and deployment of multi-tenant-aware software as a service applications[END_REF], authors propose a variability modeling technique for SOAbased multi-tenant applications, they differentiate between internal variability only visible to the developers, and external variability that is communicated to the tenants of the application. Once again, they do not propose the variability as a service which we have treated in this paper.


Kevin Lachin
email: kevin.lachin@agroparistech.fr

M Niane

M Person

J Mazet

G Delaplace

C Turchiuli

Prediction of droplets characteristic diameters and polydispersity index induced by a bifluid spraying nozzle by the means of dimensional analysis

Keywords: Dimensional analysis, bifluid nozzle, droplet size distribution modeling, atomization, machine learning modeling

A rigorous dimensional analysis (DA) was performed on bifluid spray properties • Two characteristic spray diameters and a polydispersity index were considered • Models were established with conventional equation shapes and machine-learning (ML)

• Spray properties could be accurately modelled through a DA/ML approach • Particle Size Distributions could also be predicted using this methodology

Introduction

Atomization remains essential in chemical engineering and is frequently used in many domains, functioning as a key step in processes such as, for example, fuel combustion [START_REF] Faeth | Current status of droplet and liquid combustion[END_REF] and probiotic encapsulation [START_REF] Guerin | Lactobacillus rhamnosus GG encapsulation by spray-drying: Milk proteins clotting control to produce innovative matrices[END_REF]. Atomization corresponds to the disintegration of a liquid in a surrounding gas that can be still or in motion. This disintegration occurs through a nozzle and generates a spray. This spray can then be used in its dispersed liquid form or transformed (frozen or dried) to obtain the product of interest. The dried usage is applied in spray-drying, where atomization is the first step in the formation of solid particles (e.g., for milk powder production).

Among the diversity of nozzles that can be employed, two categories are commonly used in industrial applications: pressure nozzles (where the fragmentation results from a kinetic energy dissipation inside the slots [START_REF] Lachin | Dimensional analysis modeling of spraying operation -Impact of fluid properties and pressure nozzle geometric parameters on the pressure-flow rate relationship[END_REF]) and bifluid nozzles (where atomization is air-assisted and comes from the shearing occurring at the air-liquid interface [START_REF] Petit | A dimensional analysis approach for modelling the size of droplets formed by bi-fluid atomisation[END_REF]). With both types of nozzles, two spray properties are particularly interesting in atomization: a characteristic diameter and the size distribution (PSD) of generated droplets. While the knowledge of a diameter gives general information about the spray quality (coarse or fine spray), the droplet size distribution itself is also of full relevance, especially in the case of spray-drying. The spray PSD will indeed directly impact that of the resulting powder particles as well as their behavior during drying with an impact on the final powder properties (flowability, packing properties, wettability) which are highly dependent on the PSD [START_REF] Suzuki | Effect of size distribution on tapping properties of fine powder[END_REF]. The prediction of these two spray properties is thus of great interest for the experimenter working with various fluids to efficiently pilot the atomization process to the needs.

Concerning characteristic diameter, "ready-to-use" correlations predicting diameters depending on fluid and process parameters flourish in the literature [START_REF] Petit | A dimensional analysis approach for modelling the size of droplets formed by bi-fluid atomisation[END_REF][START_REF] Mandato | Liquids' atomization with two different nozzles: Modeling of the effects of some processing and formulation conditions by dimensional analysis[END_REF][START_REF] Ballester | Drop size measurements in heavy oil sprays from pressure-swirl nozzles[END_REF][START_REF] Hewitt | Droplet size spectra produced by air-assisted atomizers[END_REF][START_REF] Chu | An experimental investigation of swirl atomizer sprays[END_REF]. While of great interest and use in the parameter range in which they are established, these correlations however fail at representing the complexity of a size distribution. Generally speaking, when determining a PSD through modeling, Population Balance Equations (PBE) are often invoked to predict time-dependent and steady-state distributions of the dispersed medium [START_REF] Somasundaran | Modeling flocculation of colloidal mineral suspensions using population balances[END_REF][START_REF] Zhang | CFD simulation of gas-liquid flow in a highpressure bubble column with a modified population balance model[END_REF]. While robust and efficient, PBE requires both kinetic information (aggregation/coalescence, rupture kernels) and computational power (especially when coupled with CFD) making them difficult to use effectively. Simple correlations capable of estimating a PSD, as it is the case for characteristic diameter, could thus be of great interest.

For such complex phenomena, dimensional analysis is an interesting tool to consider. Dimensional analysis is a methodology of prime importance in chemical engineering as it offers both a means to semi-empirically model complex phenomena for which full numerical approaches are too complex and time-consuming, and a way to qualitatively understand the prevalence of the phenomena at stake. Among some of the applications for which this approach is of high relevance, one can mention mass transfer in stirred tanks [START_REF] Hassan | A consistent dimensional analysis of gas-liquid mass transfer in an aerated stirred tank containing purely viscous fluids with shear-thinning properties[END_REF], mixing inside a novel mixing device [START_REF] Delaplace | Homogenization of liquids inside a new soft elastic reactor: Revealing mixing behavior through dimensional analysis[END_REF], or size prediction when atomizing [START_REF] Petit | A dimensional analysis approach for modelling the size of droplets formed by bi-fluid atomisation[END_REF][START_REF] Mandato | Liquids' atomization with two different nozzles: Modeling of the effects of some processing and formulation conditions by dimensional analysis[END_REF]. Recent developments in this field also proposed a method to combine design of experiments with dimensional analysis to optimize the number of experiments to be performed, with satisfactory preliminary results [START_REF] Berth | Modelling and optimisation of gasliquid mass transfer in a microporous hollow fiber membrane aerated bioreactor used to produce surfactin[END_REF].

This article focuses on bifluid nozzle atomization, and specifically on establishing new correlations predicting characteristic diameters for the experimental bifluid setup of this study. These kinds of correlations already exist in the literature and use monomial/exponential terms. Unfortunately, they often lack precision due to the shape of the equation used. We will thus overcome this obstacle by proposing a comprehensive dimensional analysis study combined with a symbolic regression software for the identification of the best fitting equation shape. The results of both approaches (monomial/exponential and symbolic regression) will be compared and discussed extensively. A methodology to quickly estimate the PSD of a spray depending on the fluid properties and operating conditions will then be proposed. This methodology relies on the liquid fragmentation theory proposed by Villermaux [15]. This theory leads to a two-parameters mathematical function describing the spray size distribution (a characteristic diameter and a polydispersity index). Dimensional analysis will be used to predict the parameters of this PSD function. The article concludes by discussing the robustness and utility of the PSD prediction.

Materials and methods

Sprayed liquids

Five mixtures of water-glycerol-ethanol were formulated based on the spraying experiments performed by Mandato et al. [START_REF] Mandato | Liquids' atomization with two different nozzles: Modeling of the effects of some processing and formulation conditions by dimensional analysis[END_REF] in order to get five model liquids with a large variety of viscosity, surface tension and density (Table 1). Glycerol (CAS Number 56-81-5, purity ≥ 97%) was purchased from VWR Chemicals (Belgium) and ethanol (CAS number 64-17-5, purity = 96%) from Sigma Aldrich (USA). Type II deionized water was also used for the mixtures.

Glycerol was mainly used in this study to change the dynamic viscosity of the solutions. Ethanol was also added in solutions 4 and 5 in order to tune the surface tension while keeping dynamic viscosities almost identical to those of solutions 2 and 3. These formulations were also chosen to ensure adequate liquid flowrate when pumping and user safety (relative to the spraying of ethanol) when performing the experiments. The results obtained with the five solutions were afterwards used to feed the dimensionless models proposed and identify the correlations (training data set). The robustness of the models was afterwards tested using a 30% w./w. reconstituted skimmed milk solution (validation data set). The milk powder used for reconstitution was bought from a local store (Cora, France). Skimmed milk is a product of particular relevance for this study, as it exhibits a Newtonian behavior [START_REF] Petit | A dimensional analysis approach for modelling the size of droplets formed by bi-fluid atomisation[END_REF] while being widely used industrially to produce milk powder with spraying nozzles. The viscosity, density, and surface tension of the liquids used (Table 1) were estimated as follows. The dynamic viscosity and Newtonian behavior of the solutions was experimentally assessed on a Anton Paar MCR 301 (Anton Paar, Austria) rheometer by varying the shear rate between 1 and 200 s -1 with a logarithmic ramp of duration 320 s inside a Couette geometry. The measurements were performed at a temperature of 20°C (spraying experiments temperature) in triplicate. Liquid density ρl was also determined three times by measuring the mass corresponding to 10 mL of solution at 20°C with a scale of precision +/-0.001 g. Ultimately, surface tension σ at equilibrium was measured in triplicate on a Teclis tensiometer (Teclis Scientific, France) using the rising drop method. The maximal coefficients of variation assessed were respectively 1.4%, 0.16%, and 0.82% for the dynamic viscosity, density, and surface tension measurements.

Spraying set-up and process conditions explored

Bifluid nozzle and related geometrical parameters

The nozzle used for the experiments was a co-current bifluid nozzle with external mixing (Minor Mobile bifluid nozzle, GEA, Germany), as presented in Figure 1. Characteristic diameters of the nozzle (dl, da,int, da,ext) were measured using the software Image J on nozzle cap snapshots taken with an optical microscope. These measurements were done three times, and the average measured values are reported in Table 2.

Pilot set-up and process operating conditions

The pilot-scale set-up used for the experiments is detailed in Figure 2: Spraying and image acquisition set-up. The sprayed solutions were firstly mechanically stirred inside a double-jacketed tank of volume V = 5 L thermally regulated at T = 20°C. The liquid flowrate inside the nozzle was set by a peristaltic pump (Masterflex L/S -EasyLoad head pump, Cole-Parmer, USA), and the mass liquid flowrate was calculated using a scale recording the mass variation through time when spraying. Air volumetric flowrate was imposed with a manual gas flowmeter (FR4500 Gas flowmeter, Key Instruments, USA) and air properties (temperature and humidity) were systematically controlled before spraying. The corresponding air density value used for the calculation was ρa = 1.204 kg.m -3 .

A vacuum system was also set up to avoid as much as possible spray visualization issues (mist accumulation inside the experiment room, droplet deposal on the camera, spurious recirculation of drops in the camera field of observation). Tests were performed before the experimental campaign and showed that this system does not disturbs the liquid flow at the measurement position.

Table 3 summarizes the process conditions explored for each of the six solutions during the experimental campaign. Three air flow rates were investigated per solution, and five liquid flowrates were swept per air flowrate considered. In total, 15 experimental points per fluid were thus obtained, leading to 75 experimental points for the training data set and 15 for the validation data set (reconstituted skimmed milk). The maximal flowrate achievable with a peristaltic pump is highly influenced by the liquid dynamic viscosity. For the solutions with the highest viscosity (solutions 3 and 5), the mass flowrate range explored was thus of a lesser magnitude than for the other four solutions.

Image acquisition and exploitation 2.3.1 Image acquisition

Images of the sprayed droplets were recorded using a laser image acquisition system designed by R&D Vision (R&D Vision, France), allowing to capture images of a maximum resolution 3296*2472 pixels with a maximal 16 Hz frequency. The acquisition camera (VA-8M, Vieworks, South Korea) is synchronized with a pulsed laser of wavelength 640 nm to illuminate the zone to photography. The laser and camera were mounted on a metallic frame whose position can be adjusted on the X, Y, and Z axes through micromotors (A-LST, Zaber Technologies, Canada), ensuring a precise and reproducible camera position.

The acquisition was performed at a distance L = 15 mm from the nozzle outlet. At this position, steadystate was experimentally observed after a few seconds. The acquisition was thus triggered after spraying for 30 seconds (to ensure steady-state) and made for another 30 seconds. Approximately 400 processed images per operating condition were recorded. Figure 3 shows examples of images obtained when spraying solution 2 in two operating conditions: air volumetric flowrate = 90 L.min -1 -liquid mass flowrate = 50 g.min -1 (left) and air volumetric flowrate = 90 L.min -1 -liquid mass flowrate = 90 g.min -1 (right).

Data treatment and statistical analysis

Each image obtained was then processed using the software HIRIS 5.1.4 (R&D Vision, France) to determine the diameter of each droplet in the observation field. Drops lying outside the depth of field (blurred) or appearing superimposed or only partially on the images were automatically eliminated using shape and contrast filters. Figure 4 shows an example of the result of image processing using HIRIS after selecting the right set of filters. For clarity concerns, only a fraction of the treated image is presented in this figure. Among all the particles in the initial picture (on the left), only the ones circled in red on the processed image (on the right) were taken into account by the software.

Figure 4: Example of HIRIS image processing

The diameter of each droplet considered by image processing was calculated from the measured surface area as the equivalent sphere surface diameter. To naturally account for the instantaneous acquisition variability in the final values used, we calculated the characteristic diameter and size distribution for each single spraying condition by gathering all the particles analyzed on the 400 images.

The subsequent calculations and optimizations were performed using the matrix calculation software MATLAB (MATLAB 9.4.0.813, MathWorks, USA). We respectively used fminsearch, integer, diff, and interp1 functions for least-squares optimization, numerical integration, differential calculation, and linear data interpolation. Monomial and exponential relations were linearized applying logarithm, allowing coefficient identification from the experimental data through a multiple linear regression based on the singular value decomposition method (function fitlm). The calculated coefficients significance was assessed through a Student test. The variance analysis was also checked to follow a F test. The machine-learning software Eureqa was used to identify relations through symbolic regression algorithms. The input variables selected for the symbolic regression were the set of experimental dimensionless numbers identified by the dimensional analysis (detailed in 3.2). All basic mathematical functions were considered, so as exponential, root-square, logarithm and power functions. Trigonometric functions were however not considered because very unlikely to seize the physics at stake. Among the relation shapes estimated by the software, we systematically kept the equation with the smallest number of coefficients as the increase in complexity in the equations identified by the software did not brought significant improvements in terms of fitness quality (R² value), while making the equation less easy to use.

Dimensional analysis

A correct and comprehensive dimensional analysis must follow some rules to be correctly established. The dimensional analysis presented in this study thus strictly follows the steps described by Delaplace et al. [START_REF] Delaplace | Dimensional analysis of food processes[END_REF] :

• Listing of the independent parameters/variables related to the target parameter; • Determination of the dimension of each parameter in terms of fundamental dimensions;

• Application of the Vashy-Buckingham theorem;

• Construction of the resulting dimensionless numbers; • Eventually, reorganization of these numbers to make appear more relevant ones.

Target variable and related parameters

When it comes to atomization, the main information the experimentalist usually needs to estimate is a characteristic diameter of the droplet size distribution. In order to detail the mathematical process employed in a generic way, a characteristic droplet diameter dchar is introduced in the subsequent sections. This characteristic diameter will be impacted by several parameters and variables:

• Material-related variables (sprayed liquid and air, at 20°C): ρl, ρa, μl, μa, σ; • Process-related variables: ul, ua, L; • Nozzle geometric parameters: dl, da,int, da,ext;

• Physical constants: g. Among the process-related variables, L (the distance between the nozzle outlet and the position where the droplet size is measured) must be taken into consideration since droplet fragmentation and coalescence phenomena may occur in the spray.

Assuming that the variables and parameters of influence are extensively listed, one can then state the existence of the following relation, where f stands for an unknown function (Equation 1):

𝑑 𝑐ℎ𝑎𝑟 = 𝑓(𝜌 𝑙 , 𝜌 𝑎 , 𝜇 𝑙 , 𝜇 𝑎 , 𝜎, 𝑢 𝑙 , 𝑢 𝑎 , 𝐿, 𝑑 𝑙 , 𝑑 𝑎,𝑖𝑛𝑡 , 𝑑 𝑎,𝑒𝑥𝑡 , 𝑔) Equation 1This relation will subsequently be simplified by making appear relevant dimensionless numbers and removing the constant parameters.

From dimensions to relevant dimensionless numbers

The above-mentioned parameters embody three fundamental dimensions: the mass dimension M, the length dimension L, and the time dimension T. The corresponding dimensional matrix is presented in Table 4.

𝒈 M 0 1 1 1 1 0 0 0 0 0 1 0 0 L 1 -3 -3 -1 0 1 1 1 1 1 -1 1 1 T 0 0 0 -1 -2 -1 -1 0 0 0 -1 0 -2
The repeated variables used in this study (µa, L, and g in the central matrix) were chosen accordingly to Petit et al. [START_REF] Petit | A dimensional analysis approach for modelling the size of droplets formed by bi-fluid atomisation[END_REF], which demonstrated the relevance of such a base while studying a similar atomization set-up. When applied, the Vashy-Buckingham theorem gives the relation between dimensionless numbers presented by Equation 2, where F stands for the mathematical function relating the causal dimensionless numbers to the target dimensionless number dchar/L.

𝑑 𝑐ℎ𝑎𝑟 𝐿 = 𝐹 ( 𝜌 𝑙 𝜇 𝑎 . 𝐿 -1.5 . 𝑔 -0.5 , 𝜌 𝑎 𝜇 𝑎 . 𝐿 -1.5 . 𝑔 -0.5 , 𝜇 𝑙 𝜇 𝑎 , 𝜎 𝜇 𝑎 . 𝐿 0.5 . 𝑔 0.5 , 𝑢 𝑙 𝐿 0.5 . 𝑔 0.5 , 𝑢 𝑎 𝐿 0.5 . 𝑔 0.5 ,

𝑑 𝑙 𝐿 , 𝑑 𝑎,𝑖𝑛𝑡 𝐿 , 𝑑 𝑎,𝑒𝑥𝑡 𝐿 ) Equation 2
These initial dimensionless numbers can afterward be rearranged to make appear more intuitive dimensionless numbers regarding the phenomena studied. In this case, the aerodynamic Weber number We (Equation 3), the Bond number Bo (Equation 4), and the Air-to-Liquid Ratio ALR (Equation 5) were chosen as of particular relevance when studying atomization: At low viscosity, the deformation of a drop inside an air stream is mostly controlled by the aerodynamic and surface tension forces. The ratio between these forces leads to the expression introduced in Equation 3 (Weber number) that compares process conditions to the physico-chemical properties of the liquid sprayed. Due to its nature, the Weber number is naturally used in a lot of studies involving droplets. The higher the We value, the more the deforming aerodynamic forces overcome the reforming surface tension forces, ultimately leading to smaller droplets [START_REF] Lefebvre | Atomization and sprays[END_REF].

𝑊𝑒 = 𝜌 𝑙 (𝑢 𝑎 -𝑢 𝑙 )
The Air-to-Liquid Ratio (ALR), also called Gas-to-Liquid Ratio, is defined as the ratio between the air and liquid mass flowrates. This dimensionless number is of typical use in the case of bifluid atomization, and frequently considered in droplet atomization models as it embodies the influence of the atomization process conditions [START_REF] Lefebvre | Atomization and sprays[END_REF][START_REF] Poozesh | Assessment of predictive models for characterizing the atomization process in a spray dryer's bi-fluid nozzle[END_REF]. The higher the ALR, the smaller the droplets generated.

The Bond number Bo is also commonly used when studying bubbles and drops as it allows their shape characterization. By comparing the gravitational forces to the surface tension ones (and thus measuring the influence of surface tension independently of the process conditions, unlike We), this number also gives an indication of the influence of the surface tension on the flow studied [START_REF] Hager | Wilfrid Noel Bond and the Bond number[END_REF].

Extensive details of the mathematical operations leading to Equation 2 can be found in a previous article [START_REF] Lachin | Dimensional analysis modeling of spraying operation -Impact of fluid properties and pressure nozzle geometric parameters on the pressure-flow rate relationship[END_REF]

Shape of the process relation and range of the numbers considered

When established, the dimensionless numbers are related through the Vashy-Buckingham theorem. However, the theorem does not state the shape of the relation between the parameters. When this relation is not known, a monomial function (that minimizes the number of constants to identify) is usually considered [START_REF] Bois | Influence of process variables on foaming ability of surfactants: Experimental study and dimensional analysis[END_REF]. This equation shape is widely found in chemical engineering studies (as for example, when studying foaming ability of surfactants in a mechanical stirrer [START_REF] Bois | Influence of process variables on foaming ability of surfactants: Experimental study and dimensional analysis[END_REF], axial dispersion in a rotary kiln [START_REF] Bongo Njeng | Effect of lifter shape and operating parameters on the flow of materials in a pilot rotary kiln: Part III. Up-scaling considerations and segregation analysis[END_REF]…). Sometimes, a combination of monomial and exponential terms can also be found in the literature [START_REF] Mandato | Liquids' atomization with two different nozzles: Modeling of the effects of some processing and formulation conditions by dimensional analysis[END_REF]. When it comes to atomization, the dimensional analysis models proposed by Mandato et al. [START_REF] Mandato | Liquids' atomization with two different nozzles: Modeling of the effects of some processing and formulation conditions by dimensional analysis[END_REF] and Petit et al. [START_REF] Petit | A dimensional analysis approach for modelling the size of droplets formed by bi-fluid atomisation[END_REF] are among the most comprehensive that currently exists. The authors respectively used monomial-exponential relation and fully monomial relation for their models. For this reason, we considered monomial and exponential terms for the conventional approach. More specifically, we used monomial terms for all the dimensionless numbers to identify characteristic diameters, except for We which influence was found to be better seized by an exponential term. For these correlations, coefficient identification was performed through multilinear regression after applying logarithm.

The overall range explored for the established dimensionless number in the conditions of this study is presented in Figure , which also indicates the range of validity of the established correlations.

Results

Based on the experimental results of the spraying trials with model solutions and reconstituted skimmed milk, different models were identified to characterize the liquid spray. Sections 4.1 and 4.2 are dedicated to the identification of the correlation presented in Equation 6 in the case of two characteristic diameters of frequent use: the median volumetric diameter dv,50 and the average number diameter dn,av. The grey dots in the subsequent figures (75 experimental conditions) represent the results obtained with the spraying of water-ethanol-glycerol mixtures (solutions 1 to 5, extensively presented in Table 1) under the process conditions introduced in Table 3. The grey dots are used to identify the equations proposed. The red dots correspond to the results obtained when spraying reconstituted milk (15 points), which are used as model validation points. Section 4.3 is dedicated to the prediction of the droplet size distribution. First, in Section 4.3.1, the atomization physical ground detailed by Villermaux is introduced. Then, using this theory, sections 4.3.2 and subsequent discuss a strategy to model a size probability distribution through dimensional analysis.

Median volumetric diameter correlations

The median volumetric diameter dv,50, also called mass median diameter (MMD) can be defined as the diameter under which half of the total volume (or mass) of the sample is accumulated. A first modelling attempt for this diameter was performed with the 75 results obtained spraying the five model solutions (trained data), taking into account all the dimensionless numbers introduced in Equation 6. However, the p-value of the coefficient associated to the density ratio 𝜌 𝑙 𝜌 𝑎 turned to be higher than 0.05, making this variable statistically insignificant. The narrow density ratio range explored here (Figure ) might explain why the model is not able to convey a clear influence of this parameter over the results. The same density ratio insignificance was found by Petit et al. [START_REF] Petit | A dimensional analysis approach for modelling the size of droplets formed by bi-fluid atomisation[END_REF] with the model established at the smallest distance from the nozzle. The authors explain that at small distances from the nozzle, physico-chemistry might not have a significant impact on the atomization process. We consequently chose to remove the density ratio from the dimensionless parameters of interest for the dv,50 modelling. Coefficient identification on the correlation led to Equation 7, with corresponding adjusted coefficient of correlation R² = 0.85 and root-mean-square error RMSE = 0.0009: 𝑑 𝑣,50 𝐿 = 0.001 * (𝐵𝑜) 0.98 * ( 𝜇 𝐿 𝜇 𝐴 ) 0.14 .(𝐴𝐿𝑅) -0.52 . 𝑒 -5.3 * 10 -4 .𝑊𝑒 Equation 7The adjustment between the experimental data and the estimated values are graphically represented on Figure 6, with dashed lines corresponding to the +-20% relative error interval. The magnitude of the coefficient associated to each dimensionless number gives some essential guidance about the phenomena' prevalence in this process. Process related numbers (ALR and We, including the air and liquid velocities) have a significant impact on the diameter calculation: the higher these velocitiesrelated numbers, the lower the droplet diameter. This result is expected, considering that the main rupture driving force in this type of nozzle is the liquid shearing through the air flow. Nevertheless, it is surprising to see that the fluid-related numbers (viscosity ratio for the liquid viscosity and Bond number for the surface tension) have such a high impact on the diameter, particularly the Bond number. Two reasons for this high impact can be invoked. The first reason is related to the distance L considered in this study, relatively close to the nozzle outlet (15 mm). At this distance and depending on the conditions, the liquid ligaments might not be fully sheared. In that case, it would not be surprising that liquid physico-chemical properties would still have a significant impact on the average droplet size. The second possibility is linked to the air volumetric flowrate values chosen. The lowest air flowrate value selected for the experiments (90 L.h -1 ) is slightly under the working recommendations of the nozzle provider. For this flowrate, results that are not fully process-driven might appear, explaining the high impact of the liquid physico-chemistry (through Bo) predicted by the model on dv,50.

Considering the range of conditions studied, the model fitness is rather satisfactory, as it reasonably predicts the diameters over the whole range of conditions considered (52 over 75 points predicted in the 20% relative error range). However, the predicted value for diameters in the range 40-60 µm is almost systematically higher than the experimental value. This deviation indicates that the mathematical equation of the model cannot perfectly account for cases in which the resulting diameter is intermediate in size. The milk validation points (red dots) confirm this hypothesis: the quality of the prediction for these experiments is similar to model solutions in this diameter range. The model depicted by Equation 7 also underpredicts the dv,50 values over 100 µm when spraying milk, failing at seizing the complexity of interactions at stake for these operating conditions. The validation data set thus clearly reveals the limits of the model.

Figure 6: Calculated vs. experimental dv,50 for the conventional model

The generation of a spray comes with three phenomena whose prevalence is highly dependent on the distance from the nozzle outlet. Following ligament formation, fine droplets are stripped from these ligaments, which then establishes coalescence/rupture equilibrium. Petit et al. [START_REF] Petit | A dimensional analysis approach for modelling the size of droplets formed by bi-fluid atomisation[END_REF] found that the higher the distance from the nozzle, the higher the impact of the physico-chemistry due to the contribution of coalescence to PSD. Whereas Petit et al. studied atomization at distances from the nozzle outlet of 5 and 14 cm, the data in this study was recorded closer to the nozzle outlet (1.5 cm). An unexpected consequence of this shorter distance is that our models consequently place significant importance on physico-chemistry. Our results thus prove that bifluid atomization is highly sensitive to the distance from the nozzle outlet, and researchers should take high precaution before applying a correlation from the literature in order to ensure reliable results.

Although the model introduced in Equation 7 is generally satisfactory, it nevertheless has some apparent flaws. One of the most problematic flaws involves how the model overestimates sizes from the validation data set between 40-60 µm and underestimates them over 100 µm. Resolving this issue required identifiying a new model using the software Eureqa (accordingly to the procedure mentioned in 2.3.1). Among the equations proposed by the software, we chose the one offering the best fitting result while keeping the number of identified constants as low as possible for practical ease of use. Equation 8 introduces the identified equation, and As was the case in Equation 7, the density ratio remains insignificant in this model. This redundancy gives confidence to the notion that such insignificance is due to physical and experimental concerns and not related to the shape of the relation. The model depicted by Equation 8 displays complex physio-chemical and process operating parameters interactions. The influence of the Bo number (gravitational forces to surface tension ones) and We number (aerodynamic forces to surface tension forces) should probably not be considered independently. These models also should include the viscous forces influence to offer a satisfactory prediction in a large range of diameters. The equation shape identified by the machine-learning approach offers an interesting representation of how physical properties can interact when spraying a fluid.

Average number diameter correlations

Under some conditions, estimating the average number diameter can also be of great interest. This diameter is usually introduced as the ratio between the moments of order 1 and 0 of the number PSD. Equation 9 presents the monomial/exponential correlation identified from our experimental results with MATLAB for dn,av.

𝑑 𝑛,𝑎𝑣 𝐿 = 2 * 10 -9 * (𝐵𝑜) 0.48 * ( 𝜌 𝐿 𝜌 𝐴 )

2.11 * (𝐴𝐿𝑅) -0.08 * 𝑒 -9.8 * 10 -5 * 𝑊𝑒 Equation 9All the dimensionless numbers were found to be statistically significant (p-values > 0.05), with the exception of the viscosity ratio The adjustment observed in Figure 8: Calculated vs. experimental dn, is excellent, with 74 of the 75 experimental points predicted within the 20% relative interval and almost all the validation points lying within this interval. The validation data also seems correctly predicted, even though a deviation starts to be visible on the largest particles. The use of the same monomial modelling structure for dn,av and dv,50 also allows some points of comparison between the prevalence of the phenomena at stake by comparing the relative contribution of the dimensionless numbers considered to the estimated diameter. The Bo number and ALR influences are significantly less pronounced in this correlation, indicating some apparent differences with the previous case in the physical phenomena impacting dn,av.

By choosing the estimation of a characteristic diameter like dn,av (related to a number distribution), the interactions between the physico-chemistry of the product and the process seem to be less complicated for simple monomial-exponential models to account for. An explanation could be that number distributions give significantly less importance to the largest particles than volume distributions, where a single large particle can have a massive impact on the distribution (as the volume is proportional to the size power three, assuming sphericity). Therefore, largest particles moderately impact characteristic diameters in the case number distributions. It would mean that complex physicochemistry/process interactions mainly impact the largest particles, explaining the deviation observed for the "middle size" particles for dv,50 models. This poor diameter prediction on midrange values could indeed be related to the fact that large median volumetric droplets hold significant importance when identifying the model, thus impacting the overall prediction quality of median volumetric diameter.

As we previously did for dv,50, a symbolic regression correlation was also searched for dn,av using the software Eureqa considering all the dimensionless numbers established by the dimensional analysis. The equation obtained is presented in Equation 10, while the graphical representation of the prediction for trained and validation data sets is shown in Figure . For this model, R² = 0.93 and the RMSE = 0.00010; 73 points over 75 are predicted in the 20% error interval. The software also found the viscosity ratio to be insignificant for the dn,av prediction, as was the case with the conventional model. As seen in Figure , the viscosity ratio interval considered is rather large; the magnitude of the values considered is thus unlikely to explain this insignificance, which is thus most likely due to the physics of the phenomena involved. The overall quality prediction for the trained data is very similar to the one obtained with Equation 9. However, clear differences are observed when considering the validation data set: the 20-30 µm points are now almost perfectly predicted, while a lesser deviation is observed for the largest particles.

When observing the shape of the relation proposed by Eureqa, it seems that We is only involved inside an interaction term somewhat similar to the one obtained with Equation 8, with Bo on the numerator and We*ALR on the denominator. This interesting fact thus gives some leads on how to accurately account for the impact of the interactions while modelling an atomization characteristic diameter.

Figure 9: Calculated vs. experimental dn,av for the model identified with symbolic regression

Based on the models obtained for dv,50 and dn,av, some general guidelines can be drawn. While the simple monomial/exponential relations were satisfactory overall, the ones obtained through a machine-learning approach were systematically more accurate while keeping the same number of parameters (5 in total) to be adjusted. The comparison between the relations also highlighted that depending on the distribution considered (number or volume), the physical phenomena at stake would not have the same impact on the final property estimation, probably due to the relative accountancy of the largest droplets of the distribution. The formulas proposed by the machine-learning approach thus highlighted some way to better seize these complex interactions from a numerical point of view. The comparison between the diameters and the methods (conventional and symbolic regression) thus offers an interesting way to make atomization modelling more complex while shedding some light on the physical interactions at stake.

Particle Size Distribution modelling 4.3.1 The physics of sprays and jets

Liquid fragmentation is a long-time studied physical phenomenon that continues to fascinate engineers and physicists Recent studies of Villermaux et al. [START_REF] Villermaux | Ligament-Mediated Spray Formation[END_REF]23,[START_REF] Villermaux | Fragmentation[END_REF] shed light on the mechanisms involved during coaxial liquid atomization, and proposed a fragmentation theory known as ligament-mediated spray theory. Their theoretical background remains used in recent atomization modelling research [START_REF] Keshavarz | Ligament Mediated Fragmentation of Viscoelastic Liquids[END_REF][START_REF] Xu | Atomization model based on the ligament-mediated spray mechanism for coaxial two-fluid air-blast nozzle[END_REF]. According to this theory, when sprayed, the liquid bulk undergoes instabilities due to shearing at the liquid/gas interface (Kelvin-Helmholtz type, then Raileigh-Taylor when the amplitude of the first waves is large enough). These Rayleigh-Taylor instabilities combined with the influence of the air flow generate ligaments (elongated liquid fractions stemming from the bulk due to the surface disturbances), which are comprised of elements called blobs.

Villermaux [START_REF] Villermaux | Ligament-Mediated Spray Formation[END_REF] showed that the distribution of these blobs follows a fragmentationcoalescence dynamics, which number probability distribution can be described with a Gamma distribution such as (Equation 11):

𝑃 (𝑥 = 𝑑 𝑑 𝑛,𝑎𝑣 ) = ℾ(𝑛, 𝑥) = [ 𝑛 𝑛 ℾ(𝑛)
] 𝑥 𝑛-1 𝑒 -𝑛𝑥 Equation 11Where 𝑑 𝑛,𝑎𝑣 is the average number diameter of the blob distribution and n is a parameter indicating the broadness of the distribution: the lower is n, the broader is the distribution (Figure 10). Villermaux successfully fitted this distribution shape to the droplets stripped from the ligaments, validating its relevance in the case of bifluid atomization. However, for some atomization conditions, a significant ligament population can also appear in the final spray distribution. In that case, the final spray distribution is a convolution of the large ligament distribution and small droplets distribution, resulting in an exponential tail displaying the same n value as the blob distribution.

Adequation of the Gamma function to the experimental droplet size distributions

Two parameters appear in Equation 11: the average number diameter of the distribution 𝑑 𝑛,𝑎𝑣 and the polydispersity index n. 𝑑 𝑛,𝑎𝑣 is also the ratio between the first-order moment and zero-order moment of the distribution, and is thus easy to estimate knowing the diameter of each droplet. After estimating 𝑑 𝑛,𝑎𝑣 , n is adjusted by fitting the numerical integral of the Gamma function to the cumulative distribution experimentally obtained through a least-squares optimization. By defining p as the number of the droplet considered (sorted by increasing size), ptot as the total number of droplets considered and xp as the normalized particle size related to the particle p, the least-squares function numerically assessed is described by Equation 12:

𝐿𝑆 = 𝑚𝑖𝑛 [∑ ( 𝑝 𝑡𝑜𝑡 -∫ [ 𝑛 𝑛 ℾ(𝑛) ] 𝑥 𝑝 𝑛-1 𝑒 -𝑛.𝑥 𝑝 𝑥 𝑝 0 . 𝑑𝑥 𝑝 ) 2 𝑝 𝑡𝑜𝑡 𝑝=1 ]
Equation 12

The results obtained for two different experimental conditions (sprayed solution 1) are represented in The agreement between the experimental distribution and the theoretical function is good for large n values (left graph, n = 13.5). While the peak height is underestimated, the Gamma function accounts well for the global shape of the distribution. More precautions should be taken when using the Gamma function for small n values, as the function clearly tends to overestimate the probability distribution in the regions near the peak while underestimating its intensity. This result is actually expected, since it has been reported elsewhere [START_REF] Keshavarz | Ligament Mediated Fragmentation of Viscoelastic Liquids[END_REF][START_REF] Marmottant | On spray formation[END_REF]. More precisely, it is mentioned that when the ligament distribution is narrowly distributed inside the droplet population (corresponding to large n values), Equation 11 is indeed suited to describe the spray distribution. Nevertheless, for cases in which air momentum is such as the ligaments exhibits more corrugation, the global size distribution tends to turn into an incomplete Gamma function. Also, the shape of the distribution can be considered relatively simple (2 fitting parameters) considering the complexity of the physical phenomena involved (air shearing, entrainment, coalescence, and rupture). 3-parameters PSD functions are usually more efficient in taking into account both the maximal height of the peak and the polydispersity in the case of atomization [START_REF] Urbán | Application of big data analysis technique on high-velocity airblast atomization: Searching for optimum probability density function[END_REF]. However, these functions could not be recommended for this study, since the primary goal is to provide the experimentalist a quick way to assess a size distribution over an extensive range of experimental conditions. A fast estimation of n is interesting to estimate the relevance of using the simple two-parameters approach proposed in this study. For small n values, the reconstitution of the distribution should be carefully managed, and more complex modelling options should be considered.

Determining a correlation for n

The physical parameters impacting the characteristic diameter will also impact the polydispersity index n, which is a dimensionless number. The dimensionless numbers quantifying the dimensionless characteristic diameter are thus also relevant for modelling n. A first attempt to identify n through a monomial equation was performed, which is shown in Equation 13. 𝑛 = 𝑎 0 . (𝐵𝑜) 𝑎

Equation 13

After removing the non-significant parameters, the coefficient of correlation between the model/experiment and the root mean-squared error were calculated and found to be R² = 0.57 and RMSE = 2.09. The monomial form is thus unable to account for the complex interaction parameters at stake in the polydispersity of the system, and is not recommended to estimate n.

The use of a machine-learning approach is especially relevant here. The obtained model can be considered to be good, as it correctly predicts the value of n for the validation data set. The shape of Equation 14 clearly exhibits strong interaction terms accounting for the cases in which the polydispersity is not fully process-driven. We can observe that the higher the value of the process parameters (mainly air velocity), the higher the value of ALR and We number and the higher the value of n, indicating a lesser polydispersity (Bo*We and 𝐵𝑜 𝐵𝑜+𝐴𝐿𝑅 ). This influence is however attenuated by the viscous forces, which can be expected to exert some resistance to the air shearing, and thus increase the polydispersity (lowering n): the term 𝑊𝑒 * 𝜇 𝐿 𝜇 𝐴 is indeed negative. These interpretations should be carefully considered, as the dimensionless numbers do not account for one single effect (as for example, surface tension is also part of the Weber number). However, this equation and these interpretations give key parameters governing the phenomena to account for when modelling polydispersity and some guidance about the way to simply consider complex interactions for atomization.

Particle size distribution prediction

The prediction quality of the PSD modelling was then tested by using the previously established models. Firstly, the prediction of n by Equation 14 was used in combination with the Gamma function (Equation 12) to predict the normalized number PSD for milk spraying (validation data set). Two results, corresponding to calculated n values respectively equal to 13.2 and 18.4 are presented in Figure 43: Experimental and predicted Gamma functions for milk3. The real number PSD were also fully reconstructed using the dn,av modelled (Equation 10) inside the Gamma function. The results obtained for the same conditions as earlier are plotted in Figure 54. The addition of this model does not have a significant impact on the global quality of the model: no significant diameter offset of the main mode is indeed observed. For the reconstruction performed, the quality of the dn,av and n models obtained with machine-learning is thus enough not to degrade the Gamma function quality of prediction. While not perfect, this attempt of reconstituting a full PSD based on simple correlations shows what can be expected from such relations and gives some solid leads for improvement.

Conclusion

The modelling of complex physical phenomena often involves high complexity in the applied mathematical approach. Atomization is no exception to that, as both process conditions and physicochemistry have a significant impact on the shearing, coalescence, and rupture phenomena involved.

Using a rigorous dimensional analysis methodology, this article explores a way to model characteristic diameters related to volume and number distributions in the case of bifluid nozzle atomization. More precisely, two models are proposed per diameter: one based on the conventional monomial/exponential function shape and the other one obtained identifying a relation between relevant dimensionless numbers using a machine-learning software. The models found with the symbolic regression approach predict the droplet diameters for the milk validation data set systematically better. However, while some clear limitations are visible for the conventional models, they can overall be qualified as good. Comparing the models' terms for the two diameters also allowed to see that, depending on the distribution type (volume or number), the phenomena involved do not have the same impact on the characteristic diameter, possibly due to the magnitude of the contribution of the largest particles.

A methodology is also proposed to estimate a number PSD, based on the fragmentation theory introduced by Villermaux [START_REF] Villermaux | Ligament-Mediated Spray Formation[END_REF]. The identification of a correlation relating the polydispersity index n to the process parameters indeed allows, through the use of a Gamma function, to predict number PSD for milk spraying. The distribution reconstruction at large n values is overall satisfactory while the PSD prediction at low n numbers is more questionable (due to the Gamma function relevance itself more than to the models accuracy). However, the value of n in itself contains very interesting information about the polydispersity that can be expected, and the need or not to consider a more complex modelling.

Figure 1 :

 1 Figure 1: Bifluid nozzle used with the pilot set-up

Figure 2 :

 2 Figure 2: Spraying and image acquisition set-up

Figure 3 :

 3 Figure 3: Example of images acquired during spraying of solution 2. Left: air flowrate = 90 L.min -1 , liquid flowrate = 50 g.min -1 ; right: air flowrate = 90 L.min -1 , liquid flowrate = 90 g.min -1

Figure 5 :

 5 Figure 5: Dimensionless number range explored

Figure 7 :

 7 Calculated vs. experimental dv,50 for the model identified with symbolic regression represents the accordance between experimental data and modelled diameters.

8 Figure 7 :

 87 Figure 7: Calculated vs. experimental dv,50 for the model identified with symbolic regression

.

  Equation 9 was thus obtained after removing the viscosity ratio for the identification. The adjusted coefficient of correlation R² and the root mean square error RMSE for this modelling were respectively R² = 0.92 and RMSE = 0.00011. The adequacy between the experimental and model values is presented in Figure , and milk points are plotted with red dots.

Figure 8 :

 8 Figure 8: Calculated vs. experimental dn,av for the conventional model

Figure 10 :

 10 Figure 10: Gamma distribution for various n values.

Figure 1 .

 1 The experimental distributions (dashed lines) were reconstituted from the cumulative experimental data by linearly interpolating the results over 1000 equidistant points relatively to size (interp1 function), then numerically differentiating the interpolated data (diff function). Modelled distributions were obtained by plotting Equation 11 with the value of n optimized through leastsquares function minimization computed by Equation12.

Figure 11 :

 11 Figure 11: Experimental distribution (dashed line) and Gamma function (solid line) for two spraying conditions with solution 1

Figure 43

 43 Figure43shows that the quality of the model for the normalized distribution is similar to the one found when assessing the relevance of using the Gamma function (Figure): the overall distribution shape is seized by the model, but the main peak intensity is underpredicted. In the case of milk, for high n numbers, both the shape of the distribution and the n prediction are thus relevant to have an indication about the expected size distribution when atomizing the solution. For low n values (not represented here), the same behavior as in Figure was observed. In these cases, while the value of n is a precious indication of polydispersity degree that can be expected, distribution reconstitution should be processed with care.

Figure 43 :

 43 Figure 43: Experimental and predicted Gamma functions for milk

Table 1 :

 1 Composition and physical properties at 20°C of the sprayed liquids

	Solution	% Water (w./w.)	% Glycerol (w./w.)	% Ethanol (w./w.)	μ (mPa.s)	ρ (kg.m -3 )	σ (mN.m -1 )
	1	100	0	0	1.1	998	65.1
	2	20	80	0	58.4	1200.6	59.8
	3	15	85	0	107.0	1202.3	58.6
	4	14.75	75.25	10	55.5	1147.6	35.8
	5	9.9	81.1	9	107.3	1162.2	36.3
	Milk (30% w./w.)	-	-	-	18.7	1105	41.9

Table 2 :

 2 Nozzle geometric parameters value with coefficient of variation over 3 measurements

		dl	da,int	da,ext
	Value (mm)	1.06	3.62	4.97
	Var coef. (%)	0.55	0.69	0.12

Table 3 :

 3 Air and sprayed liquid flowrate values chosen for atomization experiments (20°C)

	Solution	Liquid mass flowrates (g.min -1 )	Air volumetric flowrates (L.min -1 )
	1	50; 90; 130; 170; 210	90; 110; 130
	2	50; 90; 130; 170; 210	90; 110; 130
	3	50; 70; 90; 110; 130	90; 110; 130
	4	50; 90; 130; 170; 210	90; 110; 130
	5	50; 70; 90; 110; 130	90; 110; 130
	Milk (30% w./w.)	50; 90; 130; 170; 210	90; 110; 130

Table 4 :

 4 Dimension matrix D, with central matrix C (dark grey) and Residual matrix R (white)

	dchar ρl	ρa	µl	σ	ul	ua	dl	da,int da,ext µa	L

  , and the extensive number reorganization giving birth to We, Bo, and ALR is detailed in the work of Petit et al.[START_REF] Petit | A dimensional analysis approach for modelling the size of droplets formed by bi-fluid atomisation[END_REF]. The work presented in this article focuses on a single nozzle geometry. The 𝜌 𝑎 𝜇 𝑎 .𝐿 -1.5 .𝑔 -0.5 is constant as we did not consider any variation in the air properties and position of the observation zone. The relation between the ten dimensionless parameters introduced in Equation 2 thus simplifies to the relation between six parameters depicted in Equation 6:

	dimensionless numbers	𝑑 𝑙 𝐿 ,	𝑑 𝑎,𝑖𝑛𝑡 𝐿 ,	𝑑 𝑎,𝑒𝑥𝑡 𝐿	are thus constant. Also,
	𝑑 𝑐ℎ𝑎𝑟 𝐿	= 𝐹′(𝐵𝑜,	𝜇 𝑙 𝜇 𝑎	, 𝐴𝐿𝑅, 𝑊𝑒,	𝜌 𝑙 𝜌 𝑎	)	Equation 6

  1 . (

	𝜇 𝐿 𝜇 𝐴	) 𝑎 2 .(𝐴𝐿𝑅) 𝑎 3 . (𝑊𝑒) 𝑎 4 . ( 𝜌 𝐴 𝜌 𝐿	) 𝑎 5

  The equation obtained with the software Eureqa is presented in Equation 14, and the graphical representation of the adequation between the experimental results and the model is shown in Figure32: Calculated vs. experimental n for the model identified with symbolic regression2.

	𝑛 = 16.4 + 0.012 *	𝜇 𝐿 𝜇 𝐴	+ 0.027 * 𝐵𝑜 * 𝑊𝑒 -	32.2 * 𝐵𝑜 𝐵𝑜 + 𝐴𝐿𝑅	-32.87 * 𝐵𝑜 -1.4 * 10 -6 * 𝑊𝑒 *	𝜇 𝐿 𝜇 𝐴	Equation 14

Figure 32: Calculated vs. experimental n for the model identified with symbolic regression For this correlation, R²=0.87, and the RMSE is of 1.1. 67 of the 75 points of the training data are in the 20% relative error interval.

Acknowledgements

The authors would like to thank Soredab for funding the internship of Mamadou Lamine Niane.