Ola Hammoud
email: ola.hammoud@insa-lyon.fr

Denise Blanc

Maria Lupsea-Toader

Christine De Brauer

INFLUENCE OF CHEMICAL PRE-TREATMENT ON THE LEACHING BEHAVIOUR OF BOTTOM ASH

Keywords: Municipal solid waste incineration; bottom ash, mineral trace elements, chemical treatment, leaching behaviour, geochemical modelling

published or not. The documents may come

Introduction

The production of municipal solid waste has been increasing over time due to urbanization. Efficient management of this waste is required in order to reduce its volume. Currently, the most spread process for the treatment of municipal solid waste is incineration, which is a thermal method that aims to recover energy from waste. Municipal waste incineration produces two types of solid residues: bottom ash, called MSWI-BA here after, and fly ash issued from air pollution control. All incineration residues are toxic, as they contain hazardous compounds: dioxins, mineral trace elements (MTE) and soluble salts such as chlorides and sulphates [START_REF] Saikia | Pre-treatment of municipal solid waste incineration (MSWI) bottom ash for utilisation in cement mortar[END_REF]. While the valorisation of fly ash is limited because of its high content in pollutant compounds [START_REF] Tzu | Determination of the Pb, Cr, and Cd distribution patterns with various chlorine additives in the bottom ashes of a low-temperature two-stage fluidized bed incinerator by chemical sequential extraction[END_REF], the bottom ash is widely used as backfill material [START_REF] Lin | Effect of incinerator bottom-ash composition on the mechanical behavior of backfill material[END_REF][START_REF] Al-Rawas | Use of incinerator ash as replacement for cement and sand in cement mortars[END_REF]. In France, the mostly spread application of bottom ash are the road base formulations. The thresholds for such valorisation are stated in the French decree from 18/11/2011 [START_REF]Arrêté du 18 novembre 2011 relatif au recyclage en techniques routières des mâchefers d'incinération de déchets non dangereux; NOR: DEVP1131516A[END_REF]. The pollutant potential of MSWI-BA is evaluated by leaching tests performed according to the EN 12457-2 standard [START_REF]European standard. EN 12457-2: Characterization of waste -Leaching -Compliance test for leaching of granular waste materials and sludges -Part 2: One stage batch test at a liquid to solid ratio of 10 L/kg for materials with particle size below 4 mm (without or with size reduction[END_REF].

A more interesting but not developed valorisation domain is the constructionbottom ash could be integrated in construction products' formulations. In France, no regulation exists in this sense yet. In order to define reliable formulations, the material needs to be characterised at micro scale. MSWI-BA is composed of particles of variable sizes, including glass fragments, pieces of ceramic, unburned pieces and scrap. It is porous, angular with a rough texture and lighter than natural aggregates [START_REF] Singh | Effect of coal bottom ash as partial replacement of sand on properties of concrete[END_REF]. From the chemical point of view, the bottom ash has a heterogeneous nature, which differs depending on the origin of the waste and the type of the incinerators used during the process. MSWI-BA is in majority composed of silica, iron, aluminium, calcium, sulphates, chlorides etc. Besides these major compounds, bottom ash contains also MTE in small quantities, such as arsenic, copper, antimony, cadmium etc. [START_REF] Singh | Effect of coal bottom ash as partial replacement of sand on properties of concrete[END_REF][START_REF] Kim | Utilization of sieved and ground coal bottom ash powders as a coarse binder in high-strength mortar to improve workability[END_REF][START_REF] Siddique | Use of municipal solid waste ash in concrete[END_REF]. The high leachability of these trace elements from MSWI-BA can lead to problems related to human health and the environmentonce leached at contact with water, MTE can diffuse in the ground and underground water [START_REF] Nikbin | A comprehensive analytical study on the mechanical properties of concrete containing waste bottom ash as natural aggregate replacement[END_REF].

In this context, it is important considering pre-treatment methods of bottom ash before its valorisation. In France, only techniques aiming at the extraction of the soluble fraction are recognised by the authorities to open new ways of valuation as secondary raw material. As a consequence, our approach consists in chemical washings. The chemical washing of MSWI-BA has been previously tested by several authors [START_REF] Aguiar | Wet extraction of heavy metals and chloride from MSWI and straw combustion fly ashes[END_REF][START_REF] Lin | Effects of water washing on removing organic residues in bottom ashes of municipal solid waste incinerators[END_REF][START_REF] Sorlini | Recovery of MSWI and soil washing residues as concrete aggregates[END_REF][START_REF] Yang | Basic characteristics of leachate produced by various washing processes for MSWI ashes in Taiwan[END_REF]. It has been shown that very fine particles and mobile organic and inorganic contaminants (organic acids, MTE and salts) potentially leachable, can be removed by washing, or stabilised in the solid matrix, thus reducing the released quantities from the treated material. The present work serves three purposes: (i) to study the influence of different washing parameters on the leachable fraction of the main MTE, (ii) to optimize and implemente a treatment including several washing and rinsing steps, (iii) to model the washing protocol using the geochemical modelling tool PhreeqC ® , in order to better understand the geochemical composition of the MSWI-BA.

Materials and methods

About the studied bottom ash

The bottom ash used in this study was recovered from an energy enhancement unit located in the Auvergne-Rhône-Alpes region in France. Most of the waste comes from household, industry and water treatment plants. The furnace is equipped with a moving grate incinerator, which produces about 60 tons of bottom ash per day. 300 kg of MSWI-BA were recovered for this studyfor representative sampling, 10 to 20 kg bottom ash were collected daily over one month. The MSWI-BA were then dried for one week at 80 °C. Ferrous and non-ferrous metals were removed by magnetic separation and by eddy current. The samples were prepared by quartering in order to obtain homogeneous subsamples. Three particle size fractions were prepared by sieving, respectively 0-2 mm fraction (V<2), 2-5 mm fraction (V2-5) and larger than 5 mm fraction (V>5).

Leaching behaviour of bottom ash

Based on the French decree from 18/11/2011 [START_REF]Arrêté du 18 novembre 2011 relatif au recyclage en techniques routières des mâchefers d'incinération de déchets non dangereux; NOR: DEVP1131516A[END_REF], the leachability of the studied material before and after chemical treatment was evaluated, according to the European standard leaching test EN 12457-2 [START_REF]European standard. EN 12457-2: Characterization of waste -Leaching -Compliance test for leaching of granular waste materials and sludges -Part 2: One stage batch test at a liquid to solid ratio of 10 L/kg for materials with particle size below 4 mm (without or with size reduction[END_REF]. Powdered material was agitated in a rotary eluent in contact with demineralised water at a liquid-to-solid ratio (L/S) of 10 L.kg -1 during 24 h. The mixture was next filtered through a 0.45-µm cellulose acetate membrane. The concentrations of several major and trace elements in the eluates were measured by Inductively Coupled Plasma-Mass Spectrometry (ICP-MS) for metals and Ion Chromatography (IC) for the anions. All elements mentioned in the 18/11/2011 decree [START_REF]Arrêté du 18 novembre 2011 relatif au recyclage en techniques routières des mâchefers d'incinération de déchets non dangereux; NOR: DEVP1131516A[END_REF] were measured. The leaching test was performed in triplicate.

Optimisation of the chemical washing process

