domain: spi phys

text:
Benoit Varillon
email: benoit.varillon@isae-supaero.fr

Jean-Baptiste Chaudron
email: jean-baptiste.chaudron@isae-supaero.fr

David Doose
email: david.doose@onera.frcharles.lesire@onera.fr

Charles Lesire

Real-Time Polling Task: Design and Analysis

Keywords: Timing predictability, Formal methods, Modeling, Real-time Analysis, Validation

Usually, robotic systems must be designed as reactive systems because they must be able to react and adapt to their environment and also communicate with some other systems. However, the reactive behavior, because of its non-deterministic nature, is difficult to analyze and prevents designers to perform a proper real-time analysis which is usually needed for critical robotic systems. In this paper, we propose a deterministic task model to handle the reactive behavior as well as the necessary tools to analyze it and verify the respect of real-time constraints. An implementation of this model, which is used in a ROS2 patch, is also presented.

I. INTRODUCTION

The aim of autonomous systems is to evolve in an uncontrolled environment. To achieve this goal, the system needs to be able to interact with its environment, because it is continuously evolving and can't be controlled, it is important that actions are still relevant when they are executed and not due to out-of-date information. Not respecting this constraint can lead to unwanted situations with damage to the system or even to people if the system evolves alongside human. Such systems that can cause hazardous situations if they behave poorly are called critical systems. When designing such system, it is important to be able to model it using known abstractions in order to validate safety properties. For real-time systems the model needs to have a deterministic behavior, without determinism prediction of the potential executions will not be possible and real-time analysis will be difficult with classical methods. In robotic applications, the system is often subject to external requests that are not predictable, thus there is a need for a model that can represent this type of interaction where the system needs to react to unsteady arrival of messages. The best model, in terms of latency, is a reactive model where the task is activated as soon as the requests arrive. However, a fully reactive model is not deterministic.

Scheduling theory proposes multiple models and methods to analyze the timing behavior of real-time systems: a periodic model (and a simple way to analyze it) is proposed in [1], an improvement with a strictly better utilization bound is given by [2]. The multiframe model [START_REF] Mok | A multiframe model for real-time tasks[END_REF] allows more precise analysis that suites well for video encoding tasks, and has been generalized to take deadlines into account by [START_REF] Baruah | Generalized multiframe tasks[END_REF]. The recurring real-time task [START_REF] Baruah | Dynamic-and static-priority scheduling of recurring realtime tasks[END_REF] and the digraph [START_REF] Stigge | The digraph real-time task model[END_REF] models are more complex and can generalize the multiframe model, but at the cost of a higher complexity of the analysis. The method proposed in [START_REF] Ripoll | Improvement in feasibility testing for real-time tasks1[END_REF] allows determining if a system, composed of different types of tasks, is schedulable or not, only knowing the demand-bound function for each task. Another method proposed in [START_REF] Joseph | Finding response times in a real-time system[END_REF], known as Response Time Analysis (RTA), provides a way to compute the exact worst-case response time for every task of a system.

ROS is one of the most used middleware in robotics and was already studied by [START_REF] Erős | A ROS2 based communication architecture for control in collaborative and intelligent automation systems[END_REF] who shows the advantages of using ROS2 on industrial systems and [START_REF] Puck | Performance evaluation of real-time ROS2 robotic control in a time-synchronized distributed network[END_REF] who looks at its real-time capabilities. Investigation on the influence of hardware and OS on ROS2's real-time performances was done in [START_REF] Sarkar | Real-time jitter measurements under ROS2: the inverted pendulum case[END_REF]. One remaining problem with using ROS for real-time systems is its communication model based on publish/subscribe mechanism that uses reactive execution model. The method proposed in [START_REF] Casini | Response-Time Analysis of ROS 2 Processing Chains Under Reservation-Based Scheduling[END_REF] is very effective for analyzing chains of reactive tasks but assumes a provided event arrival curve, which is a big limitation. The method presented in [START_REF] Fohler | Static scheduling of pipelined periodic tasks in distributed real-time systems[END_REF] represents the system as a periodic pipeline of tasks; this approach works well when the event source that initiates the pipeline is periodic, which is an assumption we want to remove to be able to work with an unknown event source. The problem of event synchronization was addressed by [START_REF] Garcia | Schedulability analysis of distributed hard real-time systems with multiple-event synchronization[END_REF], but for very complex systems and so the analysis is also complex and not very efficent on simple systems. The communication between tasks was handled by [START_REF] Fonseca | Response time analysis of sporadic dag tasks under partitioned scheduling[END_REF], but the focus was on multi-core architectures rather than on distributed architectures where network interferences have to be taken into account. A common way to make a reactive system analyzable is to use active polling, this way the activation of tasks is limited to known instants and the system is so deterministic. The method of [START_REF] Gobillot | A Design and Analysis Methodology for Component-Based Real-Time Architectures of Autonomous Systems[END_REF] which proposes to model tasks using state machines can be used to study such tasks, but this would yield a model too complex to be verified.

In this paper, we present a method to represent reactive tasks with a deterministic execution model and provide solutions to analyze it. Section II introduces and formalizes our model for reactive tasks. Then, in section III, this model is analyzed, formulas are given for specific configurations and an algorithmic solution is presented for the general case. Section IV compares our algorithm with the resolution using an SMT solver. Finally, section V presents an implementation of our model in a real-time patch for ROS2.

II. TASK AND SYSTEM MODELS

The objective is to analyze systems composed of periodic and reactive tasks, so a model compliant with the periodic model is needed for reactive tasks. As these hypotheses are valid in most of the recent systems, a first assumption considers a preemptive scheduler with fixed priority. Another assumption is that the tasks under analysis are statically allocated among CPUs with no migration at run-time and are independent, so a CPU per CPU analysis is sufficient and can be done on many-core architectures.

A. System

In the scheduling analysis context, a system S can be considered as a set of tasks, which will each generate an ordered sequence of jobs that will be executed. A task is characterized by its priority which defines the order in which tasks will be executed, and its execution parameters, which depend on the type of task and define the possible sequences of jobs. A periodic task τ k C k , T k is described by its period T k and its worst-case execution time C k and will generate every T k unit of time a job that requires at most C k unit of CPU time, a polling task

ρ k C P k , T P k , C R k , T R k
will generate a sequence of jobs as described in II-B.

B. Polling Task

In this paragraph, the aims is to find a new model for reactive tasks that has similar characteristics in terms of reactivity, but which is deterministic and easily analyzable. The polling task model relies on active polling to ensure deterministic behavior. The state diagram of the polling model is presented in Fig. 1. At each iteration the polling state, in charge of verifying if a message has arrived, is executed for C P units of time, then if there is no message the task goes to a waiting state for the polling period T P . If there is a message, the task goes in the running state and process the callback for C C units of time and then goes in a waiting state for the execution period T R units of time. These characteristics are static. This model with two periods allows controlling individually the latency, i.e. the time between the reception of a message and its processing, and the CPU load. The polling period T P fixes the time between two polling phases and so the latency, the execution period T R limits the frequency of execution of the callback and so the CPU load generated by the task. By writing P the execution of the polling loop and R the execution of the running loop, traces of execution can be written in the form [P, R] * . Figure 2 shows the trace of an execution that could be written P 3 RP 2 RP . To make easier the writing of equations in the rest of this paper, we will write:

C R = C P + C C (1)
This allows to see the task as two independent loops. A polling task is defined by giving the characteristics of the two loops:

ρ k C P k , T P k , C R k , T R k (2) with C R > C P .
The utilization factor of each loop and the global utilization factor are written:

U P k = C P k T P k ; U R k = C R k T R k ; U k = max(U R k , U P k ) (3)
This model is a generalization of the classic periodic model, with a polling state that only returns true the execution loop executed at each iteration and the task is so periodic. An example of polling task is the handling of communication with a GNSS, if the GNSS sends messages every 50ms with a precision on the frequency of 50μs, assuming the time for computation of the message is 1ms and the time for executing the polling state is 5μs, the polling tasks that implement this could have the following characteristics 5, 25, 1000, 50000 .

C P T P C C T R

III. POLLING TASK ANALYSIS

Systems are composed of different types of tasks and even if it is possible to use the same model for every task, this will come at the cost of performance of the analysis or over-sizing the hardware part of the system. It is therefore interesting to be able to analyze them using a generic method that can be used among different types of tasks, in particular with periodic tasks which are the most often used.

A. Analysis principle

Two interesting ways of analyzing systems of tasks are the use of an utilization bound and the computation of the exact response time of each task. Utilization bounds, as introduced by Liu and Layland [1], are useful because they allow determining if a system is schedulable with a linear time but are not very precise and lead to reject systems that are schedulable. The response time analysis, although it is more complex, allows computing the exact worst-case response time of all tasks of the systems and thus determining precisely if they all meet their deadline and so if the system is schedulable or not. It is, therefore, useful to be able to analyze a polling task with utilization bounds, but also with the response time analysis. The response time analysis, presented by [START_REF] Joseph | Finding response times in a real-time system[END_REF], consists for periodic tasks in finding the fixed point of the sequence given by Eq. ( 4). The only prerequisite concerns the global utilization factor of the system U that must be lower than one.

R 0 = C i R n+1 = R n Ti • C i + j∈hp(i) Rn Tj • C j (4)
This method can be generalized for other types of tasks with Eq. ( 5), requiring only the knowledge of the release-bound function, defined as the maximum workload that a task can generate on the interval [0, t[, for each task of the system.

R 0 = ⎧ ⎨ ⎩ C i for periodic tasks max C P i , C R i for polling tasks R n+1 = rbf i (R n ) + j∈hp(i) rbf j (R n ) (5)
With rbf k the release-bound function of task τ k and hp(k) the set of tasks of equal or higher priority than τ k . Thus, to analyze systems containing polling tasks with RTA it is necessary and sufficient to provide a way of computing the release-bound function of those tasks. The rest of this section present a way to compute a linear upper bound of the rbf to allow an easy way of analyzing a system. Then, a more precise analysis of the polling task model is done. This analysis brings out a specific case and the corresponding analytical formulas to compute rbf . Next, a Sat Modulo Theory (SMT) formalization of the polling model allows handling the calculation of the rbf for the general case. Finally, a more efficient algorithm is presented.

B. Polling task request bound function

The release bound function of a polling task is a staircase function with different types of steps. The building of the curve uses two sizes of steps, one for the polling loop, and the other for the running loop. At the end of a step, i.e. the end of an iteration, the task is in its initial state and booth loops can be executed, this will generate two new possible steps that will each generate two other steps. Then the rbf is given by keeping for each t the higher step. This process of construction is shown in Fig. 3.

C. Request bound function upper bound

Because the rbf is complex to manage, it is useful to have an approximation that is easier to use during design phases when it is needed to have a quick way to find out temporal area in which the rbf will be. The loop with the higher utilization factor will drive the general shape of the rbf , the other loop will only generate interference. As shown on Fig. 3 the rbf will so follow a linear growth given by the global utilization factor of the task, U k , and the maximal interference is the given by the higher step, C R . The upper bound is so given by Eq. [START_REF] Stigge | The digraph real-time task model[END_REF].

t rbf (t) rbf k (t) C R
rbf k (t) = U k • t + C R k (6)
This linear bound allows computing the rbf with an O(1) complexity. It is a pessimistic approximation, so systems found schedulable with the bound will always be schedulable, a less pessimistic bound is given by the staircase function of Eq. ( 7), this bound is not linear but still can be computed in O(1).

rbf k (t) = U k • t + C R k (7)

D. Specific cases

When analyzing the curve, three phases are visible. The first lasts T R and has an analytic formula, the third is periodic and so easy to analyze. In some conditions, there is a more complex phase between the two previous, it comes from the interweaving of different curves and lasts a duration that is bound but not fixed. Figure 4 gives an example of the rbf of a polling task. a) Phase 1: The first phase lasts from t = 0 to t = T R the end of the running period, during this time the running loop can run only once. So the maximum CPU utilization is obtained when the polling loop runs a maximum times in [0, t[, and then the last iteration of the task is the running loop.

b) Phase 2: This phase is the most complex. The curve is the continuation of phase 1 curve inter-weaved with other curves from executions that were dominated in phase 1 but, in phase 2, overtake locally the phase 1 curve. This makes this part chaotic and prevents finding an analytic formula. In some conditions, this phase does not exist (see III-D3).

c) Phase 3: In this last phase, a pattern is repeated periodically, this pattern depends on the end of phase 2 and so no generic formula can be found. When there is no phase 2 it is phase 1 which is repeated and so a formula can be found.

The analysis of the three phases allows us to point out some particular cases for which the rbf can be described using an analytical formula and so can be computed in O(1) complexity.

1) Running loop domination: When the polling period is greater than the running period, because C R > C P the running loop dominates the polling loop, and, as it is shown in Fig. 5, any execution that contains one polling loop will have a release function lower than the release function of the 'only execution loop' execution, so the rbf is given by: t rbf (t)

