Pierre Horgue
email: horgue@imft.fr

Frédéric Augier

Michel Quintard

Marc Prat

A suitable parametrization to simulate slug flows with the Volume-Of-Fluid method

Keywords:

CFD Thin film VOF method

Diffuse-interface methods, such as the Volume-Of-Fluid method, are often used to simulate complex multiphase flows even if they require significant computation time. Moreover, it can be difficult to simulate some particular two-phase flows such as slug flows with thin liquid films. Suitable parametrization is necessary to provide accuracy and computation speed. Based on a numerical study of slug flows in capillary tubes, we show that it is not trivial to optimize the parametrization of these methods. Some simulation problems described in the literature are directly related to a poor model parametrization, such as an unsuitable discretization scheme or too large time steps. The weak influence of the mesh irregularity is also highlighted. It is shown how to capture accurately thin liquid films with reasonably low computation times.

Introduction

When cocurrent gas and liquid flows take place within a capillary tube, there is a large number of interface configurations, called flow regimes or patterns in the literature, depending on the fluid properties and flow rates which can be plotted on a flow regime map (see Fig. 1 from the work of Triplett et al. [START_REF] Triplett | Gas-liquid two-phase flow in microchannels, Part 1: Two-phase flow patterns[END_REF]). One of the most common regime for the lowest liquid flow rates is the Taylor flow, also called slug flow, characterized by the presence of elongated bubbles, with diameters close to the channel tube width, separated with liquid slugs. For a higher gas flow rate, we can observe the annular regime characterized by the flow of a liquid film along the wall while the gas phase occupies the central void space. The common point between these flow regimes is the presence of a liquid film along the wall, which can be very thin and thus pose a serious numerical problem. Taylor flows have interesting flow characteristics such as recirculation flows in the liquid slug which improve heat and mass transfer, and, for this reason, they have been widely studied experimentally and numerically.

Previous numerical studies showed that one of the main difficulties in Taylor flow modeling concerns the simulation of the thin liquid film along the wall. Most of the previous works reported that the liquid film was not captured because the grid was not sufficiently refined. More recently, Gupta et al. [START_REF] Gupta | On the cfd modelling of Taylor flow in microchannels[END_REF] showed that it was possible to capture these liquid films, even thin, provided that an ultra fine mesh grid and much computation time where used. Moreover, previous works highlighted that using diffuse-interface methods can induce some numerical errors such as velocity field fluctuations around the interfaces called "spurious currents" in the literature.

In the present Note, we simulate with the Volume-Of-Fluid method (VOF), using the many possibilities offered by FLUENT TM , two cases of slug flows with the presence of a thin liquid film along the wall: the drainage of a liquid phase [START_REF] Triplett | Gas-liquid two-phase flow in microchannels, Part 1: Two-phase flow patterns[END_REF]) and summarized in this figure by Gupta et al. [START_REF] Gupta | On the cfd modelling of Taylor flow in microchannels[END_REF]. in a vertical tube and the Taylor bubble generation in a horizontal tube. While checking the overall accuracy, a special attention is given to the liquid film modeling which have received less attention in the literature as pointed out above. Performing a parametric study on irregular meshes, VOF discretization schemes and Courant Number conditions, we first show, in the drainage case, that suitable parametrization allows to capture accurately flow properties, including the thin liquid films, while conserving small computation times and avoiding numerical errors. Applying these parameters in the Taylor bubble generation case, we show that, more generally, recommendations obtained in the previous elementary case can be extrapolated to all slug flows.

Numerical background

Slug flows with the presence of thin liquid films have been studied by different numerical method such as Volume-Of-Fluid (Taha and Cui [START_REF] Taha | Hydrodynamics of slug flows inside capillaries[END_REF], Akbar et al. [START_REF] Akbar | Simulation of Taylor flow in capillaries based on the volume-of-fluid technique[END_REF], Qian and Lawal [START_REF] Qian | Numerical study on gas and liquid slgus for Taylor flow in a t-junction microchannel[END_REF], Kumar et al. [START_REF] Kumar | Slug flow in curved microreactors: Hydrodynamic study[END_REF], Kashid et al. [START_REF] Kashid | On the hydrodynamics of liquid-liquid slug flow capillary microreactors[END_REF], Carlson et al. [START_REF] Carlson | Prediction of two-phase flow in small tubes: a systematic comparison of state-of-art cmfd codes[END_REF] and Gupta et al. [START_REF] Gupta | On the cfd modelling of Taylor flow in microchannels[END_REF]), Level-Set (Fukagata et al. [START_REF] Fukagata | Numerical simulation of gas-liquid two-phase flow and convective heat transfer in a micro tube[END_REF]), Lattice-Boltzmann (Yu et al. [START_REF] Yu | Experiment and lattice Boltzmann simulation of two-phase gas-liquid flows in microchannels[END_REF]) or Phase-field (He and Kasagi [START_REF] He | Phase-field simulation of small capillary-number two-phase flow in a microtube[END_REF]). The literature review presented here focuses on the previous works using the same method as in this study, i.e., the VOF method. Readers interested in others methods should refer to the publications mentioned above or to the review presented in the work of Gupta et al. [START_REF] Gupta | On the cfd modelling of Taylor flow in microchannels[END_REF]. Taha et al. [START_REF] Taha | Hydrodynamics of slug flows inside capillaries[END_REF] simulated Taylor flows in vertical capillaries using the VOF method implemented in the FLUENT TM software. For quite large capillary numbers (higher than 0.01), and, therefore, relatively thick films, the results on bubble velocity, bubble diameter and velocity field were in good agreement with experimental data. Akbar et al. [START_REF] Akbar | Simulation of Taylor flow in capillaries based on the volume-of-fluid technique[END_REF] demonstrated the feasibility of Taylor flow modeling with the VOF method by comparing numerical results with experimental data and correlations. Good agreement was found concerning bubble lengths and pressure drop for any capillary number. Regarding the liquid film thickness, more variable results were obtained depending on the simulation conditions. Qian and Lawal [START_REF] Qian | Numerical study on gas and liquid slgus for Taylor flow in a t-junction microchannel[END_REF] investigated Taylor flow in a T-junction with VOF method and reported that their grid was not enough refined to capture the finest liquid films on the wall. However, numerical results on bubble length and velocity showed good agreement with literature data. Kumar et al. [START_REF] Kumar | Slug flow in curved microreactors: Hydrodynamic study[END_REF] simulated for the first time Taylor flows in curved microchannels studying the effect on the bubble lengths of different parameters such as channel diameter, curvatures, inlet conditions or fluid properties. However, none of their simulations captured the liquid film on the wall, which plays an important role in this type of flow. Kashid et al. [START_REF] Kashid | On the hydrodynamics of liquid-liquid slug flow capillary microreactors[END_REF] also performed numerical simulations (Level-set and VOF) on Taylor bubble generation in a Y-junction coupled with an experimental study using the micro-PIV technique. Although they simulated a bubble train, the thin liquid film observed during their experiments could not be reproduced by any of the numerical method tested. In their work, Carlson et al. [START_REF] Carlson | Prediction of two-phase flow in small tubes: a systematic comparison of state-of-art cmfd codes[END_REF] performed simulations on the Taylor bubbles generation in a horizontal capillary tube using two diffuse-interface methods, Level-Set and Volume-Of-Fluid. They found that the VOF method implemented in FLUENT TM could not predict correctly the bubble generation in a slug flow.

More recently, Gupta et al. [START_REF] Gupta | On the cfd modelling of Taylor flow in microchannels[END_REF] performed a numerical study of Taylor flow in microchannels using the VOF method proposed in FLUENT TM . Numerical results of this study are in good agreement with the available correlations. However, they highlight the fact that an ultra fine mesh is required to capture the thin liquid film on the wall and that remain, in their simulations, some minor unphysical results such as pressure artefacts inducing spurious currents. They reduced these numerical errors by changing the way to calculate gradients in FLUENT TM (from the cell-based to the node-based methods).

These spurious currents induced by the VOF formulation were studied by Harvie and Rudman [START_REF] Harvie | An analysis of parasitic current generation in volume of fluid simulations[END_REF] on a two-dimensional bubble case. They emphasize the fact that spurious currents are correlated with the implementation of the surface tension effect in the VOF method, as explained in the next section concerning the method presentation. They showed that, depending on the situation, various parameters can influence the magnitude of these spurious currents such as the mesh refinement or the interface velocity. They especially highlights the fact that the spurious currents may have a big influence and can lead to unphysical behaviors of the simulations when the interface position is quasi-steady. This literature review emphasizes the fact that many previous studies have not captured successfully the thin liquid film when simulating Taylor flows with diffuse-interface numerical methods, except using ultra fine mesh with large computation time. This lack of liquid film reproduction can induce significant errors in the prediction of flow properties. Moreover, ultra fine mesh requires too large computation resources to be used at large scales.