Various washing parameters have been studied: the nature of the washing solution and its concentration, the washing duration and the L/S ratio. The washing efficiency has been evaluated on all granulometric fractions [START_REF] Kaufman | Etude de l'influence du lavage des Mâchefers d'Incinération de Déchets non Dangereux sur le comportement des polluants métalliques traces[END_REF]. Table 1 shows the values tested for the various parameters. After each test, the material was vacuum filtered through a 0.45-µm membrane, dried at 80 °C during 24 hours and tested according to the EN 12457-2 standard described above [START_REF]European standard. EN 12457-2: Characterization of waste -Leaching -Compliance test for leaching of granular waste materials and sludges -Part 2: One stage batch test at a liquid to solid ratio of 10 L/kg for materials with particle size below 4 mm (without or with size reduction[END_REF]. The concentrations of MTE and anions in the eluates were measured by ICP-MS and CI. The optimal set of parameters has been chosen based on the efficiency of the treatment process (see results in the next section), i.e. hydrochloric acid solution at 0.18 mol.L -1 , during 6 hours with a L/S ratio of 10 L.kg -1 (in bold in table 1).

Pre-treatment protocol

A general pre-treatment protocol, comprising several stages of pre-wash, extraction and rinsing, has been studied. This protocol includes the following steps:  Pre-wash with water: stirring during 60 min; L/S = 2 L.kg -1  Extraction with HCl (washing): 0.18M, stirring 60 min; L/S = 2 L.kg -1  First rinsing with water: stirring 60 min; L/S = 2 L.kg -1  Second rinsing with water: stirring 60 min; L/S = 2 L.kg -1 .

After each step, the eluates have been separated by centrifugation and then by filtration on a 0.45 μm cellulose acetate membrane. The concentrations of MTE and anions in the filtrates were analysed by ICP-MS and CI, respectively. Next, the MSWI-BA samples were dried at 80 °C for 24 hours and submitted to a leaching test with water according to the EN 12457-2 standard.

Geochemical modelling

All stages of the treatment protocol described above, followed by leaching tests, have been modelled using the specific geochemical modelling software PhreeqC ® version 3 [START_REF] Parkhurst | User's guide to PHREEQC (version 2) -A computer program for speciation, batchreaction, one-dimensional transport, and inverse geochemical calculations[END_REF], coupled with the LLNL database (Lawrence Livermore National Laboratory). The solid is represented by an assembly of mineral phases, called mineralogical assemblage; it is put into contact with the liquid phase, called solution. The mineralogical assemblage defined by Gonzalez [START_REF] Gonzalez | étude du comportement des polluants inorganiques traces au cours de la maturation et en fonction de la granulométrie[END_REF] for the modelling of MSWI-BA has been used as input data in this study. The calculations were made assuming thermodynamic equilibrium. All steps of the washing protocol defined above have been modelled with the corresponding parameters (L/S, solution, concentration etc.). The mineralogical assemblage resulting from one step is used as input data for the following step. The pH of the liquid phases reached at thermodynamic equilibrium and concentrations of different elements (major elements, MTE and anions) are compared to the concentrations measured experimentally. A final step simulates the leaching test performed on treated bottom ash.

Results and discussions

Leaching behaviour of bottom ash before chemical treatment

The leaching behaviour of MSWI-BA has been tested as described in the sections above, according to the EN 12457-2 standard [START_REF]European standard. EN 12457-2: Characterization of waste -Leaching -Compliance test for leaching of granular waste materials and sludges -Part 2: One stage batch test at a liquid to solid ratio of 10 L/kg for materials with particle size below 4 mm (without or with size reduction[END_REF]. Figure 1 shows the results for the three granulometric fractions studied here. The thresholds for secondary road base use [START_REF]Arrêté du 18 novembre 2011 relatif au recyclage en techniques routières des mâchefers d'incinération de déchets non dangereux; NOR: DEVP1131516A[END_REF] are satisfied for all MTE, as well as for sulphates, for all granulometric fractions. Only the chlorides exceeded the thresholds for fine and intermediate particle size.

For most elements (except antimony), the release increases with the decrease of the particle size, thus the most charged fraction is the V<2. For the majority of the elements the fine fraction (V<2) is the most charged with MTE [START_REF] Bethanis | Effect of sintering temperature on properties and leaching of incinerator bottom ash[END_REF][START_REF] Chimenos | Characterization of the bottom ash in municipal solid waste incinerator[END_REF][START_REF] Valle-Zermeño | Material characterization of the MSWI bottom ash as a function of particle size. Effects of glass recycling over time[END_REF], thus presenting the highest release. Vaporisation and condensation during the incineration process of wastes, leads to the accumulation of MTE in the fine particle size as metal oxides or combined with anions [START_REF] Valle-Zermeño | Material characterization of the MSWI bottom ash as a function of particle size. Effects of glass recycling over time[END_REF].

Less MTE are released from coarse particles, compared to the fine fraction for two main reasons: (i) they contain a high percentage of glass fragments and unburned pieces and (ii) the available contact surface is smaller than that of the fine fraction. Coarse particles might contain metal or metal alloys (e.g. copper wires, lead base alloys for soldering), metal oxide particles adherent to the surface of larger particles, or as additives used in the formulation of ceramics.

Influence of the different washing parameters on the chemical pre-treatment of bottom ash

In this section, the optimisation of various parameters of the chemical pre-treatment process, i.e. the nature of the washing solution, the concentration of the chosen washing agent, the L/S ratio and the contact time is presented. All tests have been performed in triplicate and are plotted along with the standard deviation, where this is coherent. The results are expressed in mg of leached element per kg of dry matter (DM) of bottom ash and are compared to the leaching fraction of raw MSWI-BA (V<2) and to the secondary road base thresholds.

Optimisation of the washing solution

Three different washing solutions were tested on the finest fraction of MSWI-BA (V<2) which has exhibited the highest released concentrations during the leaching tests (see figure 1). Figure 2 shows that the behaviour of elements varies depending on the washing solution. However, not all exhibit the same behaviour, e.g. arsenic, barium, molybdenum and nickel leach in higher concentrations after treatment with nitric acid, while lead and cadmium are more soluble after treatment with EDTA. HNO3 is known to be a strong oxidizing agent. On the other hand, EDTA plays the role of a complexing agent known for its ability to chelate metals [START_REF] Hu | Assessment of EDTA heap leaching of an agricultural soil highly contaminated with heavy metals[END_REF]. Concerning HCl, the chloride anions might complex with lead ions forming a water insoluble species, e.g. PbCl2 [START_REF] Huang | Leaching behavior of heavy metals with hydrochloric acid from fly ash generated in municipal waste incineration plants[END_REF]. EDTA proved to be inefficient in stabilising MTE, thus the release is higher after treatment (figure 2). Moreover, EDTA is not biodegradable and more expensive. HCl proved to be the best treatment agent, which reduces the release for most MTE. Hydrochloric acid has been selected for further tests. Fig. 2 Leachable fraction before (raw-BA) and after washing with HCl, HNO3 and EDTA; in hatched the values below or equal to the limit of quantification All three chemical agents used have interacted with the solid matrix either to solubilize the elements, thus reducing their total content), or to destabilise certain components, thus increasing their release during the following leaching test. In all cases, solubilised heavy metal ions might have been trapped on the solid matrix during the treatment process. A rinsing step proved to be necessary, in order to remove all traces of the chemical agent and other dissolved elements. The global protocol is described later.

Optimisation of the concentration of hydrochloric acid

Four different concentrations (0.18 M, 0.5 M, 1 M and 1.5 M) have been tested for the chosen washing agent, i.e. hydrochloric acid. The results are presented in figure 3 for the fine fraction (V<2). The efficiency of the washing process decreases with the increase of the concentration of acid for most elements. During the acid extraction MTE are replaced by H + ions, thus becoming soluble. Chromium and antimony present a different behaviour. They are known to be stable in this range of pH, i.e. low alkaline pH, which explains the small variations with the concentration of the acid. Concerning the behaviour of Cl -ions, their concentration increases with the increase of concentration of the washing agent. On the other hand, sulphates are soluble in acid solutions and easily leachable.

Optimisation of the L/S ratio

Three L/S ratios have been tested on the V<2 fraction using HCl 0.18 M solution during 6 hours (see figure 4). With some exceptions (chrome and antimony), the efficiency of washing increases with the decreasing of the L/S ratio. This behaviour has been also observed on fly ash by Huang et al. [START_REF] Huang | Leaching behavior of heavy metals with hydrochloric acid from fly ash generated in municipal waste incineration plants[END_REF]. The quantity in moles of acid compared to the material mass increases with the L/S ratio, which leads to a smaller alkaline pH in the eluates. The leachability of most MTE increases with the decrease of pH in this range [START_REF] Bethanis | Effect of sintering temperature on properties and leaching of incinerator bottom ash[END_REF][START_REF] Johnson | Acid neutralizing capacity of municipal waste incineration bottom ash[END_REF][START_REF] Bergfeldt | Investigation of Sintering Processes in Bottom Ash to Promote the Reuse in Civil Construction (Part 2) -Long term behaviour -Waste Materials in Construction[END_REF]. Concerning the anions, the test has been performed only at ratio L/S = 2 L.kg -1 . For economic and environmental reasons, the L/S = 2 L.kg -1 value has been selected for the global chemical treatment process of bottom ash. show small variations in the release of elements. Whatever the contact time, a reduction in the overall metal concentration in the eluates is observed relative to the initial release from raw bottom ash. Similar results have been obtained by Huang et al. [START_REF] Huang | Leaching behavior of heavy metals with hydrochloric acid from fly ash generated in municipal waste incineration plants[END_REF] on fly ashmore than 60% of MTE were released for low contact time with HCl, thus further leaching reveals unnecessary. The contact time of 1 hour has been chosen for this study.

Influence of the washing on the particle size fraction

The leaching of all considered particle size fractions before and after treatment with HCl 0.18 M solution has been tested (see figure 5). Similar results have been obtained by Gonzalez [START_REF] Gonzalez | étude du comportement des polluants inorganiques traces au cours de la maturation et en fonction de la granulométrie[END_REF] and Chen et al. [START_REF] Chen | Size distribution of metals in bottom ash of municipal solid waste incinerators[END_REF]. These results are very similar to those presented in figure 1, performed with the same experimental parameters, i.e. HCl, 0.18 M, L/S = 2 L.kg -1 , 6 hours (blue vs. red representation), thus confirming the reproducibility of the tests. Except arsenic, chrome and antimony, the general behaviour shows that for all particle size fractions the release of MTE decreases after treatment.

Fig. 5 Leachable fraction before (raw-BA) and after treatment for all particle size fractions; in hatched the values below or equal to the limit of quantification

Global chemical pre-treatment process for bottom ash

The optimisation study presented above allowed the development of a complete washing protocol, in order to meet the industrial demands (efficient, economic and environmental friendly). The most efficient extracting agent proved to be a slightly concentrated hydrochloric acid (0.2 M). Parameters such as contact time and L/S ratio presented small variations. For economic and environmental reasons (excessive production of liquid effluents to be treated later), the washing time was fixed at 1 hour and the L/S ratio at 2 L.kg -1 . The washing step is followed by two steps of rinsing with water, in order to remove the pollutants which have been previously solubilised and trapped in the solid phase. The global protocol has been detailed in the section § Materials and methods. This protocol has been next modelled; the results are presented below.

Modelling of the pre-treatment process of bottom ash

Study of the evolution of pH and the mineralogical assemblage

Figure 6 shows the evolution of the pH for each step of the treatment process. The lowest pH is observed for the extraction step, which is justified by the acidity of the extracting agent. The eluate however, has an alkaline pH (8.32), which confirms the large buffer capacity of bottom ash. The increasing of pH during the two rinsing steps and the final leaching show that the thermodynamic equilibrium is not achieved. Despite previous washing steps, the pH obtained in the leaching test remains above 10.

Concerning the theoretical evolution of the pH, the model shows that after the extraction step the pH remains constant at 11.21. This is explained by the fact that the calculations are performed with the hypothesis of thermodynamic equilibrium. Several phases are completely dissolved during the extraction step, conducting to an identical assemblage as input data for the first and the second rinsing. The pH obtained after the final leaching is slightly different (11.23), as the leaching ratio L/S has changed from 2 to 10 L.kg -1 . The 2 pH units difference between the model and the experiment for the extraction step could also be explained by the thermodynamic hypothesis used in the model. Bottom ash has shown a strong buffer capacity, taking up to a week to reach the thermodynamic equilibrium [START_REF] Gonzalez | étude du comportement des polluants inorganiques traces au cours de la maturation et en fonction de la granulométrie[END_REF]. This hypothesis being too restrictive, an improvement of the model could be obtained by modelling dissolution kinetics for several phases.

The leaching behaviour of major elements

The methodology for modelling such complex materials starts by validating the mineralogical assemblage based on the chemical behaviour of major elements present in the solid phase [START_REF] Islam | Leaching and redistribution of Cu and Pb due to simulated road runoff assessed by column leaching test, chemical analysis, and PHREEQC modeling[END_REF][START_REF] Cornelis | Antimony leaching from MSWI bottom ash: Modelling of the effect of pH and carbonation[END_REF][START_REF] Peyronnard | Study of mineralogy and leaching behavior of stabilized/solidified sludge using differential acid neutralization analysis. Part II: Use of numerical simulation as an aid tool for cementitious hydrates identification[END_REF]. Therefore, the leaching behaviour of the major elements present in the bottom ash was evaluated during the different stages of the washing process. Figure 7 shows that both experimental and modelled results have the same order of magnitude. However, depending on the element, the results are more or less convincing. Elements for which the mineral assemblage is almost complete, like calcium, aluminium and sulphur are well represented. On the other hand, the mineral assemblage for silicium is more complex, and the results are less convincing for this element.

The leaching behaviour of MTE

The chemical treatment applied in this study has proven efficient for all MTEthe release has decreased, compared to the leaching behaviour of raw bottom ash (see figure 8). Modelling has led to the same conclusions as in the case of elements, namely that well-defined elements in the assemblage, such as zinc, show good coherence between experimental and simulated leaching behaviour. On the other hand, elements such as lead or copper, are not well modelled. Further characterisation methods, such as acid neutralization capacity (ANC) leaching tests [START_REF] Afnor | A. Characterization of waste -Leaching behavior tests -Influence of pH on leaching with initial acid/base addition[END_REF] or sequential extractions [START_REF] Tessier | Sequential extraction procedure for the speciation of particulate metals[END_REF], combined with characterisation technical (XRD, SEM etc.), need to be used, in order to affine the knowledge of the mineralogical assemblage of this material.

The leaching behaviour of chlorides and sulphates

The results obtained for anions are presented in figure 9. An important quantity of chlorides and sulphates is solubilised during the pre-wash step, mainly because most of the NaCl and KCl can be removed by simple water washing [START_REF] Huang | Leaching behavior of heavy metals with hydrochloric acid from fly ash generated in municipal waste incineration plants[END_REF]. This solubilisation is well represented by modelling. However, because of the thermodynamic equilibrium hypothesis, all chlorides are solubilised in the first two steps of the treatment, which leads to zero release during the rest of the protocol. This observation is not valid for experimental data.

Conclusions

For accomplishing the valorisation of bottom ash in construction, efforts need to be made in order to stabilise this material over time. In this sense, a chemical pre-treatment has been developed in the current study. The influence of different washing parameters on the leachable fraction of the main MTE was studied. A complete washing protocol consisting of 4 steps (pre-wash, acid extraction and two rinsings) has been defined. A dynamic rinsing would be necessary in order to remove the traces of metals solubilized during the extraction with HCl and adsorbed or wedged on the solid phase. It has been observed that a low L/S ratio increases the efficiency of the pre-treatment process, reducing the production of effluents. On the other hand, the contact time does not have a great impact on the efficiency of the treatmenta short contact time is preferable. Finally, the modelling of the complete washing protocol helped at understanding the geochemical composition of this complex material. The results confirm that in order to make this model reliable, it is imperative to define a complete mineralogical assemblage of MSWI-BA.

Fig. 1

 1 Fig. 1 Leaching behaviour of the three fractions (cf. NF EN 12457-2); comparison with secondary road base thresholds according to the 18/11/2011 decree; in hatched the values below or equal to the limit of quantification

Fig. 3

 3 Fig. 3 Leachable fraction before (raw-BA) and after washing (fraction V< 2) with HCl at 0.18 M, 0.5 M, 1 M and 1.5 M; in hatched the values below or equal to the limit of quantification

Fig. 4

 4 Fig. 4 Leachable fraction before (raw-BA) and after washing of the fraction V<2 for the L/S ratios of 2, 5 and 10 L.kg -1 ; in hatched the values below or equal to the limit of quantification Optimisation of the washing time

Fig. 6

 6 Fig. 6 Evolution of the pH (experimental vs. modelling) over the various stages of the pre-treatment process

Fig. 7

 7 Fig. 7 Leachable fractions of major elements during the treatment protocol; experimental (not hatched) and modelled (hatched) results

Fig. 8

 8 Fig. 8 Evolution of the leachable fraction of MTE during the treatment protocol; experimental (not hatched) and modelled (hatched) results

Fig. 9

 9 Fig. 9 Evolution of the leachable fraction of anions during the treatment protocol; experimental (not hatched) and modelled (hatched) results

Table 1

 1 Values of the different parameters tested during the washing treatment of MSWI-BA

	Extractant (concentration)	HCl (0.18 M)	HNO3 (0.5 M)	EDTA (0.05 M)
	Washing time (h)	1	3	6	24
	HCl concentration (mol.L -1 )	0.18	0.50	1.00	1.50
	L/S ratio (L.kg -1 )	2	5	10
	Granulometric fraction (mm)	V<2	V2-5	V>5


Jianhua Fan
email: jianhua.fan@insa-rennes.fr

Franck Lominé

Mustapha Hellou

A numerical analysis of pressure drop and particle capture efficiency by rectangular fibers using LB-DE methods

Keywords: lattice Boltzmann method, discrete element method, capture efficiency, rectangular fiber, particle transport

In this work, a coupled lattice Boltzmann method (LBM) and discrete element method (DEM) are used to simulate the particle transport and deposition on rectangular fibers of clean filter. The LBM is employed to describe the fluid flow around the fibers, whereas the DEM is used to deal with the particle dynamics. The effects of the Reynolds number, the fiber aspect ratio and the arrangement of fibers (i.e., orientation angle of a fiber) on the pressure drop and capture efficiency are investigated at the initial stage of the filtration process. The quality factor, commonly used to determine the filtration performance, is also studied. The simulation results illustrate that both pressure drop and capture efficiency are dependent on the orientation angle and aspect ratio. The Reynolds number has only a slight influence on the capture efficiency but has a significant effect on the pressure drop for high aspect ratio. A good filter performance can be obtained for square fiber when the orientation angle is π/4 from the quality factor standpoint.

Introduction

Filtration is an important process for removing particles from a fluid. It is widely applied in many fields, such as water treatment, reactors for catalysis, particle transport in porous media and air purification. Thus, the optimization of the filters is a crucial challenge and should cover fluid and transported particles at the same time. Filtration by fibrous filter has gained great attention among researchers and industries due to its low cost and easy manufacture [START_REF] Zhu | Inertial impaction-dominated fibrous filtration with rectangular or cylindrical fibers[END_REF][START_REF] Wang | The influence of fiber geometry and orientation angle on filtration performance[END_REF]. Some mathematical models have been proposed and experimental investigations have been carried out to predict and improve pressure drop and capture efficiency, since they are two important parameters to evaluate the performance of a filter [START_REF] Wang | Filtration of aerosol particles by elliptical fibers: A numerical study[END_REF]. Fiber filtration involves interaction between transported particles, the fiber and the fluid. Most initial articles on the subject were concerned with circular fibers and used stream function developments to solve fluid velocity field. For example, Kuwabara [START_REF] Kuwabara | The forces experienced by randomly distributed parallel circular cylinders or spheres in a viscous flow at small reynolds numbers[END_REF] gave theoretical formula of drag force by solving the velocity field of viscous Stokes flow around a circle. Kuwabara considered a circular outer boundary of the domain. This hypothesis can be valid for solitary fiber or for a weak density of fibers, but fails when the gap between the fibers is smaller than a given threshold. Indeed, in such a case, the outer boundary does not cover all the fluid domain. Lee and Liu [START_REF] Lee | Theoretical study of aerosol filtration by fibrous filters[END_REF] based on a boundary layer approach using the Kuwabara flow field combine the diffusion and interception mechanisms together to drive an expression for predicting the single fiber efficiency. Schweers et al. [START_REF] Schweers | Experimental investigation of particle collection on single fibres of different configurations[END_REF] experimentally determined the collection efficiency of single circular fiber and two circular fibers in different arrangements. Augusto et al. [START_REF] Augusto | Microfiber filter performance prediction using a lattice Boltzmann method[END_REF] used the bounce-back scheme and slip condition to simulate the pressure drop through filters with the OpenLB library. A good agreement is observed by comparing with experimental and numerical results. Some studies include, but are not limited to, pressure drop and capture efficiency of Liu and Wang [START_REF] Liu | Pressure drop and interception efficiency of multifiber filters[END_REF], capture efficiencies due to three mechanisms (Brownian diffusion, interception and inertial impaction) of Stechkina et al. [START_REF] Stechkina | Studies in fibrous aerosol filters iv. calculation of aerosol deposition in model filters in the range of maximum penetration[END_REF], numerical simulation work of Fan et al. [START_REF] Fan | Numerical study of particle capture efficiency in granular filter[END_REF], Rabiee et al [START_REF] Rabiee | Investigation of the characteristics of particulate flows through fibrous filters using the lattice boltzmann method[END_REF], Wang et al. [START_REF] Wang | Numerical simulation of particle capture process of fibrous filters using lattice boltzmann two-phase flow model[END_REF] as well as experimental work of Lee and Liu [START_REF] Lee | Theoretical study of aerosol filtration by fibrous filters[END_REF].

However, with the development of filter manufacturing technology, researchers recently showed interest for the study of filtration by non-circular fibers [START_REF] Wang | Filtration of aerosol particles by elliptical fibers: A numerical study[END_REF][START_REF] Jafari | Investigation of particle dispersion and deposition in a channel with a square cylinder obstruction using the lattice boltzmann method[END_REF][START_REF] Hosseini | On the importance of fibers' cross-sectional shape for air filters operating in the slip flow regime[END_REF][START_REF] Kirsh | Stokes flow and deposition of aerosol nanoparticles in model filters composed of elliptic fibers[END_REF]. In recent years, several researches have been dedicated to square or rectangular fibers. The rectangular fiber is one of the non-circular fibers which have larger specific surface area compared to the circular cross-section fiber with the same volume fraction. Due to this, it is easy to perform higher collection efficiency of particles and to increase mechanical strength, porosity, particle loading capacity for fibers [START_REF] Huang | Numerical study of pressure drop and diffusional collection efficiency of several typical noncircular fibers in filtration[END_REF]. Fardi and Liu [START_REF] Fardi | Flow field and pressure drop of filters with rectangular fibers[END_REF] investigated the pressure drop with staggered array of rectangular fiber in a viscous fluid flow. Then, they studied the capture efficiency of rectangular fiber dominated by the Brownian diffusion and interception mechanisms [START_REF] Fardi | Efficiency of fibrous filters with rectangular fibers[END_REF]. Wang [START_REF] Wang | Stokes flow through an array of rectangular fibers[END_REF], Ouyang and Liu [START_REF] Ouyang | Analytical solution of flow field and pressure drop for filters with rectangular fibers[END_REF] solved the biharmonic equation of 2D Stokes flow across an array of rectangular fibers with different boundary conditions. Wang [START_REF] Wang | Stokes flow through an array of rectangular fibers[END_REF] used a rectangular boundary around each fiber, while Ouyang and Liu [START_REF] Ouyang | Analytical solution of flow field and pressure drop for filters with rectangular fibers[END_REF] used a circular boundary like Kuwabara did. Wang [START_REF] Wang | Stokes flow through an array of rectangular fibers[END_REF] solved the Stokes flow across an array of rectangular fibers and proposed an empirical formula to calculate the drag of a single square fiber. Ouyang and Liu [START_REF] Ouyang | Analytical solution of flow field and pressure drop for filters with rectangular fibers[END_REF] calculated pressure drop of filter with rectangular fiber perpendicular to the flow for various fiber aspect ratios and packing density. Zhu et al. [START_REF] Zhu | Inertial impaction-dominated fibrous filtration with rectangular or cylindrical fibers[END_REF] numerically and analytically investigated the impaction-dominated filtration process with rectangular fiber for particulate size much smaller than the one of fiber. The influence of fiber aspect ratio, filter solid fraction, particle size and Stokes number on the collection efficiency of a rectangular fiber has been determined for impaction-dominated filtration. Chen et al. [START_REF] Chen | Numerical simulation of aerosol collection in filters with staggered parallel rectangular fibres[END_REF] used the control volume method to analyze the effects of fiber aspect ratio, filter volume fraction, size of particles and Reynolds number on the collection efficiency with staggered parallel rectangular fibers. Hosseini and Tafreshi [START_REF] Hosseini | On the importance of fibers' cross-sectional shape for air filters operating in the slip flow regime[END_REF] investigated the effect of the fiber's cross-section shape on pressure drop and collection efficiency. They found that the fiber geometry is more important for pressure drop prediction rather than capture efficiency estimation. Most of these studies on rectangular fibers were conducted for a fiber Reynolds number below unity. In this condition, they put in evidence that the drag force applied on rectangular fiber is independent of the fiber Reynolds number.

In order to investigate the filtration process and assess filter performance, it is essential to acquire an accurate fluid velocity field. This can be achieved with numerical simulations. For example, the lattice Boltzmann method (LBM), as an efficient approach to simulate fluid flow, has been widely used for such studies [START_REF] Jafari | Investigation of particle dispersion and deposition in a channel with a square cylinder obstruction using the lattice boltzmann method[END_REF][START_REF] Filippova | Lattice-boltzmann simulation of gas-particle flow in filters[END_REF][START_REF] Przekop | Lattice-boltzmann approach for description of the structure of deposited particulate matter in fibrous filters[END_REF][START_REF] Lantermann | Particle monte carlo and latticeboltzmann methods for simulations of gas-particle flows[END_REF][START_REF] Wang | Simulation of filtration process for multi-fiber filter using the Lattice-Boltzmann two-phase flow model[END_REF][START_REF] Lin | An early stage of aerosol particle transport in flows past periodic arrays of clear staggered obstructions: A computational study[END_REF][START_REF] Ansari | Numerical study of solid particles motion and deposition in a filter with regular and irregular arrangement of blocks with using lattice boltzmann method[END_REF]. This method utilizes the Boltzmann equation to describe the fluid phase from a mesoscopic perspective. It is different from conventional numerical methods, which solve the incompressible Navier-Stokes equation [START_REF] Deng | Topology optimization of unsteady incompressible Navier-Stokes flows[END_REF][START_REF] Deng | Topology optimization of steady Navier-Stokes flow with body force[END_REF]. In the LBM, the fluid phase is considered as fictitious fluid particles which can collide and propagate on a lattice. The macroscopic properties (i.e., pressure, fluid density and velocity) can be obtained through the density distribution functions of fluid particles. Compared with conventional methods, LBM demonstrates significant advantages due to its parallel computation and capability to deal with complex boundary conditions [START_REF] Wang | Numerical simulation of particle capture process of fibrous filters using lattice boltzmann two-phase flow model[END_REF][START_REF] De Rosis | Lattice Boltzmann analysis of fluid-structure interaction with moving boundaries[END_REF][START_REF] Ernst | A lattice Boltzmann method for simulating transport and agglomeration of resolved particles[END_REF]. As presented in the work of Trunk et al. [START_REF] Trunk | Inertial dilute particulate fluid flow simulations with an Euler-Euler lattice Boltzmann method[END_REF], the fluid and particles detached from the walls of bifurcation are treated using LBM. Jafari et al. [START_REF] Jafari | Investigation of particle dispersion and deposition in a channel with a square cylinder obstruction using the lattice boltzmann method[END_REF] computed fluid flow past a square in a channel using the LBM for small particle (diameter less or equal to 10 µm) and Reynolds number larger than 100. In their work, particles transport and deposition are treated with a Lagrangian description, where the presence of particles does not influence the fluid velocity field. Ansari et al. [START_REF] Ansari | Numerical study of solid particles motion and deposition in a filter with regular and irregular arrangement of blocks with using lattice boltzmann method[END_REF] also employed the LBM to calculate the laminar flow field in regular and irregular structure filter to study the transport and deposition of fine particles. They also used a Lagrangian description neglecting the particle influence on fluid due to the small size of particles involved. Concerning the geometry, it can be noted that all these studies are limited to a unique square obstacle in a channel flow.

According to solid particle volume fraction φ p [START_REF] Elghobashi | Particle-laden turbulent flows: direct simulation and closure models[END_REF][START_REF] Elghobashi | On predicting particle-laden turbulent flows[END_REF], for value of φ p ≤ 10 -6 , the effect of particle on the fluid can be neglected; thus, the one-way coupling could be sufficient. In the second regime 10 -6 < φ p ≤ 10 -3 , the mutual interaction between the fluid and particle need to be considered: two-way coupling. As increased in the volume fraction of particle φ p > 10 -3 , both the fluid-particle interaction and particleparticle interaction should be included, four-way coupling is needed. In many fluid-particle simulations, the fluid and solid phases are generally dealt separately, particle is usually viewed as point without volume [START_REF] Wang | Numerical simulation of particle capture process of fibrous filters using lattice boltzmann two-phase flow model[END_REF], the effect of particles on fluid velocity field, and interparticle interactions are often ignored [START_REF] Jafari | Investigation of particle dispersion and deposition in a channel with a square cylinder obstruction using the lattice boltzmann method[END_REF][START_REF] Ansari | Numerical study of solid particles motion and deposition in a filter with regular and irregular arrangement of blocks with using lattice boltzmann method[END_REF][START_REF] Elghobashi | On predicting particle-laden turbulent flows[END_REF]. However, in this work, we consider two-way coupling leading to the action of the fluid on solid particles and reciprocally, though there is a small amount of interparticle interactions. As a first step toward high particle concentrations, a more accurate resolution of the fluid velocity field in the presence of particles is required. This can be achieved through the lattice Boltzmann method (LBM) and the discrete element method (DEM). Interparticle interactions can be simulated with the widely used DEM [START_REF] Cundall | Discrete numerical model for granular assemblies[END_REF][START_REF] Lominé | Dispersion of particles by spontaneous interparticle percolation through unconsolidated porous media[END_REF][START_REF] Lominé | Transit time during the interparticle percolation process[END_REF]. Therefore, the coupled methods, LBM for the fluid flow and DEM for interparticle interactions, provide an accurate and predictive tool for simulating fluid-particle problems, as the one treated in this paper [START_REF] Feng | Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows: Computational issues[END_REF][START_REF] Lominé | Modeling of fluid-solid interaction in granular media with coupled lattice boltzmann/discrete element methods: application to piping erosion[END_REF].

The originality of the present work is to use this coupled LB-DE methods to simulate the transport and deposition of solid particles with geometric boundary in filtration process. Firstly, the hybrid method handles a simple clean filtration problem and demonstrates its ability to reproduce and extend results from previous studies. Then, simulating fluid and particle flow in fibrous filter with rectangular fiber offers a fundamental insights into the physical phenomena in particle-fluid systems. Finally, the pressure drop, capture efficiency and quality factor of periodic arrangement of rectangular fibers with different orientation angles and aspect ratios are simulated at various Reynolds numbers.

In the following section, the numerical method is presented. More specifically, each method is presented separately; then, the coupling between them is presented in turn. In the next section, results of the simu-lation of deposition of particles on a rectangular fiber for different aspect ratios and different orientation angles are presented through the analysis of pressure drop, capture efficiency and quality factor.

2 Numerical method

Lattice Boltzmann method for fluid flow

The LB method LBM is based on gas kinetic theory and is applicable to calculate the flow field. The classical D2Q9 model proposed by D'Humières et al. [START_REF] D'humières | Lattice gas models for 3d hydrodynamics[END_REF] is adopted. In the D2Q9 model, a computational domain is discretized using a square lattice with nine discrete velocities e i (shown in Fig. 1). The fictitious fluid particles of LBM are constrained to be located at a node and move only along links with nine directions i, and collide with each other at the nodes. As presented in Fig. 1, discrete velocities in the D2Q9 model are defined as:

e i =              (0, 0) (i = 0) C cos π(i -1) 2 , sin π(i -1) 2 (i = 1, ....4) √ 2C cos π(2i -9) 4 , sin π(2i -9) 4 (i = 5, ....8) (1)
in which C means the lattice speed given by:

C = h/dt (2)
h is the lattice spacing and dt is the time step. In the LBM, the fundamental notion is particle distribution function f i representing the probability density of fluid particles with velocity e i along the ith direction at position x and time t. The macroscopic properties of the fluid such as mass density ρ and momentum density velocity ρu at (x, t) can be determined through moments of f i :

ρ = 8 i=0 f i (3) u = 1 ρ 8 i=0 f i e i (4)
The principle of LBM consists of two sequential substeps: collision and propagation. Collision substep:

f i (x, t + ) = f i (x, t) + ∆ i (5)
Propagation substep:

f i (x + e i dt, t + dt) = f i (x, t + ) (6)
The term f i (x, t + ) in Eqs. ( 5) and ( 6) is called post-collision term; t + is an intermediate time t < t + < t+dt.

For any grid node x, x + e i dt is its nearest node along the direction i. ∆ i is the collision operator. By introducing a simplified and linearized Bhatnagar-Gross-Krook (called BGK) collision operator [START_REF] Bhatnagar | A model for collision processes in gases. I. small amplitude processes in charged and neutral one-component systems[END_REF], the evolution of the density distribution functions is obtained [START_REF] Qian | Lattice bgk models for navierstokes equation[END_REF]:

f i (x + e i dt, t + dt) = f i (x, t) - 1 τ [f i (x, t) -f eq i (x, t)] (7)
where τ is a dimensionless relaxation time. f eq i is the ith equilibrium distribution functions. It is expressed as

f eq i (x, t) = κ i ρ 1 + 3 C 2 e i • u + 9 2C 4 (e i • u) 2 - 3 2C 2 u • u (8)
where κ i is the weight coefficient related to the D2Q9 model,

κ 0 = 4/9, κ 1 , 2 , 3 , 4 = 1/9, κ 5 , 6 , 7 , 8 = 1/36.
The kinematic viscosity ν and time step dt are related to the lattice velocity C, dimensionless relaxation time τ and lattice spacing h. The equations for calculating kinematic viscosity and time step are given by:

ν = 1 3 τ - 1 2 Ch (9) dt = 1 3ν τ - 1 2 h 2 (10)
One of the main assumptions of the LBM approach used in this work consists in considering a small compressibility of the fluid. Hence, its pressure is calculated by the following state law:

p = C 2 s ρ (11)
where C s is the sound speed, defined as

C s = C/ √ 3.
The value of the Mach number is limited to M ≤ 0.1, which ensures the compressibility error M 2 < 1%.

As macroscopic parameters, velocity and pressure are derived from the density distribution functions f i , they cannot be imposed directly on the system boundaries. The boundary condition method proposed by Zou and He [START_REF] Zou | Lattice bgk models for navier-stokes equation[END_REF] is applied on the system. At the interface of solid and fluid, bounce-back on the link between nodes [START_REF] Ladd | Numerical simulations of particulate suspensions via a discretized boltzmann equation. part 1. theoretical foundation[END_REF][START_REF] Behrend | Solid-fluid boundaries in particle suspension simulations via the lattice boltzmann method[END_REF] is considered to ensure the no-slip boundary condition on the solid surface.

Discrete element method for particle interactions

The discrete element method (DEM) introduced by Cundall and Strack [START_REF] Cundall | Discrete numerical model for granular assemblies[END_REF] is employed to deal with interactions between solid particles in contact. In the DEM, particles are allowed to slightly overlap among each other during their contact. When two particles a and b located at positions r a and r b contact each other, the overlap is

χ n = R a + R b -(r a -r b ) • n,
where R a and R b are particle radius and n is a unit vector normal to the contact line in 2D (or contact plane in 3D). However, seeking of interactions between particles is the most CPU time-consuming part of the method. In order to avoid this consumption, a doubly linkedlist method [START_REF] Lominé | Dispersion of particles by spontaneous interparticle percolation through unconsolidated porous media[END_REF] is implemented to optimize interaction search. We used "linear spring dashpot" (LSD) model [START_REF] Lominé | Transit time during the interparticle percolation process[END_REF][START_REF] Schwager | Coefficient of restitution and linear-dashpot model revisited[END_REF] to calculate the normal force F n and an extension method proposed by Cundall and Strack [START_REF] Cundall | Discrete numerical model for granular assemblies[END_REF] to compute the tangential force F s respectively. The expressions are given by:

F n = K n χ n + γ n χn ( 12
)
F s = F s + ∆F s ( 13
)
where K n is the normal stiffness coefficient of the interaction and γ n is a damping coefficient. χn is the time derivative of χ n . In Eq. ( 13), F s , computed with a shear force increment ∆F s , is initialized to zero when the contact is formed. The ∆F s can be calculated by accumulating the path of moving contact point between the two particles in contact as follows:

∆F s = -K s ∆χ s ( 14
)
where K s represents the interaction shear stiffness and ∆χ s is the increment of the tangential relative displacement of the particles contact. Tangential sliding occurs if

|F s | > min(λ a , λ b )|F n | (15)
where λ a and λ b are particle friction coefficients.

For each solid particle, we use the contact law described above to calculate the total contact force F c which is equal to the summation of contact force over all interacting partners. At each time step of the DE method, the new position and velocity of each solid particle are deduced from the explicit time integration of Newton's second law:

m p du p dt = F c + F h (16)
where m p and u p are, respectively, particle mass and velocity. F h is the hydrodynamic force which will be expressed in Sect. 2.3. Torques applied to the solid particle are computed in the same way by integrating the equation

T c + T h = J d p dt ( 17
)
where T c is the torque resulting from the inter-particle forces, T h the hydrodynamic torque (see Sect. 2.3) and J the moment of inertia of particles and p the angular velocity. The time integration is performed using Verlet's algorithm [START_REF] Verlet | Computer experiments on classical fluids i. thermodynamical properties of lennard-jones molecules[END_REF].

LB-DE coupling

The hydrodynamic force and torque determined from LBM are added to contact forces in the DEM. Then, position and velocity of each solid particle are calculated. Thus, a subcycling scheme is needed in the coupled LB-DE procedure. The DEM time step dt DE is generally smaller than the LBM time step dt. As a consequence, a number of DEM loops should be included in a single LBM cycle, which reduced a new value dt DE = dt/n, where n is an integer greater or equal to dt DE /dt. n substeps of DEM are performed during one computational step of LBM. Generally, the subcycling scheme works well for small n [START_REF] Feng | Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows: Computational issues[END_REF]; however, further investigations for large n are required to understand its behavior. As discussed in Ref. [START_REF] Lominé | Modeling of fluid-solid interaction in granular media with coupled lattice boltzmann/discrete element methods: application to piping erosion[END_REF], to not alter LB-DE coupling, we choose n values smaller or equal to 3. Solid particles are viewed as objects with geometric boundaries, and they are mapped on LBM lattice. During the propagation phase of the LBM, a density distribution function migrates from one lattice node to its nearest node along the direction i. The nodes belonging to an obstacle which are the closest ones to its external surface are called "Solid Boundary" (SB) nodes. All the fluid nodes located outside the obstacle, with at least one direction pointing toward a SB node, are called "Fluid Boundary" (FB) nodes. Then, between a SB node and a FB node, we introduce a boundary link. Here we use the "mid-plane bounce-back" boundary condition on the boundary link to describe the transfer of momentum between a SB and a FB node [START_REF] Ladd | Numerical simulations of particulate suspensions via a discretized boltzmann equation. part 1. theoretical foundation[END_REF]. At the end of propagation step, the fluid particles reflect along a boundary link σi (pointing from fluid to solid nodes) at mid-distance. Reflected distribution functions are given by [START_REF] Ladd | Numerical simulations of particulate suspensions via a discretized boltzmann equation. part 1. theoretical foundation[END_REF]:

f -σi (x F B , t + dt) = f σi (x F B , t + ) -2α i u s .e i f σi (x SB , t + dt) = f -σi (x SB , t + ) + 2α i u s .e i (18)
where -σi simply denotes the opposite direction of σi. x SB = x F B + e σi dt, α i = 3κ i ρ/C 2 . u s is the boundary velocity at the middle of the boundary link σ. It can be expressed as u s = u p + p ∧ r c , where u p and p are the translational and angular velocities, respectively, at the center of mass of the solid obstacle. r c = x F B + e σi dt/2 -x c is the vector joining the center of mass x c of the solid obstacle to the middle of the boundary link.

We can determine the fluid force and torque applied on an obstacle by calculating the momentum exchange over all the boundary links of the obstacle [START_REF] Lominé | Modeling of fluid-solid interaction in granular media with coupled lattice boltzmann/discrete element methods: application to piping erosion[END_REF]:

F h t + 1 2 dt = σi 2 Ω dt f σi (x, t + ) -α i u s .e i e σi (19)
T h t + 1 2 dt = σi T σ x, t + 1 2 dt ( 20
)
where Ω is the cell lattice surface in 2D (or cell lattice volume for the 3D case). To be time consistent, hydrodynamic forces and torques applied on obstacles are obtained by an average over two consecutive LBM time steps:

F h (t) = 1 2 F h t + 1 2 dt + F h t - 1 2 dt ( 21
)
T h (t) = 1 2 T h t + 1 2 dt + T h t - 1 2 dt ( 22
)
During the DEM loops, the hydrodynamic forces and torques applied on particles remain unchanged between two LBM consecutive time steps. Then, Eqs. ( 16) and ( 17) can be solved. New positions and velocities of the solid particles are obtained.

3 Numerical simulation for filtration process of filter A schematic view of a periodic cell of the filter is illustrated in Fig. 2. The fiber, represented by its rectangular cross section, is located in the center of a fluid domain of width H and length L. This domain is also the computational one. The fiber can be oriented with angle θ relatively to x-axis. w and b are, respectively, the width and height of the fiber. A constant velocity U 0 is applied at the inlet of the computational domain. At the outlet, we have imposed ∂u x (x, 0)/∂y = 0 and ∂u y (x, 0)/∂y = 0, where u x and u y are components of velocity vector U 0 , respectively, along x and y directions. Both left and right sides of computational domain are periodic boundary conditions for both fluid and particles. After establishment of steady state, particles are randomly injected into the domain on a line at the inlet (see Fig. 2). Since the particle size used here is d p = 30 µm, the Brownian diffusion is relatively weak, and the Brownian force on a particle is far less than the drag force exerted by the fluid [START_REF] Wang | Numerical simulation of particle capture process of fibrous filters using lattice boltzmann two-phase flow model[END_REF]. For this reason, the Brownian diffusion is not considered in this work. The initial velocity of each injected particle is equal to the fluid velocity at this location.

The filtration process usually includes two stages: initial and aging stages [START_REF] Duval | Two-dimensional lattice-boltzmann model of hydrosol depth filtration[END_REF]. At initial stage, the captured particles inside the filter is relatively small. The influence of deposited particles on the properties of the fiber is negligible, fiber shape does not change with time, and filter is regarded as a clean one. Thus, it is reasonable to assume that particles are removed from the flow once they are captured by the fiber at this stage of filtration [START_REF] Wang | Simulation of filtration process for multi-fiber filter using the Lattice-Boltzmann two-phase flow model[END_REF][START_REF] Lin | An early stage of aerosol particle transport in flows past periodic arrays of clear staggered obstructions: A computational study[END_REF]. However, at the aging stage, a large number of particles will deposit on the surface of the fiber, which change the shape of fiber and lead to the progressive clogging of filter and the increase of the pressure drop. In this study, we focused on the initial stage of the filtration process. Moreover, particles that are not captured by the fiber are removed when they escape from the outlet. The capture efficiency can be obtained as follows: where N inject and N escape represent, respectively, the number of solid particles injected at the inlet and removed from the outlet. The grid resolution for LBM, different values of d p /h = 8, 10, 12, 14, 16 are examined. The drag coefficient is compared with the result obtained from commercial software Comsol. It is found that good agreement can be observed for d p /h ≥ 10 and the number of nodes does not affect the drag coefficient. In this work, we fix grid resolution of LBM here d p /h ≈ 16, which is good enough to compute a reliable flow field and hydrodynamic force. Other parameters used in the simulations are presented in Table 1.

η = (N inject -N escape )/N inject (23)
In the following part, we investigate the effects of Reynolds number Re of the fiber, aspect ratio w/b and orientation angle of the fiber θ on the pressure drop and capture efficiency of the filter. The Reynolds number is defined as Re = U 0 D h /ν, where D h is the hydraulic diameter of the fiber equal to D h = 4S/P = 2wb/(w + b), S is cross-section area and P is the perimeter. Four different values of aspect ratio (w/b = 1, 2, 3, 4) are used in the simulations while keeping D h constant. Five orientation angles (θ = 0, π/8, π/4, 3π/8, π/2) are studied for each aspect ratio. The solid volume fraction of the fibers is computed as β = bw/(LH).

Pressure drop for clean fiber

The pressure drop ∆p between inlet and outlet of the periodic cell, which is related to drag force F D on a unit length of single fiber, is a key parameter to evaluate filtration performance [START_REF] Liu | Pressure drop and interception efficiency of multifiber filters[END_REF]. The determination of the dimensionless drag force F * = F D /(µU 0 ) permits to calculate the pressure drop, where µ is the dynamic viscosity, with the following relation:

∆p L = F D β bw (24)
Ouyang and Liu [START_REF] Ouyang | Analytical solution of flow field and pressure drop for filters with rectangular fibers[END_REF] determined the dimensionless drag force for a rectangular fiber with θ = 0 in Stokes flow regime. It is worth to note that the volume fraction has also an influence on the dimensionless drag force. In order to compare with these results and to validate our model, we have plotted in Fig. 3 the dependency of the dimensionless drag force, obtained with our program, on the volume fraction of fiber for different aspect ratios w/b and Re = 0.2 and more generally Re < 1 [START_REF] Ouyang | Analytical solution of flow field and pressure drop for filters with rectangular fibers[END_REF]. Figure 3 shows a good agreement and demonstrates that our model can be used to simulate the pressure drop of a filter with rectangular fiber for low Reynolds number.

In addition to the volume fraction, the aspect ratio, orientation angle of rectangular fiber as well as Reynolds number influence the pressure drop. In previous works, many researchers generally choose the Reynolds number < 1 [START_REF] Rabiee | Investigation of the characteristics of particulate flows through fibrous filters using the lattice boltzmann method[END_REF][START_REF] Huang | Numerical study of pressure drop and diffusional collection efficiency of several typical noncircular fibers in filtration[END_REF]. As Re < 1, Liu and Wang [START_REF] Liu | Pressure drop and interception efficiency of multifiber filters[END_REF] found that the dimensionless drag force is independent of Reynolds number for circular fiber. Furthermore, Fardi and Liu [START_REF] Fardi | Flow field and pressure drop of filters with rectangular fibers[END_REF] concluded also that F D is relatively independent of Re for rectangular fiber. In this study, we focus on the study of the relationship between dimensionless force and Reynolds number for the case of Re > 1. To vary the Reynolds number, we have chosen to vary the fluid velocity at the inlet.

Figure 4 shows the dimensionless drag force F * as a function of Reynolds number, varying from 2 to 10, with four typical aspect ratios and θ = 0. As presented in Fig. 4, F * is influenced by Reynolds number and increased with it in the region of 2 ≤ Re ≤ 10. When Re increases, the growth rate of F * becomes gradually enhanced for aspect ratio (w/b) from 1 to 4. Thus, Reynolds number has significant effect on the pressure drop for large value of aspect ratio. We know that a lower pressure drop across a filter is more favorable. Here at a fixed Reynolds number, it is obtained for rectangular fiber with a small aspect ratio (i.e., w/b = 1).

The streamlines around rectangular fiber with aspect ratio of w/b = 2 and different orientation angles are presented in Fig. 5. It is obviously observed that the orientation angle of fiber modifies the distribution of the flow field, which furthermore influences the pressure drop. The streamlines near the fiber change remarkably as the orientation angle decreases. Figure 8 shows the relationship between dimensionless drag force and orientation angle with different aspect ratios. Complementary to this figure, Fig. 6 shows the evolution of the windward area l c , defined as:

l c = w cos θ + b sin θ ( 25
)
and represented in Fig. 7. For l c , we vary w and b in order to keep D h fixed in the simulation. The maximum value of l c for each aspect ratio is also plotted. Figure 8 shows the dimensionless drag force varying with θ.

It illustrates that when the aspect ratio w/b = 1 (i.e., square), F * is larger for θ = π/4 than for the other angles due to the largest windward area in comparison with the other angle values.

The maximum windward area l c is arctan b/w. However, F * decreases with θ for another aspect ratios w/b = 2, 3, 4, the larger w/b is, the faster the F * decrease is. It can also be noted that F * at θ = π/2 is lower than that at θ = 0 for w/b = 2, 3, 4. As observed in Fig. 6, the reason is that the windward area of θ = π/2 is much smaller than that of θ = 0. Therefore, as θ increases, the windward area of the rectangular fiber decreases, which allows the fluid to pass through the filter more easily. It can also be noted, in Fig. 6 that the values of l c are quite similar for θ = 3π/8 and π/2, leading to quite similar values of drag force as observed in Fig. 6. Based on this analysis, it is concluded that a larger θ with a smaller windward area of rectangular fiber is more stable to obtain a lower pressure drop. Figure 9 presents the same simulations results of Fig. 8 in terms of cross section l c for all aspect ratios investigated. It illustrates that drag force increases with l c . In this analysis, we present the relationship between F * and θ (or l c ) in the case where b ≤ w. However, for the case of b > w (i.e., w/b = 1/2, 1/3, 1/4), symmetrical results about θ = π/4 with those of Fig. 8 can be obtained, which are not demonstrated in this work.

Capture efficiency

Capture efficiency, defined in Eq. ( 23), is another important parameter used to assess the filter performance. The capture efficiency obtained using LB-DE methods has been validated by comparison of empirical correlation for single circular fiber in our previous work [START_REF] Fan | Modelling particle capture efficiency with lattice boltzmann method[END_REF]. Good agreement is observed between the results of LB-DE methods and empirical correlation. Thus, the LB-DE methods developed in this work can be used to simulate the capture efficiency. In the present simulations, the total number of solid particles is 500 and 10 particles are randomly injected every 0.00016 s at the whole inlet, as illustrated in Fig. 2, once the flow field reaches a steady state. The velocity and trajectory of solid particles depend obviously on the fluid velocity field, whose regime is characterized by Reynolds number. Thus, the relationship between capture efficiency and Reynolds number is sought in Fig. 10. It is shown that the capture efficiency varies slightly with the Reynolds number from 2 to 10 for a fixed aspect ratio and orientation θ = 0. The reason is that particles are mainly dominated by the inertial impaction, which means that particles deviate from the flow streamlines at first, and after several time steps, they do not follow them any more. For Re = 10, particles move straight ahead from inlet to outlet since they have a strong inertial force. For this case, the Reynolds number does not have a significant effect on capture efficiency. In the region of 2 ≤ Re ≤ 10, the rectangular fiber with large aspect ratio (i.e., large w) is more stable to obtain a high capture efficiency due to its large windward area. The effect of orientation angle θ on the capture efficiency is analyzed. Different orientation angles corresponding to different cross-section area of fiber have significant influence on the capture efficiency. In Fig. 11, capture efficiency as a function of orientation angle is presented for Re = 8. We can easily recognize that the capture efficiency of square fiber (b/w = 1/1) at θ = π/4 is higher than for other orientation angles, since l c is maximum for this orientation angle. In addition, the value of capture efficiency of square fiber is symmetrical about θ = π/4. For the aspect ratios of w/b = 2, 3, 4, the highest capture efficiency is obtained for θ = arctan b/w. It is also observed that the capture efficiency for different aspect ratios is very close when the orientation angle is 3π/8 due to similar values of l c for this orientation angle, as observed in Fig. 6. More generally, the similarity between Figs. 6 and 11 implies that the capture efficiency is mainly governed by the cross section l c .

Figure 12 shows, for Re = 8, the capture efficiency η versus the windward area l c for all aspect ratios w/b investigated in this paper and for different values of orientation angle. Figure 12 illustrates that the larger the cross-section area is, the higher a capture efficiency can be obtained. The capture efficiency increases almost linearly with the cross-section area, which results from the strong inertial force. The particle Stokes number is St = ρ p d 2 p U 0 /(18µD h ), where ρ p is the particle density. With our data, the Stokes number is O F + i i c i n i y i η i . i 1 i 5

Fig. 12 Capture efficiency as a function of the cross-section area of the rectangular fiber for Re=8 and different orientation angles. .

Figure 13 presents the same results for different Reynolds number. The same behavior can be observed, with a small dependency on Re, also clearly indicating that the capture efficiency is mainly governed by the windward area l c in this flow regime.

Figure 14 depicts an instantaneous image of the velocity field of the fluid and particle distribution for St = 8.33 with different orientation angles θ. It provides an obvious evidence that the particle transport affects the flow field. The link between fluid streamlines presented in Fig. 5 and particle positions presented in Fig. 14 shows that the particles do not follow the streamlines downstream of the fiber. As explained previously, this is due to the inertial regime characterized by large Stokes number as in our simulation. Additionally, since the particle volume is considered in the simulation, a small amount of particles interactions are observed in Fig. 14, though the number of particles within the fluid domain is not huge.

Quality factor

A best-performing filter first of all not only ensures a high capture efficiency η but also maintains a low pressure drop ∆p. In order to judge the performance of fibrous filter, the definition of a quality factor (or the figure of merit) is an indicator of the ratio of penetration rate P (P = 1 -η) to the pressure drop. The quality factor QF is expressed as [START_REF] Wang | The influence of fiber geometry and orientation angle on filtration performance[END_REF][START_REF] Wang | Filtration of aerosol particles by elliptical fibers: A numerical study[END_REF][START_REF] Kuo | Filtration and loading characteristics of granular bed filters[END_REF]:

QF = -ln(P ) ∆p (26)
From the filter quality factor standpoint, a larger value of the quality factor QF indicates a better quality of the fibrous filter. Figure 15 depicts the relationship between quality factor and Reynolds number with different aspect ratios. It can be observed that the quality factor decreases rapidly when Re varies from 2 to 4. Then, it gradually changes as Re varies from 4 to 10. Indeed, the pressure drop increases faster than η, which leads to a decline in the quality factor. As a reference, simulation results for a circular fiber with diameter D f = D h are also presented in Fig. 15. The performance of the filter for w/b = 1 and w/b = 2 is approximately equal in the region of 2 ≤ Re ≤ 10. They have better quality than the fiber with the aspect ratio of w/b = 3 and circular fiber. They also show great advantages compared with the case of w/b = 4. Finally, we explore the effect of orientation angle and aspect ratio for Re = 8. Figure 16 shows the effect of orientation angle on the quality factor. Figure 16 illustrates that the quality factor for a rectangular fiber increases with the orientation angle at first and then decreases. We find that a higher quality factor can be obtained when θ = π/4 for aspect ratio w/b = 1. Moreover, for θ = 3π/8, the values of the quality factors are quite similar for the different aspect ratio investigated due to the fact that the dimensionless drag force and capture efficiency are not very dependent on w/b for this value of θ. The previous analysis of dependency of drag force of capture efficiency on cross section l c shows clearly that both of them increase with l c . In terms of the quality factor, there exists an optimal value of l c /H (around 0.3), as illustrated in Fig. 17, that maximizes the quality factor QF . In Fig. 18, we have plotted the QF as a function of l c /H. QF is strongly influenced by the variation of l c . However, it is good for low Reynolds number. Although an increasing Reynolds number induces performance drop, a uniform performance can be obtained.

Conclusion

In this study, the filtration processes of fibrous filters with rectangular fibers have been simulated by coupled lattice Boltzmann and discrete element methods. The effects of Reynolds number, aspect ratio and orientation angle on filtration performance parameters, pressure drop and capture efficiency, are investigated for a rectangular fiber during the initial stage of filtration. The quality factor is also used to characterize the filtration performance in order to obtain a high capture efficiency and low pressure drop. The following conclusions can be drawn based on the presented simulation results. First, the dimensionless drag force is dependent not only on the orientation angle but also on the Reynolds number when 2 ≤ Re ≤ 10. For large aspect ratio (i.e., w/b = 4), the dimensionless drag force increases rapidly with Re. Due to the variation of the windward area l c of the rectangular fiber, the dimensionless drag force decreased greatly with the orientation angle at this aspect ratio. Indeed, it has been illustrated that the drag force increase significantly with l c .

Second, the capture efficiency is dependent on both orientation angle and aspect ratio during the filtration process. A larger aspect ratio (i.e., w/b = 2, 3, 4) at θ = 0, π/8 generally leads to a higher capture efficiency. For a square fiber (w/b = 1) at θ = π/4, the capture efficiency is higher than for of other orientation angles due to its larger cross-section area. Indeed, we found that for Re > 1, the capture efficiency increases almost linearly with the windward area. However, it can be noted that the capture efficiency increases slightly with Re. Finally, in terms of quality factor, it is valuable to measure the filtration performance concerning the capture efficiency and pressure drop. For these different aspect ratios, the quality factor increases with orientation angle and reaches its maximum at a certain value and then decreases. The largest quality factor is obtained for square fiber at θ = π/4. Filter performance, through the quality factor, has also been analyzed in terms of windward area. This study shows that there exists an optimal value l c that guarantees the minimization of pressure drop and maximization of capture efficiency.

Generally, through the analysis of pressure drop, capture efficiency and quality factor, the current study provides a contribution for filter design and optimization. Future work will extend the current investigation to more complex studies such as filtration processes with particle loading on fibers implying fiber shape changes. For this, the LB-DE method validated and used in this study will be a valuable tool.

Fig. 1

 1 Fig. 1 Discrete velocity in the D2Q9 model

Fig. 2

 2 Fig. 2 Schematic view of the computational domain for simulating the flow field around a rectangular fiber

Fig. 3 Fig. 4

 34 Fig. 3 Dimensionless drag force as a function of volume fraction, LB simulation results are presented and compared with the work of Ouyang and Liu [20].

Fig. 5

 5 Fig. 5 Streamlines around rectangular fiber with an aspect ratio of w/b = 2 and orientation angles. a) θ = 0, b) θ = π/8, c) θ = π/4, d) θ = 3π/8, e) θ = π/2. Re = 8 for all configurations.

