M Bongrand

Ch Bourgeois

D Breton

R Breier

E Birdsall

V B Brudanin

Alimardon V Rakhimov
email: alimardon@jinr.ru

A S Barabash

A Basharina-Freshville

S Blot

H Burešova

J Busto

S Calvez

M Cascella

C Cerna

J P Cesar

E Chauveau

A Chopra

G Claverie

S De Capua

F Delalee

D Duchesneau

V G Egorov

G Eurin

J J Evans

L Fajt

D V Filosofov

R Flack

X Garrido

H Gomez

B Guillon

P Guzowski

R Hodák

K Holý

A Huber

C Hugon

A Jeremie

S Jullian

D V Karaivanov

M Kauer

A A Klimenko

O I Kochetov

S I Konovalov

V Kovalenko

K Lang

Y Lemière

T Le Noblet

Z Liptak

X R Liu

P Loaiza

G Lutter

J Maalmi

M Macko

F Mamedov

C Marquet

F Mauger

A Minotti

A A Mirsagatova

N A Mirzayev

I Moreau

B Morgan

J Mott

I B Nemchenok

M Nomachi

F Nova

H Ohsumi

G Oliviero

R B Pahlka

J R Pater

V Palušová

F Perrot

F Piquemal

P Povinec

P Pridal

Y A Ramachers

A Rebii

A Remoto

B Richards

J S Ricol

E Rukhadze

N I Rukhadze

R Saakyan

I I Sadikov

R Salazar

X Sarazin

J Sedgbeer

Yu A Shitov

F Šimkovic

L Simard

A Smetana

K Smolek

A A Smolnikov

S Snow

S Söldner-Rembold

B Soulé

M Špavorova

V I Umatov

J Thomas

S Torre

C Vilela

D Waters

: Lal

: Cppm

: Cenbg

V Timkin

V I Tretyak

I Štekl

F A Tashimova

G Warot

M Zampaolo

published or not. The documents may come     Development of methods for the preparation of radiopure 82 Se sources for the SuperNEMO neutrinoless double-beta decay experiment Alimardon V.

centrifugation, washing, grinding, and drying), and 82 Se foil production. The conditioning stage, during which highly dispersed elemental selenium is obtained by the reduction of purified selenious acid (H 2 SeO 3 ) with sulfur dioxide (SO 2 ) represents the crucial step in the preparation of radiopure 82 Se samples. The natural selenium (600 g) was first produced in this procedure in order to refine the method. The technique developed was then used to produce 2.5 kg of radiopure enriched selenium ( 82 Se). The produced 82 Se samples were wrapped in polyethylene (12 μm thick) and radionuclides present in the sample were analyzed with the BiPo-3 detector. The radiopurity of the plastic materials (chromatographic column material and polypropylene chemical vessels), which were used at all stages, was determined by instrumental neutron activation analysis. The radiopurity of the 82 Se foils was checked by measurements with the BiPo-3 spectrometer, which confirmed the high purity of the final product. The measured contamination level for 208 Tl was 8-54 μBq/kg, and for 214 Bi the detection limit of 600 μBq/kg has been reached.

Keywords: Selenium-82, ion exchange chromatography, purification, mass spectrometry, neutron activation analysis, SuperNEMO, LSM Modane.

Introduction

Selenium-82 is one of the most suitable isotopes for studying double-beta decay (ββ decay) and searching for its neutrinoless mode (0νββ). Interest in studies of 0νββ decay has increased significantly in recent years since the discovery of neutrino oscillations. NEMO-3 was one of the leading underground experiments for ββ decay studies in different isotopes ( 48 Ca, 82 Se, 96 Zr, 100 Mo, 116 Cd, 130 Te, 150 Nd) with a total mass of about 10 kg. The experiment used a unique tracking-calorimetric technique to measure the full topology of ββ decay and obtained a number of important ββ decay results [START_REF] Arnold | Technical design and performance of the NEMO 3 detector[END_REF][START_REF] Arnold | Vanushin[END_REF][START_REF] Arnold | [END_REF]. The SuperNEMO project extrapolates the successful NEMO-3 technique to a nextgeneration of 0νββ-experiments, planning to reach a sensitivity t 1/2 0νββ ~ 10 26 years with ~100 kg of enriched isotope ( 82 Se, 150 Nd) [START_REF] Perrot | Status of the SuperNEMO Demonstrator, 38 th International Conference on High Energy Physics (ICHEP 2016)[END_REF].

Use of ultrapure materials with respect to impurities of 40 K, 232 Th, 238 U is a key factor in the success of low-background underground experiments. Therefore, all materials used (e. g. construction materials, detector shielding, photomultipliers, cables, etc.) are carefully selected after screening for radioactive impurities. Inductively coupled plasma mass spectrometry (ICP-MS) [START_REF] Karandashev | Analysis of germanium and germanium dioxide samples by mass-spectrometry and atomic emission spectroscopy[END_REF], instrumental neutron activation analysis (INAA) [START_REF] Rakhimov | Neutron activation analysis of polyethylene from neutron shield of EDELWEISS experiment[END_REF][START_REF] Leonard | Systematic study of trace radioactive impurities in candidate construction materials for EXO-200[END_REF], and ultra-low background HPGe and BiPo-3 detectors have frequently been used [START_REF] Leonard | Systematic study of trace radioactive impurities in candidate construction materials for EXO-200[END_REF][START_REF] Budjás | Gamma-ray spectrometry of ultra-low levels of radioactivity within the material screening program for the GERDA experiment[END_REF][START_REF] Argyriades | Vorobel[END_REF][START_REF] Barabash | [END_REF] for the determination of the radiopurity of materials.

The most critical factor is the ultra-high radiopurity of the ββ source foils. The SuperNEMO spectrometer is designed in such a way that all radio-decays in source foils are registered with maximal efficiency. Unfortunately, this equally applies to the useful ββ signal and the background from radioactive impurities in foils. This background in the ROI can degrade critically the sensitivity of the experiment. In the SuperNEMO detector the ββ-sources, i. e. the enriched selenium ( 82 Se), are wrapped in 12 μm-thick polyethylene. The radiopurity of the NEMO-3 experiment's 82 Se foils was measured precisely: 208 Tl -0.39 ± 0.01 mBq/kg, 214 Bi -1.5 ± 0.04 mBq/kg (1σ statistical errors are given) [START_REF] Arnold | [END_REF]. A significant reduction (down to ~ μBq/kg level) in the radioactive contamination of the 82 Se was one of the main tasks of the SuperNEMO R&D.

The work presented here concentrates on the development of the radiochemical purification and conditioning) of enriched selenium ( 82 Se) samples produced in a "standard" state, and also on radioanalytical studies of radiopure (conditioned) 82 Se produced for the preparation of new 82 Se foils. For the "standard" state we mean samples with certain characteristics: isotopically enriched, chemically pure, radiochemically pure, in a physico-chemical state of elemental selenium, and with characteristic particle size of the powder of up to 5 μm.

This paper is a continuation of our previous work [START_REF] Rakhimov | Purification of selenium from thorium, uranium, radium, actinium and potassium impurities for the low background measurements[END_REF] concerning chromatographic separation of a small quantity (100 g) of selenium and preparation of selenium oxide (SeO 2 ) as the final product. In this paper we present the developed method for production of radiopure 82 Se in large quantities (1.5 kg and 1 kg batches from JINR Dubna and France, respectively) in the appropriate form (elemental selenium powder) required for the production of thin 82 Se foils. Therefore, special attention has been paid to a number of new aspects of the technology: -Creation of a new chromatographic system with a new larger ion-exchange column optimized for obtaining kilogram quantities of highly radiopure selenium; -Development of a technique and construction of a facility (system, equipment) for reduction of purified selenious acid to an elemental state; -Preparation of the final product by decantation, centrifugation, washing, grinding and drying.

Particular attention has been paid to materials from which chemical vessels and instruments are made. These materials should be as pure as possible with respect to radioactive impurities ( 40 K, decay products of 232 Th and 238 U) and sufficiently resistant to chemical reagents. It should be noted that it is extremely difficult to find chemical vessels from materials which meet the requirements of modern low-background experiments. The radiopurity of plastic materials of the chromatographic column (PEEK -chromatographic column material, Teflon -column filter) and polypropylene chemical vessels (glasses, bottles, centrifuge tubes) was determined by instrumental neutron activation analysis (INAA). The radiopurity of the final 82 Se wrapped in the foils was measured by the BiPo-3 spectrometer [START_REF] Barabash | [END_REF] (operating in the Canfranc underground laboratory in Spain).

The distribution coefficients of a series of lanthanides ( 88 Y, 139 Ce, 143 Pm, 167 Tm, 169 Yb, 173 Lu), 60 Co and 137 Cs (a potassium analog) on the DOWEX 50W × 8 resin, used to produce the high radiopurity selenium, were also determined.

Experiment

Chemical reagents, materials and equipment