T P T P C P C R T R Fig. 5. Domination of the running loop rbf (t) = t T R • C R iff T P ≥ T R (8)
2) Phase 1: This phase can always be described with an analytical formula. The running loop can run only once, so every execution has the form P n R, and the execution with maximum released CPU time is the one with the most polling loops and then one running loop, rbf is given by 9.

rbf (t) = ⎧ ⎨ ⎩ 0 if t = 0 t T P -1 • C P + C R iff t ∈]0, T R ] (9)
3) Simple harmonic case: When the running period is a multiple of the polling period every execution that did not take part in maximizing the rbf in phase 1 will overlap a curve that leaded to maximizing the rbf in phase 1, so no new interference is generated after T R and phase 2 does not exist. When the utilization factor of the polling loop is greater than the one of the running loop, the polling loop will dominate the execution, and every execution that uses the running loop will not leads to maximizing the rbf anymore. The execution that maximizes the released CPU time is the one using the most polling loops and then uses one running loop. The following formulas represent this "simple harmonic case":

T R = k • T P with k ∈ N + and U P ≥ U R (10) rbf (t) = ⎧ ⎨ ⎩ 0 iff t = 0 t T P -1 • C P + C R iff t > 0 (11)
4) Complex harmonic case: If the utilization factor of the running loop is greater than the one of the polling loop, the formula is more complex. In that case, the global shape of the rbf is given by the running loop and the pattern of phase 1 is repeated on every step of the running loop. The formula is built by shifting right the equation of phase 1 by the number of running periods in [0, t[ and adding the same number of C R :

T R = k • T P with k ∈ N + and U R ≥ U P ( 12
)
rbf (t) = ⎧ ⎨ ⎩ 0 iff t = 0 m T P -1 • C P + t T R + 1 • C R iff t > 0 (13)
with

m = t - t T R • T R (14)

E. Optimization problem

For the general case, phase 2 exists and there is no analytical formula. However, the calculation of the rbf can be written as an optimization problem that consists, for every t, of finding the number of iterations of each loop that maximize the released CPU time. An execution is made of a series of loops, each one could be the running or the polling loop. The amount of released CPU time on an interval is the sum of the execution time of each loop that occurred during this interval, every executions path with the same number of each loop will request the same amount of CPU time regardless of the order of the loops. By writing r(t) the amount of request CPU time on [0, t[, i the number of execution loops and j the number of polling loops in [0, t[ we have:

r = i • C R + j • C P (15)
To find the release bound function, we need to maximize r(t). For t which are a linear combination of the periods, it means to find (i, j) that respects [START_REF] Gobillot | A Design and Analysis Methodology for Component-Based Real-Time Architectures of Autonomous Systems[END_REF] and maximizes [START_REF] Fonseca | Response time analysis of sporadic dag tasks under partitioned scheduling[END_REF].

t = i • T R + j • T P (16)
For values that are not a linear combination of the periods, we have to transform ( 16) into (17) but this condition does not take into account the last loop, which has not finished its execution yet. This loop could be either the polling or the execution loop but to maximize r(t) we should use the maximum released time which is C R , so (15) becomes (18):

x < t x = i • T R + j • T P (17) r = i • C R + j • C P + C R (18)
Figure 6 shows a description of this optimization problem using SMT-lib21 [START_REF] Barrett | The SMT-LIB Standard: Version 2.6[END_REF]. This language was chosen because it is both sufficiently expressive and easily readable for humans. To resolve this problem, the Z32 [START_REF] De Moura | Z3: an efficient SMT solver[END_REF] solver which is known for its efficiency was used. It is important to notice that other solvers or modelization languages could be used.

This SMT model was used to prove the earlier formulas, by writing a SMT modelization of the formulas3 , like in Fig. 7 for the bound, it is possible to ask the solver to find a counter example with (assert (> r b)), if the result is unsat then there is no counter example and the formulas is proven. The formulas of phase 1 and the bound (Eqs. (6, 9)) was proved with any other condition, the formulas of the running domination and the harmonic cases (Eqs. [START_REF] Stigge | The digraph real-time task model[END_REF][START_REF] Joseph | Finding response times in a real-time system[END_REF][START_REF] Sarkar | Real-time jitter measurements under ROS2: the inverted pendulum case[END_REF][START_REF] Fohler | Static scheduling of pipelined periodic tasks in distributed real-time systems[END_REF]) were proved for bounded intervals.

; ----------inputs --------- ;------------------------------(define-const Cr Int (+ Cp Cc)) (declare-const i Int) (declare-const j Int) (define-const x Int (+ ( * i Tp) ( * j Tr))) (define-const r Int (+ ( * i Cp) ( * j Cr) Cr)) ;-----(assert (< x t)) (assert (< 0 i)) (assert (< 0 j)) (assert (< 0 r)) ;----------optimize ----------(maximize r) (check-sat); returns sat; rbf(100) = 19 (get-value (r)) ; 19 (get-value (i)) ; 1 (get-value (j)) ; 5 (get-value (x)) ; 96 Fig. 6. SMT modelization and rbf computation for a polling task

F. Algorithmic solution

The SMT modeling of the problem introduced in the previous section finds an exact solution to the request-bound function of the polling tasks. However, this solution suffers from an efficiency problem. Indeed, its computation time is not negligible, and this problem increases with t. Moreover, this calculation is performed many times during the complete analysis of the system.

In this section, we present an algorithmic solution to the problem of computing the request-bound function of a polling task. This solution is more efficient than the previous one because it takes into account the structure of the rbf described represents the lasso structure, which is composed of a set of nodes N = {n 0 , . Fig. 8. Lasso representation

1) Computation principle:

The first node contains the first execution of the polling loop and the first execution of the running loop (n 0 = {(x 0 , y 0 ), (x 1 , y 1 )}).

The principle is as follows: we build the lasso structure step by step, following the rbf by moving from one point of interest to another. The next point of interest is the next (minimum t axis) end of the segments. At each iteration, all the new possible executions (polling and running) are added and filtered if necessary. When a new node is added we check if the same node has already been computed; in this case, the loop is detected and the lasso computation is complete, otherwise, the computation goes on. Algorithm 1 describes the computation of the lasso path. The centering function (Algorithm 2) computes the new edge while taking into account the new point of interest.

Segment drop: It is important to notice that it is not useful to add all the segments. Indeed, some segments are useless for the computation of the request bound function because it exists other segments that dominate the new one. Figure 9 highlight two cases in which the segment (x j , y j ) is dominated by (x i , y i ). Because the rbf is built recursively by adding the same two new segments to every point of interest, the curve originated from a point P 1 is the same as the curve originated from another point P 2 shifted by ---→ P 2 P 1 , so if a segment ends later and lower than another segment the worst curve originated from the first will never overtake the one from the second.

Adding a new segment (node.add) can also lead to the dropping of dominated segments already in the current node. Table I details all the dropping conditions. It is also important to notice that thanks to the dominated segments dropping, the next point of interest (next center) is unique.

TABLE I SEGMENT DROP CONDITIONS

x i < x j x i = x j x i > x j y i < y j drop (x,y i ) drop (x i , y i ) y i = y j drop (x j , y j ) drop one drop (x i , y i ) y i > y j drop (x j , y j ) drop (x j , y j )
Computation example: Figure 10 highlights the first steps of the lasso structure computation in a simple example. Fig. 10a represent the initial node n 0 = {(x 0 , y 0 ), (x 1 , y 1 )} with x 0 = T P , y 0 = C P , x 1 = T R , and y 1 = C R . Then in this example, the next point of interest is the end of the polling loop (x 0 , y 0 ). Fig. 10b show the new node centered with according to the first edge e 1 = (x 0 , y 0 ) in which the first segment is consumed by the centering function and the second is translated (x 1 = x 1x 0 , y 1 = y 1y 0 ). The second node is computed (Fig. 10c) and then the third one (Fig. 10d).

2) Request bound function from Lasso: Once the lasso structure is computed, the calculation is straightforward. Let's notice X(i) (resp. Y (i)) the sum of the dx (resp. dy) of the edges to reach node n i :

X(i) = j∈1..i dx j (19)
This notation will be used with S (resp. L) to represent the sum of the dx of edges between nodes of the sequence (resp.

t x i y x j t x i y i
x j y j Fig. 9. Segment domination. Dashed segments represent the next step if the corresponding segment is kept

t x 0 x 1 y 0 y 1 next center (a) Initial node (n 0 ) t x 1 y 1 (b) Initial node centered t x 1 y 1 x 2 y 2 x 3 y 3 (c) Second node (n 1 ) t x " 1 y " 1 x 3 y 3 x 4 y 4 x 5 y 5 (d) Third node (n 2 )
Fig. 10. Lasso computation example. Blue and green segments are the two new segments added at the current step, black segments are remaining from previous step. loop). Let's I(t) be the node index of the node useful to compute rbf (t):

I(t) = {max i∈0..N i | X(i) < t} (20)
The request bound function has two different formulas: one if no loop is needed to reach t, and another one if some loops are needed.

No loop: No loop is needed if the last node "is after" t:

t < X(N ) + dx N +1 (21)
In this case, we just have to find the maximum y coordinate of the segments, before t, of the useful node:

rbf (t) = Y (I(t)) + max (x,y)∈ni y | X(I(t)) + x < t (22) With loop: With loops the principle is to represent t = X(S) + k • X(L) + t with k = t-X(S) X(L)
and t = tk • X(L)-X(S). Then the release function is defined as follows:

rbf (t) = Y (S) + k • Y (L) + Y (L(t )) (23)

G. Analysis conclusions

Different ways of computing the rbf have been presented to suit the different situations. When a conservative approximation is sufficient a linear bound is given by Eq. [START_REF] Stigge | The digraph real-time task model[END_REF]. When an exact value of the rbf is wanted, if the situation corresponds to one of the cases in Tab. II the corresponding analytical formula can be used. In the general case the most efficient method is to use the lasso algorithm presented in section III-F. Eq. ( 8)

∀t ≤ T R Eq. ( 9)

T R = k • T P U R ≤ U P
Eq. ( 11)

U P ≤ U R Eq. ( 13
)
IV. ALGORITHMS COMPARISON Different methods to compute the release bound function for polling tasks are presented in this paper. This section presents the results of a benchmark that compares these methods in terms of performance. The aim of these algorithms is to be used in a more global algorithm, like the RTA, to compute the worst-case response time of a polling task in a complete system. In this context, the algorithm will be instantiated with a task, then the rbf will be computed for multiple values of t. In order to compare the performance of algorithms, we need to emulate these real conditions, so the benchmark will consist in creating a set of tasks, then for each algorithm, we measure for each task the time for instantiating the algorithm and compute the value of the rbf (t) for different t. To mimic real systems, tasks were generated using random values for their parameters, worst-case execution times (wcet) are generated using a uniform distribution on [0, 1000] and the periods with a uniform distribution on [wcet, 100000].

Results

Tab. III shows the time for analyzing a system in function of the number of access to the rbf value for systems of different sizes. This is visible that the lasso algorithm is more efficient than solving the optimization problem with Z3. The experiments done shows that the lasso algorithm is more than 10 5 faster than the SMT solver. This can be explain by the way this two methods work. The optimization problem call the maximize function of the SMT solver every time a value of the rbf is needed, this results in solving multiple optimization problem which are each complex to solve. On the other hand the lasso algorithm compute the lasso of Fig. 8 once, then accessing a value can be done by simply reading the lasso in a constant time.

V. POLLING TASKS IMPLEMENTATION IN ROS2

Because ROS2 suffers from the impossibility to be analyzed since it widely uses reactive tasks, we develop an extension, Corail 4 , which aims to solve this problem by providing an implementation using polling tasks.

A. ROS2

ROS is an open-source middleware dedicated to robotic applications and distributed systems which is based on the publish/subscribe paradigm. Elements of the system are called nodes and have to register to topics as a publisher or subscriber to be able to communicate. The version 2 was designed to be more real-time than the version 1. The use of DDS greatly improved the real-time capacity of the communication but the execution part still lack of real-time capabilities. The execution in ROS2 is handled by executors; the two default executors provided allow neither the preemption of tasks nor to order them by priority. Furthermore reactive model is used for subscription which prevent real-time analysis.

B. Corail

Corail is a ROS2 package that aims to make ROS2 more suitable for real-time systems. Corail relies on the concept of executors of ROS2 and provides a new real-time executor that ensures the respect of the assumptions needed for real-time analysis. This executor allocates a POSIX thread for each task of the system, which makes the task preemptive and allows using the POSIX scheduler to handle priority between tasks. Moreover, the Corail executor uses the polling task model to implement subscriptions and services. The use of POSIX and the polling task model enable Corail to ensure the respect of the constraints needed for real-time analysis.

C. Comparison

To compare ROS2 and Corail, two systems were implemented using both plain rclcpp and Corail. The first presented in Tab. Va is made of two periodic tasks, the second, presented in Tab. Vb is made of one periodic task and one subscription. Trace of the executions were made using LTTng 5 , Figs. 11 and12 show the traces of the implementation using only ROS 2. We can observe that there is no synchronization at the beginning of the run-time and executions of τ 1 when τ 2 is running are missed and that there is no preemption. This leads to systems that do not respect their deadlines while the two systems are found schedulable using the scheduling theory. The traces from the Corail implementations, detailed in Figs. 13 and14, show executions that follow the hypotheses (preemption, synchronization, priority) and so the prevision of the scheduling theory. This is due to the fact that the execution follow the hypothesis needed for the scheduling theory. Moreover the use of the polling task model allows to compute the exact worst case response time of the reactive task without any knowledge of the event source. In this paper, we address a problem that is regularly faced in our robotic experiments. Indeed, the embedded systems that we implement are reactive critical real-time systems. In this sense, they must be able to react to their environment but also communicate with each other while being analyzable from a real-time point of view. That is to say, it is essential to be able to show the respect of the deadlines of the different tasks. The difficulty lies in the fact that we have no guarantee either of communication times or of the dates of incoming events. To address this issue, we have proposed a method based on three pillars: a new task model that allows taking into account event polling, an analysis method that allows studying the respect of deadlines, and an implementation corresponding to this execution model (including a ROS2 patch).

We have given different solutions for the analysis method depending on the possible system configurations. Some have an analytical formula. For the general case, we have proposed a model in SMT that allows solving this problem and a specific algorithm. We explained the efficiency gain of the algorithmic version with the help of benchmarks. Finally, we presented our implementation in ROS2 for this new task model.

Our future work will first focus on the experimentation and deployment of these developments. Then we want to extend this execution model to take into account communication time because in certain configurations the networks we use allow it.

ACKNOWLEDGMENT

Fig. 1 .tFig. 2 .

 12 Fig. 1. Polling task state diagram. Round states are execution that lasts the held value, diamond states wait to the beginning date of the curent execution plus the held value.

Fig. 3 .

 3 Fig. 3. Polling task request bound function with polling steps in green, running steps in blue and linear upper bound in red.

3 Fig. 4 .

 34 Fig. 4. Curve structure. Only steps that contribute to the rbf are shown.

(Fig. 7 .

 7 Fig. 7. SMT modelization of the bound formula

Algorithm 1 edges := ∅ 3 nodes := ∅ 4 node 6 loop 7 prev

 13467 Lasso computation algorithm 1 function COMPUTE LASSO 2 := [(T P , C P ), (T R , C R )] 5 nodes.push(node) := nodes.last() 8 edge, node := prev.center() 9 edges.push(edge) 10 node.add((T P , C P )) 11 node.add((T P , C P )) 12 id = same id(nodes, node) 13 if id then 14 sequence.nodes := nodes[1..id] 15 sequence.edges := edges[1..id] 16 cycle.nodes := nodes[id..len(nodes)] 17 cycle.edges := edges[id..len(nodes)] 18 return (sequence, cyle) 19 else 20 nodes.push(node) Algorithm 2 Centering algorithm (used in Alg. 1) 1 function NODE::CENTER 2 first := self.f irst() 3 edge := (first.x, first.y) 4 new node := ∅ 5 for seg in self do 6 new node.add((seg.xedge.x, seg.yedge.y)) 7 return (edge, new node)

Fig. 11 .

 11 Fig. 11. Execution trace of system 1 with ROS 2

Fig. 12 .

 12 Fig. 12. Execution trace of system 2 with ROS 2

  This work was supported by the Defense Innovation Agency (AID) of the French Ministry of Defense (research project CONCORDE N°2019 65 0090004707501) REFERENCES [1] C. L. Liu and J. Layland, "Scheduling algorithms for multiprogramming in a hard-real-time environment," Journal of the Association for Computing Machinery, vol. 20, no. 1, p. 46-61, Jan 1973. [2] E. Bini, G. C. Buttazzo, and G. M. Butazzo, "Rate monotonic analysis: The hyperbolic bound." EEE Transactions on Computers, vol. 52, no. 7, pp. 933-942, 2003.

Fig. 13 .

 13 Fig. 13. Execution trace of system 1 with Corail, hatched areas correspond to preemption.

Fig. 14 .

 14 Fig. 14. Execution trace of system 2 with Corail, hatched areas correspond to preemption.

  ., n N } and edges E = {e 1 , .., e N +1 }. The first part of nodes represents the sequence noted S, and the second part is the loop noted L. Each node represents a part of the query function. A node contains a set of segments (n 1 = {(x 1 , y 1 ), . . . (x n , y n )}). Each segment represents either the execution of the polling loop or the execution of a running loop. Each edge is a shift on the axes of the request bound function (e 1 = (dx 1 , dy 1 )). This shift is a centering of the next node on the next point of interest. The notion of centering will be detailed later.

							e 5
	n0	e 1	n1	e 2	n2	e 3	n3	e 4	n4
			S				L

http://smtlib.cs.uiowa.edu/

https://github.com/Z3Prover/z3

https://gitlab.com/corail1/smt proof in III-B. Indeed, the latter has a form of "lasso path": it starts with a sequence composed of phase 1 then phase 2; then a loop formed by phase 3. The algorithm we propose calculates an efficient representation of this "lasso path". Once this computation is done, the cost of computing a query for a given t becomes negligible, which makes the method very efficient for the complete analysis of the system. We will show this efficiency in the next section. Figure8

https://corail1.gitlab.io


Wala Aloulou
email: walaaloulou6@gmail.com

Hajer Aloulou
email: hajer.aloulou89@yahoo.fr

Mohamed Romdhani

Joelle Duplay
email: jduplay@unistra.fr

Raja Ben Amar

Evaluation of TiO 2 /smectite nanoparticles as an alternative low-cost adsorbent for chromium removal from industrial wastewater

Keywords: Adsorption, Chromium, TiO 2 /smectite, Nanoparticles, Electroplating industry wastewater, Modeling

In this work, chromium removal efficiency from electroplating industry wastewater was studied by adsorption using TiO 2 /smectite nanoparticles (NP) synthesized by modification of natural smectite by colloidal route. The effect on chromium adsorption of adsorbent dose, pH, contact time, temperature and initial concentration was then determined in batch system. Chromium concentration can be reduced to 24 mg/L (70%) under the experimental condition (pH = 2.9, adsorbent dose = 1.6 g/L, contact time = 75 min, T = 298 K and C 0 = 80 mg/L) when initial chromium concentration of 80 mg/L is employed. The chromium adsorption on NP was described by the Langmuir isotherm and the maximum chromium adsorption capacity was found as 35 mg/g. Kinetics data were best described by the pseudo-second-order model. The thermodynamic studies proved that the adsorption was exothermic and spontaneous.

Introduction

Nowadays, the increased industrialization generates multiple metallic pollutants and causes a serious threat on the environment (air, water and soil) due to the toxicity generated by those substances [START_REF] Cuhorka | Removal of micropollutants from water by commercially available nanofiltration membranes[END_REF]. Therefore, the treatment of the industrial effluents before their discharge into water bodies is required to minimize the pollution of the water resource [START_REF] Damtie | Removal of fluoride in membranebased water and wastewater treatment technologies: performance review[END_REF].

In particular, heavy metals including cobalt (Co), cadmium (Cd), chromium (Cr), lead (Pb) and mercury (Hg), even in trace amounts, cause harmful effects on aquatic fauna such as respiratory irritation and hepatic cell destruction [START_REF] Banni | Monitoring pollution in Tunisian coasts: application of a classification scale based on biochemical markers[END_REF][START_REF] Barhoumi | Cadmium bioaccumulation in three benthic fish species, Slaria basilisca, Zosterisessor ophicephalus, Solea vulgaris collected from the Gulf of Gabes in Tunisia[END_REF][START_REF] Kessabi | Involvement of Cd bioacumulation in spinal deformities occurrence in natural populations of Mediterranean Killifish[END_REF]. In aquatic environment near industrialized coastline areas, these heavy metals can affect fish, mollusks, tuna [START_REF] Chaffai | Physico-chemical froms of storge of metals (Cd, Cu and Zn) and metallothionein-like proteins in gills and liver of marine fish from Tunisia coast: ecotoxicological consequences[END_REF][START_REF] Kessabi | Evaluation of environmental impact on natural populations of the Mediterranean Killifish Aphanius fasciatus by quantitative RNA biomarkers[END_REF]. Heavy metals are found in wastewater of various industries such as surface treatment time allergy and electroplating industries [START_REF] Panayotova | Purification and reuse of heavy metals containing wastewaters from electroplating plants[END_REF]. Especially, chromium is a toxic contaminant, even in very low concentrations. Chromium in its hexavalent form is the predominant form of chromium species in most water streams and surface water [START_REF] Hashim | Remediation technologies for heavy metal contaminated ground water[END_REF]. Cr(VI) is well known for its health-related issues in humans including carcinogenic and mutagenic risks [START_REF] Igwe | A bioseparation process for removing heavy metals from waste water using biosorbents[END_REF][START_REF] Park | Removal of chromium by activated carbon fibers plated with copper metal[END_REF]. Due to very high solubility and mobility of chromium(VI) ions in ecosystems, chromium(VI) ions are considered the most toxic when compared to other forms of chromium salts [START_REF] Muthukrishnan | Effect of pH on rejection of hexavalent chromium by nanofiltration[END_REF]. Moreover, Cr(VI) is included in the list of the US Environmental Protection Agency (EPA) as one of the most toxic heavy metals requiring prioritized control. According to the World Health Organization (WHO) standards, the acceptable level of chromium ion in drinking water is 0.05 mg/L for Cr(VI) and 0.1 mg/L for total Cr [START_REF] Mahmoud | Remediation of Cr(VI) via combined self-redaction and adsorption by chemically modified carbon sorbents[END_REF].

The conventional techniques for the removal of chromium(VI) wastes include precipitation, adsorption, membrane processes, oxidation and ion exchange [START_REF] Alemu | Removal of chromium(VI) from aqueous solution using vesicular basalt: a potential low cost wastewater treatment system[END_REF][START_REF] Aloulou | Developing of titania-smectite nanocomposites UF membrane over zeolite based ceramic support[END_REF]. However, most of these methods results in partial removal of metal ions, low selectivity, high consumption of reagents, high operational cost and production of secondary pollutants [START_REF] Wu | Solvent extraction of scandium(III), yttrium(III), lanthanum(III) and gadolinium(III) using Cyanex 302 in heptane from hydrochloric acid solutions[END_REF][START_REF] Bayuo | Adsorptive removal of chromium(VI) from aqueous solution unto groundnut shell[END_REF]. Especially, adsorption has distinct advantages which include selectivity for specific metal, reusability of biomaterial, short operation time and low operating cost [START_REF] Rahmani | Bioremoval of lead by use of waste activated sludge[END_REF][START_REF] Shah | Sorption isotherms and kinetics of chromium uptake from wastewater using natural sorbent material[END_REF].

Bayuo et al. [START_REF] Bayuo | Adsorptive removal of chromium(VI) from aqueous solution unto groundnut shell[END_REF] found that the groundnut shell can be an efficient alternative adsorbent material to remove chromium(VI) from aqueous solution. The influence of operating conditions such as contact time, pH, adsorbent dose, initial metal concentrations and temperature was studied. Labied et al. [START_REF] Labied | Adsorption of hexavalent chromium by activated carbon obtained from a waste lignocellulosic material (Ziziphus jujuba cores): kinetic, equilibrium, and thermodynamic study[END_REF] reported that 62.08% of hexavalent chromium removal was achieved by activated carbon obtained from a waste of lignocelluloses material (Ziziphus jujuba cores).

In the present work, separation performance of adsorption process during the removal of chromium solution by TiO 2 /smectite nanoparticles (NP) was evaluated. The effect of operating parameters such as dose of NP, pH, contact time, temperature and initial concentration on the adsorption capacity of NP during the chromium removal was investigated. Furthermore, adsorption kinetic and isotherm studies were undertaken to determine the adsorbate removal rate and the maximum adsorption capacity. Thermodynamic properties, such as enthalpy, entropy and the free energy were also estimated to provide insight into the adsorption mechanism during the adsorption removal of the chromium. These NP are used in different applications such as membrane separation [START_REF] Aloulou | Developing of titania-smectite nanocomposites UF membrane over zeolite based ceramic support[END_REF][START_REF] Aloulou | Synthesis and characterization of clay-based ultrafiltration membranes supported on natural zeolite for removal of heavy metals from wastewater[END_REF][START_REF] Aloulou | Clay-based ultrafiltration membranes supTiO 2 and nanocomposite clay materials over zeolite support for oily wastewater purification and heavy metals removal[END_REF][START_REF] Aloulou | Characterization of an asymmetric ultrafiltration membrane prepared from TiO 2 -smectite nanocomposites doped with commercial TiO 2 and its application to the treatment of textile wastewater[END_REF], hybrid system and adsorption [START_REF] Romdhani | Performance studies of indigo dye removal using TiO 2 modified clay and zeolite ultrafiltration membrane hybrid system[END_REF] to treat different effluents. High performances of these NP encourage their use in the elimination of heavy metals by adsorption.

Experimental procedures

Materials

The nanoparticles (NP) of TiO 2 /smectite were synthesized in a previous work via sol-gel method [START_REF] Rahmani | Bioremoval of lead by use of waste activated sludge[END_REF]. The NP powder was obtained by incorporating titanium(IV) isopropoxide with in the organo-modified smectite (Fig. 1). The resulting (NP) have a homogeneous structure with a uniform distribution of the TiO 2 nanoparticles and a mean particle size in the range of 8-12 nm.

Characterization of NP

The Fourier-transform infrared spectroscopy (FTIR) of purified smectite, cetyltrimethylammonium bromide modified smectite (CTAB/Sm) and NP were obtained by using a PerkinElmer spectrometer to observe the surface functional groups. The zeta-potentials of the samples were determined using a zeta sizer. The measurements were conducted in wide pH range (pH was adjusted with HCl and NaOH aqueous solutions) in order to identify the surface charge of the NP. The average pore sizes, surface areas, and pore volumes of NP were determined using the Brunauer-Emmett-Teller (BET) method and Micromeritics ASAP 2020 apparatus.

Industrial effluent

The industrial wastewater used in this study is produced by a Tunisian electroplating industry. The characterization of the wastewater with initial pH of 2.9 presents chromium concentration of 80 mg/L. Physicochemical parameters of raw and treated effluents were measured according to the standard methods suggested by American Public Health Association. pH measurement was done by a pH-meter (ISTEK pH-220L, Japan). The content of heavy metals was measured by Atomic Absorption Spectroscopy (AAS) (PerkinElmer AAnalyst 200).

Adsorption in a batch system

Adsorption experiments were performed in batch system. The effect of the operational parameters, adsorbent dosage (0.2-2 g/L), contact time (5-120 min) and temperature (298-303 K) were carried out. All the uptake experiments were conducted using 50 mL of the test solution in 100 mL Erlenmeyer flask on a magnetic shaker at 200 rpm. After equilibrium, the samples were filtered and the filtrate was then analyzed. The percentage of chromium adsorption was calculated as follows by Eq. ( 1):

Adsorption % § © ¨• ¹ ¸u C C C i f i 100 (1)
The adsorption capacity of chromium compounds by NP was determined by Eq. ( 2):

q C C m V i f § © ¨• ¹ ¸u ( 2
)
where C i and C f are the initial and final (or equilibrium) chromium concentrations, respectively. V is the volume of the chromium solutions (L), and m is the weight of the adsorbent (g).

Results and discussion

NP characterization

From our previous work [START_REF] Aloulou | Developing of titania-smectite nanocomposites UF membrane over zeolite based ceramic support[END_REF][START_REF] Aloulou | Synthesis and characterization of clay-based ultrafiltration membranes supported on natural zeolite for removal of heavy metals from wastewater[END_REF], the FTIR data shows that the smectite particles are strongly coated with TiO 2 nanoparticles. The point of zero charge (pH pzc ) of the synthesized NP was found around 1.7 or 1.8. Therefore, the NP has a negative surface charge in the pH range of 1.7-12. The NP size was determined by BET analysis. The surface area, average pore size and pore volumes of NP were 92 m 2 /g, 8-12 nm and 0.19 cm 3 /g, respectively. These values are more important than the results of the purified smectite that did not exceed 72 m 2 /g [START_REF] Aloulou | Developing of titania-smectite nanocomposites UF membrane over zeolite based ceramic support[END_REF]. The higher surface area presents more available active sites [START_REF] Tan | Role of biochar surface characteristics in the adsorption of aromatic compounds: pore structure and functional groups[END_REF]. As a result, raising the surface area can be a workable tactic for further enhancing adsorption on the NP.

Adsorption in a batch system 3.2.1. Effect of the NP dose

The effect of the dose of NP on the adsorption of chromium was evaluated by changing the dosage in the range from 0.2 to 2 g/L while the pH and the initial chromium concentration were kept constant at 2.9 and 80 mg/L respectively (Fig. 2). The results showed that the adsorption efficiency increased with the increase of the adsorbent dose. This indicates that the higher the NP dosage, the higher the number of active adsorption sites. However, the removal of the chromium remained unchanged when the adsorbent dose is over 1.6 g/L. This behavior assumes that the adsorption reaction can reach a dynamic equilibrium, and that the maximum chromium removal of 70%, was achieved for a minimum dosage value of 1.6 g/L.

Effect of pH

The performance of Cr(VI) removal is strongly influenced by pH. This parameter is related to the existing of ionic species of Cr(VI) in the solution and also to the surface charge of the adsorbent. The Cr(VI) form shows distinct ionic species. It often takes the forms of HCrO 4 -and Cr 2 O 7 2-when the pH is between 1 and 6.5. CrO 4 2-is the predominant form when the pH is higher than 6.5 (Fig. 3a) [START_REF] Ramirez | Removal of Cr(VI) from an aqueous solution using an activated carbon obtained from teakwood sawdust: kinetics, equilibrium, and density functional theory calculations[END_REF]. Therefore, the study of the effect of this parameter is very useful to study the mechanism and the interaction adsorbent-Cr(VI). The Cr(VI) adsorption (%) was investigated in the pH range of 2.9-10 as shown in Fig. 3b. The adsorption (%) was strongly decreased from 70% to 12.7% when the pH increased from 2.9 to 10. It is clear that the adsorption was high at low pH (2.9). This behavior can be explained by the electrostatic attractions between the positively charged sites formed under acidic pH condition (H + predominate) by the protonation of the NP and the predominate ionic species of Cr(VI) (HCrO 4 -) [START_REF] Jia | Efficient adsorption to hexavalent chromium by iron oxalate modified D301: characterization, performance and mechanisms[END_REF]. However, the increasing of the pH caused a reduction in the protonation degree of the NP and consequently the NP ability to interact with Cr(VI) reduced [START_REF] Ramirez | Removal of Cr(VI) from an aqueous solution using an activated carbon obtained from teakwood sawdust: kinetics, equilibrium, and density functional theory calculations[END_REF][START_REF] Jia | Efficient adsorption to hexavalent chromium by iron oxalate modified D301: characterization, performance and mechanisms[END_REF]. Besides, the -OH in the NP can act as powerful electron donors and it can cause a direct reduction for Cr(VI) to Cr(III). In this, the negative charge surface of the NP can interact with the Cr(III) species and the hydrolyzed form such as Cr(OH) 2 + , Cr 2 (OH) 2 4+ , and Cr(OH) 2 + [START_REF] Liu | Effects of Cu(II) on the adsorption behaviors of Cr(III) and Cr(VI) onto kaolin[END_REF]. When the pH was over 5, the precipitation of Cr(OH) 3 increased as a result of increasing OH -ions. This precipitation was not uptake by the acidic functional groups of the NP [START_REF] Islam | Recent innovative research on chromium(VI) adsorption mechanism[END_REF]. As a result, the adsorption (%) was strongly reduced at higher pH. Therefore, in this study, the initial industrial wastewater pH of 2.9 was considered for all the experiments.

Effect of contact time

The adsorption efficiency of the chromium onto NP as a function of contact time is illustrated in Fig. 4. It is noticed that the adsorption increases with contact time and reaches equilibrium within 75 min showing a maximum efficiency of chromium retention of about 70%. Therefore, 75 min was selected as the period of optimum contact time for further experiments.

Effect of the temperature

The adsorption of chromium on NP was studied by varying the temperature from 298 to 308 K at the optimal adsorbent dose of 1.6 g/L during 75 min as contact time. According to Fig. 5, chromium removal decreased from 70% to 61% when the temperature increased from 298 to 308 K, leading to a decrease in the adsorption capacity. The increase of the temperature allows the chromium to desorbs from the adsorbent surface to the solution, or to damage the active binding sites in the adsorbent which explains the decrease in the retention of chromium [START_REF] Bayuo | Adsorptive removal of chromium(VI) from aqueous solution unto groundnut shell[END_REF]. As a result, the optimum temperature was selected as 298 K for further adsorption experiments.

Effect of initial Cr(VI) concentration

The effect of initial Cr(VI) concentration (20-120 mg/L) on the efficiency of the NP adsorbent was investigated and the adsorption (%) was determined as shown in Fig. 6. The higher efficiency of adsorption (%) of Cr(VI) ions was observed at lower concentrations. Specifically, a total retention was observed for initial concentration between 20 and 40 mg/L. This indicates that for any Cr(VI) solution with a concentration below 40 mg/L, the NP will retain almost totally the Cr(VI) ions. Beyond this value, the adsorption (%) was decreased from 87.3 (at 60 mg/L) to 39 (at 120 mg/L). This behavior might be explained by the large accessible active sites of NP available at low concentration, allowing a total removal of ions [START_REF] Masindi | Facile synthesis of halloysite-bentonite clay/magnesite nanocomposite and its application for the removal of chromium ions: adsorption and precipitation process[END_REF]. Consequently, 40 mg/L was considered as the highest possible concentration when a dose of NP of 1.6 g/L is used.

Adsorption kinetics

Adsorption kinetic models are utilized to interpret the experimental data in order to identify the mechanism controlling the adsorption process. The equilibrium data were calculated using kinetic models such as pseudo-firstorder, pseudo-second-order and intra particle diffusion model. The pseudo-first-order model can be expressed as follows (Eq. ( 3) [START_REF] Lagergren | Zur Theorie De r Sogenannten, AdsorptionGeloster Stoffe, Kunglia a[END_REF]):

ln ln q q q k t e t e 1 (3)
where q t and q e are the adsorption capacity (mg/g) at time t (min) and at equilibrium respectively and k 1 , the rate constant of the model (min -1 ). The adsorption rate constants k 1 can be obtained experimentally from the plot of ln(q e -q t ) against t (Fig. 7a). Experimental data were also determined by the pseudo-second-order kinetic model described by Eq. (4): where k 2 (g/mg•min) is the rate constant of the model and q e is the maximum adsorption capacity (mg/g) [START_REF] Ys | Study of the sorption of divalent metal ions onto peet[END_REF]. The values of k 2 and q e were determined from the slopes and intercepts of plots of t/q t against t (Fig. 7b).

t q k q t q t e t § © ¨• ¹ ¸ § © ¨• ¹ 1
The intra particle diffusion model is usually used for identifying the adsorption mechanism for design purposes. The equation related to this model is expressed by Eq. ( 5) [START_REF] Balarak | The use of low cost adsorbent (Canola residues) for the adsorption of methylene blue from aqueous solution: isotherm, kinetic and thermodynamic studies[END_REF]:

q k t C t 3 0 5 .
(

) 5
where q t is the amount of fluoride ion adsorbed (mg/g) at time t (min), k 3 is the intraparticle diffusion rate constant (mg/g•min), which can be calculated from the slope of the linear plots of q t vs. t 0.5 and C is the intercept (mg/g) (Fig. 7c).

All constants of the different models are summarized in Table 1. It is clear that the determination coefficient (R 2 ) obtained from pseudo-second-order model is higher to that found from the pseudo-first-order and intra particle diffusion model. In addition, the value of q e (q e,cal ) determined from the model perfectly fitted with the experimental values of q e (q e,exp ), showing that the removal process chromium by adsorption on NP follows the pseudo-second-order kinetic model.

Adsorption isotherm models

The determination of the adsorption isotherms is important for describing the adsorption mechanism related to the interaction of the chromium ions on the NP adsorbent surface. The determination of the adsorption isotherm corresponds to the distribution of molecules between the liquid phase and the solid phase when the adsorption process reaches an equilibrium state.

In the current work, two sorption isotherm models, Langmuir [START_REF] Langmuir | The adsorption of gases on plane surface of glass, mica, and platinum[END_REF] and Freundlich [START_REF] Freundlich | Over the adsorption in solution[END_REF] were studied.

Langmuir model is characterized by the unique coating layer on sorption surface. It can be defined according to the following Eq. ( 6):

C q q K C q e e L e m 1 max (6)
The separation factor is expressed by Eq. ( 7): where q e the equilibrium chromium concentration on the adsorbent (mg/g), C e the equilibrium chromium concentration in the solution (mg/L), q max the monolayer adsorption capacity of the adsorbent (mg/g) and K L the Langmuir adsorption constant (L/mg). The values of Langmuir constants q max and K L were obtained from the slope and intercept of the linear plot of C e /q e against C e (Fig. 8a). The value of R L indicates that when R L is between 0 and 1 the adsorption is favorable, while R L > 1 represents unfavorable adsorption, and R L = 1 represents linear adsorption, while for R L = 0, the procedure of adsorption is irreversible. The Freundlich isotherm is based on the adsorption on the heterogeneous surface. The linear form of this isotherm can be expressed by Eq. ( 8) [START_REF] Freundlich | Over the adsorption in solution[END_REF]:

log log log q K n C e F e § © ¨• ¹ 1 (8)
where K F an empirical constant (mg/g), indicative of the adsorption capacity and 1/n an empirical parameter, indicative of the adsorption intensity. The Freundlich isotherm constants K F and 1/n were obtained from the slopes and intercepts of the linear plot of logq e vs. logC e (Fig. 8b). Table 2 summarizes the adsorption constants determined according to Langmuir and Freundlich models. It is to notice that the linear equation of Langmuir model presents higher determination constant (R 2 = 0.997) than for the Freundlich model (R 2 = 0.985). Therefore, Langmuir model is chosen as the best model to describe the adsorption of chromium on NP. This model suggests that the retention of chromium implies the coverage of a monolayer on the NP surface.

Adsorption thermodynamics

Temperature is an important parameter affecting the adsorption capacity of sorbents and the transport/kinetic process of the chromium adsorption. Three thermodynamic parameters, including Gibbs free energy change (ΔG°), standard enthalpy change (ΔH°) and standard entropy change (ΔS°) of chromium adsorption on NP were determined by using Eqs. ( 9) and ( 10) [START_REF] Wang | Adsorption of basic dyes on activated carbon prepared from Polygonum orientale Linn: equilibrium, kinetic and thermodynamic studies[END_REF]:

' q G RT K d ln (9) ln K S R H RT d ' q § © ¨• ¹ ¸ ' q § © ¨• ¹ ¸ ( 10
)
where R is the universal gas constant (8.314 J/mol•K), T is the absolute temperature (K) and K d is the adsorption equilibrium constant. ΔH° and ΔS° values can be determined from the slope and intercept of the Van't Hoff equation by plotting lnK d against 1/T (Fig. 9). The thermodynamic parameters are summarized in Table 3. The negative values of ΔG° at different temperatures show the spontaneous nature of the adsorption and the thermodynamically feasible process. The decrease in ΔG° values with the increase in temperature indicates the possibility of the decrease of the chromium retention at higher temperatures.

The negative ΔH° (-18.75 kJ/mol) shows that the adsorption process is exothermic. In addition, the adsorption type, physical or chemical, can be determined from the magnitude of enthalpy change. It is conventional that the adsorption is physical when the magnitude of ΔH° is less than 84 kJ/ mol. However chemical adsorption is ranging from 84 to 420 kJ/mol [START_REF] Errais | Efficient anionic dye adsorption on natural untreated clay: kinetic study and thermodynamic parameters[END_REF]. Therefore, in this study the determination of the ΔH° value shows a physical adsorption involving weak interaction forces such as van der Waals. The ΔS° parameter was found to be -53.19 J/mol•K. This negative value suggests reduction in the randomness at the solid/ solution interface during the adsorption process [START_REF] Ra | Equilibrium, thermodynamic and kinetic studies on biosorption of Pb(II) and Cd(II) from aqueous solution by macro fungus (Lactarius scrobiculatus)[END_REF].

Adsorption mechanism

In general, the pH affects the degree of ionization and specifications of the adsorbate as well as the surface charge of the adsorbent during adsorption process [START_REF] Bayuo | Adsorptive removal of chromium(VI) from aqueous solution unto groundnut shell[END_REF][START_REF] Babel | Cr(VI) removal from synthetic wastewater using coconut shell charcoal and commercial activated carbon modified with oxidizing agents and/or chitosan[END_REF]. So the effect of pH directly impacts the electrostatic interactions between the adsorbate and adsorbent's surface.

For the industrial wastewater, at pH 2.9 a large number of H + ions exist in the solution medium and the dominant form of Cr(VI) is HCrO 4 -. Consequently, the surface protonation of NP leads to the formation of positively charged sites. These protons interact with chromium atoms by an electrostatic attraction phenomenon so that HCrO 4 is retained on NP [START_REF] Bayuo | Adsorptive removal of chromium(VI) from aqueous solution unto groundnut shell[END_REF][START_REF] Labied | Adsorption of hexavalent chromium by activated carbon obtained from a waste lignocellulosic material (Ziziphus jujuba cores): kinetic, equilibrium, and thermodynamic study[END_REF].

Discussion

The performances of NP for chromium removal were compared with others nanoparticles reported in the literature (Table 4). Zou et al. [START_REF] Zou | Adsorptive removal of Cr(III) and Fe(III) from aqueous solution by chitosan/attapulgite composites: equilibrium, thermodynamics and kinetics[END_REF] synthesized and tested chitosan/attapulgite composites for the removal of Cr(III) from aqueous solution. The maximum adsorption capacity was 65.37 mg/g with adsorbent dosage of 0.2 g/L at pH 5 and a temperature of 45°C. The chromium adsorption was described by Langmuir isotherm and intraparticle diffusion model. Wang et al. [START_REF] Wang | Adsorptive removal of Cr(VI) from aqueous solutions with an effective adsorbent: cross-linked chitosan/montmorillonite nanocomposites in the presence of hydroxy-aluminum oligomericcations[END_REF] prepared chitosan-Al-pillared montmorillonite nanocomposite to remove Cr(VI). They reported that the maximum adsorption capacity of chromium was 15.68 mg/g for pH 6.38 and ambient temperature. On the other hand, the kinetic results indicate that the adsorption of chromium follows the pseudo-second-order model and that the adsorption isotherm is described by the Langmuir model.

Alemu et al. [START_REF] Alemu | Removal of chromium(VI) from aqueous solution using vesicular basalt: a potential low cost wastewater treatment system[END_REF] used vesicular basalt volcanic rock to remove Cr(VI). It was reported that the maximum removal achieved at pH 2 was of 79.20 mg/kg at an initial concentration of 5.0 mg/L and adsorbent dosage of 50 g/L. The study of adsorption isotherms shows in this case that the pseudo-second-order kinetic model and the Freundlich isotherm model best describe the adsorption process [START_REF] Alemu | Removal of chromium(VI) from aqueous solution using vesicular basalt: a potential low cost wastewater treatment system[END_REF].

Another strategy was developed by Kumar et al. [START_REF] Kumar | Application of cellulose-clay composite biosorbent toward the effective adsorption and removal of chromium from industrial wastewater[END_REF] by synthesizing Na-montmorillonite/cellulose nanocomposite as adsorbent. The application of this nanocomposite for the removal of Cr(VI) from industrial wastewater showed an adsorption capacity of chromium of 22.3 mg/g at the pH range 3.8-5.5. The removal kinetics was found to follow the Freundlich Pseudo-first-order [START_REF] Romdhani | Performance studies of indigo dye removal using TiO 2 modified clay and zeolite ultrafiltration membrane hybrid system[END_REF] pseudo-second-order and the maximum amount adsorbed was best simulated by the Langmuir isotherm model [START_REF] Kumar | Application of cellulose-clay composite biosorbent toward the effective adsorption and removal of chromium from industrial wastewater[END_REF].

Conclusion

In this work, the potential of NP nanoparticles for the removal of chromium from aqueous solutions was studied. The effect on chromium removal efficiency by NP of the different parameters, such as adsorbent dose, pH, contact time, temperature and initial concentration, was investigated. It was found that the adsorption mechanism followed the Langmuir isotherm model. The monolayer adsorption capacity of NP was of 35 mg/g, adsorbent dose = 1.6 g/L, pH = 2.9, contact time = 75 min, T = 298 K and C 0 = 80 mg/L. The kinetic study of the equilibrium data showed that the removal of chromium ions followed best the pseudosecond-order model. The thermodynamic study indicated the feasibility, exothermic and spontaneous adsorption of chromium onto NP at 298-308 K. Compared to other nanocomposites, NP is a promising material for chromium removal from industrial wastewaters.

Fig. 1 .

 1 Fig. 1. Preparation of nanoparticles (NP) nanostructured materials based on sol-gel method.

Fig. 2 .

 2 Fig. 2. Effect of adsorbent dose (C 0 = 80 mg/L; V = 50 mL; pH = 2.9; t = 75 min; T = 298 K).

Fig. 3 .

 3 Fig. 3. Distribution of Cr(VI) species in water as a function of pH level and Cr(VI) concentration (a) and effect of pH (C 0 = 80 mg/L; V = 50 mL; t = 75 min; T = 298 K) (b).

Fig. 4 .Fig. 5 .

 45 Fig. 4. Effect of contact time (C 0 = 80 mg/L; m = 1.6 g; V = 50 mL; pH = 2.9; T = 298 K).

Fig. 6 .

 6 Fig. 6. Effect of initial Cr(VI) concentration (m = 1.6 g; V = 50 mL; pH = 2.9; t = 75 min; T = 298 K).

Fig. 7 .

 7 Fig. 7. Pseudo-first-order (a), pseudo-second-order (b) and intraparticle diffusion (c) kinetic plot for the removal of chromium by NP.

Fig. 9 .

 9 Fig. 9. Plot of lnK d vs. 1/T for estimation of thermodynamic parameters for adsorption of removal chromium by NP.

Fig. 8 .

 8 Fig. 8. Isotherm plots for adsorption of chromium onto NP (m = 1.6 g; V = 50 mL; t = 75 min; C 0 = 80 mg/L; T = 298 K: Langmuir (a) and Freundlich (b).

Table 1

 1 Calculated kinetic parameters for pseudo-first-order, pseudo-second-order and intraparticle diffusion for the removal of chromium by NP

	Temperature	q e,exp	Pseudo-first-order	Pseudo-second-order		Intraparticle diffusion
	(K)	(mg/g)	k 1 (min -1 )	q e (mg/g)	R 2	k 2 (g/mg•min)	q e (mg/g)	R 2	k 3 (mg/g•min)	C	R 2
	298	35	0.02	3.896	0.945	0.325	34.48	0.999	0.541	29.24	0.953

Table 2 Isotherm

 2

		Langmuir	Freundlich
	R 2	0.997	R 2	0.985
	q m (mg/g)	21.73	1/n	0.54
	K L (L/mg)	0.103	K F (mg/g)	192
	R L	0.108	-	-

model parameters for removal of chromium by NP (m = 1.6 g; V = 50 mL; pH = 2.9; t = 75 min; C 0 = 80 mg/L; T = 298 K)

Table 3

 3

	Thermodynamic parameters for chromium adsorption onto NP
	T (K)	ΔG° (kJ/mol)	ΔH° (kJ/mol)	ΔS° (J/mol•K)
	298	-2.93
	303	-2.501	-18.75	-53.19
	308	-2.409

Table 4

 4 Water treatment for removal Cr(VI) from different inorganic contaminants over clay-polymers nanocomposites (NCs) adsorbents

	Adsorbent	Temperature	pH	Adsorption	Type of	Isotherm	Kinetics model	References
		(°C)		capacity	solution	model
	Nanoparticles (NP)	25	2.9	35 (mg/g)	Industrial	Langmuir Pseudo-second-order This study
					wastewater
	Chitosan/attapulgite	45	5	65.37 (mg/g) Aqueous	Langmuir Intraparticle diffusion [37]
					solution
	Chitosan-Al-pillared	25	6.38	15.68 (mg/g) Aqueous	Langmuir Pseudo-second-order [38]
	montmorillonite				solution
	Na-montmorillonite/	-	3.8-5.5 22.3 (mg/g)	Aqueous	Langmuir Pseudo-second-order [39]
	cellulose				solution
	Vesicular basalt	25	2	79.20 mg/kg Aqueous	Freundlich Pseudo-second-order [14]
	volcanic rock				solution
	Activated carbon	25	2	72.46 mg/L	Aqueous	Langmuir Intraparticle diffusion [23]
					solution
	Oxalate modified weak	20	5	201.30 mg/L Aqueous
	basic resin (IO@D301)				solution

Acknowledgments

This paper is supported by the PRIMA program under grant agreement No2024 -TRUST project. The PRIMA program is supported by the European Union. This work was also supported by the Tunisian Ministry of Higher Education and Scientific Research (DGRST) TR in the frame of the SETPROPER research project (grant n°13-043, ERANETMED_WATER program).

Disclosure statement

No potential conflict of interest was reported by the authors. S. Khemakhem, Elaboration and characterization of ceramic microfiltration membranes from natural zeolite: application to the treatment of cuttlefish effluents, Desal. Water Treat., 95 (2017) 9-17.

Author Query

AQ1

Kindly provide the Manufacturer Location for the company "PerkinElmer spectrometer", "Micromeritics ASAP 2020".


Damien Leduc
email: damien.leduc@univ-lehavre.fr

Bruno Morvan

Pascal Pareige

Jean-Louis Izbicki

Measurement of the effects of rough surfaces on Lamb waves propagation

Keywords: 43.35.Cg, 43.35.Zc Lamb waves, rough surface, surface topography

This work deals with the sensitivity to the plate roughness of Lamb waves. An experimental study is performed involving an air-coupling transducer system. Signal processing allows us to extract the Lamb waves characteristics: phase velocity and attenuation. Plate surface topographies are obtained by means of an optical surface profiler. The acoustic characteristics and the surface topographies are finally linked.

Introduction

The use of ultrasonic waves to control adhesive bonding is widespread nowadays (see e.g. [START_REF] Guyott | The non-destructive testing of adhesively bonded structure: a review[END_REF]). In most cases of an adhesive joint, the substrate surface is treated before application of the adhesive. The goal is to increase the interfacial energy in the contact region. There are many different chemical or mechanical treatments; most of them consist in increasing the substrate surface roughness. Indeed, the surface roughness increases the area of physical contact and improves the adhesion strength. Estimation of the joint quality by means of acoustic methods cannot be performed without understanding the roughness effect on the ultrasonic waves propagation. In a bond structure, the plate rough side is located at the interface and there is of course no access to this surface for direct evaluation with surface waves (i.e. Rayleigh wave type). That is why Lamb waves seem accurate for future study. While these waves will be generated from the outer side of the plate, they may provide us with the bonding inner surface.

Many studies deal with roughness influence on ultrasonic waves reflection at a fluidsolid interface ( [START_REF] Nagy | Surface roughness induced attenuation of reflected and transmitted ultrasonic waves[END_REF], [START_REF] Boyd | Forward and specular scattering from a rough surface: theory and experiment[END_REF]). Few of them focus on guided waves propagation ( [START_REF] Nagy | Non destructive evaluation of adhesive joints by guided waves[END_REF], [START_REF] Kundu | A Lamb wave scanning approach for the mapping of defects in [0/90] titanium matrix composites[END_REF]).

Brekhovskikh [START_REF] Brekhovskikh | Propagation of surface Rayleigh waves along the uneven boundary of an elastic body[END_REF] has formulated an analysis that predicts the attenuation of Rayleigh waves travelling on a periodic rough surface over a particular frequency band. Chimenti and Lobkis [START_REF] Chimenti | The effect of rough surfaces on guided waves in plates[END_REF] show that guided waves damping increases with surface roughness and frequency for an immersed plate. However roughness has no influence on the phase velocity of Lamb waves.

This work is an experimental study involving air-coupling measurement of Lamb waves on rough plates. Signal analysis allows us to extract, at a given frequency, the real ( ' x k ) and imaginary ( ''

x k ) components of the wave vector of the Lamb waves propagating along the x direction. The real part '

x k of the wave number gives the wave phase velocity whereas its imaginary part ''

x k gives the amplitude attenuation. We will study these parameters evolution versus the roughness parameter.

In a first part, the surface topography of many samples is studied. Then acoustic measurements are reported and discussed in the next section. Finally a brief conclusion is given.

Surface topography

Characterization of the studied samples

The studied samples are glass plates with 200 mm sides and a thickness d = 5 mm. The measured acoustics parameters are the glass density: 2490 kg m -3 and the shear and longitudinal velocities: cT = 3485 m s -1 and cL = 5825 m s -1 . These plates are treated on one side only in order to obtain a rough surface except one that is considered as a reference plate. The other surface is roughless and exactly the same for all samples: this surface stands for a reference surface. Two quite different techniques are used to create surface roughness: sanding and shot blasting. We may observe the studied surface topography by means of an optical surface profiler. The analysed surface is about 25 mm 2 .The resolution is 5 µm in x and y directions and 10 nm in z direction.

Quantitative surface evaluation by optical profiler

We characterize the roughness by use of the usual statistic parameters Ra, Rq and Rt [START_REF]ASME B46.1-1995 Surface texture, surface roughness, waviness and lay[END_REF].

The Ra is called the roughness average; it is the arithmetic average of the absolute values of the measured height deviations. Height deviations Zij are measured from the mean surface. Ra mathematical expression is:

11 1 MN a ij ji RZ MN == =  (1)
where M and N are the number of data points in each direction of the analysed surface.

Rq, also called the root mean square (RMS) roughness, is obtained by squaring each Zij value over the evaluation area, then taking the square root of the mean:

2 11 1 MN q ij ji RZ MN == =  (2)
Rt is the maximum peak-to-valley height; it is the vertical distance between the highest and lowest points of the surface. The statistic parameters of the four studied samples are reported in table 1. These mean values are calculated from six acquisitions on each plate providing from different zones. As shown in table 1, roughness values cover a large range of Rq.

Acoustic measurements

Experimental set-up

The experimental set-up is reported in figure 1. A pulse generator delivers a very short pulse voltage (about 300 V during 300 ns) to a Panametrics emitting piezo-composite transducer. Lamb waves are generated by the wedge method. The receiver transducer is a Sofratest air-coupling piezo-electric transducer. The air-coupling transducer detects the re-radiated signal. The emitter transducer remains unmoved whereas the receiver is translated along the propagation direction of Lamb waves. The use of a receiver aircoupling transducer is very interesting for such a study. Experiments are repeatable because we avoid all the problems due to the contact of the transducer with the sample (to maintain unchanged the thickness of coupling gel and the pressure on the transducer during the measurement). Moreover, air-coupling measurements enable motorized translation and automatic data acquisition.

The two transducers are set on the non-rough side. The displacement amplitudes are collected from x = 10 mm to x = 90 mm by 0.1 mm step (origin x = 0 corresponds to the wedge position). For each position of the air-coupling transducer, a 200 µs long signal is acquired on 10000 points. In order to improve the signal to noise ratio, we performed an average of 1000 successive shots.

Signal processing

Lamb waves '

x k evolution for different surface roughness

The previous study for reference sample is now applied on the three plates with increasing roughness. We notice that, for the roughest sample, no mode is excited above x k values of the Lamb waves for each studied sample and reported them in figure 3a. At a given frequency, the experimental '

x k for each sample are close to one another. Figure 3b shows the results for the S1 mode: the shift between the theoretical and the experimental '

x k increases with the roughness.

The relative difference defined by:

' ' ' '' ( ) ( ) ( ) ( ) x x x xx k k rough k ref k ref k ref - = (3)
between rough and reference samples is reported in table 2. This difference does not significantly change from a mode to one another, in other words it is independent to frequency shifts. It increases however with roughness even if its variation is not linear.

Indeed the Rq coefficient for the shot blasting plate is 5 times higher than the sanding plate (table 1) whereas the relative difference is almost the same. On the contrary the Rq value between the shot blasting and the strong shot blasting plates is twice greater and the relative difference increases from 2 % to 4.6 %. For the two smoothest samples, the relative difference variation is in the range of detectability. For the roughest plate, the value of the relative difference is outside the experimental error domain and it can be concluded that, for a large value of roughness, an effect on the value of '

x k exists.

Nevertheless the value of '

x k does not look like a good criterion to discriminate different roughness. In the next paragraph, we will study the effect of roughness on Lamb waves attenuation.

Comparison of Lamb waves attenuation

In order to isolate one mode for its further characterisation, a spatial Short Time Fourier Transform (STFT) is achieved. First a temporal Fast Fourier Transform of the signal s(x,t) is performed providing a ( , ) t s x f signal in the spatial-frequency domain.

Then, the frequency corresponding to the studied mode is isolated using a spatial window sliding in the entire spatial domain (80 mm). A spatial FFT is computed for each position and provides a , ( , ) xt s k f signal. Thus the amplitude attenuation versus distance between the emitting and receiving transducer is visualized in the (x,k) space.

After the STFT treatment, the amplitude for each mode is obtained and plotted in figure 4 and 5 versus propagation distance. On the figures, the level one corresponds to initial position. For a given mode, this normalization allows us to compare the attenuation for different plates.

Firstly we focus on a particular mode and we study the evolution of the wave amplitude versus the propagation distance x for different roughness (figure 4). This evolution is almost linear when the roughness is weak (Rq < 6 µm). Attenuation exists for the reference sample. This attenuation is close to the sanded sample one. It is mainly due to the diffraction phenomena occurring by the limited size of the emitter transducer. It has to be noticed that a small part of the energy is also re-radiated in the surrounding medium.

This attenuation can be qualified as "apparent". With frequency increasing, the effect of the diffraction diminishes nevertheless this attenuation still exists for the rough samples.

For a given frequency, the attenuation increases with increasing roughness (figure 4a, 4b and 4c). This fact has been experimentally verified for each of the considered mode.

Secondly the amplitude evolution of a selected roughness for different modes is shown in figure 5. We note that the higher frequency the higher the attenuation. By example, for a propagation distance of 70 mm, the attenuation varies from 3.5 dB (at 5.9 MHz mm) to 6.9 dB (at 13.1MHz mm) for the shot blasted plate. The sensitivity to the roughness is then stronger in high frequency range.

As a partial conclusion, contrary to the '

x k value, the attenuation seems to be a good parameter in order to discriminate between different roughness.

Conclusion

After determining first the surface topography, secondly the effect of roughness of ultrasonic waves, we join the two parts of our work. Whereas the '

x k value is only slightly affected by the Rq parameter as discussed in section 3.2.1, the attenuation is selected to be an efficient acoustic parameter. Thus we plot the attenuation versus Rq statistic parameter for a given Lamb mode (figure 6). This representation links explicitly acoustic parameters and surface topography. Obviously a larger study involving more samples with a larger range of roughness would be necessary to obtain calibration curves.

Future works will concern adhesive interfaces in order to estimate the influence of the roughness on the bond quality.

Table caption

Table 1: Roughness statistic parameters for the studied samples.

Table 2: '

x k relative difference values between rough and reference samples.

For the strong shot blasted sample, the S2 mode is too attenuated. R a (µm) R q (µm) R t (µm) a.

Figure caption

Figure 2

 2 Figure2is the result of a spatial followed by a temporal Fast Fourier Transform[START_REF] Alleyne | A two-dimensional Fourier transform method for the measurement of propagating multimode signals[END_REF] of

  MHz mm = . This is probably due to the attenuation of the high frequencies modes propagating in rough plates. At high frequencies, the wavelengths of the Lamb modes become commensurate with the size of the roughness. A multi-scattering process can occur and consequently the propagation of the Lamb modes cannot be detected. That is why the ' x k comparison is performed only on the three first modes A1, S1 and A2. The experimental values of ' x k are the maximum values of signals in (k,f) representation. In that way, we have evaluated '

Figure 1 :

 1 Figure 1: Experimental set-up.

Figure 2 :

 2 Figure 2: Representation of the signal in the dual space (k,f) for reference sample.The lines correspond to the theoretical dispersion curves: solid lines for antisymmetric modes and dotted lines for symmetric modes. Experimental data are in grey level. They practically coincide with the theoretical data.

Figure 3 :

 3 Figure 3: Comparison of the ' x k for different samples: a. For given mode, the ' x k value are almost unchanged whatever the sample. b. Zoom on the S1 mode ' x k values.

Figure 4 :

 4 Figure 4: Evolution of Lamb waves attenuation versus propagation distance: a. for the A1 mode, b. for the S1 mode, c. for the A2 mode.

Figure 5 :

 5 Figure 5: Evolution of Lamb waves attenuation for the shot blasted plate at several frequencies.

Figure 6 :

 6 Figure 6: Evolution of the acoustic normalized amplitude versus Rq surface parameter for different propagation distances (x = 30, 50 or 70 mm): a. A1 mode, b. S1 mode, c. A2 mode.

  Fig.1

  Fig.3

Acknowledgement

The authors would like to thank Pascal Rembert for Matlab help, Liot Techni-verres for providing samples and CRITT Analyses and Surfaces for the surface topography.


Ivan Fedioun
email: ivan.fedioun@univ-orleans.fr

Two-fluids plane or axisymmetric Poiseuille flow at optimum A (quasi)analytical study

Context of the study

This work was initiated in the context of a final project by engineering students at Polytech Orléans, devoted to the optimization of the suction function of a medical endoscope. The manufacturer observed that the flow rate of aspirated body fluid has a maximum when an air inlet is present in the circuit. However, if the air intlet is too large, the performance of the device drops. Liquid-gas two-phase flows are very complex1 , and it is not intended here to address the problem in its entirety. We propose a simple analytical study based on the axisymmetric or plane Poiseuille flow model, which allows to explain to some extent the observed phenomenon.

Classical Poiseuille flow : reminder

The plane or axisymmetric Poiseuille flow is one of the very few analytical solutions of Navier-Stokes equations, a classical exercice for any student in fluid mechanics2 . Under the asumptions of a steady and parallel flow of incompressible fluid (ρ =cste), the problem to be solved is (2)

dp dx = µ∆u = K < 0 (
• In the axisymmetric case (0 ≤ r ≤ R), the general solution of (1) is

u(r) = K µ r 2 4 + A ln r + B (3)
The finite velocity on the axis (r = 0) imposes A = 0, and the no-slip condition u(±R) = 0 gives B = -R 2 /4. The solution is then

u(r) = U 1 - r R 2 (4) U = - K µ R 2 4 > 0 (5)
where U = 2U ave is the maximum velocity along the axis of the pipe. The volume flow rate is

Q POIS axi = πR 2 U ave = -2π K µ R 4 16 (6)
• In the plane case (-H ≤ y ≤ H), the general solution of (1) is

u(y) = K µ y 2 2 + Ay + B (7)
The maximum velocity in the symmetry plane (y = 0) imposes A = 0, and the no-slip condition u(±H) = 0 gives B = -H 2 /2. The solution is then

u(y) = U 1 - y H 2 (8) U = - K µ H 2 2 > 0 (9)
where U = 3 2 U ave is the velocity in the symmetry plane. The volume flow rate (per m) is

Q POIS plane = 2HU ave = - 2 3 K µ H 3 (10)
3 Axisymmetric two-fluids Poiseuille flow

General solution

Let us consider a two-fluids axisymmetric configuration, fluid 1 being in the center of the pipe, and fluid 2 in an annular configuration around fluid 1 (figure 1). Fluid 1 is therefore not constrained by the no-slip condition at the wall. The radius 0 ≤ R 12 ≤ R is at the interface of the two fluids assumed non-miscible.

In a cross section of the pipe, both fluids are at the same pressure. The problem reads : and the general solution is

dp dx = µ 1 r d dr r du 1 dr = µ 2 r d dr r du 2 dr = K < 0 (11)
u α (r) = K µ α r 2 4 + A α ln r + B α ; α = 1, 2 (12)
Boundary conditions are :

r = 0 : du 1 dr = 0 (13) r = R 12 : u 1 (R 12 ) = u 2 (R 12 ) = U 12 (14) µ 1 du 1 dr = µ 2 du 2 dr (15) r = ±R : u 2 (±R) = 0 (16)
We deduce the constants :

A 1 = A 2 = 0 ; B 1 = - R 2 -R 2 12 4η - R 2 12 4 ; B 2 = - R 2 4 where η = µ 2 µ 1 > 0 (17)
is the first physical parameter of the problem. If η < 1, the most viscous fluid (the "liquid") is in the center, side 1, and vice versa. Velocity profiles are :

u 1 (r) = - K µ 1 R 2 12 -r 2 4 + R 2 -R 2 12 4η ; 0 ≤ r ≤ R 12 (18) u 2 (r) = - K µ 2 R 2 -r 2 4 ; R 12 ≤ r ≤ R (19)
The volume flow rates of fluid 1 and fluid 2, functions of R 12 , are obtained from

Q 1 (R 12 ) = 2π R 12 0 u 1 (r) rdr ; Q 2 (R 12 ) = 2π R R 12 u 2 (r) rdr (20) Q 1 (R 12 ) = -2π K µ 1 R 4 12 16 + R 2 12 8η R 2 -R 2 12 (21) Q 2 (R 12 ) = -2π K µ 2 R 4 + R 4 12 16 - R 2 R 2 12 8 (22)
In the case η = 1, one can check that

Q 1 + Q 2 = Q POIS axi
where Q POIS axi is the flow rate (6) of the classical axisymmetric Poiseuille flow, and if

R 12 = R (resp. R 12 = 0), Q 1 = Q POIS 1 axi (resp. Q 2 = Q POIS 2 axi ).

Search for the optimum

The problem is, given the viscosity ratio η, to find the radius R 12 -if it exists-that maximizes the flow rate of the most viscous fluid, i.e. the flow rate

Q 1 if η < 1 or the flow rate Q 2 if η > 1.
The flow rate Q 2 (22) does not depend on η. Therefore, the search for optimum only makes sense for η ≤ 1, i.e. for Q 1 . Unfortunately, the most physically probable case is η >> 1 : a liquid film on the wall drained by the faster and less viscous gas in the center.

Let us define :

r = R r ; 0 ≤ r ≤ 1 (23) R 12 = R ρ ; 0 ≤ ρ ≤ 1 (24) Q α (R 12 ) = Q POIS α Q α (ρ) ; α = 1, 2 (25)
u α (r) = U 1 u α (r) (26)
The reference velocity must be the same for both flows. We choose the maximum velocity U 1 (5) of the most viscous internal fluid. The flow rates are normalized by those of the respective Poiseuille flows (6). It gives :

u 1 (r) = ρ 2 -r 2 + 1 η (1 -ρ 2 ) ; 0 ≤ r ≤ ρ < 1 (27) u 2 (r) = 1 -r 2 /η ; 0 < ρ ≤ r ≤ 1 (28) Q 1 (ρ) = 2 η ρ 2 (1 -ρ 2 ) + ρ 4 (29) Q 2 (ρ) = 1 -2ρ 2 + ρ 4 (30)
One can see that the second physical, or rather geometrical, parameter on which the problem depends is

ρ = R 12 R (31)
from which the volume fractions of both fluids are straightforward :

X 1 = ρ 2 ; X 2 = 1 -ρ 2 (32)
The normalized flow rates Q 1 (ρ) and Q 2 (ρ) are displayed in figure 2. Indeed, the flow rate Q 1 (ρ) shows a maximum that moves from ρ = 1 to ρ ≈ 0.7 when η decreases. This finding is consistent with the experimental obvervations described in the introduction.

In order to find the maximum flow rate, one writes dQ 1 dρ = 0, which leads to

ρ = 0 ; Q 1 = 0 trivial minimum 1 -(2 -η)ρ 2 = 0 ; ρOPT = 1 √ 2 -η Q 1 maximum sought (33) with limit 1/ √ 2 ≈ 0.707 • • • when η → 0.
The flow rates at optimum are then

Q OPT 1 = 1 η(2 -η) ; Q 2 = 1 -η 2 -η 2 (34) Q1 U K=1 0.5 0.2 0.1 Q2 Figure 2: Normalized flow rates Q 1 (ρ) (29)
and Q 2 (ρ) (30). It can be seen that a viscosity ratio of 4 (η ≤ 0.25) more than doubles the flow rate of the inner most viscous fluid. All these results are analytical.

Velocity profiles

Velocity profiles are shown in figure 4. One can see that the velocity of the "liquid" can reach much higher values than that of the classical Poiseuille flow, by effect of entrainment by the annular gas flow. The maximum velocity along the axis is

u 1 (0) = 1 η(2 -η) = Q OPT 1
and can theoretically reach infinity in the limit η → 0. The problem reads :

dp dx = µ 1 d 2 u 1 dy 2 = µ 2 d 2 u 2 dy 2 = K < 0 (35)
and the general solution is

u α (y) = K µ α y 2 2 + A α y + B α ; α = 1, 2 (36)
Boundary conditions are :

y = -H : u 1 (-H) = 0 (37) y = H 12 : u 1 (H 12 ) = u 2 (H 12 ) = U 12 (38) µ 1 du 1 dy = µ 2 du 2 dy (39) y = +H : u 2 (+H) = 0 (40)
We deduce the constants :

A 1 = A 2 = A 2 = (η -1)(H 2 -H 2 12 ) 2 [(η + 1)H + (η -1)H 12 ] (41)
B 1 = H 2 (A -H) ; B 2 = - H 2 (A + H)
Velocity profiles are :

u 1 (y) = - K µ 1 H 2 -y 2 2 - A 2 (H + y) ; -H ≤ y ≤ H 12 (42) u 2 (y) = - K µ 2 H 2 -y 2 2 + A 2 (H -y) ; H 12 ≤ y ≤ H (43)
The volume flow rates of fluid 1 and fluid 2, functions of H 12 , are obtained from

Q 1 (H 12 ) = H 12 -H u 1 (y)dy ; Q 2 (H 12 ) = H H 12 u 2 (y)dy (44) Q 1 (H 12 ) = - K µ 1 H 2 (H + H 12 ) 2 - H 3 + H 3 12 6 - A 2 H(H + H 12 ) - H 2 -H 2 12 2 (45) Q 2 (H 12 ) = - K µ 2 H 2 (H + H 12 ) 2 - H 3 -H 3 12 6 + A 2 H(H -H 12 ) - H 2 -H 2 12 2 (46)

Search for the optimum

Again, given a value 0 < η ≤ 1, we want to find the height H 12 that maximizes the flow rate Q 1 of the most viscous fluid (the problem is symmetric).

Let us define :

y = H y ; -1 ≤ y ≤ 1 (47)
H 12 = H ρ ; -1 ≤ ρ ≤ 1 (48) Q α (H 12 ) = Q POIS α Q α (ρ) ; α = 1, 2 (49)
u α (y) = U 1 u α (y) (50)
The reference velocity is the maximum velocity U 1 (9) of the most viscous fluid. Flow rates are normalized by those of the respective Poiseuille flows (10). It gives :

u 1 (y) = 1 -y 2 -A(1 + y) ; -1 ≤ y ≤ ρ < 1 (51) u 2 (r) = 1 -y 2 + A(1 -y) /η ; -1 < ρ ≤ y ≤ 1 (52) Q 1 (ρ) = 3 4 1 + ρ - 1 3 (1 + ρ 3 ) - A 2 (1 + ρ) 2 (53) Q 2 (ρ) = 3 4 1 -ρ - 1 3 (1 -ρ 3 ) + A 2 (1 -ρ) 2 (54)
with

A(ρ) = (η -1)(1 -ρ 2 ) η + 1 + (η -1)ρ (55)
The volume fractions of both fluids are

X 1 = 1 2 (1 + ρ) ; X 2 = 1 2 (1 -ρ) (56)
The function A(ρ) with parameter η, from which the solutions differ compared to the classical Poiseuille flow, is shown in figure 6. Its limit when η → 0 is A η→0 (ρ) = -(1 + ρ), for -1 ≤ ρ < 1 and A η→0 (1) = 0. The normalized flow rates Q 1 (ρ) and Q 2 (ρ) are displayed in figure 7. Flow rate Q 1 (ρ) shows a maximum that moves as η decreases, but non-monotonically unlike the axisymmetric case (figures 2 and 3).

The optimum flow rate 1 is obtained from dQ 1 dρ = 0, that leads to

ρ = -1 ; Q 1 = 0 trivial minimum 1 -ρ -A - A ′ 2 (1 + ρ) = 0 = f 1 (ρ) ; Q 1 maximum sought (57)
where A ′ (ρ) is the derivative of ( 55) with respect to ρ. Symmetrically, we get for flow rate 2 :

ρ = +1 ; Q 2 = 0 trivial minimum 1 + ρ + A - A ′ 2 (1 -ρ) = 0 = f 2 (ρ) ; Q 2 maximum sought ( 58
)
The optimization function f 1 (ρ) is shown in figure 8 for different values of η. Its roots ρ OPT 1 (η) are not analytical. They are computed with the Newton-Raphson method (absolute convergence at 10 -12 ). They are plotted in figure 9 and some numerical values are given in table 2 in appendix. It is remarkable in figure 9 that at the specific point where ρ OPT 1 is minimum, the optimum flow rate is

Q OPT 1 ≈ √ 2 and ρ OPT 1 ≈ √ 2/2 .
If we accept these values, we can get after some algebra the theoretical corresponding viscosity ratio η and flow rate Q 2 :

η = 1 5 47 -24 √ 2 3 + 4 √ 2 ≈ 0.30170023 ; Q 2 = 3 2 3 -2 √ 2 3 + 2 √ 2 ≈ 0.04415588
With this theoretical value of η, the Newton-Raphson method converges to ρ OPT

1 = 0.70721664 • • • ≈ √ 2/2 within 1.1 10 -4 (-0.015%), Q OPT 1
≈ √ 2 within 3.8 10 -8 (-0.0000027%), and Q 2 = 0.044128 • • • , which is again very close to the theoretical value. It was not possible to further investigate this issue, neither analytically nor numerically, but this specific point reminds us the asymptotic limit √ 2/2 of the optimal radius (33) of the axisymmetric case. Asymptotic behaviors of different variables for limit values of η are gathered in table 1.

η 0 (-1 ≤ ρ < 1) 1 +∞ (-1 < ρ ≤ 1) A -(1 + ρ) 0 1 -ρ A ′ -1 0 -1 f 1 (ρ) 2 + (1 + ρ)/2 1 -ρ (1 + ρ)/2 f 2 (ρ) (1 -ρ)/2 1 + ρ 2 + (1 -ρ)/2 ρ OPT 1 1 1 - Q OPT 1 4 1 -
Table 1: Limit behaviors of solutions.

Velocity profiles

Velocity profiles are shown in figure 10. As in the axisymmetric case, one can see that the "liquid", drained by the less viscous gas, reaches higher velocities than that of the classical Poiseuille flow. The maximum velocity is obtained for η → 0 :

Conclusions and perspectives

This analytical study is of course very schematic and does not pretend to represent all the complexity of two-phase flows in pipes. However, it has allowed us to identify a mechanism that can explain the experimental observations. To complete it, one could introduce the surface tension forces at the liquid-gas interface. The axisymmetric case is unfortunately optimal for η < 1 only (liquid in the center, gas around). Optimization of the η > 1 case (gas in the center, liquid at the walls) would have been more interesting for industrial applications. However, the flow rate of the central liquid can be theoretically infinite in the limit η → 0.

The plane case is no more analytical. The maximum reachable flow rate is 4 times that of the classical Poiseuille flow in the limit η → 0. The specific point ρ

OPT 1 = √ 2/2, Q OPT 1 = √ 2 is striking, but is not an exact optimum.
The two-phase "slug" type of flow, where long bubbles of gas are drained at the center of the pipe, then collapse and split again, could be modeled alternating η > 1 ↔ η < 1 with an intermittency factor to be designed more or less empirically, in order to reproduce the experimental optimum. The question arises of course as how to achieve this optimum for a given viscosity ratio η, i.e. how to concretely realize the air intake in the pipe (parietal, axial...) and how to calibrate its flow rate, if not by the practitioner's feeling.

Appendix

η ρ OPT 1 Q OPT 1 Q OPT 2 X 1 X 2 0 1 4 0 1 0 0.

Figure 1 :

 1 Figure 1: Two-fluids axisymmetric Poiseuille flow configuration.

Figure 3 :

 3 Figure 3: Optimum normalized radius (33), flow rates (34) and volume fractions.

Figure 4 :

 4 Figure 4: Velocity profiles that maximize the volume flow rate Q 1 (η < 1).

Figure 5 :

 5 Figure 5: Two-fluids plane Poiseuille flow configuration.

Figure 6 :

 6 Figure 6: Function A, eq.(55).

Figure 7 :

 7 Figure 7: Flow rates Q 1 (ρ) and Q 2 (ρ).

Figure 8 :

 8 Figure 8: Optimization function (57).

Figure 9 :

 9 Figure 9: Optimum normalized height, volume flow rates and volume fractions.

Figure 10 :

 10 Figure 10: Velocity profiles that maximize the volume flow rate Q 1 (η < 1).

Figure 11 :

 11 Figure 11: VOF ANSYS/Fluent simulation of a plane (2D) "slug" flow (Delavault A. / Maberi-Grodet M., 2023)

Table 2 :

 2 Numerical roots of (57), and corresponding optimum solutions. In bold, numerical solution at the remarkable point, very close to ρ OPT

	0001	0.988607239	3.898045090	0.000001856	0.99430362	0.00569638
	0.0005	0.974943430	3.777271089	0.000019820	0.987471714	0.012528286
	0.001	0.988607239	3.898045090	0.000054198	0.994303619	0.005696380
	0.005	0.925693963	3.354808198	0.000526320	0.962846981	0.037153019
	0.01	0.898894644	3.132626433	0.001341430	0.949447322	0.050552678
	0.02	0.864583746	2.854801700	0.003280031	0.932291873	0.067708127
	0.05	0.807820625	2.405606500	0.009748029	0.903910312	0.096089688
	0.1	0.759984915	2.021779694	0.020006793	0.879992457	0.120007543
	0.2	0.718661678	1.631679570	0.035428058	0.859330839	0.140669161
	0.25	0.710431818	1.511426731	0.040530850	0.855215909	0.144784091
	0.3	0.707250263	1.417070415	0.044035615	0.853625132	0.146374868
	0.30170023 0.707216642 1.414213601 0.044128111 0.853608321 0.146391679
	0.35	0.708055811	1.340840836	0.046059677	0.854027905	0.145972095
	0.4	0.712221090	1.278007399	0.046714672	0.856110545	0.143889455
	0.45	0.719368769	1.225491404	0.046106818	0.859684385	0.140315615
	0.5	0.729281773	1.181191183	0.044340018	0.864640887	0.135359113
	0.55	0.741856425	1.143621618	0.041520381	0.870928213	0.129071787
	0.6	0.757077754	1.111707824	0.037761984	0.878538877	0.121461123
	0.65	0.775007768	1.084660856	0.033194269	0.887503884	0.112496116
	0.7	0.795782520	1.061900134	0.027971887	0.897891260	0.102108740
	0.75	0.819616365	1.043004186	0.022288332	0.909808182	0.090191818
	0.8	0.846813471	1.027679729	0.016395568	0.923406736	0.076593264
	0.85	0.877788146	1.015743714	0.010633270	0.938894073	0.061105927
	0.9	0.913097411	1.007115719	0.005473940	0.956548705	0.043451295
	0.95	0.953492312	1.001820172	0.001595138	0.976746156	0.023253844
	1	1	1		0			1	0
			1	=	√	2/2, Q	OPT 1	=	√ 2.

See for example : Jean-Michel Fitremann, Ecoulements diphasiques gaz-liquides, Les techniques de l'Ingénieur Ref. A722 V1, 1983

Schlichting H. Boundary-Layer Theory, 6th Ed., Mc Graw Hill, 1968