Fig. 6 7Fig. 8 Fig. 9

 689 Fig. 6 Windward area with orientation angle θ.

Fig. 10

 10 Fig. 10 Capture efficiency as a function of Reynolds number for θ = 0.

2. 08 ≤Fig. 11

 0811 Fig. 11 Capture efficiency as a function of orientation angle for Re = 8.

Fig. 13 Fig. 14

 1314 Fig. 13 Capture efficiency as a function of cross-section area of rectangular fiber for different Re and different orientation angles.

Fig. 15

 15 Fig. 15 Quality factor as a function of Re for θ = 0.

Fig. 16

 16 Fig. 16 Quality factor as a function of orientation angle for Re = 8.

Fig. 17 Fig. 18

 1718 Fig. 17 Quality factor as a function of lc for Re = 8.

Table 1

 1 Model parameters used in the simulations

	Channel length, L	0.0017	(m)
	Channel height, H	0.0009	(m)
	Fluid velocity, U 0	0.157, 0.314, 0.471, 0.628, 0.758	(m/s)
	Fluid kinetic viscosity, ν	1.57 × 10 -5	(m 2 /s)
	Fluid density, ρ	1.2	(kg/m 3 )
	Relaxation time, τ	1	(-)
	Time step, dt	3.45 × 10 -8	(s)