The following reagents, materials and equipment were used: bidistilled water; nitric acid (qualification "osch 18-4", GOST 11125-84); sulfur dioxide -SO 2 (GOST 2918-79, purity -99.9 %); isopropyl alcohol (qualification "osch", STP TU COMR

Production of radiopure selenium (purification of selenium)

The complete purification procedure consists of four main stages:

-Preparation of the initial solution of selenious acid (H 2 SeO 3 ); -Cation-exchange chromatographic separation (purification); All processes were conducted in a specially prepared clean room equipped with air supply, exhaust and recirculation ventilation systems, air intake devices, air distributors with final filters, air regulating devices, sensor equipment and automation elements, etc. The level of dust particles in the room was measured to be ~100 000 (≥0.5 μm) per m 3 and zero in the laminar box (measurement by the AERO-TRAK Handheld Particle Counter TSI 9303-01) where the chromatography stage was carried out.

Preparation of initial solution of H 2 SeO 3

The selenium was divided into portions for purification. The volume of one portion of initial selenious acid solution was 5 L, with selenium concentration of 0.5 mol/L. In the case of natural selenium ( nat Se), the solution was prepared by dissolving 200 g of selenium in concentrated nitric acid, evaporating to SeO 2 , followed by dissolution in water. In the case of enriched selenium (99.88 % of 82 Se), 285 g of selenium oxide ( 82 SeO 2 ) was dissolved in water (equivalent to 200 g of selenium, taking into account mass fraction corrections). The obtained solutions were then filtered.

Purification of selenium

The basic aspects of chromatographic purification of selenium for low-background measurements are described in our previous article [START_REF] Rakhimov | Purification of selenium from thorium, uranium, radium, actinium and potassium impurities for the low background measurements[END_REF]. Figure 1 shows the purification system with a reverse flow used in this work. The main changes from the previously used system relate to the components and the overall organization of the system. The unit is specially adapted for multiple repetitions of purification cycles for separate small portions in order to purify, in total, kilogram quantities of selenium.

A chromatographic column (item 10, Figure 1) is filled with cation exchange resin DOWEX 50W × 8 (200-400 mesh). The purification system uses PTFE pumping tubes and connection hoses, and the flow of liquid is ensured by a peristaltic pump. Particular attention is paid to the selection of three-way valves, which must be sufficiently resistant to increased pressure (~5 atm) and to HNO 3 (5 mol/L). Before starting the work, the resin loaded in the column was thoroughly washed (reverse direction) with HNO 3 (2 L, 5 mol/L) and water (0.5 L).

Each purification cycle of a single portion of selenium (200 g) lasted about 3 h and included the following operations, presented in Figure 1: 1. Passing of unpurified selenious acid (5 L, 0.5 mol/L) solution through the column at a rate 40 mL/min in the forward direction (vessel 1→column→vessel 4). 2. Washing the column with HNO 3 (1 L, 0.1 mol/L) in the forward direction (vessel 2→column→vessel 4). The solutions from vessels 1 and 2 were eluted from the column and collected in a plastic vessel 4 for the purified H 2 SeO 3 (6 L). 3. Elution of impurities adsorbed on the resin in the reverse direction (vessel 6→column→vessel 3) with HNO 3 (6 L, 5 mol/L). 4. Rinsing the resin from the HNO 3 (5 mol/L) residuals in the reverse direction (vessel 5→column→vessel 3) with HNO 3 (1 L, 0.1 mol/L). These solutions from vessels 6 and 5 were eluted through the column and collected in a separate vessel 3 (7 L waste (impurities) solution).

After the last operation, the purification system is ready to process the next portion of selenium.

Reduction of selenious acid to the elemental state

In order to obtain elemental selenium the purified selenious acid (vessel 4, Figure 1) was reduced with sulfur dioxide (SO 2 ) according to the reaction

2 3 2 2 2 4 H SeO 2 SO H O Se 2 H SO . + + → + (1)
Figure 2 shows the scheme of the reduction system used. A purified solution of H 2 SeO 3 (vessel 4, Figure 1) was placed in vessel 6 (6 L; 5 L of original + 1 L of washing solution) and saturated through the surface with sulfur dioxide.

To increase the absorption of SO 2 , the solution was agitated with a magnetic stirrer. The flow of SO 2 was regulated (valve 2 and rotameter 3, Figure 2) to provide almost complete absorption of SO 2 with a small excess of gas released by bubbles in vessel 8. The average flow of SO 2 was ≈10 L/h. The reduction of one portion took 6-7 h. The saturation of the solution through the surface is more convenient than bubbling due to the good solubility of SO 2 , otherwise the tube is quickly clogged with reduced selenium. The use of gas (SO 2 ) in the recovery process is preferable for our purposes to using other reagents (liquid, solid) since it minimizes radioactive contamination of the final product.

Conditioning

The reduced elemental selenium was separated from the residual solution (Figure 2, vessel 6) by decanting and centrifugation, and then was collected in 750 mL centrifuge tubes. About 200 g of selenium (a single processed portion) was stored in each tube (~300 mL of the wet material). Even small residues of sulfuric acid formed in the process of selenium reduction create serious difficulties in the process of making the selenium source. Therefore, the selenium was further washed, first 10 times with distilled water, and then 5 times with isopropyl alcohol. For this purpose, water (or isopropyl alcohol) was added to the centrifuge tube with a single selenium portion until the tube was completely filled (750 mL mark), then gently blended with a Teflon stick, decanted and centrifuged.

During the reduction, elemental selenium was recovered in the form of tiny particles (<10 nm) which, over time, especially during centrifugation, coalesce into large crumbs and flakes of millimeter portions. They are not suitable for the production of thin 82 Se foils. Therefore, a final milling of the selenium was carried out with a homogenizer in isopropyl alcohol until a homogeneous suspension was obtained (~5 min, ~500 rpm). After grinding, the characteristic size of selenium particles was of the order of μm, which is suitable for the production of thin foils (average thickness of 300 μm). The use of surfactant isopropyl alcohol in the final stages of washing and grinding also reduces the adhesion of selenium. In addition, because of its good volatility, isopropyl alcohol facilitates the drying process.

Selenium drying was carried out for a day in a vacuum cabinet, without heating, at an average vacuum of ~5 mbar. For each selenium portion (200 g), all washing, grinding and drying processes were carried out in the same centrifuge tube.

Magnetic stirrer

During conditioning, we tried to reduce the speed and time of centrifugation (typical value ~1300 rpm for 2 min), and, if possible, to replace this procedure with settling the solution; this reduces the formation of coarsegrained selenium and reduces the grinding time. The grinding process can lead to contamination of selenium by the attrition products of the homogenizer nozzles made of stainless steel (Fe, Cr, Ni), which in turn can contain radioactive impurities.

It should be emphasized that reduction and washing leads to the removal of thorium, uranium, radium, actinium and potassium as they remain in solution.

A total of 600 g of natural selenium was purified and conditioned in order to test and refine the method. The resin in the column was then renewed and used to purify 2.5 kg of enriched selenium ( 82 Se).

Measurement of the distribution coefficients of radionuclides

The following radionuclides were used as indicators: 60 169 Yb, and 173 Lu were obtained by irradiating a tantalum target with 660 MeV protons [START_REF] Lebedev | Poluchenie radioaktivnix preparatov redkozemelnix elementov i istochnikov izlucheniya dlya yadernospectroskopicheskix issledovaniy[END_REF]. The 60 Co and 137 Cs radionuclides were of commercial origin.

To determine the distribution coefficients (K D ), the static adsorption method was used. Certain volumes of selenium acid solution, H 2 O, and 10 μL of the prepared radionuclide cocktail were added to a sample of 50 mg resin (DOWEX 50W × 8) and placed in a 1.5 mL centrifugation tube ("Eppendorf", Germany) so that the volume of each sample was 1 mL. The samples were then thoroughly mixed and left to stand for 24 h to achieve equilibrium. Then, centrifugation was carried out and the radioactivity of the solution was measured with a HPGe γ-spectrometer. The distribution coefficients were calculated using the equation:

1 g 5 0 mg D 1 mL 50 L (resin) ( resin) K , (solution) (solution) A A A A µ = = (2)
where A 50 mg (resin) is the radioactivity (Bq) of 50 mg of resin, and A 50 μL (solution) is the radioactivity (Bq) of 50 μL of the solution. To reduce uncertainties in the determination of these values, samples were weighed immediately after preparation and after fluid extraction.

Instrumental neutron activation analysis of plastic materials

Determination of thorium and uranium contamination in plastic materials of the chromatographic column (PEEK -column material, 2 Teflon samples -construction material and filter for the column) and polypropylene chemical vessels (glasses, bottles, centrifuge tubes) was carried out by INAA. Prepared samples of the plastic materials together with standard reference samples were irradiated at the IBR-2 reactor (JINR, Dubna) in a flux of epithermal neutrons of 2.2 × 10 11 cm -2 s -1 . The standard materials IAEA 433 (marine sediment) and BCR 667 (estuarine sediment) were used as reference samples. The radioactivity of the samples and standards was measured with a HPGe detector GC5519 (Canberra, USA) with relative registration efficiency of 55 %. The collection and processing of gamma-spectra were carried out using Genie-2000 software [START_REF]Genie 2000 Software Family[END_REF].

Measurements with a BiPo-3 spectrometer

The BiPo-3 detector is a low radioactivity detector dedicated to the measurement of 208 Tl and 214 Bi in thin materials. The detector, running in the Canfranc Underground Laboratory since 2013, was initially developed to measure the radiopurity of the ββ-source foils for the SuperNEMO experiment. A comprehensive description of the design and performance of the detector, as well as the results of the radioactivity measurements of the first 82 Se foils prepared for the SuperNEMO experiment, can be found in [START_REF] Argyriades | Vorobel[END_REF]. The BiPo-3 detector represents at present the most sensitive radiometrics detector for radiopurity measurements of samples wrapped in foils. The obtained detection limits were 2 μBq/kg (90 % CL) for 208 Tl and 140 μBq/kg (90 % CL) for 214 Bi after 6 months of measurement.

Source foil production

The prepared source foils from the radiopure (conditioned) natural selenium (600 g) and enriched 82 Se (2.5 kg) were made in LAPP (Annecy, France); details of the foil production can be found in [START_REF] Jeremie | The SuperNEMO ββ source production[END_REF]. The surface density of the foil is ~55 mg/cm 2 . The foils were carefully packed, transported by car from Annecy to the Canfranc underground laboratory, unpacked and installed in the BiPo-3 detector. It is necessary to emphasize once again that all materials, substances, and tools used in the foil preparation procedures have been carefully screened using low-background HPGe spectrometers and the BiPo-3 (foil) detector in order to select samples with the best radiopurity. All the details of the foil production protocol will be presented in a dedicated publication.

Results and discussion

The radiopurity of the final product (at a given radiopurity of the initial selenium) depends on the following factors: -The characteristics of the chromatographic process (distribution coefficients, parameters of the column, the flow of solutions, etc.). In such chromatographic separation systems, it is usually possible to achieve very high purification coefficients. Our experience shows that the obtained purification coefficients on thorium and uranium in the chromatographic system exceed 10 5 . -Radiopurity of the liquid chemical reagents ( bi-distilled water, nitric acid and isopropyl alcohol). -Radiopurity of plastic materials of the chromatographic column and polypropylene chemical vessels. -General organization of the procedure; for example, one needs to minimize the number of operations with higher risk of contamination of the final product.

Chemical yield of selenium

The chemical yield of the selenium recovery in the chromatographic purification procedure is more than 99 %. The main losses occur during conditioning. The residue is collected for subsequent extraction of 82 Se. The losses during chromatography are negligible -less than 0.1 %, according to the ICP-AES/MS analysis of the (waste (impurities) solution (vessel 3 in Figure 1).

Distribution coefficients of radionuclides

The distribution coefficients (K D ) of lanthanides ( 88 Y, 139 Ce, 143 Pm, 167 Tm, 169 Yb, 173 Lu) and technogenic radionuclides -60 Co and 137 Cs (a potassium analog) in the DOWEX 50W × 8 -H 2 SeO 3 system were determined experimentally (Table 1). It can be seen that the K D values of rare-earth elements with decreasing concentration of H 2 SeO 3 (from 4 to 0.5 mol/L) exceed 10 4 .

The K D values of cobalt and caesium are lower, but with a decrease in the concentration of selenious acid to a working concentration of 0.5 mol/L, the coefficients increase significantly, exceeding the critical values with respect to U. Thus, when the solution passes in the forward direction in a volume exceeding 100 times the volume of the column, the radionuclides and their analogs do not enter the purified solution.

Radiopurity of plastic materials

INAA was used to determine the radiopurity of the plastic materials used: the chromatographic column and polypropylene chemical vessels, as explained in the Section 2.4. Table 2 presents the INAA results for the plastic materials that are in direct contact with solutions: PEEK (chromatographic column material), Teflon (filter for the column and construction material), polypropylene vessels (glasses, bottles, centrifuge tubes) and polyethylene terephthalate (often used for making chemical vessels). The last two materials were not used in the setup but were measured for completeness and comparison.

The most radiopure plastic materials are usually Teflon and high-pressure polyethylene, although it should be noted that there can be exceptions to this rule. Typically, low-density polyethylene and polypropylene have approximately the same level of contamination in thorium and uranium. Polypropylene is still more mechanically resistant than polyethylene, and therefore in many cases its use is preferable. The obtained results show that the materials used are reliable for the whole of our purification procedure. Of course, it is better to use the purest material (if available), however, the contamination levels in the selected plastic materials do not exceed ~10 -8 g/g for uranium. From our point of view, such levels are not critical for vessels and materials that are used repeatedly, especially for a chromatographic column. The surface of such materials will be cleaned during multiple cycles and will be in balance with the solutions used in the purification process. The radiopurity of the materials will be reduced to the radiopurity of the water and chemicals used in the purification process. As a result, we used chemical vessels made of polypropylene (with the exception of pump hoses from teflon). According to Table 2, these were the purest products available to us from the relevant materials, and their use should not substantially affect the radiopurity of the resulting selenium samples. It should be noted that INAA is, because of good sensitivity and the ability to obtain results relatively quickly, one of the most successful methods for analysis of thorium, uranium and potassium in plastic materials.

Purity of liquid chemical reagents

Data on the purity of bi-distilled water and nitric acid can be found in our previous work [START_REF] Rakhimov | Purification of selenium from thorium, uranium, radium, actinium and potassium impurities for the low background measurements[END_REF]. The bi-distilled water does not contain thorium and uranium (detection limits were <4.6 × 10 -13 g/g and <1.7 × 10 -13 g/g, respectively). ICP-AES/MS analysis of isopropyl alcohol (1 L of isopropyl alcohol was evaporated in a Teflon vessel at a temperature of 50 °C; the residue remaining in the pot was analyzed with ICP-AES/MS) showed that the content of thorium and uranium was 6.7 × 10 -14 g/g and 1.4 × 10 -13 g/g, respectively. As can be seen, the content of Th and U is very low, therefore it should not affect the purity of the selenium samples produced.

Measurement of selenium foils with the BiPo-3 spectrometer

Unfortunately, the ICP-AES/MS analyses were not sensitive enough to determine K, Th and U levels in the initial and final (purified) selenium samples (both the natural and the enriched). Only detection limits for the K, Th and U were obtained. For thorium and uranium detection limits obtained by ICP-MS were 8 × 10 -9 g/g and 3 × 10 -10 g/g for the natural Se and for the enriched selenium ( 82 Se) samples, respectively. Further work is in progress to analyze these samples by ICP-MS and accelerator mass spectrometry. The obtained detection limit for K by ICP-AES was 10 -7 g/g. The source foils containing ~150 g of natural selenium ( nat Se) and 2.5 kg of purified (conditioned) 82 Se were measured, however, with the BiPo-3 spectrometer (Table 3).

Measurements (with interruptions) took place between December 2016 and August 2017, with a total exposure of ~130 days. The results obtained (Table 3) show lower 208 Tl ( 232 Th) and 214 Bi ( 238 U) levels for enriched selenium ( 82 Se) than for the natural ( nat Se) sample which was the first sample used for optimization of the method. This shows the suitability of the technique. The radionuclide levels in 82 Se (purified and conditioned in JINR Dubna) foils which will be used in the SuperNEMO experiment, are lower by factors of about 10 ( 208 Tl) and 3 ( 214 Bi) than in the 82 Se foil used in the NEMO-3 experiment [START_REF] Arnold | [END_REF].

Conclusions

The methods for the preparation of radiopure 82 Se sources for the SuperNEMO neutrinoless double-beta decay experiment have been developed. The chromatographic separation, reduction and conditioning have been used consistently to produce 2.5 kg (from total 6.23 kg) of highly dispersed, high radiopurity (conditioned) 82 Se. Grinding and washing were used additionally to obtain the highly dispersed powder (particle size ~1 μm) of elemental selenium as the final product, required for the production of thin source foils.

Key features of the proposed technology should be noted: -use of ion-exchange chromatography for elution of impurities with reverse method; -minimum number of chemical reagents (bidistilled water, nitric acid and isopropyl alcohol) used; -use of sulfur dioxide (SO 2 ) as a reducing agent for obtaining of elemental selenium; -high chemical yield (>99 %), unique radionuclide purity and high dispersion of final product.

The productivity of technology proposed is 300 g of 82 Se per day, which makes it possible to obtain tens and hundreds of kilograms of purified (conditioned) selenium. Particular attention was paid to the selection of pure materials, liquid chemical reagents (water, nitric acid and isopropyl alcohol) and plastic vessels which were used at all stages. One needs also to highlight the methodological studies of the various plastic materials radiopurity (ultralow contaminations of thorium and uranium) performed by INAA method, which is found to be the best for these purposes. It has been shown, it is preferable to use products from Teflon and, if necessary, polypropylene.

The radiopurity of the final 82 Se foils has been measured by BiPo-3 spectrometer. The levels of impurities were found to be 8-54 μBq/kg level for 208 Tl (the 232 Th chain), and ≤600 μBq/kg for 214 Bi (the 238 U chain). This corresponds to an improvement in the radiopurity of the selenium foils by about a factor of 10 for 208 Tl and at least by about a factor of 3 for 214 Bi in comparison with the NEMO-3 experiment.

The 82 Se source foils produced by described techniques were approved for using in SuperNEMO, an experiment aimed to search for neutrinoless double-beta decay of 82 Se. The foils have been installed in the Demonstrator (the first module of the experiment) and data-taking will start at the end of 2019.

Figure 2 :

 2 Figure 2:The Scheme for the reduction system of selenious acid (H 2 SeO 3 ) to an elemental state (Se). 1 -a cylinder with SO 2 , 2 -a reducer, 3 -a control valve, 4 -a rotameter, 5-8 -polypropylene vessels (10 L), 9 -a magnetic stirrer.

  Reduction of purified selenious acid to an elemental state (Se); -Conditioning (decantation, centrifugation, washing, grinding and drying) of the final product.

Brought to you by | Imperial College London Authenticated Download Date | 10/11/19 12:45 PM -

Table 1 :

 1 Distribution coefficients (K D ) of lanthanides, as well as cobalt and cesium on the cation -exchange resin DOWEX 50W × 8 at various concentrations (C) of selenious acid in solutions.

	C(H

2 SeO 3 ), mol/L Distribution coefficients (K D ) of radionuclides 60

		Co	88 Y	137	Cs	139 Ce
	4	400	>10 4	55	>10 4
	3	1350	>10 4	160	>10 4
	1	6700	>10 4	970	>10 4
	0.5	~10 4	>10 4	1400	>10 4
		143	Pm	167 Tm	169	Yb	173 Lu
	4	>10 4	>10 4	~10 4	~10 4
	3	>10 4	>10 4	>10 4	>10 4
	1	>10 4	>10 4	>10 4	>10 4
	0.5	>10 4	>10 4	>10 4	>10 4
		Brought to you by | Imperial College London
					Authenticated
			Download Date | 10/11/19 12:45 PM

Table 2 :

 2 Radiopurity of various plastic samples measured by INAA.

	Samples	Mass, g	Irradiation/measurement	Th content, g/g	U content, g/g
			time, h
	PEEK	0.49	93/0.5-1.5	(3.8 ± 0.7) × 10 -8	(4.3 ± 0.2) × 10 -8
	Teflon (filter for column)	0.69	93/0.5-1.5	(4.9 ± 0.8) × 10 -8	(4.8 ± 0.2) × 10 -8
	Teflon -material for construction	1.10	71/0.5-7	<3 × 10 -10	<9 × 10 -11
	Glass, PP	0.37	93/0.5-1.5	(5.5 ± 1.0) × 10 -8	(1.3 ± 0.1) × 10 -7
	Bottle 5 L, PP	0.17	71/0.5-10	<2.9 × 10 -9	<1.4 × 10 -9
	Bottle 10 L, PP	0.19	71/0.5-10	<2.6 × 10 -9	<1.1 × 10 -9
	Centrifuge tube, PP	2.76	71/0.5-10	(6.9 ± 0.4) × 10 -9	(2.7 ± 0.2) × 10 -9
	Polyethylene terephthalate	0.18	71/0.5-10	<2.4 × 10 -8	<2.6 × 10 -8
	PEEK, Polyether ether ketone; PP, polypropylene.

Table 3 :

 3 Results of208 Tl and 214 Bi in the purified82 Se samples wrapped in polyethylene foils measured by the BiPo-3 detector (90 % CL).

	Origin of selenium foils	The specific activity, μBq/kg
		228 Th ( 232 Th)/ 208 Tl	226 Ra/ 214 Bi
	High radiopurity	19-336/7-121	<870/<870
	(conditioned) nat Se
	High radiopurity	22-150/8-54	<600/<600
	(conditioned) 82 Se "Dubna"
	High radiopurity	<294/<106	<1370/<1370
	(conditioned) 82 Se "France"
	Brought to you by | Imperial College London
		Authenticated
	Download Date | 10/11/19 12:45 PM

Acknowledgments: We acknowledge support by the staff of the Modane Underground Laboratory and by the MEYS of the Czech Republic (grant number EF16_013/0001733), APVV in Slovakia, CNRS/IN2P3 and Labex Enigmass (A. Remoto) in France, RFBR in Russia (Project No. 16-52-16018 NCNIL a and No. 19-52-16002 NCNIL a), STFC in the UK and NSF in the USA.


Arnesh Palanisamy

Marco Ramaioli

Paul Menut

Artemio Plana-Fattori

Denis Flick
email: denis.flick@agroparistech.fr

Multiscale modelling of flow, heat transfer and swelling during thermo-mechanical treatment of starch suspensions

Keywords: Starch Granule, Swelling kinetics, DEM, CFD, Modelling

come

Introduction

Starch is one of the most abundantly available naturally occurring biopolymers on earth.

It is predominantly found in grains and tubers and is the source of 70-80 % of the calories consumed by humans worldwide in various forms of food products [START_REF] Bertolini | Trends in starch applications. Starches: Characterization, Properties, and Applications[END_REF]. Starch also finds many applications apart from food in industries such as bio-medical, packaging, paper, textile and pharmaceutical industries. Starch in its native form occurs as granules which consists of two polysaccharides, amylose and amylopectin. Both amylose and amylopectin have α -glucose as the monomer but differ in the way they are bonded. Amylose has α(1-4) bonds and is linear whereas amylopectin has α(1-6) bonds and is highly branched.

The granule structure, molecular weight distribution, degree of amylopectin branching and the ratio of amylose to amylopectin varies according to the botanical origin of the starch granule. This, in turn, can have drastic consequences in the physical and chemical properties of the granules such as rheology and gelatinization behaviour [START_REF] Banks | Studies on starches of high amylose content. part 17. a review of current concepts[END_REF][START_REF] Hoover | Composition, molecular structure, and physicochemical properties of tuber and root starches: a review[END_REF][START_REF] Srichuwong | Starches from different botanical sources ii: Contribution of starch structure to swelling and pasting properties[END_REF][START_REF] Visser | Some physicochemical properties of amylose-free potato starch[END_REF]. Some kinds of starch are chemically modified to have enhanced structural and functional properties of the starch granules [START_REF] Debnath | Chemically modified starches and their applications in pharmacy[END_REF][START_REF] Alcázar-Alay | Physicochemical properties, modifications and applications of starches from different botanical sources[END_REF]. By 2025, the world starch production is projected to be around 156 million metric tons, including native and modified starches (Reportlinker.com, 2020). As a result starch has been in the past, a subject of perennial research interest for researchers from diverse backgrounds and will be in the near future as well.

Gelatinization is the process in which starch granules undergo an irreversible transformation when subjected to thermomechanical treatment in the presence of water. During this process, starch granules swell by imbibing water and the overall viscosity of the suspension increases. It is also sometimes referred to as starch pasting in literature. Starch granules are semi-crystalline with concentric rings of crystalline and amorphous phases [START_REF] Ratnayake | Starch gelatinization[END_REF]. On the microscopic scale, the process corresponds to the melting of crystallites, loss of birefringence, leaching of amylose into the continuous phase and sometimes rupture of the granular structure [START_REF] Atwell | The terminology and methodology associated with basic starch phenomena[END_REF][START_REF] Biliarderis | Structure and phase transitions of starch in food system[END_REF]. The evolution of viscos-ity during gelatinization is also influenced by the process parameters such as concentration, shear rate, thermal history, pH and presence of other ingredients [START_REF] Thomas | Eagan press handbook series: Starches. st. paul[END_REF][START_REF] Malumba | The effect of heating rates on functional properties of wheat and potato starch-water systems[END_REF]. For high amylopectin cereal starches, such as waxy maize starch, the granules tend to hydrate with ease, swell rapidly and rupture to a great extent which results in loss of paste viscosity. As starch is used as a viscosity enhancer in many food applications such as in ketchup and also in non-food applications such as enhanced oil recovery in oil drilling, the rupture of the granules which leads to loss in suspension viscosity should be avoided, this can be avoided by chemically modifying the starch.

Coupling starch swelling, flow and heat transfer

It is of particular interest to understand the coupled evolution of temperature, starch gelatinization (notably swelling) and rheology in heat exchangers. Pioneer work was done in this field by [START_REF] Yang | Numerical study of parameters affecting broken heating curve[END_REF] and more recently by [START_REF] Plana-Fattori | An eulerian-lagrangian approach for coupling fluid flow, heat transfer and liquid food product transformation[END_REF] and [START_REF] Plana-Fattori | A deterministic approach for predicting the transformation of starch suspensions in tubular heat exchangers[END_REF]. In these studies, the rheology of partially gelatinized suspension is related to the temperature evolution directly [START_REF] Yang | Numerical study of parameters affecting broken heating curve[END_REF] or via a swelling kinetic model where the parameters of this model are estimated empirically. Then fluid flow and heat transfer are studied after solving the governing equations in interaction with the rheology/swelling model. The simulations were in relatively good agreement with some experimental measurement. However, these studies had several limitations: the size distribution of the granules was ignored (only average diameter evolution was considered), the variability from a granule to another in terms of swelling was ignored, the granules followed exactly the fluid flow, inside a fluid parcel (on a given trajectory) all the granules were assumed to have the same behaviour. It is known that particles migrate towards a region of lower shear rate, lower volume fraction and lower viscosity [START_REF] Shauly | Shear-induced particle resuspension in settling polydisperse concentrated suspension[END_REF]. It is also known, that even without a gradient of these variables, there is shear-induced self-diffusion of the particles due to elastic collisions and lubrication interaction [START_REF] Foss | Self-diffusion in sheared suspensions by dynamic simulation[END_REF]. This is well documented in literature [START_REF] Eckstein | Self-diffusion of particles in shear flow of a suspension[END_REF][START_REF] Leighton | The shear-induced migration of particles in concentrated suspensions[END_REF][START_REF] Breedveld | The measurement of the shear-induced particle and fluid tracer diffusivities in concentrated suspensions by a novel method[END_REF][START_REF] Phillips | A constitutive equation for concentrated suspensions that accounts for shear-induced particle migration[END_REF]. Shear induced self diffusivity is proportional to the shear rate and to the squared particle radius. The proportionality parameter de-pends on the particle volume fraction (φ) and different correlations were proposed (linear or quadratic). [START_REF] Eckstein | Self-diffusion of particles in shear flow of a suspension[END_REF] found that this parameters increases almost linearly from 0 to 0.02 when the volume fraction rises from 0 to 0.2. This general trend was confirmed by [START_REF] Breedveld | The measurement of the shear-induced particle and fluid tracer diffusivities in concentrated suspensions by a novel method[END_REF] who compared different studies. Inertial effect can also be involved depending on the particle Reynolds number. Also, particle migratory self-diffusion can induce an increase of momentum and heat transfer, in other words, this phenomenon can increase the apparent viscosity and conductivity.

The main question we would like to answer in this study is, are these various granular phenomena of appreciable importance in the case of starch suspension subjected to fluid flow and heat transfer? To answer this question, an approach combining Computational Fluid Dynamics (CFD) and the Discrete Element Method (DEM) is here proposed, similar to one implemented in [START_REF] Qian | Cfd-dem simulation of the filtration performance for fibrous media based on the mimic structure[END_REF] for simulating granular filters. One specificity is that the suspension evolves because of the starch granule swelling. We apply this modelling approach to the case of a small gap Couette rheometer during fast heating of the external cylinder.

In this case, we also have shear and heat transfer with momentum and energy exchanges between fluid and particles, possible particle migration and self-diffusion (random walk) but the fluid flow is simpler than in a heat exchanger.

Starch swelling modelling

Many models exist in the literature to predict the average degree of swelling depending on the thermal history of a given starch suspension [START_REF] Lagarrigue | Swelling kinetics of waxy maize and maize starches at high temperatures and heating rates[END_REF][START_REF] Malumba | The swelling behaviour of wheat starch granules during isothermal and non-isothermal treatments[END_REF][START_REF] Plana-Fattori | A deterministic approach for predicting the transformation of starch suspensions in tubular heat exchangers[END_REF]. More recently, Desam et al. (2018aDesam et al. ( ,b, 2020) ) proposed a mechanistic model based on the Flory-Rehner theories, for the starch swelling of various grades. However, this deals at the population level assuming same model parameters for all granules and hence, does not take into account the inter-granule variability. The intergranule variability has biological origins and the alternating crystalline and amorphous phase of the starch granule is bound to slightly vary even among starch granules from the same plant. Thus, modelling the average behavior may not be sufficient to understand some mechanisms such as granule migration during thermal processing in presence of high velocity and temperature gradients. Therefore, in our previous work we have described the granules swelling kinetics at the individual granule scale using a semi-mechanistic model [START_REF] Palanisamy | Kinetic modelling of individual starch granules swelling[END_REF] for modified waxy maize starch. According to the microscopic observation of starch granules during heating in the presence of water, we observed variability in swelling kinetics from granule to granule. Each granule has a specific delay time associated with the onset of swelling. Also, the granules which began to swell at a later instant tended to swell slower compared to the granules which began to swell earlier for a given time-temperature profile. Below a minimal temperature, no swelling was observed. Experiments showed that high temperatures resulted in lower delay time and faster increase in granule diameter. These phenomena were cast mathematically in a model [START_REF] Palanisamy | Kinetic modelling of individual starch granules swelling[END_REF]. [START_REF] Nayouf | Effect of thermomechanical treatment on the rheological properties of crosslinked waxy corn starch[END_REF] indicated that the volume fraction of granule can be estimated if the mean granule size can be accurately determined. As discussed in the next section, volume fraction of starch is a critical parameter used in the closure relations for predicting viscosity in mixture models. This inter granular level kinetic model of starch granule swelling is coupled here, with a CFD-DEM model for prediction of fluid flow, heat transfer and granule trajectories.

Starch suspension rheology

As discussed in the previous sections, starch has many different varieties and is used as a viscosity enhancer in various fields. Rheology measurements are generally carried out at a constant temperature and after complete gelatinization. However, having the evolution of rheology during gelatinization is quintessential for robust process design and can be of interest to design control strategies as well. This is done notably with the Brabender farinograph method to characterize different starches before processing. This method involves measuring the torque required to move the impeller at a predetermined rpm while measuring the temperature of the sample independently throughout the gelatinization process. However, this method does not allow measuring directly the usual rheological material functions but gives only an indication of the viscosity in arbitrary Brabender units [START_REF] Eliasson | Starch in food: Structure, function and applications[END_REF]. To carry out Couette rheometry of starch suspensions during heating, at a rate comparable to the industrial processes (around 1 • C/s), [START_REF] Lagarrigue | Modélisation de l'évolution du comportement à l'écoulement de suspensions d'amidon soumises à des cinétiques rapides de température et de cisaillment[END_REF] developed a high heating/cooling rate rheometer. She took into account (via a thermal model) the thermal gradient in the rheometer gap and was able to propose a model predicting the rheology of the starch suspension as a function of its thermal history.

The properties of starch vary based on its botanical origin and processing it has undergone. Hence, suspension rheology also varies widely, it can be either shear thinning or shear thickening depending on the variety and concentrations [START_REF] Plana-Fattori | A deterministic approach for predicting the transformation of starch suspensions in tubular heat exchangers[END_REF][START_REF] Crawford | Shear thickening of corn starch suspensions: Does concentration matter[END_REF]. Origins of such varied behaviour of starch suspensions can be deciphered by leveraging the literature existing on viscosity predictions of broad class of materials known as non-brownian particulate suspensions.

Researchers over the years have developed two methods of modelling the suspension rheology namely, two-fluid model and mixture model [START_REF] Vollebregt | Suspension flow modelling in particle migration and microfiltration[END_REF]. In a twofluid model the conservation equations are written for both particle phase and the fluid phase with the interaction terms between the phases. In the mixture approach, the suspension is modelled as an equivalent single phase. Both two-fluid model and mixture model regard the liquid and solid phases as interpenetrating continua.

The two-fluid model is obtained by averaging the equations of motion over the liquid and the particle phases. This results in a set of equations with a closure problem as it contains a number of indeterminate terms which are not related to the averaged variables but associated with the details in the microscopic length scale. Multiple methods exist for coming up with this closure problem. Two of which are effective medium theory and fully resolved CFD-DEM simulation. The latter is an Eulerian-Lagrangian method which is computationally expensive. Effective medium theory approach is favourable for thermodynamic and computational point of view. The indeterminate terms of the two-fluid model are represented by the effective stress tensors of the phases and interaction force exchanged by the phases. The expressions for the effective stress tensors and fluid-particle interaction forces are important as they indicate the origins of these terms.

According to [START_REF] Jamshidi | On the closure problem of the effective stress in the eulerian-eulerian and mixture modeling approaches for the simulation of liquid-particle suspensions[END_REF], the liquid effective stress tensor is the sum of four terms. The first is the product of the liquid viscosity(µ), the liquid volume fraction (1-φ) and the mean liquid shear rate ( γ). The fourth term is similar to Reynolds stress and becomes important in turbulent case. The second and third terms together are called the 'particleparticle stress'. The 'particle-presence stress' takes notably into account that the presence of particles distort the fluid streamlines and generate larger local velocity gradients, which increase viscous dissipation. This occurs even if the particles follow the liquid i.e. with zero slip velocity with the surrounding fluid. For dense mixtures, it considers also lubrication films between a significant number of particles in which viscous dissipation occurs.

For non-Brownian, neutrally buoyant, spherical particles dispersed in a Newtonian fluid in laminar flow and low volume fractions (typically φ < 0.05), [START_REF] Jamshidi | On the closure problem of the effective stress in the eulerian-eulerian and mixture modeling approaches for the simulation of liquid-particle suspensions[END_REF] recovered the expression obtained by [START_REF] Einstein | Eine neue bestimmung der moleküldimensionen[END_REF] for the part of the 'particle-presence stress' due to streamline distortion (in this derivation the lubrication interactions are neglected). [START_REF] Jamshidi | On the closure problem of the effective stress in the eulerian-eulerian and mixture modeling approaches for the simulation of liquid-particle suspensions[END_REF] further claim that the approach often used in two-fluid modelling for the fluid/fluid interaction, which considers only the first term of the liquid effective stress tensor is incorrect and should be avoided notably because it predicts a decrease in effective viscosity as the solid concentration rises.

For dense suspensions and soft particle pastes, besides fluid/fluid interactions, a critical point is how to treat the particle/fluid and particle/particle interactions [START_REF] Vollebregt | Suspension flow modelling in particle migration and microfiltration[END_REF][START_REF] Van Der .Sman | Effective temperature for sheared suspensions: A route towards closures for migration in bidisperse suspension[END_REF]. In particular, shear induced diffusion is a phenomenon of main importance. It has been shown that particle-particle interaction leads to osmotic pressure of the suspension as derived by [START_REF] Seth | A micromechanical model to predict the flow of soft particle glasses[END_REF]. These authors studied through simulations the variation of osmotic pressure and suspension shear modulus as a function of the particle properties. It was shown that gradient in osmotic pressure of a suspension contributes to the shear induced particle migration [START_REF] Vollebregt | Suspension flow modelling in particle migration and microfiltration[END_REF]. Further [START_REF] Van Der .Sman | Effective temperature for sheared suspensions: A route towards closures for migration in bidisperse suspension[END_REF] went on to define an effective granular temperature based on osmotic pressure for sheared suspensions.

In the mixture approach, the suspension is considered like a single phase with an apparent viscosity also called suspension viscosity. The latter includes not only the effect of streamlines distortion due to the simple presence of the particles but also the particle/particle and fluid/particle interactions. Analytical and empirical correlations were proposed (e.g. [START_REF] Batchelor | The determination of the bulk stress in a suspension of spherical particles to order c2[END_REF]). The four most important factors that dictate a suspension viscosity are (a) suspending medium's viscosity, (b) granule volume fraction, (c) granule deformability, and (d) imposed shear rate [START_REF] Chen | A dynamic model of crosslinked corn starch granules swelling during thermal processing[END_REF]. Granules deformability (granule's Young's modulus) is of particular importance for large volume fractions when particle-particle interactions becomes predominant [START_REF] Li | Methodology to predict the time-dependent storage modulus of starch suspensions during heating[END_REF][START_REF] Evans | Concentration dependence of the linear elastic behaviour of model microgel dispersions[END_REF].

The governing momentum equation for the mixture can also be obtained by adding the two momentum equations of the two-fluid approach, this allows also developing approximate expressions for the suspension viscosity [START_REF] Jamshidi | On the closure problem of the effective stress in the eulerian-eulerian and mixture modeling approaches for the simulation of liquid-particle suspensions[END_REF].

For concentrated suspensions, when fluid/fluid interactions becomes negligible compared to contact and lubrication forces, DEM simulations can also predict the rheology including shear-thickening [START_REF] Ness | Shear thickening regimes of dense non-brownian suspensions[END_REF].

Despite these efforts to understand how rheology evolution is connected to the transformation of the starch suspension during heating, there is no unequivocal single answer to explain or model all the experimentally observed phenomena.

Modelling

Starch swelling modelling

The previously developed individual starch granule swelling model [START_REF] Palanisamy | Kinetic modelling of individual starch granules swelling[END_REF] is presented briefly here. In this study we use chemically stabilized cross-linked waxy maize starch (C*Tex 06205), provided by Cargil (Baupte, France). This has a D 10 greater than 8 µm and thus is non-colloidal in nature.

For a constant reference temperature T ref , the onset delay time t onset follows an exponential distribution for which the cumulative distribution function is: 1 -e -tonset/τ ref , where

τ ref is the characteristic latency time at T ref .
To traduce this, a parameter α called 'difficulty of swelling' was introduced. For each granule, α is selected at random from uniform distribution in the range [0;1]. This way, for constant temperature T ref , the delay time t onset follows:

t onset = -τ ref • ln(1 -α) (1)
Regarding the swelling onset, equation 1 arises from an experiment carried out a constant temperature of 62 • C during one hour [START_REF] Palanisamy | Kinetic modelling of individual starch granules swelling[END_REF].

For variable temperature, the time delay in onset of swelling is given by solving equations 2 and 3, which means that the delay in onset of swelling is higher at lower temperatures that at higher temperatures.

ln(1 -α) = - tonset 0 dt τ (2) τ = max 0, T ref -T min T -T min τ ref (3)
Dimensionless diameter S is given by

S(t) = D(t) -D i D f -D i (4)
where D(t) is the diameter at an instant t, D i is the initial diameter of the granule, D f is the final diameter of the granule.

After the delay time t onset is overcome, the swelling kinetics is given by:

dS dt = K(1 -S) (5)
with

K = K ref (1 -α) 0.5 T -T min T ref -T min 2 (6)
where K ref is the rate constant at T ref for a granule which has no difficulty to swell (α = 0).

For chemically stabilized cross-linked waxy maize starch (C*Tex 06205) provided by Cargill (Baupte, France) the following parameters were estimated:

T min = 60 • C, τ ref = 190 s and K ref = 0.0023s -1 at T ref = 62 • C.
From the above equations, one can describe the swelling kinetics of these starch granules.

The following assumptions regarding the starch population were made. Initial diameters and swelling ratio are from a random normal distribution: D i ⊂ N (14.77µm, 4.44µm),

D f /D i ⊂ N (2.34, 0.33).
The details and experimental basis of this swelling model for individual starch granules of cross-linked waxy maize starch were presented in previous work [START_REF] Deslandes | Estimation of individual starch granule swelling under hydro-thermal treatment[END_REF][START_REF] Palanisamy | Kinetic modelling of individual starch granules swelling[END_REF]. The model was validated by microscopic observation for different temperature ramps (up to 90 • C) and also by comparing size distribution obtained by Malvern mastersizer after different non-linear heat treatments. These model parameters should be re-calibrated for other waxy cross-linked starches and the model should be modified for other types of starches and is not applicable if the granule ruptures.

CFD-DEM Modelling

A two-phase modelling approach is used as by [START_REF] Jamshidi | On the closure problem of the effective stress in the eulerian-eulerian and mixture modeling approaches for the simulation of liquid-particle suspensions[END_REF] except that the solid phase (starch granules) is not treated as a continuous phase but as a discrete one: each particle is characterized by its centre position and velocity. This approach is similar to the one used by [START_REF] Qian | Cfd-dem simulation of the filtration performance for fibrous media based on the mimic structure[END_REF].

An Eulerian approach is used for solving flow and heat transfer of the continuous fluid (in this case pure water). In each representative elementary volume (REV) there are a number of solid particles (in this case starch granules).

For the fluid phase, in a general case, the momentum equation holds:

(

1 -φ)ρ f ∂ - → v f ∂t + (1 -φ)ρ f - → v f . - → ∇ - → v f = - - → ∇. S + (1 -φ)ρ f - → g + --→ f p→f (7)
where --→ f p→f is the sum of the forces exerted on the fluid by the particles per unit volume (notably drag forces), and S is the fluid effective stress tensor.

For non-Brownian, neutrally bouyant, spherical particles dispersed in a Newtonian fluid in laminar flow and low volume (typically φ < 0.05), [START_REF] Jamshidi | On the closure problem of the effective stress in the eulerian-eulerian and mixture modeling approaches for the simulation of liquid-particle suspensions[END_REF] obtained the following approximation.

S = p Ī -(1 + 2.5φ)µ γ (8)
where γ is twice the deformation rate tensor. we can use the deformation rate tensor of the fluid because particle velocity rapidly tends to the surrounding fluid velocity (gradient of the mean particle velocity is close to the fluid velocity gradient).

In the gap of a Couette rheometer only the variations in radial directions are considered; taking a small gap, i.e. (R e -R i ) << R e , a Cartesian approximation can be used. The radial position in the gap is referred as x (unit vector -→ u x ). On the internal wall x = 0, on the external wall x = e. The flow direction is referred as y.

For the fluid, only momentum balance in y direction is considered and only the shear stress τ xy is involved. Fluid velocity is assumed to be in y direction (unit vector -→ u y ) its magnitude is referred as v f ; hence:

(1 -φ)ρ f ∂v f ∂t + ∂τ xy ∂x = f y,p→f (9)
where f y,p→f is the sum of the forces exerted on the fluid by the particles located in the REV per unit volume in y direction. Since there is no variation in pressure nor in fluid velocity in y direction, the convective term of momentum and the pressure gradient term vanish.

In this equation τ xy accounts for the fluid/fluid interaction. Therefore, this term should take into account the increased viscous dissipation due the distortion of the streamlines caused by the simple presence of the particles. According to equation 8Jamshidi et al.

(2019), for dilute suspensions, it can be approximated by Einstein's correlation. Unfortunately, as per our knowledge, there is no available expression for higher volume fraction, which only integrates the additional stress due to streamlines distortion. As a case in point, the correlation of [START_REF] Batchelor | The determination of the bulk stress in a suspension of spherical particles to order c2[END_REF] applies for higher volume fractions but includes lubrication interactions as well which are here directly calculated in the DEM part of the model. Thus, we intend to use [START_REF] Einstein | Eine neue bestimmung der moleküldimensionen[END_REF] expression for the effective viscosity even though it certainly underestimate the additional stress due to streamline distortion for intermediate solid volume fractions of 0.05 < φ < 0.3. At least, this relation (equation 11) takes into account that higher the particle volume fraction, the more the streamline are distorted, which increases viscous dissipation.

τ xy = -µ e (T f , φ) ∂v f ∂x (10) µ e (T, φ) = µ(T f ) • (1 + 2.5φ) (11)
The energy balance for the fluid is :

(

1 -φ)ρ f C pf ∂T f ∂t + ∂q ∂x = Q p→f + Q diss ( 12
)
where

q = -λ ∂T f ∂x (13)
where Q p→f is the sum of the heat fluxes from the particles to the fluid and Q diss is the viscous dissipation of mechanical energy into heat due to the viscous forces.

Q diss = µ e (T f , φ) ∂v f ∂x 2 (14)
The external wall has zero velocity and imposed temperature evolution, say

x = e : v f = 0 and T f = T e(t) .

The internal wall has a given velocity, this corresponds to a strain controlled rheometer (constant global shear rate):

x = 0 : v f = v i = γ.e
The internal wall is also characterized by a given thermal inertia per unit surface (

m iw C p.iw A iw
):

x = 0 :

m iw C p.iw A iw dT iw dt = λ ∂T f ∂x
Periodic boundary conditions apply in y and z direction.

The trajectory of each particle of index i is given by :

m i d 2-→ x i dt 2 = --→ F f →i + j --→ F j→i (15) with --→ F f →i = 6πµr i (v f . - → u y -- → v i ) (Stokes viscous drag force).
Here, v f is the fluid velocity interpolated at particle centre position -→ x i (piecewise linear velocity profile is considered) and r i is the radius of granule i.

The maximal observed particle/fluid relative velocity was about 1 mm/s in our case, the maximal granule diameter is around 50 µm, and the lowest liquid viscosity is 0.32 mPa.s.

Therefore, Reynolds number was always lower than 0.15, inertial effect is negligible and Stokes law is applicable. Therefore we can write that

F p→f = - i∈REV --→ F f →i V REV . - → u y (16)
Lubrication forces in the normal direction are expressed by the following expression. Only the first-order normal force is considered, a more detailed lubrication interaction [START_REF] Kim | Microhydrodynamics: principles and selected applications[END_REF] will be integrated in further work.

--→

F j→i = 6πµr i β 2 (1 + β) 3 r i + r j h ij.ef f (( - → v i -- → v j ). -→ n ij ) -→ n ij (17) -→ x ij = - → x j -- → x i d ij = | -→ x ij | -→ n ij = -→ x ij /d ij h ij = d ij -(r i + r j )
-→ n ij is the centre to centre unit vector and h ij is the surface to surface distance.

h ij.ef f = max(h ij , 0.005(r i + r j )) β = r j /r i
No gravity, no Brownian or repulsive forces were considered.

If two particles come in contact (h ij < 0) an instantaneous elastic rebound is assumed.

Friction between particles are assumed negligible, this assumption should be reconsidered, especially for high volume fraction, since starch granules cannot be considered as smooth spheres. They exchange momentum -→ P in the normal direction -→ n ij given by:

--→ P j→i = 2m i m j m i + m j (( - → v j -- → v i ). -→ n ij ) -→ n ij ( 18
)
m i δ - → v i = --→ P j→i ; m j δ - → v j = --→ P i→j ( 19
)
where δ -→ v i and δ -→ v j are the difference in velocities before and after the collision for particles i and j which are applied during next time step.

Elastic rebound was assumed on the internal and external walls (x = 0 and x = e).

Periodic boundary conditions were assumed in y and z direction for both the particle and fluid phases. [START_REF] Seth | Elastic properties of soft particle pastes[END_REF] have used a similar approach of elastohydrodynamic lubrication for microgels. However, they have used Hertzian contact model when particles touch each other in place of the instantaneous shocks in the model presented here. In their case with a high concentration suspension, a detailed contact model was essential because the volume fraction was higher than the random close packing fraction. In our case, the maximum volume fraction is about 30 %, much lower than the random close packing value of 68 % which was determined for our poly-disperse suspension using the methodology given in [START_REF] Farr | Close packing density of polydisperse hard spheres[END_REF]. Thus, the current approach (instantaneous elastic collision) seems a valid approximation.

The temperature evolution of each particle is given by:

m i C pp dT i dt = Q f →i + Q ld (20) with Q f →i = 4πr 2 i h(T f -T i ), h = N u λ 2r i and N u ≈ 2.
Here, Q ld is the dissipation term for the lubrication forces and drag. The same is calculated using the following expression:

Q ld = 1 2 j [ --→ F j→i • ( - → v i -- → v j )] + --→ F f →i • (v f - → u y -- → v i ) (21)
Initially the granules are randomly distributed, the fluid and particle velocity varies linearly in the gap and the temperature is uniform and equal to the minimal swelling temperature (60 • C)

Elsewhere researchers have often used fabric tensor to quantify anisotropy [START_REF] Ness | Shear thickening regimes of dense non-brownian suspensions[END_REF][START_REF] Barreto | Quantifying the evolution of soil fabric under different stress paths[END_REF]. However, fabric tensor is typically used for the cases of quasi-static systems and we have a dynamic system with almost no true contacts. This is why we consider 'quasi-contact' when the centre-centre distance between particle i and j is less than 1.2 (r i + r j ) where r i and r j are the radius of granules i and j. There exists two definitions of fabric tensor and multiple definitions for the fabric deviator. Here, we follow the definitions as reported in [START_REF] Barreto | Quantifying the evolution of soil fabric under different stress paths[END_REF]. They are as follows

Ω ij = 1 N c Nc k=1 n k i n k j ( 22
)
where N c is the total number of quasi-contacts and n k i denotes the component of the normalized centre to centre vector in direction i for the k-th quasi-contact.

Thus, we have the fabric tensor which is a 3×3 matrix. The eigenvectors give the dominant directions of these quasi-contacts. The eigenvalues

(Ω 1 , Ω 2 < Ω 1 , Ω 3 < Ω 2 )
are called major, intermediate and minor magnitudes of the fabric tensor. In an isotropic

situation Ω 1 = Ω 2 = Ω 3 . The deviator fabric is defined by the following equation as reported in [START_REF] Barreto | Quantifying the evolution of soil fabric under different stress paths[END_REF]. The higher its value the more the configuration is anisotropic.

Ω d = 1 √ 2 [(Ω 1 -Ω 2 ) 2 + (Ω 2 -Ω 3 ) 2 + (Ω 1 -Ω 3 ) 2 ] 0.5 (23)

Numerical solving

The simulated domain was of 1 mm in x-direction (R e -R i of the Couette rheometer)

and 200 µm in the two other directions. Particle centroid method is used for coupling the DEM to the CFD. As the name sug-gested the entire particle is assigned to a particular cell based on the location of the centroid of the particle. This method is known to be accurate only when the fluid cell size is significantly larger than the diameter of the particle. Thus, the fluid equations were solved with the finite volume method with N = 5 cells and initially 100 granules in each cell as a compromise. Also, if N is too low, the fluid velocity and temperature profile are not accurately represented and if N is too high, the number of granules in each cell becomes too low to be representative. To ensure periodic condition in y and z direction, for each particle at position (x,y,z) inside the domain, 4 mirror particles (same radius and velocity)

were considered at (x, y + ∆y, z), (x, y -∆y, z), (x, y, z + ∆z)and(x, y, z -∆z).

First, for each granule, its initial radius (r 0 ), its maximum swelling ratio (r max /r 0 ) and its difficulty of swelling (α) are chosen from the distributions as defined in section 2.1. Then its position is randomly chosen to avoid any overlap. This is ensured by checking for overlap at every step and if overlap occurs a new position is chosen until there is no overlap. This methodology works as the initial volume fraction of the granules (φ) is very low.

At each time step, the swelling of each granule is calculated (according to the temperature history of the granule and difficulty of the swelling parameter) in order to re-evaluate its radius and mass. All the different forces and heat fluxes are then estimated. Finally, the velocities, temperatures (fluid and particles) and particle positions are actualized by an explicit first order scheme. Fluid velocity is assumed to be in y direction. Particle positions and velocities are calculated in 3D.

The maximum time step is dictated by the characteristic Stokes time (ρr 2 /µ); ∆t = 10 -5 s was chosen (smallest granules having radius of about 3 µm). The code was implemented in Matlab (MATLAB, 2018). For the fluid (Eulerian part), the velocity and the temperature profiles were considered as piecewise linear and an explicit first order scheme was used.

The presented results correspond mainly to one simulation of N = 5 cubic cells with a side length of 0.2 mm each. In each cell, 100 granules were initially randomly located. This allows comparing what happens with the granules depending on their initial position in the gap. The average shear rate was of 10 s -1 (internal wall velocity 10 mm/s). The precise results depend on the initial position and radius of the granules, four sets of random initial conditions were simulated. The results were comparable and the reported analyze(trends, dominant phenomena, etc..) applies also for the other runs. A comparison between two runs for the evolution of some key variables is presented and the variability of some parameter between the four runs is also presented.

The simulation duration was of 1 min, during which the external wall temperature increased linearly from 60 to 90 • C. The mean volume fraction was initially 2.6 % and 30.0 % at the end of the simulation. If all the granules were completely swollen the volume fraction would have been 32.8 %, so on average, 91 % of the possible swelling is achieved after 1 min of heating.

The liquid properties are those of pure water:

C p = 4180J.K -1 .kg -1 λ f = 0.6W.m -1 .K -1 µ = 2.414 10 -5 10 (247.8/(T f +133)) Pa.s where T f is in • C.
The starch gelatinization enthalpy was neglected. In fact, the enthalpy of gelatinization is less than 10 kJ per kg of dry starch (Desam et al., 2018a;[START_REF] Li | Methodology to predict the time-dependent storage modulus of starch suspensions during heating[END_REF], which corresponds to 0.3 kJ per kg of suspension. This is negligible compared to the sensible heat increase from 60 to 90 • C, which is about 125 kJ per kg of suspension.

Results and Discussion

m p C pp dT i dt = hS(T f -T i ).
A characteristic time appears as follows: t th = m p C pp /(h.S), where S is the surface area of the particle and h is obtained considering N u = 2. This results in a t th which is less s. In our case, ρ f C pf ∆T /∆t is about 2 10 6 W/m 3 .

Figure 3 presents at 9 different time instants, the fluid and particles velocity profiles (v y vs x). Cumulative particle density is superimposed, i.e. the fraction of particles whose centre is at a position lower than x. For a uniform distribution, this line would correspond to a straight line. From this we can see that near the wall there is a region where particles are absent and this region increases slightly as the swelling increases. This is arising due to the wall effects as there exists a minimum distance of particle radius between the centre of the particle and the wall (around 20 µm).

At a given time, since the curvature is ignored, the shear stress is almost constant throughout the gap, i.e. versus x in a quasi steady state approximation. This is not exactly true since the particles are growing and moving during the time, but it should apply on average. Therefore, if we assume that the shear stress is almost the same across the rheometer gap, at the beginning of the simulation, one could expect a higher velocity gradient near the heated wall where the viscosity of the continuous phase is lower. This effect is observed at t = 15 s, the local shear rate varies monotonically between 9.8 s -1 to 10.2 s -1 from internal to an external wall. However, this variation is very small because of the limited temperature difference. At the later instances when granules swell first near the external wall, the opposite effect is expected because the suspension viscosity will be higher due to the increased volume fraction of the granules. A lower shear rate is effectively observed near the heated wall (10.7 s -1 ) than near the internal wall (12.2 s -1 ) at t = 45 s. But in the same time, the particle density decreases near the wall. This explains that the velocity gradient is higher than the average one (10 s -1 ) near the wall and reaches its minimum value (8 s -1 ) near the centre of the gap. This phenomenon of depletion of particles near the wall and subsequent higher shear rate is discussed in literature Coussot (2016). In addition, the interaction between a particle and the smooth wall generates no friction contrary to the interaction between particles. This highlights that the boundary effects are not completely negligible even for a low ratio of particle diameter to gap width, which is about 4 % in our case. In future works the impact of wall lubrication forces and friction forces can be interesting to investigate. It is also observed that the particle velocity fluctuations increase with time. This is because the granules swell and occupy more volume and they interact more and more in a random manner. The granule swelling begins significantly after about 30 s. As we expected a higher degree of swelling is observed near the external heated wall (right side of the graphs) but the difference through the gap appears small, perhaps because of the particle movements. It should be noted that the swelling model was developed by observing a very dilute suspension without any motion or collision and here, we employ the same in presence of particle interactions, these interactions could modify the swelling kinetics. This should be investigated further.

In the figure at 60 s, we see a small grouping of particles.

Figure 5 shows some parameters evolution as the simulation progresses. Here we show two runs of the simulation with different initial properties(radius, swelling difficulty) and positions of the starch granules (reference case at left). Figure 5a wall when its temperature reaches around 70 • C. Similar to the temperature, there is also a time delay for granule swelling from the external to the internal wall. Figure 5c presents the evolution of the mean particle volume fraction φ over all the gap, the internal cell and the outer cell. These time series are directly related to those of granules radius. Near the heated external wall, the volume fraction reaches a plateau whereas the mean value over the gap and the value near the internal (colder) wall are still increasing at the end of heating time of 1 min. Figure 5d shows the evolution of the wall shear stress. Shear stress first decreases slightly because when the temperature is increased, liquid viscosity decreases. But when the volume fraction reaches about 10 %, the shear stress increases significantly. The swelling effect then becomes predominant compared to liquid viscosity decrease. This stress is mainly due to the viscous friction between the fluid layers (represented here by the 5 cells). This friction is proportional to the effective viscosity which increases with the particles volume fraction. Another contribution of global stress is due to the movement of the particle relative to the fluid because of the lubrication forces and the collisions. For example, one particle located near the internal wall with a high velocity in the y-direction can move toward a more external zone, and then via the Stokes forces, it will communicate some momentum (in the y-direction) to this new fluid layer. Alternatively, it can also interact with slower particles of this layer. For a given run, the shear stress on the internal and external wall are almost equal at a given time. Random collisions can explain these fluctuations and the small differences between internal and external values. Nevertheless, averagely constant shear stress across the gap seems ensured.

From figure 5, it can be observed that the temperature evolution is the same for the two runs. The mean radius and volume fraction evolutions over the whole gap are slightly different between the two runs. This can be explained by a small difference in initial volume fraction (run 1: 2.6 % ; run 2: 2.9 %) as well as differences in swelling difficulty and swelling ratio which are randomly chosen. When we observe the internal or external cells, there are more fluctuations because of the random walk of the particles. The most important difference between the two runs is observed for the shear stress (τ xy ). Initially the time series are almost the same, but when swelling becomes significant (t > 30 s) the shear stress becomes higher for run 2. This can be explained by a slightly higher volume fraction (at the end of simulation run1: 30 %, run2: 32 %) this leads to a higher effective viscosity but also to an increase of particle/particle interactions (collisions and lubrication).

$ average for 55 s < t < 60 s 1 presents the variability of some parameters between the four runs (some of them will be presented later). The variability is mainly linked to the range of initial volume fraction. For example run 3 has the lowest initial value of φ and logically the shear stress and the random walk diffusivity (D rw ) are also the lowest. Variability could be reduced by increasing box sizes in y and z direction. To compare easier different runs, the particle properties could be fixed and only the initial position could be varied. To obtain precise quantitative values, average over large number of simulations could be taken. Nevertheless the trends are the same and the phenomena analyzed of run 1 applies also to the other runs.

φ inital (%) φ f inal (%) T iwf inal ( • C) τ xyf inal $ (mP a) D $ rwf inal (µm 2 /s) Ω $ d.f inal d nn.f inal (
In the following only results for run 1 are presented.

To illustrate the movement of the particles inside the gap, Figure 6d presents the time evolution of the x-displacement of 20 particles (initially 4 particles in each of the 5 cells).

At the beginning of the simulation, the granules are small, they interact only weakly by lubrication forces; random walk of low amplitude is observed. After about 40 s, when the volume fraction reaches about 10 % (see Fig 5c), the grain grain interactions induce a random walk movement which becomes more significant and the displacement sometimes exceed 100 µm, which is not negligible compared to the gap width. Figure 6a presents the temperature evolution of these particles. In the beginning, the particle temperatures increase smoothly but when a random walk becomes significant, small fluctuations appear as can be observed in the zoom view in figure 6b and6c. This means that the particle moves randomly between layers of different temperatures (more or less close to the heating wall). When a particle moves towards a cold layer, heat is transferred from this particle to the fluid. This participates in heat transfer from the external to the internal wall. Figure 6e presents the radius evolution of these 20 particles. Initially, the four tracked particles in a cell have the same radius and the difficulty to swell as the four particles in the other cells to see the influence of the initial position in the gap of identical granules. As expected, the particles close to the heating wall begin to swell earlier and faster.

Figure 7a presents the variance of displacement (σ 2 ) in x-direction of all the particles respectively from 35 to 40 s (φ ≈ 10 %) and from 55 to 60 s (φ ≈ 30 %). For one-dimensional random walk it is expected that σ 2 ≈ 2D rw t where D rw is the self-diffusivity of the particles. Different formulas are given in [START_REF] Leighton | The shear-induced migration of particles in concentrated suspensions[END_REF], in a first approach the findings of [START_REF] Eckstein | Self-diffusion of particles in shear flow of a suspension[END_REF] can be used.

D rw r 2 γ ≈ φ 10 (24)
In our case γ = 10 s -1 ; when φ ≈ 10 % most of the granules are unswollen, mean radius is about 7 µm and when φ ≈ 30 % mean radius is about 17 µm, based on this formula the self-diffusivity is expected to be around 4.9 10 -12 m 2 /s and 8.7 10 -11 m 2 /s respectively. This is close to our estimates of D rw , which is 5 10 -12 and 1 10 -10 m 2 /s for these periods respectively (values obtained by linear regression). This is much higher than the Brownian diffusivity of the particles due to the thermal agitation (collisions of the water molecules on the granules) D B = k B T 6πµr ≈ 5 10 -14 m 2 s -1 .

The self-diffusivity due to the random walk of the particles is notably a function of the volume fraction (as shown here) and of the shear rate as shown in the literature [START_REF] Leighton | The shear-induced migration of particles in concentrated suspensions[END_REF]. Its value can be interesting to analyze in other situations, like heat treatment of starch suspension. For example, for a tubular heat exchanger (tube radius R t , residence time t r ), it can reasonably be assumed that the particles follow the fluid flow only

if √ D rw t r << R t .
As explained earlier, the random walk of the particles participates to heat transfer. It is expected that the apparent conductivity is increased by about φρ p C pp D rw . The order of magnitude of this additional conductivity (about 10 -4 W.m -1 .K -1 ) is however low compared to the ordinary conductivity of water (about 0.6 W.m -1 .K -1 ). When this kind of heat transfer is switched off (Nu = 0) there is almost no difference (< 0.15 • C) in the heating rate of the internal wall.

To identify the part of the friction due to particle random walk, Figure 7b At the end of the simulation(φ ≈ 0.30), the two contributions are of the same order.

This highlights the importance of particle random walk in the rheology evolution of a starch suspension.

As shown previously, swelling begins near the heated wall. In addition to the random walk, a global particle migration from the zone of highest volume fraction (near the heated wall) towards zones of lower volume fraction (near the internal wall) is expected.

Figure 7c shows the mean displacement in the x-direction of the granules initially located in the cell near the heated wall vs time. The phenomenon starts effectively with the swelling but it remains small: only about 15 µm. Part of this displacement is simply due to the fact that the centre of the particles is at least one radius from the wall. This is on average 7 µm initially and 17 µm once the granules are swollen. Another part is due to diffusion in the direction opposite to the volume fraction gradient, but this gradient remains weak (Fig5c). It would be interesting in further work to relate particle migration to gradients in osmotic pressure as proposed by [START_REF] Vollebregt | Suspension flow modelling in particle migration and microfiltration[END_REF] but this is not easy in our case due to the heterogeneity across the gap. Thus, we can see that the granule migration towards the internal walls is significant only when the granules begin to swell. The temperature gradient and velocity gradients, which exist initially, do not lead to any significant diffusion unless the granules begin to swell.

The figure 7d shows the number of quasi contacts defined by the criteria that centre to centre distance between particles is less than 1.2×(r i +r j ). Initially when the starch granules are in an unswollen state the number of quasi-contacts are too low to determine the fabric tensor. When this number exceeded 100, the fabric deviator was found to be less than 0.1.

At the end of the simulation, the fabric tensor is as follows

Ω ij,run1 =     
0.3549 -0.0009 -0.0027 -0.0009 0.3151 -0.0025 -0.0027 -0.0025 0.3300

     (25)
and the fabric deviator, Ω d,run1 = 0.060.

One can see very faint anisotropy based on these metrics in the system as the nondiagonal elements are very smaller than the unity. For perfect isotropy, all the non-diagonal elements are equal to zero.

The quadratic mean velocity components of the particles (relative to the fluid) in x, y and z-direction at the end of the simulation are respectively 0.097, 0.096 and 0.050 mm/s.

Particle agitation is anisotropic, this can be explained by the particular direction of shear (in the velocity gradient only the term ∂vy ∂x = 0). This finding is coherent with the results of [START_REF] Breedveld | The measurement of the shear-induced particle and fluid tracer diffusivities in concentrated suspensions by a novel method[END_REF] which found that the self diffusivity was 1.7 higher in the velocity gradient direction than in the other direction. A characteristic distance before direction change (like mean free path before the collision) can be obtained by dividing 2D rw by the quadratic mean velocity. The order of magnitude is of some microns at the end of the simulation. Figure 8b shows the cumulative distribution of the centre to the centre distance to the nearest neighbour at the beginning and the end of the simulation. The centre to centre average distance increases from 25.7 µm (initial random positioning) to 31.2 µm at the end of the simulation. The particles tend to increase their mutual distance, by the swelling effect and by the lubrication forces. However, the final value is far from the maximal possible distance (hexagonal arrangement) : 2 1/6 (V box /I) 1/3 = 48 µm in our case. It was mentioned previously that the particles velocity fluctuations are not isotropic. Indeed, alignment in ydirection could reduce friction (diminish dissipation) as expected from the minimal entropy production principle proposed by Prigogine [START_REF] Prigogine | Time, structure, and fluctuations[END_REF]. From the visual observation of Figure 4, such trends are not obvious even at the final instant, two zones orientated in the y-direction of lower particle density seem to be present. Therefore, to quantify possible heterogeneity/anisotropy of particles repartition in addition to the fabric tensor discussed previously, the position of the nearest neighbour of all the particles were analyzed. The root mean squared average distance to the closest neighbour of each particle in x, y and z directions at the end of the simulation are respectively 19.1, 17.5 and 18.8 µm. This indicates a small trend toward alignment in the y-direction (lower y-direction distance, higher x-direction distance). Isothermal simulations at different shear rates (not presented here) indicate that the alignment trend increases with the shear rate which could imply shear thinning. This has to be confirmed by further work. Figure 9c shows the distribution of the surface to the surface distance to the closest neighbour. The mean distance decreases from 11.1 to 4.0 µm as expected since the volume fraction increases from 2.6 to 30 %. For about 35 % of the particles, the closet neighbour is at less than 1µm at the end of the simulation.

Conclusion

Starch granules swelling inside a heated Couette rheometer is simulated using a CFD-DEM model with a kinetic model for granule swelling. One can see, as expected that the 29 sample is not uniform across the rheometer gap and the same has been quantified through various parameters. The velocity profile is not always linear except at the start of the simulation. The particle temperature is very close to the temperature of the nearby fluid always. It is also showed that random walk of the particles (self diffusivity) arising due to the lubrication force is significant. There is a small migration of particles from the outer walls to the inner cells. The transport of heat by particles is negligible compared to heat transfer by the fluid.

A small anisotropy exists in the direction of shear for the position of the particles in the system. This study was limited to a volume fraction of about 30 % whereas volume fraction up to 70 % are common in processing. For higher volume fraction, thermal effects are still expected to be negligible. However, the particle-particle and particle-fluid interaction will dominate the fluid-fluid interaction. The lubrication/drag interactions and contact frictions will gain prominence; thus, warranting a pure DEM approach by neglecting the fluid phase at φ close to φ rcp . This model should be enhanced, to account for the shearthinning behaviour which is experimentally observed for gelatinized starch suspensions of intermediate solid volume fraction. Further work could integrate full lubrication terms along with angular velocity components as well as lubrication forces between the walls and particles. In conclusion, this work gives insights into the important phenomena that one needs to consider during thermo-mechanical treatment of starch suspensions.

Acknowledgments

This research has been carried out with funding from European Union as part of EU RESEARCH FRAMEWORK PROGRAMME: H2020 / Marie Sk lodowska-Curie Actions ITN MATHEGRAM [813202]. We acknowledge and thank the organizers of the conference Dr. Ruud Vandersman and Dr. Duynhoven John Van. We thank the participants of the conference for their comments, suggestions and inputs. [START_REF] Palanisamy | Multiscale modelling of flow, heat transfer and transformation during thermal treatment of starch suspensions[END_REF] is the conference presentation corresponding to this article.

Figure 1 :

 1 Figure 1: Simulation domain

Figure

  Figure 2a presents at 5 different instants the temperature profile of the fluid vs position in the gap. Similarly, Figure 2b shows the temperature profile of the 500 granules and of surrounding fluid vs position in the gap at a time of 30 s from the start of simulation. The temperature of the particles is always very close to the local fluid temperature as shown in Figure 2b. This is always the case irrespective of the time. This can be explained by the very short characteristic time to reach thermal equilibrium. If we consider the energy balance of one particle (taking into account only the heat exchange with the fluid):

Figure 2 :

 2 Figure 2: Temperature of granules and fluid

Figure 3 :

 3 Figure 3: velocity of granules and fluid

Figure 4 :

 4 Figure 4: xy cross section of the simulation domain

Figure 5 :

 5 Figure 5: Time evolution of suspension parameters at different regions of the simulation domain

Figure 6 :

 6 Figure 6: Sample granule properties evolution during simulation

Figure 7 :

 7 Figure 7: Granule population properties evolution with time

  presents the evolution of the dissipation due to viscous friction between layers of fluid and the dissipation due to lubrication and drag forces (per unit volume: W.m -3 ). The total dissipation is the product of the wall shear stress and the global shear rate. So this figure illustrates how the shear stress and similarly the mixture apparent viscosity evolve and which part is due to fluid/fluid interaction and particle/particle/fluid interactions. The fluid/fluid interaction first decreases because the liquid viscosity decreases when temperature rises. Then it increases because the solid volume fraction increases which enhances the distortion of the fluid streamlines. The particle/particle/fluid interaction (lubrication and drag) is first very low (t < 30 s) except at the simulation starting. Initially, due to the random localisation, some particles are very close to another, this leads to some collisions. Then the lubrication forces tend these particles to move away from each other. When the volume fraction becomes significant (t > 30 s) and the particles moves randomly across the gap the particle/particle interaction increases significantly the friction, or in other terms the mixture apparent viscosity.

Figure 8 :

 8 Figure 8: cumulative distributions of distance from neighbors

Table 1 :

 1 Variability between runs




  Titre de la présentation Date / information / nom de l'auteur Working hypothesis Is multimodal MRI able to detect tumor microenvironment changes during therapy ? Context: TH-302, an innovative treatment NIH -National cancer institute definition: A hypoxia-activated prodrug […] introduces intra-and inter-strand DNA crosslinks in nearby cells; the crosslinks inhibit both DNA replication and cell division, and may lead to apoptosis of cells in the tumor. The inactive form of the prodrug is stable under normoxic conditions, which may limit systemic toxicity. TH-302  Phase III clinical trials (x2)  TH-302 is effective in osteosarcoma  No information about chondrosarcoma  Liapis V. et al., 2015, Cancer Lett. Harms J. et al., 2019, Cells Paratibial Swarm rats follow-up during 5 weeks with TH-302 → (n=35)  Treatment: • 50mg/kg • i.p injection once a day for 5 days followed by 2 days of rest • d+11 to d+22 Aigner T. et al., 2002, Cancer  Vidal A. et al., 2015, Appl. Radiat. Isot.  Ling W. et al., 2008, PNAS

Acknowledgments

Thank you for your attention Thanks to SFRMBM 2021 organization team

Any questions ?