Numerical method

According to the VOF model, dynamics of two immiscible fluids is governed by a single set of Navier-Stokes equations.

Fluid properties (viscosity µ, density ρ) in each cell are calculated using the volume fraction, φ, which represents the proportion of each fluid in a cell, using the following relationships:

ρ = φρ liquid + (1 -φ)ρ gas (1) µ = φµ liquid + (1 -φ)µ gas (2)
Like other Eulerian methods, the mesh is fixed and interfaces can be located by interpolation of the volume fraction. The Navier-Stokes system of equations contains an additional term representing the effect of surface tension (Brackbill et al. [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF]):

F = σ (∇ • n)∇φ (3)
where σ is the surface tension between the two fluids and n is the unit normal to the interface. This representation of the surface tension has been identified as responsible of the spurious currents along interfaces by Harvie and Rudman [START_REF] Harvie | An analysis of parasitic current generation in volume of fluid simulations[END_REF] due to the discrete representation of the domain. To illustrate their argument, they used the case of a static bubble. In a "ideal" VOF simulation, i.e., a continuous domain, the unit normal n to the interface is perfectly oriented to the center of the bubble. By solving the Navier-Stokes system, the VOF method simulates a zero velocity field with an opposite pressure gradient equivalent to the surface force. In a "real" VOF simulation, i.e., a domain discretized in cells, the volume fraction varies by step which implies that the unit normal n varies slightly and is not perfectly oriented towards the bubble center.

The calculated force vector is not exact and the amplitude of the error depends on the variation of φ within the cell. The error induced results in local recirculations of the velocity field and pressure artefacts which are relatively important in the case of static interfaces.

The volume fraction φ is transported using the VOF equation:

∂φ ∂t + V • ∇φ = 0 ( 4
)
where V is the fluid velocity. Usually, for unsteady-state problems, the VOF equation is solved explicitly and is decoupled from pressure and velocity computations. Variable time steps are controlled through the Courant number Co defined as follows:

Co = V fluid t x ( 5
)
The Courant number condition may be different for the Navier-Stokes or VOF equations. Standard FLUENT TM parameters set the Courant number condition to 0.25 for the VOF equation. In the following parametric study, we discuss the Courant number condition for the Navier-Stokes system. Since this system is solved implicitly, it represents most of the computing time and, therefore, the best gain prospect.

Another important parameter in the VOF method is the scheme used for the discretization of the convective term in Eq. ( 4). As the fluids are immiscible, the interface is represented by a step in φ. If we use classical "Upwind" or "2nd

Order Upwind" schemes, the numerical diffusion will lead to a strong diffusion of the interface. Several strategies have been designed to overcome this numerical diffusion using specific schemes, most of them being implemented in FLUENT TM (hence our choice of this software for this study). The ability of a discretization scheme to transport and conserve a sharp interface is called the scheme compressibility. We described briefly below the four discretization schemes proposed by FLUENT TM , ranked from the least to the most compressive: the QUICK scheme [START_REF] Leonard | ULTRA-SHARP nonoscillatory convection schemes of high-speed steady multidimensional flow[END_REF], High-Resolution schemes (HRIC [START_REF] Muzaferija | A two-fluid Navier-Stokes solver to simulate water entry[END_REF] and CICSAM [START_REF] Ubbink | A method for capturing sharp fluid interfaces on arbitrary meshes[END_REF]) and the Geometric Reconstruction scheme [START_REF] Youngs | Time-dependent multi-material flow with large distortion[END_REF].

The QUICK scheme is based on a weighted average of 2nd order upwind and central interpolation. The central interpolation tends to reduce the numerical diffusion but we shall see in this paper that it remains the less compressive scheme. Based on the work of Leonard [START_REF] Leonard | The ultimate conservative difference scheme applied to unsteady one-dimensional advection[END_REF], high-resolution schemes have in common the alternate use of a downwind-type scheme (signal preserving and interface compression but unstable) and an upwind-type scheme (stable but diffusive). Some authors (see Ubbink and Issa [START_REF] Ubbink | A method for capturing sharp fluid interfaces on arbitrary meshes[END_REF]) reported that using a downwind-type scheme can bring numerical errors such as the alignment of the interface on the mesh grid. Improvements have been brought to high-resolution schemes by reducing the use of the downwind-type scheme depending on the interface orientation. The main difference between the HRIC and the CICSAM schemes concerns the downwind-type scheme. In the CICSAM method, the compressibility of the downwind-type scheme increases when the local Courant number decreases (due to the downwind-type scheme called HYPER-C, see Ubbink and Issa [START_REF] Ubbink | A method for capturing sharp fluid interfaces on arbitrary meshes[END_REF]). HRIC scheme is independent of the local Courant number and its compressibility is equivalent to the CICSAM scheme with a local Courant number condition equal to 0.5. Recommended maximum Courant number for VOF iterations is 0.25. As a result, the CICSAM scheme in our simulations is always more compressive than the HRIC scheme. A previous study performed by Waclawczyk and Koronowicz [START_REF] Waclawczyk | Comparison of CICSAM and HRIC high-resolution schemes for interface capturing[END_REF] confirms that CICSAM preserves the shape of an interface more accurately than the HRIC scheme.

The last method called Geo-Reconstruct in FLUENT TM is based on the interface reconstruction concept and represents the most compressive solution. In a first step, the interface is transported using standard schemes which can be quite diffusive. In a second step, the numerical diffusion is canceled by interface compression using a piecewise linear interpolation (PLIC). This method is the most common method and is used in all the previous studies presented above, except in the work of Carlson et al. [START_REF] Carlson | Prediction of two-phase flow in small tubes: a systematic comparison of state-of-art cmfd codes[END_REF] who used the CICSAM scheme.

Parametric study on the drainage of a liquid phase in a vertical tube

Properties of slug flows were studied by Bretherton [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF] who investigated the motion of a long bubble into a liquid-filled tube of diameter D and showed that the liquid film thickness, e, on the wall in a tube depends on the capillary number Ca (defined as the ratio of the viscous forces to interfacial tension, Ca = µV /σ , where µ is the liquid viscosity, V the drainage average velocity and σ the two-fluid interfacial tension), according to the relationship: 2e/D = aCa 2/3 [START_REF] Kumar | Slug flow in curved microreactors: Hydrodynamic study[END_REF] valid for Ca → 0 with a = 1.338 for a horizontal tube. From a study of the drainage of a liquid phase in a vertical tube, Lasseux [START_REF] Lasseux | Caractérisation expérimentale, analytique et numérique d'un film dynamique lors du drainage d'un capillaire[END_REF] extended this relationship for non-zero Bond numbers using a semi-analytical solution of the problem and then for a wide range of capillary numbers with direct numerical simulations. These simulations were performed using a Lagrangian-method, the Boundary Element Method, which discretizes explicitly the interface and follows the nodes displacements. This type of method has the advantage of following accurately the interface but it is quite difficult to simulate changes in the number of interfaces, and, also, it produces filled matrices in the solution procedure, which leads to heavy computations.

A vertical capillary tube of diameter 2 mm and length 10D is used for the simulations. Geometry and boundary conditions are reported in Fig. 2 and physical conditions are reported in Table 1. The capillary number is modified by varying the liquid phase viscosity.

Before discussing the influence of the different simulation parameters, we first perform a comparison as regards the computation of the liquid film thickness during drainage (cf. Fig. 3). VOF results fit well to the reference results obtained by Lasseux [START_REF] Lasseux | Caractérisation expérimentale, analytique et numérique d'un film dynamique lors du drainage d'un capillaire[END_REF] with the Boundary Integral Element method, thus showing that the VOF method can simulate accurately a slug flow for a wide range of capillary numbers including thin liquid films, at least when the simulation is correctly parametrized as we shall see below. As expected, we can observe that the numerical results are in good agreement with the theory for sufficient low capillary numbers (Ca < 0.1) and diverge for higher Ca. As can be seen in Fig. 4, the velocity and pressure profiles do not present any pressure artefacts or spurious currents. How were these results obtained? Does the choice of the various numerical parameters influence significantly the results? A parametric study was done to evaluate the impact of the various choices.

The influence of the irregularity of the mesh on the liquid film thickness is first tested using four different meshes refined along the tube wall (cf. Fig. 5). Meshes have the same number of cells and, therefore, computation times for the different simulations are nearly the same. Using the HRIC scheme and a Courant number set to 0.25, our results, which are summarized in Table 2, show that irregular meshes have a small influence on the simulation results for different capillary numbers provided that the liquid film covers at least two cells. Having more than two cells in the liquid film does not seem Dimensionless liquid film thickness (normalized by the tube radius) for different meshes and capillary numbers (HRIC, Co = 0.25) compared with the reference solution (Lasseux [21]).