Independent Fausto Galetto
email: fausto.galetto@gmail.com

Proof of the Syracuse_Collatz Conjecture

Keywords: Quality Methods, Numerical Methods, Hailstone Conjecture, SPQR

We analyse the conjecture and we provide a proof, using the Reliability Integral Theory.

Introduction

A problem posed by L. Collatz in 1937, also called the 3x+1 mapping [Collatz Algorithm], states that the system of the two difference equations [START_REF] Deming | Out of the Crisis[END_REF] given the initial condition y 0 = any integer positive number, arrives after some (n is a number not known in advance) "continued" iterations at the value y n =1. The system (1) can be reduced to a non-linear difference equation, as the following one (1bis)

The numbers y k+1 of the sequence provided by the previous (Collatz) equations are sometimes named hailstone numbers. The Collatz algorithm has been tested and found to always reach 1 for all numbers up to 5.48 times 10 18 . (Oliveira e Silva 2008) It is considered a very difficult problem to be solved, in spite of its very simple definition; they say that Erdős commented that "mathematics is not yet ready for such problems".

The proposed solution

Let's consider an infinite-dimensional matrix (with infinite rows and columns); rows and columns are indexed by the natural numbers 1, 2, 3, 4, ..., n, n+1, ...; every a ij entry is 0, except ........ [START_REF] Deming | The new economics for industry, government, education[END_REF] where the indexes i and j are given by (1). In every row there is only one number 1, while all the others are 0. Therefore any row is a probability vector. This idea is important because we can consider the matrix as a transition probability matrix, with the same properties of the matrices used for modelling the Markov Chains. To take advantage of the Reliability Integral Theory (RIT) we arrange the rows and columns in the following order: 1, 2, 4, 3, 5, 6, ....., n, n+1, ............: rows and columns 3 and 4 are exchanged as the reader can see in the figure 1, where it is presented a part of the infinite matrix, with its entries. The a ij =1 entries mean that the Markov Process makes a transition from the state i to the state j, with probability 1. In the figure 1, the 3 by 3 matrix with rows and columns indexed by the numbers 1, 2, 4, is highlighted due to its importance:  when the system is in the state 1, the next transition is to state 4: 1  4  when the system is in the state 2, the next transition is to state 1: 2  1  when the system is in the state 4, the next transition is to state 2: 4  2

All this means that when the system enters one of those 3 states [START_REF] Deming | Out of the Crisis[END_REF][START_REF] Deming | The new economics for industry, government, education[END_REF][START_REF] Galetto | Numero dei guasti di un sistema e determinazione di un modello reale atto a rappresentarlo[END_REF] it never leaves out of them, the system (or the process) circulates in the set [START_REF] Deming | Out of the Crisis[END_REF][START_REF] Deming | The new economics for industry, government, education[END_REF][START_REF] Galetto | Numero dei guasti di un sistema e determinazione di un modello reale atto a rappresentarlo[END_REF] forever. It is a "periodic process". This matrix is named P and its entries are indicated p ij (to remember "probability").

state 1 2 4 3 5 6 7 8 9 10 11 12 13 14 15 16 17 18 [START_REF] Galetto | Qualità. Alcuni metodi statistici da Manager[END_REF] 20 ……………………………… 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 4 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 3 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 ……………………………… 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 ……………………………… 6 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 8 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 10 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 12 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 13 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 14 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 15 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 16 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 17 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 18 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 19 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… 20 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 ……………………………… … 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… … 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 ……………………………… … 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… … 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 ……………………………… … 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… … 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 ……………………………… … 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… … 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 ……………………………… … 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… … 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 ……………………………… … 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ……………………………… … 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 ………………………………