Ca

Lasseux [START_REF] Lasseux | Caractérisation expérimentale, analytique et numérique d'un film dynamique lors du drainage d'un capillaire[END_REF] Regular mesh Mesh Fig. 8. Alignment of the interface on the mesh when VOF scheme is too compressive, for example with the Geo-Reconstruct scheme (top) compared with the HRIC scheme (bottom).

to lead to a better accuracy. Our results show that it is possible to capture accurately thin films (less than 1% of the radius) while conserving small computation time by locally refining the mesh. In a second phase, a parametric study on the VOF schemes and the Courant number was performed. In the case of a poor model parametrization, spurious currents (in parallel with pressure artefact) can occur around interfaces, which was previously observed by Gupta et al. [START_REF] Gupta | On the cfd modelling of Taylor flow in microchannels[END_REF] and Harvie and Rudman [START_REF] Harvie | An analysis of parasitic current generation in volume of fluid simulations[END_REF]. For all the tested VOF schemes, we observed that, when Co > 0.25, these spurious currents occur around the meniscus curvature and increase with the Courant number Co (an example is shown in Fig. 6). These pressure artefacts and spurious currents also appear at the triple point (intersection of the interface with the wall) and move along the liquid film (cf. Fig. 7) when the liquid film becomes thin (2e < 0.10D) for the two most compressive methods, Geo-Reconstruct and CICSAM. As in the study of Harvie and Rudman [START_REF] Harvie | An analysis of parasitic current generation in volume of fluid simulations[END_REF], we observe, in this area, a curvature of the interface on a small number of cells. We can presume that the compressibility of these two schemes, which reduces the thickness of the interface, increases the intensity of the interfacial force on a small number of cells (i.e., the ∇φ term in the Navier-Stokes equations), and, therefore, the local spurious effects. Moreover, the link of these spurious effects to interface compression can be directly emphasized using the CICSAM scheme properties. By reducing the Courant number for the same simulations, we increase the compressibility of the CICSAM scheme (see Ubbink and Issa [START_REF] Ubbink | A method for capturing sharp fluid interfaces on arbitrary meshes[END_REF]) and we can observe that it simultaneously increases the amplitude of the spurious currents generated. The interface compression can also lead to the alignment of the interface on the mesh grid that was previously described by Ubbink and Issa [START_REF] Ubbink | A method for capturing sharp fluid interfaces on arbitrary meshes[END_REF]. Then, the interface thickens by step along the tube wall as it can be seen in Fig. 8.

We also performed simulations using different discretization schemes which are available in Fluent TM ("1st Order Upwind", "2nd Order Upwind" and "3rd Order MUSCL"). As the influence of these schemes was very weak compared to the previously cited parameters, we did not incorporate these parameters to this parametric study. The interface thickness (in numbers of cells) at the front and along the film for different schemes and Courant numbers are displayed in Table 3. The CICSAM and Geo-Reconstruct methods conserve sharp interfaces but, as seen before, there is no valid Courant number to obtain consistent velocity and pressure fields. The QUICK and HRIC schemes are both able to simulate correctly the properties (velocity profile, liquid film thickness) of the slug flow. However, as shown in Fig. 9, the QUICK scheme is more diffusive that the HRIC scheme and it can lead to interface disappearance when the flow becomes more complex.

As a conclusion to this parametric study, the HRIC scheme, with Co = 0.25, seems to be the optimal combination of parameters to simulate accurately slug flows with thin liquid films while avoiding the numerical errors described above. This optimal combination of parameters have been used to perform the comparison related in Fig. 3.

So far, we have studied a simple flow geometry, with the emphasis on the film accurate simulation. The existence of a well established reference helped us to draw practical conclusions in terms of suitable numerical parameters. What happens in the case of more complex situations, for which there is a little a priori knowledge of the obtained configuration?

Validation case: Taylor bubble generation

Taylor bubble generation in a horizontal capillary tube is a good test case since films and meniscii are present in an unsteady-state flow. Like in the study of Gupta et al. [START_REF] Gupta | On the cfd modelling of Taylor flow in microchannels[END_REF], we inject air and water (superficial velocities are respectively 0.245 and 0.255 m/s) in an initially liquid filled microchannel with diameter, D, of 0.5 mm and length of 10D as depicted on Fig. 10. The conditions of the simulations, summarized in Table 4, are the same as Gupta et al. [START_REF] Gupta | On the cfd modelling of Taylor flow in microchannels[END_REF] except for the VOF scheme and the mesh grid used. As the expected dimensionless film thickness is around 0.02D, we use a sufficiently refined irregular mesh with a width ratio of 1.10, i.e., between 2 and 3 cells included in the liquid film. Plots of the liquid volume fraction during the simulation at different times can be seen in Fig. 11 and results of simulations are reported in Table 5.

Concerning the film thickness, bubble velocity and gas hold-up, our simulations fit pretty well with reference values. We can note that bubble and slug lengths are quite different between the two simulations even if the gas hold-up is nearly the [START_REF] Kashid | On the hydrodynamics of liquid-liquid slug flow capillary microreactors[END_REF] where Re is the Reynolds number and Ls = L slug /D is the dimensionless slug length. Pressure drop can be calculated by the following equation:

P L = f 1 2 ρ L U 2 4 D β L ( 8
)
where β L is the liquid hold up and U the sum of gas and liquid superficial velocities. Difference between correlated and computed pressure drops is around 20% for the two types of simulations. As can be seen from Fig. 12, the pressure field simulated with the HRIC scheme does not present any artefact in the tail or the head of the Taylor bubble. We can note that the creation of bubbles induces a small diffusion at the tail of the bubble (cf. Fig. 11 and Fig. 12) contrary to the Geo-Reconstruct method which always conserves a sharp interface. This interface diffusion gradually disappears thanks to the scheme compressibility. In terms of computation time (cf. Table 4), capturing the liquid film with the Geo-Reconstruct method requires ultra fine mesh and several days of calculations [START_REF] Gupta | On the cfd modelling of Taylor flow in microchannels[END_REF]. Using a suitable irregular mesh with the HRIC scheme reduces the computation time to a couple of hours for the same physical conditions. As a test, we also performed a qualitative comparison by changing the simulation conditions. Carlson et al. [START_REF] Carlson | Prediction of two-phase flow in small tubes: a systematic comparison of state-of-art cmfd codes[END_REF], reported that the VOF method (with the CICSAM scheme) failed to simulate the generation of a regular bubble train. Similar simulations were carried out in the same conditions using the HRIC and CICSAM schemes. Similar results as in the work of Carlson et al. [START_REF] Carlson | Prediction of two-phase flow in small tubes: a systematic comparison of state-of-art cmfd codes[END_REF] were observed concerning the CICSAM scheme, i.e., the random generation of bubble of different sizes, while the use of the HRIC scheme allows simulating correctly the flow phenomenology (see Fig. 13). However, we can also see in Fig. 13 that the interface dissipation becomes more important, which suggests a limit to the HRIC scheme.

Conclusions

We emphasized in this Note that simulating slug flows with the VOF method needs suitable parameterization to provide computation speed and accuracy. We mainly focus on capturing thin liquid films on the walls which is the most difficult aspect of the problem. We first validated quantitatively VOF simulations in the case of the classical Bretherton dynamic film and conclusions provided by the parametric study highlight some aspects which must be considered to simulate slug flows with the VOF method:

-irregular meshes have a small influence on computed results provided that the liquid film is covered by at least two cells, which allows refining locally the mesh and gain much computation time; -Courant number must be chosen sufficiently small to avoid pressure artefacts and spurious currents around interfaces; -using the most compressive scheme induces numerical errors and does not represent the best choice for simulating thin liquid films; less compressive schemes allow using coarser meshes while maintaining good accuracy and avoiding the different numerical errors described above.

In the specific case of slug flows in capillaries, the HRIC scheme with a Courant number of 0.25 and an irregular mesh depending on the expected liquid film is the best suited to simulate accurately the properties of slug flows. However, the process of bubble creation induces interface dissipation which is obviously more important for the least compressive schemes. This observation suggests that the most compressive schemes will be probably more suited to simulate flow regimes with many interface changes such as the bubbly flow.

Fig. 1 .

 1 Fig. 1. Regime flow map obtained from experiments (Triplett et al. [1]) and summarized in this figure by Gupta et al. [2].

Fig. 2 .

 2 Fig. 2. Geometry and boundary condition for the drainage simulations.

Fig. 3 .

 3 Fig. 3. Dimensionless liquid film thickness as a function of capillary number during the drainage of a liquid phase in a vertical tube. Comparison between VOF simulation and a reference solution (Lasseux [21]).

Fig. 4 .

 4 Fig. 4. Simulation of the drainage of a liquid phase (HRIC, Co = 0.25, Ca = 0.005), (a) volume fraction, (b) velocity profile, (c) pressure distribution.