Figure 1. The transition matrix P [only a part is shown]

Our system has a total number of states N=: we indicate a S the set of all the states. S can be partitioned into two disjoints sets S 1 , the set 1, 2, 4 [of the periodic process] of the N 1 =3 states, and S 2 , the set 3, 5, 6, ..., n, n+1, .... of the N 2 = states, S=S 1 S 2 with S 1 S 2 =: N= N 1 + N 2 . See the figure 2.

We partition the matrix P into four submatrices P 11 , P 12 , P 21 , P 22 , related to the transitions ij (i S 1 , j S 1 ) within the set S 1 =1, 2, 4, related to the transitions ij (i S 1 , j S 2 ) from S 1 =1, 2, 4 to S 2 =3, 5, 6, ..., n, n+1, ...., related to the transitions ij (i S 2 , j S 1 ) from S 2 =3, 5, 6, ..., n, n+1, .... to S 1 =1, 2, 4 and related to the transitions ij (i S 2 , j S 2 ) within the set S 2 =3, 5, 6, ..., n, n+1, ..... You can see part of the submatrices in the formula (3); only the left upper submatrix is written in full, as one can see the set S 1 . The entries of such submatrices are related to the figure 1. The periodic process" circulating in the set S 1 =1, 2, 4 is ruled by the submatrix P 11 .