Fig. 5 .

 5 Fig. 5. Mesh used for VOF simulations refined from the tube wall (left) to the center (right) with different width ratio between cells (a. 1.00 (regular), b. 1.05, c. 1.10, d. 1.20).

Fig. 6 .

 6 Fig. 6. Simulation of drainage with spurious currents along the meniscus curvature (HRIC, Co = 1.00, Ca = 0.03), (a) volume fraction, (b) velocity profile, (c) pressure distribution.

Fig. 7 .

 7 Fig. 7. Simulation of drainage with spurious currents along the liquid film (CICSAM, Co = 0.05, Ca = 0.03), (a) volume fraction, (b) velocity profile, (c) pressure distribution.

Fig. 9 .

 9 Fig. 9. A comparison of the interface dissipation from the less to the most compressive scheme (Co = 0.25, Ca = 0.03), (a) QUICK, (b) HRIC, (c) CICSAM, (d) Geo-Reconstruct.

Fig. 10 .

 10 Fig. 10. Geometry and boundary conditions for the Taylor bubbles simulations.

Fig. 11 .

 11 Fig. 11. VOF simulations of the process of bubble generation in a microchannel using the HRIC scheme at different times (3.6 ms, 7.0 ms and 8.5 ms).

Fig. 12 .

 12 Fig. 12. Volume fraction and Pressure distribution around a typical Taylor bubble simulated with the VOF method (HRIC, Co = 0.25).

Fig. 13 .

 13 Fig. 13. A qualitative comparison of VOF simulations of bubble generation using HRIC scheme (top) and CICSAM scheme (bottom).

Table 1

 1 Simulation conditions for the drainage of a liquid phase in a vertical tube.

	Drainage velocity	0.007 m s -1
	Surface tension	0.002
	Liquid phase density	210 kg m -3
	Bond number	1.03

Table 2

 2

Table 3

 3 Interface thickness at the front/film (in number of cells).

	Courant number	0.05	0.25	0.50	1.00
	QUICK	6.86/3.40	5.76/3.33	9.68/3.40	14.02/3.73
	HRIC	4.32/3.33	3.34/2.80	3.96/2.73	4.58/2.80
	CICSAM	2.56/1.93	2.61/2.73	2.34/2.47	2.54/2.60
	Geo-Reconstruct	1.98/1.92	2.05/1.93	1.93/1.92	1.86/1.94

Table 4

 4 Parameters of the two numerical simulations studied.

	Parameters	VOF Geo-Reconstruct [2]	VOF HRIC
	Pressure-velocity scheme	QUICK	QUICK
	VOF scheme	Geo-Reconstruct	HRIC
	Courant number	0.25	0.25
	Number of mesh cells	2 × 10 3 to 4.5 × 10 5	8 × 10 3 (irregular mesh)
	Computing time	1 hour to 13 days	2-3 hours

Table 5

 5 Comparisons of numerical results (with reference values).

	Flow properties	VOF Geo-Reconstruct [2]	VOF HRIC
	Film thickness (µm)	12 (11)	12 (11)
	Bubble velocity (m s -1 )	0.55 (0.54)	0.56 (0.54)
	Bubble length	2.0D	1.62D
	Slug length	1.3D	1.07D
	Gas hold-up	0.455 (0.41)	0.414 (0.41)
	Pressure drop (Pa m -1 )	77 700 (94 400)	137 000 (112 000)


R Lozi
email: lozi@unice.fr

Chaotic mathematical circuitry

Keywords:

come

Introduction

Since the seminal work of E. N. Lorenz [START_REF] Lorenz | Deterministic nonperiodic ow[END_REF] in 1963, who discovered by accident the rst chaotic strange attractor, chaotic dynamical systems have been fully studied. Fifty years after, only few chaotic attractors involving di¤erential equations remain actively explored. Among them Chuas attractor is nowadays incredibly used, because both of its realizations: electronic circuit or system of di¤erential equations can be combined for multiple purpose. [START_REF] Duan | Multiinput and multi-output nonlinear systems interconnected Chuas ciruits[END_REF] Following the rst studies applying such combinations to crypted transmission, our aim is to build an analog of paradigm of electronic circuitry, which is the design of electronic circuit: the paradigm of chaotic mathematical circuitry, in order to improve the performance of well known chaotic attractors for application purpose (cryptography, generic algorithms in optimization, control,...). An electronic circuit is composed of individual electronic components, such as resistors, transistors, capacitors, inductors and diodes, connected by conductive wires through which electric current can ow. The combination of components and wires allows various simple and complex operations to be performed. In the same way a mathematical circuit is composed of individual components we introduce (generators, couplers, samplers, mixers, reducers and cascaders,...) connected through streams of data. The combination of such mathematical components leads to several news applications such as improving the performance of well known chaotic attractors (Belykh, [START_REF] Loskutov | Stabilization of chaotic oscillations in systems with a hyperbolic-type attractor[END_REF] Lorenz, Rössler, [START_REF] Rössler | Chaotic behavior in simple reaction system[END_REF] ...) for application purpose (chaotic cryptography, genetic algorithms in optimization, control,...).

In Sec. 2 we recall a rst historic example of chaotic mathematical circuitry: the cascading of two identical receivers in Chuas circuit. In Sec. 3, we introduce the new paradigm of chaotic mathematical circuitry, some examples of applications of which are given in Sec. 4.

A rst historic example of chaotic

mathematical circuitry: the cascading of two identical receivers in Chuas circuit

Chuas circuit

In october 1983, visiting T. Matsumoto at Waseda University, L. O. Chua invented an electronic circuit (Fig. 1(a) and (b)) mimicking directly on an oscilloscope screen a chaotic signal (Fig. 1(c)).

Only two autonomous systems of ordinary di¤er-ential equations were generally accepted then as being chaotic, the Lorenz equations [START_REF] Lorenz | Deterministic nonperiodic ow[END_REF] and the Rössler equations. [START_REF] Rössler | Chaotic behavior in simple reaction system[END_REF] The nonlinearity in both systems is a the product function of two variables which is very di¢cult to build in electronic circuit. L. O. Chua [START_REF] Chua | The genesis of Chuas Circuit[END_REF] says, The fault lies on the dearth of a critical nonlinear IC component with a near-ideal characteristic and a su¢ciently large dynamic range; namely, the analog multiplier. Unfortunately, this component was the key to building an autonomous chaotic circuit in 1983. He adds, Suddenly [in the evening of the precise day he attended an unsuccesful presentation of an electronic circuit realization of the Lorenz Equations] it dawned upon me that since the main mechanism which gives rise to chaos, in both the Lorenz and the Rössler Equations, is the presence of at least two unstable equilibrium points -3 for the Lorenz Equations and 2 for the Rössler Equationsit seems only prudent to design a simpler and more robust circuit [than that built by Matsumotos team] having these attributes. Having identied this alternative approach and strategy, it becomes a simple exercise in elementary nonlinear circuit theory to enumerate systematically all such circuit candidates, of which there were only 8 of them, and then to systematically eliminate those that, for one reason or another, can not be chaotic. The Chuas equations he laid down

<

:

: x = B (y 0 x 0 f (x)) ; : y = x 0 y + z; : z = 0Cy; (1)
where