Figure 2. The two disjoint sets

The partitioned submatrices are [START_REF] Galetto | Nuovi sviluppi nel calcolo dei parametri affidabilistici dei sistemi[END_REF] and can be written simply as where P11 and P22 are square matrices, The rectangular submatrix in the upper right corner has only 0 entries. The n-th power P n of the matrix P provides the transition probabilities, between the states i  j, in n steps:

: is the conditional probability of entering the state j at the nth step, given the initial state i. We have

States partitioned and transitions

S1 S2

The submatrix V n is related to the transitions to the set S 1 in n steps. Now we consider the theory given in the Fausto Galetto book Reliability and Maintenance, Scientific Methods, Practical Approach, Volume ONE.

There, two vectors z 1 and z 2 are defined, as follows,  z 1 is the vector of the (steady state) probabilities of entering into one of the states 1, 2, 4 [S 1 ], when there is a transition S 2 =3, 5, 6, ..., n, n+1, ....  S 1 =1, 2, 4.

 z 2 is the vector of the (steady state) probabilities of entering into one of the states 3, 5, 6, ..., n, n+1, .... [S 2 ], when there is a transition from S 1 =1, 2, 4  S 2 =3, 5, 6, ..., n, n+1, ..... z 1 is by definition a three-dimensional row vector [0, 0, 1] related to the set S 1 ; see the figure 2. The system enters into the set S 1 only through the state 4. After entering S 1 the probability of being in the states 1, 2, 4 is given by a vector =[1/3, 1/3, 1/3] solution of the relationship =P 11 . z 2 is by definition an infinite dimensional row null vector [0, 0, 0, ...., 0, ...] related to the set S 2 ; see the figure 2. The system never enters into the set S 2 from S 1 .

This means that

 IF a person choose a number y 0  S 2 =3, 5, 6, ..., n, n+1, ....  and applies the rules  then the Markov system enters the state y 0  S 2 for the first time,  it makes all the needed transitions within the set S 2  until it goes into the set S 1  and it remains there, with probability vector : periodic process By making the state 1 an absorbing state the system stays in the state1, forever.

The hailstone numbers y k+1 of the sequence provided by the Collatz equations end into the number 1, whichever is the starting number y 0 .