f (x) = bx + 1 2 (a 0 b) [jx + 1j 0 jx 0 1j] ; (2
were soon numerically analyzed by T. Matsumoto. [START_REF] Matsumoto | A chaotic attractor from Chuas circuit[END_REF] The nonlinear characteristic which is in fact only piecewise linear allows some exact computations. Henceforth, L. O. Chua et al. [START_REF] Chua | The Double Scroll Family[END_REF] proved in the same breath that the mechanism of chaos exists in this attractor.

Secure communication via chaotic synchronization

Synchronization

The synchronization of two Chuas circuit was studied experimentaly height years later in 1992, [START_REF] Chua | Experimental chaos synchronization in Chuas circuit[END_REF] soon followed by its application to crypted transmission. The rst laboratory demonstration of a secure communication system which uses a chaotic signal for masking purposes [START_REF] Oppenheim | Signal processing in the context of chaotic signals[END_REF] and which exploits the chaotic synchronization techniques to recover the signal was reported in 1992. [START_REF] Kocarev | Experimental demonstration of secure communication via chaotic synchronization[END_REF] While the transmitter in this system is a direct implementation of the method proposed in oppenheim et al., [START_REF] Oppenheim | Signal processing in the context of chaotic signals[END_REF] the receiver di¤ers from their computer simulation approach in that it actually contains two subsystems of the chaotic transmitter (Chuas circuit in that case).

Single chaotic synchronization

The mathematical translation of the dynamics of the circuit used by Kocarev et al., [START_REF] Kocarev | Experimental demonstration of secure communication via chaotic synchronization[END_REF] for the experimental demonstration of secure communication is as follows: the basic building block is a Chuas circuit, the dy-namics of which is given by the Chuas equations ( 1) and ( 2). Here x(t) is used as noise-like masking signal. Let s(t) be an information-bearing signal. The transmitted signal is r(t) = x(t) + s(t), where the power level of s(t) is assumed to be signicantly lower than that of x(t), in order to have the signal e¤ectively hidden. The receiver consists of two subsytems.

The rst one is driven by the transmitted signal r(t):

:

y 1 = r (t) 0 y 1 + z 1 ;
:

z 1 = 0Cy 1 : (4)
The second subsytem is driven by the signal y 1 (t):

:

x 2 = B (y 1 (t) 0 x 2 0 f (x 2 )) : (5)
Then s(t) is recovered as

s 2 (t) = r (t) 0 x 2 (t) s (t) : (6)
Actually the dynamics of the experimental set-up (see Fig. 2) is described by C 2 . The triangular symbols are OpAmps which decoupled the systems, acting as the signal drive elements, from. 1 5 However, as long as we do not need z 2 (t) to recover s 2 (t), we continue to use Eq. 6 instead of Eq. 7 in the following improved system.

: x 2 = B (y 1 (t) 0 x 2 0 f (x 2 )) ; : z 2 = 0Cy 1 (t) : (7)

Cascade chaotic synchronization

Ten years after the invention of his ubiquitous realworld example of a chaotic system, leading some to declare it a paradigm for chaos, and few months Fig. 4. Electronic circuit implementation of the two stage receiver consisting of two identical copies of the circuit given in Fig. 3.

after the experimental demonstration of secure communication using the properties of this chaotic generator, Professor L. O. Chua was visiting us for one month, in May 1993, at the University of Nice.

In both implementations -electronic circuit realization (Fig. 3) or computer simulation (Fig. 2)-of the circuit used by Kocarev et al., [START_REF] Kocarev | Experimental demonstration of secure communication via chaotic synchronization[END_REF] there is an inevitable error introduced by the signal s(t). Then the claim was to enhance the performance of the chaotic masking technique by improving the convergence of the recovered signal s 2 (t) towards the informationbearing signal s(t). Having in mind the knowledge of relaxation methods used in numerical analysis (in numerical mathematics, relaxation methods are iterative methods for solving systems of equations, including nonlinear systems), we proposed to iterate the process of recovering the signal, cascading a second identical receiver to the rst one, i.e. introducing a second system of equations comparable to Eqs. 4-5 driven by x 2 (t) instead of r(t); as displayed in Fig. 5. The second receiver also consists of two subsytems. The rst one is driven by the signal x 2 (t) which is assumed to be more synchronized to x(t) than the transmitted signal r(t):

: y 3 = x 2 (t) 0 y 3 + z 3 : z 3 = 0Cy 3 (8)
The second subsystem of the second receiver is then driven by the signal y 3 (t) from Eq. 8:

: x 4 = B (y 3 (t) 0 x 4 0 f (x 4 )) (9)
Then s(t) is recovered as

s 4 (t) = r (t) 0 x 4 (t) s (t) (10)
In practice we simply make two copies of the receiver as shown in Fig. 4. By identifying the symbols (V C1 ; V C2 ; I l ) in Chuas circuit (see Figs. 2 and 3) with (x; y; z), the electronic circuit implementation in Fig. 4 can be translated into the block diagram shown in Fig. 5. Although no two electronic circuits can be made perfectly identical in practice, this ideal situation can be approached with the help of the integrated circuit technology demonstrated in Delgado-Restituto & Rodriguez-Vasquez. [START_REF] Delgado-Restituto | A CMOS monolithic Chuas circuit[END_REF] By fabricating several identical Chuas circuits on the same silicon chip, the resulting circuits are almost clones of each other. This technique has the additional security adavantage in that even if someone else has discovered the parameters (B; C) used in the system, integrating it into another silicon chip invariably introduces discrepancies due to the di¤erent processing parameters from di¤erent silicon foundries. We have shown by computer experiments that by connecting two identical receivers, a signicant amount of noise can be reduced, thereby allowing the recovery of a much higher quality signal. 18

Numerical experiments

Due to the limited extend of this chapter, we give only a sketch of results of numerical experiments. Assuming that the input (information-bearing) signal s(t) is a single tone (sine wave) of amplitude k 1:

s(t) = k sin(!t) with k > 0: (11)
Therefore s(t) has a power level signicantly lower than that of x(t).

Let us dene the errors

ke s (k; !)k 2 = ks 2 (t) 0 s(t)k 2 (12) ke c (k; !)k 2 = ks 4 (t) 0 s(t)k 2 (13)
where

kf (t)k 2 , lim T 1 T " Z T 0 f 2 (t) dt # 1 2 (14)
denotes the quadratic norm of f (t).

Extensive computer studies show that both ke s (k; !)k 2 = ks(t)k 2 and ke c (k; !)k 2 = ks(t)k 2 are independent of k, and that both norms ke s (k; !)k 2 and ke c (k; !)k 2 are decreasing in accordance to a power law with an exponent greater than 2 when ! is increasing. Moreover ke c (k; !)k 2 is always less than ke s (k; !)k 2 . This shows that the cascade chaotic synchronization technique o¤ers a good improvment over the single-stage chaotic synchronization results.

Mathematical circuitry

Our aim is to build an analog of paradigm of electronic circuitry, which is the design of electronic circuit: the paradigm of mathematical circuitry and especially chaotic mathematical circuitry in order to improve the performance of well known chaotic attractors (Belykh, [START_REF] Loskutov | Stabilization of chaotic oscillations in systems with a hyperbolic-type attractor[END_REF] Lorenz, [START_REF] Lorenz | Deterministic nonperiodic ow[END_REF] Rössler, 28 ...) for application purpose (chaotic cryptography, genetic algorithms in optimization, control, emergence of ramdoness from chaos,...).

An electronic circuit is composed of individual electronic components, such as capacitors, diodes, inductors, resistors, transistors and connected by conductive wires through which electric current can ow. Recently, in 2009, three more components discovered by L. O. Chua et al., [START_REF] Ventra | Circuit elements with memory: memristors, memca-pacitors, and meminductors[END_REF] have been added: memristors, memcapacitors and meminductors (Fig. 6). The combination of components and wires allows various simple and complex operations to be performed: signals can be amplied, computations can be accomplished, and data can be moved from one place to another. In the same way a mathematical circuit is composed of individual components we introduce now: generators, couplers, samplers, mixers, reducers and cascaders, connected by streams of data. The combination of such mathematical components leads to several performing news applications, such as Chaotic Pseudo Random Number Generators (CPRNG), as we will show in the rest of this article.

Generator

Analog circuits are very commonly represented in schematic diagrams, in which wires are shown as lines, and each component has a unique symbol. We present in this section the symbols we design in order to draw mathematical schematic diagrams. The rst symbols we describe, generator symbols, are, from a mathematical point of view, equivalent to a battery or a current generator in electronic circuit. However we consider that they generate a numerical signal (in one or more dimensions) rather than a voltage or an intensity variation (nonetheless, a voltage or intensity variation can be considered as a physical signal which can be discretized). This signal can be either continuous, as in Chuas circuit, or discrete as in Hénon map. [START_REF] Hénon | A Two-dimensional mapping with a strange attractor[END_REF] In the expanded symbol of continuous generator (see Fig. 7(a)), the solid line arrows coming out from the generator represent the three components of the signal x (t) = (x (t) ; y (t) ; z (t)) (see Eq. 1), the dashed line arrow points at which stands for the parameter value dened by Eq. 3, and the dot line arrow points at x 0 = x (0), the given initial value of the signal. If there is no ambiguity on the nature of the generator used, the symbol can be simplied as in Fig. 7(b).

We need also to design chaotic circuitry for discrete signal. To this aim two generators can be considered: in dimension 2, the Hénon map (see Figs.

which is associated to the dynamical system,

x n+1 = y n + 1 0 ax 2 n y n+1 = bx n ( 17
)
and the logistic map

f r (x) = rx (1 0 x) (18)
linked to the one dimensional system which is chaotic when r = 4 (see Fig. 7(e)).

x n+1 = rx n (1 0 x n ) (19)
Thereafter, another 1-dimensional chaotic generator, the symmetric tent map, will be, also, represented by the same symbol.

Coupler

Two years after the experimental study of the synchronization of two Chuas circuit, in 1994, experimental observation of hyperchaotic attractors in open and closed chain of Chuas circuit was reported. [START_REF] Kapitaniak | Experimental hyperchaos in coupled Chuas circuits[END_REF] The layout of the ve identical coupled Chuas circuit forming a ring is displayed on Fig. 8, the state equations of this circuit are as follows

8 > > > > > > > > > > > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > > > > > > > > > > > : C 1 dv (1) C 1 dt = G v (1) C2 0 v (1) C1 0 f v (1) C1 ; C 2 dv (1) C 2 dt = G v (1) C1 0 v (1) C2 + i (1) L + K 1 v (2) C2 0 v (1) C2 ; L di (1) L dt = 0v (1) C2 ; C 1 dv (2) C 1 dt = G v (2) C2 0 v (2) C1 0 f v (2) C1 ; C 2 dv (2) C 2 dt = G v (2) C1 0 v (2) C2 + i (2) L + K 2 v (3) C2 0 v (2) C2 ; L di (2) L dt = 0v (2)
C2 ;

1 1 1 1 1 1 1 1 1 C 1 dv (5) C 1 dt = G v (5) C2 0 v (5) C1 0 f v (5) C1 ; C 2 dv (5) C 2 dt = G v (5) C1 0 v (5) C2 + i (5)
L + K 5 v (1) C2 0 v (5) C2 ; L di (2) L dt = 0v (5)
C2 :

Each Chuas circuit (see Fig. 1 (a)) contains three linear energy-storage elements (an inductor and two capacitors), a linear resitor, and a single nonlinear resistor, namely, Chuas diode with three segment linear characteristic dened by

f (v R ) = m 0 v R + 1 2 (m 1 0 m 0 ) [jv R + B p j 0 jv R 0 B p j] (21)
where the slopes in the inner and the outer regions are m 0 and m 1 , respectively, and 6B p denotes the breakpoints (see Fig. 1 (b)). Equation 21is equivalent to Eq. 2. By identifying the symbols

V (i) C1 ; V (i) C2 ; I (i) l in
each Chuas circuit with 0 x i ; y i ; z i 1 , the state equations of the circuit can be translated into the di¤erential equations (22) and the electronic circuit symbolized by the mathematical circuit of Fig. 9.

8 > > > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > > > : : x 1 = B 0 y 1 0 x 1 0 f 0 x 1 11 ;
:

y 1 = x 1 0 y 1 + z 1 + k 1 0 y 2 0 y 1 1 ; : z 1 = 0Cy 1 ; : x 2 = B 0 y 2 0 x 2 0 f 0 x 2 11 ;
:

y 2 = x 2 0 y 2 + z 2 + k 2 0 y 3 0 y 2 1 ; : z 2 = 0Cy 2 ; 1 1 1 1 1 1 1 1 1 : x 5 = B 0 y 5 0 x 5 0 f 0 x 5 11 ;
: y 5 = x 5 0 y 5 + z 5 + k 5 0 y 1 0 y 5 1 ;

:

z 5 = 0Cy 5 : (22)
In this gure the double arrows symbolize the coupling k i 0 y i+1 0 y i 1 of one Chuas circuit to the next one. In order to represent the coupling between mathematical equation, depending on the nature of the coupling, we can use both symbols: the ring coupler corresponding to the coupling of one generator to the next one (Fig. 10 It has been shown, a few years ago [START_REF] Lozi | Giga-periodic orbits for weakly coupled tent and logistic discretized maps, Modern Mathematical Models, Methods and Algorithms for Real World Systems[END_REF] that the ultra-weak coupling of several logistic maps (Eq. 18) or symmetric tent maps

f a (x) = 1 0 a jxj ( 23
)
x n+1 = f a (x n ) ( 24
)
when a = 2; allows the production of long series of chaotic numbers equally distributed over the interval J = [01; 1] 4. When a dynamical system is realized on a computer using oating point numbers, the computation is of discretization, where nite arithmetic replaces continuum state space. For chaotic dynamical systems, the discretizations often have collapsing e¤ects to a xed point or to short cycles. [START_REF] Lozi | Can we trust in numerical computation of chaotic solutions of Dynamical systems ?[END_REF] Instead, the ultra-weak coupling of logistic or symmetric tent map 8 < :

x 1 n+1 = (1 0 2" 1 ) f 0 x 1 n 1 + " 1 f 0 x 2 n 1 + " 1 f 0 x 3 n 1 x 2 n+1 = " 2 f 0 x 1 n 1 + (1 0 2" 2 ) f 0 x 2 n 1 + " 2 f 0 x 3 n 1 x 3 n+1 = " 3 f 0 x 1 n 1 + " 3 f 0 x 2 n 1 + (1 0 2" 3 ) f 0 x 3 n 1 (25
) symbolized by Fig. 11(a) gives rise to sterling model of generator of chaotic numbers with a uniform distribution of these numbers of the interval [01; 1]. Ultra-weak coupling means " i 2 2 10 015 ; 10 07 3 for computations using double precision numbers. [START_REF] Lozi | Complexity leads to randomness in chaotic systems Mathematics in Science and Technology: Mathematical Methods, Models and Algorithms in Science and Technology[END_REF] The periodic solutions of such generator have period far greater than ten billions. Of course, more than three 1-dimensional generators can be coupled under the same process, enhancing the ergodic properties of the generator.

The considered system of the p-coupled dynamical systems is with

X n+1 = F (X n ) = A: 0 f (X n ) 1 (26)
X n = 0 B @ x 1 n . . . x p n 1 C A f (X n ) = 0 B @ f (x 1 n ) . . . f (x p n ) 1 C A (27)
and the coupling matrix A: F is a map of J p = [01; 1] p 4 p into itself. The mathematical circuit of this system is displayed on Fig. 11(b).

A = 0 B B B B B B B B B B B B B @
" 1;1 = 1 0 j=p P j=2 " 1;j " 1;2 1 1 1 " 1;p01 " 1;p " 2;1 " 2;2 = 1 0 j=p P j=1;j$ =2 " 2;j 1 1 1 " 2;p01 " 2;p . . . . . . . . . . . . . . . . . . . . . . . .

" p;1 1 1 1 1 1 1 " p;p01 " p;p = 1 0 j=p01 P j=1 " p;j 1 C C C C C C C C C C C C C A (28)
It is noteworthy that these families of very weakly coupled maps are more powerful than the usual formulas used to generate chaotic sequences mainly because only additions and multiplications are used in the computation process; no division being required. Moreover the computations are done using oating point or double precision numbers, allowing the use of the powerful Floating Point Unit (FPU) of the modern microprocessors. In addition, a large part of the computations can be parallelized taking advantage of the multicore microprocessors which are common in laptop computers, nowadays.

We display in Table 1 the discrepancies in quadratic norm (Eq. 14) between the distribution of the iterated values x 1 n and the Lebesgue measure vs. the number of iterates N iter , for 2, 3 and 4-coupled symmetric tent maps. Computations are done using double precision numbers ( 14 0 15 digits), " i;j = i" 1 , j = 1; 4, " 1 = 10 014 , initial values: x 1 0 = 0:330000013113, x 2 0 = 0:338756413113, x 3 0 = 0:331353442113, x 4 0 = 0:333213583113. Remark: the ring coupler of Fig. 10(a) corresponds to an hollow matrix, in which only one diagonal and one other coe¢cient are not empty. The full coupler of the same gure stands for a matrix more lled with non vanishing coe¢cients as in Eq. 28.

Sampler

However chaotic numbers are not pseudo-random numbers because the plot of the couples of any component (x l n ; x l n+1 ) of iterated points (X n ; X n+1 ) in the corresponding phase plane reveals the map f used as one-dimensional dynamical systems to generate them via Eq. 26. Nevertheless we have recently introduced a family of enhanced Chaotic Pseudo Random Number Generators (CPRNG) in order to compute very fast long series of pseudorandom numbers with desktop computer 21 ;22 ;23 and its properties have been analyzed. 10 ;11 ;12 This family is based on the previous ultra-weak coupling which is improved in order to conceal the chaotic genuine function. In order to hide f of Eq. 26, in the phase space (x l n ; x l n+1 ); two mechanisms are used. The pivotal idea of the rst one mechanism is to sample chaotically the sequence (x l 0 ; x l 1 ; x l 2 ; ::: ; x l n ; x l n+1 ; :::) generated by the l-th component x l , selecting x l n every time the value x m n of the m-th component x m , is strictly greater (or smaller) than a threshold T 2 J, with l 6 = m, for 1 l; m p.

That is to say, to extract the subsequence (x l n (0) ; x l n (1) ; x l n (2) ; ::: ; x l n (q) ; x l n (q+1) ; :::) denoted here (x 0 ; x 1 ; x 2 ; ::: ; x q ; x q+1 ; :::) from the original one, in the following way:

Given

1 l; m p, l 6 = m ( n (01) = 01 x q = x l n (q) ; with n (q) = min r N fr > n (q01) j x m r > T g (29)
the chaotic mixing process. For sake of simplicity we have only displayed a circuit with three 1-dimensional generators. However the mixing process runs better when more generators are coupled. We display in Table 2 the discrepancies in quadratic norm (Eq. 14) between the autocorrelation distribution of the iterated values and the Lebesgue measure vs. the number of subsampled or mixed iterates N iter , for 4-coupled symmetric tent maps. Sam-pling: the rst component x 1 is sampled by x 4 for the threshold value 0:998. Mixing: the three components x 1 , x 2 , x 3 are mixed and sampled by x 4 for the threshold values T 1 = 0:998, T 2 = 0:9987, T 3 = 0:9994. Computations are done using double precision numbers ( 14 0 15 digits), " i;j = i" 1 , j = 1; 4, " 1 = 10 014 , initial values: x 1 0 = 0:330, x 2 0 = 0:3387564, x 3 0 = 0:50492331, x 4 0 = 0:0. We can say that the design of mathematical circuit including couplers, samplers or mixers allows the emergence of complexity in chaotic systems which leads to randomness. [START_REF] Lozi | Complexity leads to randomness in chaotic systems Mathematics in Science and Technology: Mathematical Methods, Models and Algorithms in Science and Technology[END_REF]

T

Reducer

The accelerated development of modern data transaction applications such as telecommunications requires encoding techniques with higher standards of security.

Classically, these encoding sequences are obtained using Pseudo Random Number Generators (PRNG). Since the seminal work presented in Sec. 2, an e¢cient alternative, the chaotic-based generators (CPRNG) are used to achieve even higher demanding encryption standards. Indeed, the chaotic systems exhibit a plethora of properties which make them suitable to meet the above requirements.

The advantage to use chaotic systems lies in their extreme sensitivity to small parameter and initial conditions variations: in this way, as many di¤erent chaotic carriers as wanted can be generated. However, the appropriate selection of a chaotic map that satises cryptographic applications requirements is a huge problem. It has to be emphasized that all chaotic maps are not applicable, because the chaotic gen-erator -which is deterministic-has to satisfy at the same time the criteria for closeness to random signals. Therefore many practical problems arise, from the choice of the chaotic generator and its parameters, to the chaotic properties verication after the quantisation. Ideally, for cryptographic applications and higher security, an everywhere dense chaotic attractor is required, so all chaotic signal samples will appear with the same probability. We have shown that highly e¢cient discrete-time chaotic generators can be obtained from coupling, sampling, and mixing quite simple models such as the logisic map or the symmetric tent map. We have built the corresponding mathematical circuit. Moreover there exist other combinations of such 1-dimensional chaotic attractor.

We introduce now, as an example, another combination which can directly provides random number without sampling or mixing, although it is possible to combine these processes with it.

The idea leading to this system is to conne on J p = [01; 1] p 4 p , considered as a torus, a ring of p0coupled symmetric tent map (or logistic map). [START_REF] Espinel | Dynamical and Statistical Analysis of a New Lozi Function for Random Numbers Generation[END_REF] Consider the equation

> > > <

> > > :

x 1 n+1 = 1 0 2 C C x 1 n C C + k 1 x 2 n x 2 n+1 = 1 0 2 C C x 2 n C C + k 2 x 3 n . . . x p n+1 = 1 0 2 jx p n j + k p x 1 n (31)
where the parameters k i 2 f01; 1g. In order to conne the variables x i n+1 on J p , we do, for every iteration the transform

if x i n+1 < 01; add 2 (32) if x i n+1 > 1; substract 2 (33)
We desing a new symbol: the reducer (on the right hand side of Fig. 14) in order to give a schematic representation of the projection of the variable on the torus J p . For sake of simplicity we have only displayed a circuit with three 1-dimensional generators. However this new pseudo-random number generator works better when more generators are coupled.

To evaluate the random properties of these generators, a set of statistical based test known as NIST test developed by the National Institute of Standards and Technology have been used.

The random properties validation of a 4dimensional system has been carried out. Addition-ally, the chaotic carrier output needs to be quantised and binarised (0 and 1) in order to be validated as being random using NIST tests. Therefore, di¤erent methods of binarisation (converting real signals to binary ones) have been implemented and compared.

A rst 10bit binarisation has been applied to the system (31) with n = 4, output: the results showed to be highly sensitive to the type of binarisation. Eventually, after testing several di¤erent methods, a 320bit binarisation has been chosen as being the most suitable solution. Because the system is conned to the p0dimensional torus, 31 bits are assigned to represent the decimal part, and 1 bit to the sign. To illustrate the results, the NIST tests for the four dimensional system (31) with parameters k i = (01) i+1 are shown in Fig. 15. The chosen conditions are:

Length of the original sequence: 10 8 bits, length of bit string: 10 6 , quantity of bit strings: 100. The output of the system has been arbitrary chosen as being: y = x 4 n . Furthermore, as the results show their independence from the initial conditions, every bit string in this test is the resulting sequence of a di¤erent randomly chosen initial condition.

The criterion for a successful test is that the p0value has to be superior to the signicance level (0:01 for this case). For the present model, all tests were successful thus the sequence can be accepted as being random. As we said above, in order to improve even more the random properties of that random signal, two possible strategies are possible: increasing the system order or under-sampling the output signal, which is possible with the circuit of Fig. 14 in which a sampler is added on hand right side.

Cascader

Finally, turning back to the problem of chaotic masking, via the cascading of Chuas attractor, the last symbol we design in order to schematise mathematically the cascading method is the cascader displayed on Fig. 16(a). In the limited extend of this chapter, we give only two examples of the engineery of mathematical circuits we have introduced in the previous Section. The rst one is an improvement of the cascading of two identical receivers of Chuas circuit. Albeit this improvement is rather from a numerical point of view than a practical one, it is given in order to illustrate more in deep the combination of circuit elements. It is possible to combine two cascading receiver as in Fig. 16(b).

The equations of this circuit are : y 5 = x 4 (t) 0 y 5 + z 5

:

z 5 = 0Cy 5 (34)
:

x 6 = B (y 5 (t) 0 x 6 0 f (x 6 )) (35)
:

x 8 = B (y 7 (t) 0 x 8 0 f (x 8 )) (37)
Then s(t) is recovered as

s 8 (t) = r (t) 0 x 8 (t) s (t) (38)
Numerical experiments show an improvement of the results. 1 ;19

Noise-resisting cryptography

The second example which highligt the complex combination of circuit element belongs to the new ourishing eld of chaotic based cryptography. We proposed last year, 2 a novel noise-resisting ciphering method resorting to a chaotic multi-stream pseudorandom number generator (denoted Cms-PRNG) described below. This Cms-PRNG co-generates an arbitrarily large number of uncorrelated chaotic sequences. These cogenerated sequences are actually used in several steps of the ciphering process. Noisy transmission conditions are considered, with realistic assumptions. The e¢ciency of the proposed method for ciphering and deciphering is illustrated through numerical simulations based on a Cms-PRNG involving ten coupled chaotic sequences.

The CPRNG decribed in Fig. 14, can be improved in order to generate uncorrelated sequences of pseudo-random numbers, possessing a large number of keys. This is simply obtained by adding a coupler as a keyer as in the circuit of Fig. 17 corresponding to Eq. 39 for 4-streams, or Eq. 40 for p0streams.

> > <

> > : The ring coupling which is expressed as a diagonal matrix in Eq. 31 is completed with numerous other non vanishing coe¢cients of the matrix used as secret keys in Eq. 40. Then the originality of the noise-resisting ciphering method introduced which uses this Cms-PRNG is twofold. First a novel ciphering method is proposed aimed at resisting to a noisy transmission channel. The main idea is to establish, between the transmitter and the receiver, a correspondence between the alphabet constituting the plain text and some intervals dening a partition of [01; 1]. Some realistic assumption about the noise boundedness allows to restrict the bounds of the aforementioned intervals in order to precisely re-sist to the e¤ects of the noise. An extra scrambling resorting to a co-generated chaotic sequence enhances the ciphering process. Then a new chaotic substitution method is developed: considering a chaotic carrier, belonging to the set of cogenerated and coupled pseudo-random chaotic sequences, the idea is to randomly/chaotically (in fact, this is determined by a second pseudo-random chaotic sequence) replace some elements of the carrier by a ciphered element (a letter here) of the message. At the receiver end, a copy of the Cms-PRNG, with the same parameters (hence we deal with a symmetrical ciphering method) allows to generate the necessary chaotic sequences and therefore to retrieve the initial message.

x 1 n+1 = 1 0 2 C C x 1 n C C + k 1 (1 0 " 1;3 0 " 1;4 ) x 2 n + " 1;3 x 3 n + " 1;4 x 4 n x 2 n+1 = 1 0 2 C C x 2 n C C + k 2 (1 0 " 2;4 0 " 2;1 ) x 3 n + " 2;4 x 4 n + " 2;1 x 1 n x 3 n+1 = 1 0 2 C C x 3 n C C + k 3 (1 0 " 3;1 0 " 3;2 ) x 4 n + " 3;1 x 1 n + " 3;2 x 2 n x 4 n+1 = 1 0 2 C C x 4 n C C + k 4 (1 0 " 4;2 0 " 4;3 ) x 1 n + " 4;2 x 2 n + " 4;3 x 3 n (39) 8 > > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > > : x 1 n+1 = 1 0 2 C C x 1 n C C + k 1 1 0 p P j=3 " 1;j ! x 2 n + p P j=3 " 1;j x j n ! . . . x m n+1 = 1 0 2 jx m n j + k m 1 0 p P j=1;j$ =m;m+1 " m;j ! x m+1 n + p P j=1;j$ =m;m+1 " m;j x j n ! . . . x p01 n+1 = 1 0 2 C C x p01 n C C + k p01 1 0 p02 P j=1 " p01;j ! x p n + p02 P j=1 " p01;j x j n ! x p n+1 = 1 0 2 jx p n j + k p 1 0 p02 P j=1 " p;j ! x 1 n + p02 P j=1 " p;j x j n ! (40)
This process can be summarized in both circuits of Figs. 18, 19. Due again to the limited extend of this chapter, we cannot expand these gures in order to show the constituting symbols in each oval shaped region of the circuit. The originality of the method lies in the use of a chaotic pseudo-random number generator: several co-generated sequences can be used at di¤erent steps of the ciphering process, since they present the strong property of being uncorrelated. Each letter of the initial alphabet of the plain text is encoded as a subinterval of [01; 1]. The bounds of each interval are dened in function of the known bound of the additive noise. A pseudo-random sequence is used to enhance the complexity of the ciphering. The transmission consists of a substitution technique inside a chaotic carrier, depending on another cogenerated sequence. The e¢ciency of the proposed scheme is illustrated on some numerical simulations. [START_REF] Cherrier | Noise-resisting ciphering based on a chaotic multi-stream pseudorandom number generator[END_REF] As further work, some studies should be performed of several sets of unknown parameters, since with the considered Cms-PRNG with 10 states, the number of possible parameters amounts to 90 (the " i;j and the k i ).

Conclusion

Following the worldwide tradition of use of Chuas circuits for various purposes, we have introduced the paradigm of chaotic mathematical circuitry which shows some similarity to the paradigm of electronic circuitry -the design of electronic circuits. This new paradigm allows, as an example, the building of new chaotic and random number generators. In this beginning of the third Millenium, the old tradition of design of electronic circuits is drastically revolutionized by the introduction by L. O. Chua of new components with memory, namely memristor, memcapacitor and meminductor. These devices are common at the nanoscale and their combination in circuits open up new functionnalities in electronics. Apart from the obvious use of these elements in nonvolatile memories, several applications can be already envisioned, especially in neuromorphic devices to stimulate learning, adaptative and spontaneous behavior. [START_REF] Ventra | Circuit elements with memory: memristors, memca-pacitors, and meminductors[END_REF] Alongside to this electronic circuits revolution, the new theory of mathematical circuits allows many new applications in chaotic cryptography, genetic algorithms in optimization and in control,... Due to the versatility of the new components we introduce, the combined operation of these chaotic mathematical circuits is still largely unexplored. We hope our work will motivate experimental and theoretical investigations in this direction. For the 75th birthday of Professor L. O. Chua, one can highlight that his seminal work since more than forty years is fruitful exceeding his wide-ranging eld of research. In this chapter (as suggested in the title of the anniversary book) we followed him from chaos to memristor and we went beyond.

Fig. 1 .

 1 Fig. 1. (a) Realization of Chuas circuit from. 5 (b) Three-segment piecewise-linear v-i characteristic of nonlinear resistor in Chuas circuit. (c) Chua attractor.

Fig. 2 .

 2 Fig. 2. Experimental set up. Block diagram of the system. It contains one Chuas circuit and two partial Chuas circuits, that is, subsystems #1 et #2 of Fig.3, from.1 5

Fig. 3 .

 3 Fig. 3. Practical realization of the receiver. The rst subsystem is a partial Chuas circuit consisting of the 0 v C 2 ; i L 1 -subsytem driven by the transmitted signal r (t). The second subsystem is a partial Chuas circuit consisting of the 0 v C 1 1 -subsytem driven by the transmitted signal v

Fig. 5 .

 5 Fig. 5. Block diagram of the electronic circuit implemented in Fig. 4.

Fig. 6 .

 6 Fig. 6. Symbols of the three recently discovered circuit elements with memory.

  7(c) and (d))H a;b :x y = y + 1 0 ax 2 bx (15)with a = 1:4; b = 0:3;

Fig. 7 .

 7 Fig. 7. Generators. (a) Continuous generator (Chuas circuit, expanded symbol). (b) Simplied continuous generator. (c) Discrete generator (Hénon map, expanded symbol). (d) Simplied discrete generator. (e) 1-dimensional logistic generator.

Fig. 8 .

 8 Fig.8. Five identical coupled Chuas circuits forming a ring, from.1 4

  (a)), and the full coupler when the coupling involves more connections between the couplers (Fig 10(b)). Remark 1.1. In the rest of this article, we use solid

Fig. 9 .

 9 Fig. 9. Numerical circuit corresponding to electronic circuit of Fig. 8.

Fig. 10 .

 10 Fig. 10. (a) Ring coupler. (b) Full coupler.

Fig. 11 .

 11 Fig. 11. (a) Circuit of ultraweak coupling of three 1-dimensional chaotic maps. (b) Circuit of ultraweak coupling of p 1dimensional chaotic maps.

Fig. 12 .

 12 Fig. 12. Circuit of enhanced Chaotic Pseudo Random Number Generators (CPRNG) based on chaotic sampling.

Fig. 14 .

 14 Fig. 14. Reducer for the circuit of Eq. 31, with n = 3.

Fig. 15 .

 15 Fig. 15. Successful result of the NIST tests.

Fig. 16 .

 16 Fig. 16. (a) Cascader symbol. (b) Two cascading receivers combined.

Fig. 17 .

 17 Fig. 17. Chaotic multistream PRNG (Cms-PNRG).

Fig. 18

 18 Fig. 18. Transmitter

The sequence (x 0 ; x 1 ; x 2 ; ::: ; x q ; x q+1 ; :::) is then the sequence of chaotic pseudo-random numbers.

The mathematical formula (29) can be best understood in algorithmic way. The pseudo-code, for computing iterates of (29) corresponding to N iterates of Eq. 26 is: X 0 = (x 1 0 ; x 2 0 ; :::; x p01 0 ; x p 0 ) = seed n = 0; q = 0; do { while n < N do { while (x m n T ) compute (x 1 n ; x 2 n ; :::; x p01 n ; x p n ); n++} compute (x 1 n ; x 2 n ; :::; x p01 n ; x p n ); then n(q) = n; x q = x 1 n(q) ; n++; q++} This chaotic sampling is possible due to the independence of each component of the iterated points X n vs. the others. [START_REF] Lozi | New Enhanced Chaotic Number Generators[END_REF] We introduce the symbol on the right hand side of Fig. 12 in order to give a schematic representation of this sampling (also called subsampling) process.

Mixer

A second mechanism can improve the unpredictability of the pseudo-random sequence generated as above, using synergistically all the components of the vector X n , instead of two. Given p 0 1 thresholds

and the corresponding partition of

simple mechanism is based on the chaotic mixing of the p 0 1 sequences (x 1 0 ; x 1 1 ; x 1 2 ; ::: ; x 1 n ; x 1 n+1 ; :::), (x 2 0 ; x 2 1 ; x 2 2 ; ::: ; x 2 n ; x 2 n+1 ; :::), ... , (x p01 0 ; x p01 1 ; x p01 2 ; ::: ; x p01 n ; x p01 n+1 ; :::), ... using the last one (x p 0 ; x p 1 ; x p 2 ; ::: ; x p n ; x p n+1 ; :::) in order to distribute the iterated points with respect to this given partition dening the subsequence (x 0 ; x 1 ; x 2 ; ::: ; x q ; x q+1 ; :::) by 8 < :

The pseudo-code, for computing the iterates of Eq. 30 corresponding to N iterates of Eq. 26 is X 0 = (x 1 0 ; x 2 0 ; :::; x p01 0 ; x p 0 ) = seed n = 0; q = 0 ; do { while n < N do {while (x p n 2 J 0 ) compute (x 1 n ; x 2 n ; :::; x p01 n ; x p n ); n++} compute (x 1 n ; x 2 n ; :::; x p01 n ; x p n ) let k be such that x p n 2 J k then n(q) = n; x q = x k n(q) ; n++; q++} We introduce the symbol on the right hand side of Fig. 13 in order to give a schematic representation of