Tung Le Duc
email: tung.le-duc@g2elab.grenoble-inp.fr

Olivier Chadebec

Jean-Michel Guichon

Gérard Meunier

T Le-Duc

J-M Guichon

COUPLING BETWEEN PEEC METHOD AND AN INTEGRO-DIFFERENTIAL APPROACH FOR SOLVING ELECTROMAGNETICS PROBLEMS

published or not. The documents may come

Coupling between partial element equivalent circuit method and an integro-differential approach for solving electromagnetics problems

Introduction

In order to analyse eddy-current in electromagnetic devices, an original formulation is developed. It is particularly well-adapted to model devices that include thin electromagnetic shells (electromagnetic shieldings) and complicated conductor systems.

For this kind of application, the direct use of differential methods, such as the finite element method (FEM), is often difficult and may be sometimes inefficient because of the mesh generation difficulty. Indeed, in this method, mesh refinement around conductors is needed in order to take into account the large variation of the magnetic field close to conductors. To generate this mesh can be a hard task. Moreover, even if this mesh is satisfactorily achieved, such modelling lead to an huge numbers of degrees of freedom and so is limited by memory needs and CPU time requirements if modelling of 3D problems with a large number of conductors is required [START_REF] Guérin | A Shell Element for Computing 3D Eddy Current -Application to Transformers[END_REF], [START_REF] Biro | A Finite Element Formulation for Eddy Current Carrying Ferrromagnetic Thin Sheets[END_REF], [START_REF] Koch | Modeling Thin Conductive Sheets Using Shell Elements in Magnetoquasistatic Field Simulations[END_REF], [START_REF] Sabariego | h-and a-formulations for the time-domain modelling of thin electromagnetic shells[END_REF].

The PEEC method (Partial Element Equivalent Circuit) is mainly used for the modelling of complex electrical interconnections and can be applied to a large range of devices where the air region is dominant [START_REF] Hoer | Exact Inductance Equations for Rectangular Conductors with Applications to More Complicated Geometries[END_REF], [START_REF] Ruehli | Equivalent circuit models for three dimensional multiconductor systems[END_REF], [START_REF] Clavel | Electrical Modeling of Transformer Connecting Bars[END_REF], [START_REF] Antonini | PEEC capacitance extraction of 3-D interconnects[END_REF]. In this method, only conductors are described and meshed, so it generally requires less memory and CPU time than the FEM method. However, the classical PEEC method does not enable the 3D modelling of thin electromagnetic regions.

On the other hand, some integral formulations (usually called integro-differential one) have been developed in order to analyze eddy current in conducting shells [START_REF] Kalimov | Dynamic Processes in Laminated Magnets: Simulation and Comparison with Experimental Results[END_REF], [START_REF] Salon | An Intergro-Differential approach to eddy currents in thin plates[END_REF], [START_REF] Tsuboi | Eddy Current and Deflection Analyses of a Thin Plate in Time-Changing Magnetic Field[END_REF]. Like the PEEC method, it does not require the meshing of the air region.

In this paper, an original IDM/PEEC coupling is presented. The modelling of thin electromagnetic shells will be achieved thanks to the integro-differential method (Part II). The PEEC method allows us to model the contributions of inductors fed (Part III).

The coupling of those formulations, in order to compute the interactions between conductors and electromagnetic shells, will be developed in part IV. Finally, a numerical example will be presented in part V. Results obtained with this formulation and the Finite Element Method (FEM) will be compared in this last part.

2 Integro-differential approach for modelling thin electromagnetic shells

Low frequency problems are focused (i.e. the skin depth is high in comparison with the thickness of the shell) so the assumption of a uniform current distribution across the thickness of the shell is made. Moreover, eddy currents C J are assumed to flow tangentially to the shell.

Let us now consider a thin electromagnetic shell placed in an inductor field H0.

Equations governing the problem in shells are:

B E    j Curl (1) C Curl J H  (2)
E J   C (3) ) ( 0 0 M H H B       r ( 4
)
where  is the conductivity, C J is the eddy current distribution in the shell, r is the relative permeability of the material and M is the magnetization. In our application, we are mainly interested in shielding effects, so the shell will be placed in a quite low magnetic field. So, the permeability will be considered as uniform and the problem as linear.

The Ampere-Maxwell equation (2) implies the existence of an electric vector potential T as:

C Curl J T  (5)
As the current is tangential and uniform across the thickness, T is normal to the shell and can be consider as scalar (i.e. its normal component).

Equation ( 1) can be rewritten as follows [START_REF] Salon | An Intergro-Differential approach to eddy currents in thin plates[END_REF]:

n r j Curl Curl H T      0 )) ( 1 ( (6)
where n H is the normal magnetic field. Considering the relation linear between the magnetization M and the magnetic field intensity H, equation ( 6) can be rewritten as:

n r r j Curl Curl M T 1 )) ( 1 ( 0        (7)
These equations have to be discretized. The easiest way is to mesh the thin shell into n triangular prisms elements associated with a constant magnetization (0-order shape functions for M). The unknown current potential variables T is approximated with linear functions (1-order shape functions so uniform value of current J per element).

Let p denote the number of node. A Galerkine projection method of ( 7) is then applied in order to get the following matrix system:

0 M B T A     (8)
where unknowns are T a ( 1  p ) vector and M a (

1 3  n ) vector. Matrix A is a ( p p  ) matrix and B is a ( n p 3  ) matrix.
The total magnetic field H is the sum of H0, the source field, Hm the reaction of the magnetization and Hc the reaction of the eddy current flowing in the thin plate. Let Vm and  denote the volume of magnetic material and its middle surface respectively. An integral equation linking the local field to the magnetization and the electric vector potential can be written [START_REF] Kalimov | Dynamic Processes in Laminated Magnets: Simulation and Comparison with Experimental Results[END_REF]:

            d r Curl e dVm r r P P Vm 3 3 5 ) ( 4 ) ) ( 3 ( 4 1 ) ( ) ( r T M r r M H H 0 ( 9
)
where P is a point located in the domain region, r is the vector linking the integration point to P and e is thickness of shell. Thanks to a point matching approach of this material law at the centre of each element, we get a linear matrix system: [START_REF] Salon | An Intergro-Differential approach to eddy currents in thin plates[END_REF] where d I represents the identity matrix; h0 is the source fields at each barycentre, F is a (3n×3n) matrix and C is a (3n×p) matrix.

  0 ) 1 /( h T C M F I        r d
Equations ( 8) and [START_REF] Salon | An Intergro-Differential approach to eddy currents in thin plates[END_REF] are solved together to obtain the magnetization M of each element and the current potential T of each node. This method is a hybridisation between a differential equation and an integral one; this is why it is usually called an integro-differential method.

Inductive PEEC method

Let us consider m volume conductors fed with voltage sources placed in a surrounding air region without any magnetic materials. The well-known PEEC method is particularly reliable to solve this kind of problem. It is based on the determination of partial voltage generated on each conductor by electromagnetic sources. To compute these voltages, volume integration on the conductor of the magnetic vector potential created by all the others conductor is provided. Let us assume that the current density in each conductor is uniform. The expression of the magnetic vector potential j A created by conductor j is [5], [6]:

    j V j j j j j dV r S I P u A 4 ) ( 0 ( 11
)
where r is the distance between the integration point and the point P ;

          k j k V k V j k j k j j V k k j k kj dV dV r S S I dV S u u u A 4 1 0 (12)
The mutual inductance between two conductors is defined as:

        k j V k V j k j k j j kj kj dV dV r S S I m u u 1 4 0 (13)
Thus, the voltage appearing on the conductor k is given by:

       m j kj j k k cond k m I j R I U 1 _ ( 14
)
where  is the angular frequency, k R is the resistance of k -th conductor. By applying equation ( 14) to all conductors, we get a matrix system known as impedance matrix system:

U I Z   ( 15
)
where Z is a (m×m) matrix, I is a (m×1) vector and U is a (m×1) vector.

Coupling principle

Let us now consider a system composed of m conductors with unknown flowing currents I and thin magnetic conductive shell. The thin shell is still meshed into n triangular prism elements associated with a constant magnetization and a constant value of current per element. As in part two, the interactions between each element in the thin shell are calculated thanks the integro-differential method. The PEEC method is used to model interconnection between conductors. It remains to write explicitly the couplings matrix between both methods.

Influence of current conductor in thin electromagnetic shell

The term h0 of equation ( 10) can be computed by using Biot and Savart law. The magnetic field created by m conductors is given by:

j m j V j j dV r I P j       1 3 0 4 ) ( r u H (16)
with same notations as in [START_REF] Tsuboi | Eddy Current and Deflection Analyses of a Thin Plate in Time-Changing Magnetic Field[END_REF]. We write equation ( 16) at barycentre of each thin shell element to obtain the following matrix block:

I BS h   0 ( 17
)
where BS is a (3n×m) matrix, I is a (m×1) vector. We substitute (17) into ( 10) and obtain:

0 I BS T C M F I          ) ) 1 /( ( r d (18)

Influence of magnetization and eddy currents of thin shell on conductors

Like in the standard PEEC method, we have to integrate the magnetic vector potential i A created by each thin shell elements i on the conductor k. The expression of this magnetic vector potential is:

           i i i i Vm i i i d r Curl e dVm r P ) ( 4 4 ) ( 0 3 0 T r M A (19)
with same notations as in [START_REF] Kalimov | Dynamic Processes in Laminated Magnets: Simulation and Comparison with Experimental Results[END_REF]. The magnetic vector potentials of n shell elements create a new voltage on the conductor k:

k k i n i i V i n i Vm k shell k dV d r Curl e dVm r S j U i k i u T r M i                  ) ( 4 1 1 3 0 _ (20)
Thus, the total voltage appearing on the conductor k is given by:

) ( 1 1 1 _ _                  n i p i ki i ki i m j kj j k k shell k cond k k mc T mf M m I j R I U U U (21)
where ki mf is defined as the mutual inductance between the magnetization Mi of shell element and the conductor k; ki mc is the mutual inductance between the current potential Ti of the shell node and the conductor k.

If we write this equation for all conductors, we get a new impedance matrix system:

T L M L I Z U       C f (22)
where f L is a (m×3n) matrix and C L is a (n×p) mutual inductance matrix between thin shell elements and conductors.

Finally, equations ( 8), ( 18) and ( 22) are brought together in the global (3n+p+m) × (3n+p+m) following matrix systems:

                                 U 0 0 I T M Z L L 0 A B BS C MM C f (23)
where MM is the matrix (

F I    ) 1 /( r d
)( see equation [START_REF] Salon | An Intergro-Differential approach to eddy currents in thin plates[END_REF]).

In (23), I and U are the current and voltage values of each branch of the conductors system. In order to solve this equation system with a reduced number of degrees of freedom, we use the Kirchhoff's mesh rule. Consequently, current and voltage values become associated to each of independent circuit mesh. We get then a lighter system and after solving it with a direct linear solver (LU decomposition), we obtain magnetizations, current potentials in thin shell and currents in conductors.

Numerical Example

Let us consider a conducting ferromagnetic plate with the conductivity  =9.6E+6 S/m and relative permeability r =50 (reversible linear permeability representative of a classical shielding effect). The thickness of the plate e=1.5mm and the dimensions are 200x200 mm. This plate is placed close to a conductor coil which is fed with a voltage source (1V, 50Hz) (Figure 1).

FIGURE 1 HERE

Our coupling is compared with results given by FEM in Flux software [12]. We focus on the computed eddy current distribution and Joule losses in the plate (Figure 2,3).

Whereas the convergence is quickly reached with our coupling, the problem needs a very fine mesh to be accurately solved with FEM because of the large variation of the field close to the conductor [START_REF] Tran | Coupling PEEC-Finite Element Method for Solving Electromagnetic Problem[END_REF].

We can notice small differences between Joule losses computed values. If we consider that FEM method as our reference, the IDM/PEEC coupling leads to an error of 0.8% for the Joule losses in plate (Figure 3). Results provided by our coupling are very encouraging, the convergence being reached with a very few number of elements (around 2000). Of course, the obtained matrix is fully dense, but the computation time is divided per ten in comparison with FEM. Let us notice that this problem of fully dense matrix and associated memory need could be fixed with the use of a compression technique like « Fast Multipole Method » (FMM) [START_REF] Buchau | Parallelization of a Fast Multipole Boundary Element Method with Cluster OpenMP[END_REF] for instance.

Conclusions

In this paper, we have presented a coupling between PEEC and an integro-differential method in order to introduce the thin electromagnetic shell in PEEC method. Through a simple example, we have validated this formulation. The result shows that the coupling is much faster and more accurate than the FEM method on this simple geometry. However, our formulation can only be applied for problem where the skin depth is high in comparison with the thickness of the shell ( e  

). In further work, formulations enabling smaller skin depth will be investigated.

j I is the current in the conductor j ; j S

 j is the cross-section of the current flow; j V is the volume and j u is the unit vector whose direction is that of conventional current j I .The flux created by j I through the conductor k is:

  FIGURE 2 HERE

UFigure 1 :Figure 2 :

 12 Figure 1: Geometry of the numerical example

Figure 3 :

 3 Figure 3: Joule losses in thin plate (W) calculated by IDM/PEEC and FEM method in


Indra Taufiq

Ryohei Rukmana1

Gabriela Yasukuni1

Rachel Moran2

Gilles Méallet-Renault2

Tadashi Clavier3

Misato Kunieda4

Ohtani4

Taku Demura4

Yoichiroh Hosokawa1 Affiliations

Direct observation of nanoparticle diffusion in cytoplasm of single plant cells realized by photoinjection with femtosecond laser amplifier

Diffusion is an important process for molecular transport inside plant cells. Recent advancement in plant physiological study demands verification of the diffusion process at the single cell level. In this work, the real-time intracellular diffusion of nanoparticles in the cytoplasm of single plant cells was realized using photoinjection with femtosecond laser amplifier. The diffusion behavior was analyzed by estimating the diffusion coefficient in cytoplasm. In addition, the effect of particle size to the photoinjection efficiency and diffusion was evaluated. Surprisingly, the intercellular diffusion of nanoparticles between single plant cells was also directly observed.

Particle diffusion in cytoplasm is an essential process for plant cells to function properly. The diffusion of functional molecules, such as proteins, nucleic acids, hormones, and nutrients, among others, through cell cytoplasm to the intended position in cells is important to be ascertained in the study of plant physiology. [1][2][3][4] Recent development shows an increasing interest for the direct observation of the diffusion process at the single cell level, especially to understand molecular transport, biological phase separation, and cell-to-cell communication.5-8) However, there has been a lack of available method to directly observe the diffusion in cytoplasm of single plant cells, particularly for large molecules.

Recently, the photoinjection of large molecules into single intact plant cells using femtosecond (fs) laser amplifier has been reported.9,10) The macromolecules diffused through a pore on the cell wall and membrane into cytoplasm. Considering this breakthrough, the method offers a way to monitor real-time diffusion of macromolecules in cell cytoplasm. In this work, the intracellular diffusion of nanoparticles in single plant cells was directly observed using photoinjection with NIR fs laser amplifier coupled with a confocal laser scanning microscope. The diffusion behavior was assessed by estimating the diffusion coefficient in cytoplasm, and the effect of particle size was briefly evaluated. In addition, the direct observation of intercellular diffusion between cells was also demonstrated.

As a representative plant cell sample, we used tobacco BY-2 (TBY-2) cell. Wild-type TBY-2 cells were cultured as described previously.9-12) Two-day old cultured cells were used for the photoporation experiments. 0.5M mannitol was added into the final culture medium to control the osmotic pressure. To partially degrade the cell wall, the cells were incubated in the cellulase and pectolyase enzyme solution in a dark place at 27 °C for 10 min.9-11,13) After these cell preparations, 100 μl of the prepared cell suspensions were placed in a glass-bottom dish (j = 12 mm, Iwaki). All the experiments were conducted within 6 h to maintain a healthy cell condition.

Femtosecond laser pulses from a regeneratively amplified fs Ti:Sapphire laser system (Spectra-Physics, Solstice-Ref-MT5W, 800 nm, 150 fs) were led into a laser-scanning confocal microscope (Olympus, IX71FVSF-2) and focused on a TBY-2 cell through a 100x oil-immersion objective lens (Olympus; Plan N, NA = 1.25). A single laser pulse was obtained using a mechanical shutter with a gate time of 1/ 125 s from pulse trains of 125 Hz. The laser pulse energy through the objective lens was measured with a laser power meter (Ophir, Nova Display-Rohr). A sample was put on a motorized microscope stage (Sigma Koki, E-65GR).

Ultrabright green fluorescent nanoparticles (G-FNPs) made of copolymer of styrene and fluorescent BDPMA (borondipyrromethene (BODIPY) methacrylate) were used for sensitive visualization of nanoparticle injection into TBY-2 cells. The volume of G-FNPs is about 20 times and 400 times larger than average genes (2 MDa) and proteins (30 kDa), respectively.14,15) The G-FNP preparation was described in previous reports. 16,17) The photoinjection and diffusion of G-FNPs were observed by confocal fluorescence imaging. A diode-pumped solid-state (DPSS) laser (Spectra-Physics, PC14763) at 488 nm wavelength was applied for the excitation. 1 μl of the G-FNP solution was added to 100 μl of the TBY-2 cell suspension in the glass-bottom dish, and the sample was put on the microscope stage. All fluorescence images were captured with the same acquisition parameters.

The representative results of photoinjection of the G-FNPs into single TBY-2 cell and subsequent diffusion of the fluorescent nanoparticles in the cell cytoplasm are shown in Fig. 1. The prepared G-FNPs were 80 ± 2.0 nm in diameter (d) with concentration of 7.3 ± 0.34 × 108 particles ml -1 . A single 20 nJ fs laser pulse was focused at the contact point of cell wall and membrane. After the single shot irradiation, intracellular diffusion of the G-FNPs was evaluated based on the differential fluorescence intensity (ΔI) as a function of distance from the laser focal point. The ΔI value was obtained by subtracting the mean intensity value before photoporation (Ibefore) from that of after photoporation (Iafter). ΔI = Iafter -Ibefore. (1) The temporal evolution of fluorescence intensity is shown in Fig. 2. Fluorescence intensity increased immediately after laser irradiation at the laser focal point (0 μm) and reached to half maximum at 45 s. At the distances of 10 μm and 20 μm from the laser focal point, the fluorescence intensity increased slowly from 45 s and 60 s after laser irradiation and reached to half maximum at 120 s and 150 s, respectively. The ΔI value was the highest at the laser focal point. Indeed, the area nearby laser focusing position was the brightest part, as shown in Fig. 1. It was suggested that the fluorescent nanoparticles were accumulated at the laser focal point before diffusion in the cell cytoplasm.

The orange, green, and purple circles in (d) indicate the points of fluorescence intensity measurement for the evaluation of intracellular diffusion, approximately at 0 μm, 10 μm, and 20 μm from the laser focal point, respectively. The fluorescence intensity of whole cytoplasm area, surrounded by blue and brown dashed lines in (d) for the target cell and its adjacent cell, respectively, was measured for the evaluation of intercellular diffusion.

Fig. 2. Temporal evolution of fluorescence intensity in the cytoplasm for enzyme-and-mannitol treated cells after single 20 nJ pulse fs laser irradiation at various distances from the laser focal point (0 μm, 10 μm, and 20 μm). The graphs are averaged for 8 cells and normalized by the ΔI at 240 s. The error bar indicates standard deviation.

The diffusion coefficient of the G-FNPs in the cytoplasm was estimated with this equation. 18,19)

MSD = 2nDt. (2)
MSD is the mean squared displacement (|Δx| 2 ), with Δx as the distance between each measurement point; n is the dimension of particle movement (3); D is the diffusion coefficient; and t is the time differential between each point to reach the half maximum intensity. The calculated diffusion coefficients of the G-FNPs in cytoplasm were about 0.2 μm2 s-1 and 0.6 μm2 s for 0-10 μm (D0-10) and 10-20 μm (D10-20), respectively. The D0-10 value was lower than the D10-20 value, confirming that the initial nanoparticle diffusion in cytoplasm, from 0 to 10 μm, was slower than the later diffusion, from 10 to 20 μm. Beside nanoparticle accumulation at the laser focal point, the different value of the diffusion coefficient implied that there might be a biological phase separation in the cytoplasm. Meanwhile, the nanoparticle diffusion coefficient in water (Dw), calculated with Stokes-Einstein equation,20,21) was about 6.1 μm2 s-1, suggesting that the nanoparticle diffusion in cytoplasm was much slower than in water.

The effect of particle size to the photoinjection efficiency was also examined, as summarized in Fig. 3. Beside the GFNPs (d = 80 nm), other objects with different particle sizes were used: fluorescein isothiocyanate (FITC) conjugated with 20 kDa dextran (FITC-20k), d = 3.2 nm; FITC conjugated with 2 MDa dextran (FITC-2M), d = 26.7 nm; and red fluorescent nanoparticles (R-FNPs), d = 110 nm. The diameter of dextran molecules was calculated with this formula.22,23)

SR[nm] = 0.81x(MW[kDa]) 0.46 . (3)
SR is the Stokes radius in nm and MW is the molecular weight in kDa. Smaller particles were injected more efficiently into cell cytoplasm.24) The diffusion through a similar sized pore became slower with increasing particle size.25) The maximum size of injected particles into enzyme-and-mannitol treated cells after irradiation with single 20 nJ pulse from fs laser amplifier was estimated between 80 nm and 110 nm in diameter. Interestingly, in Fig. 1, we found that the nanoparticles also diffused into the adjacent cell attached to the target cell. To ascertain this, the differential fluorescence intensities before and after laser irradiation in the target cell and adjacent cell were also measured as shown on Fig. 4. In this case of intercellular diffusion, the whole cytoplasm area of the target and adjacent cells was chosen as a region of interest. It turned out that the fluorescence intensity also increased in the adjacent cell, even though it was lower than that of the target cell. There are many tunnels that connect the cytoplasm between plant cells, called plasmodesmata. Plasmodesmata typically have a diameter of about 50 to 60 nm in tobacco cells.26) The plasmodesmata could be closed or enlarged in response to cell conditions, such as turgor pressure. In this case, turgor pressure was decreased and the plasmodesmata might be enlarged.27) In addition, the enzyme treatment might have some effects in the enlargement of plasmodesmata, although it still needs to be verified.28) Nevertheless, the gradual diffusion of 80 nm nanoparticles through this tunnel into the adjacent cell was observed. In previous studies, the movement of nanoparticles in plant organs is analyzed indirectly measuring the concentration of nanoparticles in leaves and stems after submerging roots or spraying leaves with the nanoparticle suspension. [29][30][31] In contrast, the intercellular diffusion of nanoparticles between single plant cells was directly observed in this work, offering a real-time analysis, which is very promising for the study of plant physiology, particularly to examine molecular transport and cell-to-cell communication.

In conclusion, the direct observation of intracellular and intercellular diffusion of nanoparticles in the cytoplasm of single plant cells was realized by photoinjection utilizing fs laser amplifier. The ability to observe and evaluate such processes is a major development in the study of plant physiology, especially at the single cell level.

Fig. 1 .

 1 Fig. 1. Confocal fluorescence images of enzyme-and-mannitol treated TBY-2 cells before (0 s) and after single 20 nJ pulse fs laser irradiation (60, 120, 180, and 240 s) under the presence of G-FNPs. The top left image (a) shows the transmission image before laser irradiation. The white and yellow dashed lines in (b) indicate the outlines of a target single cell and its adjacent cell, respectively. The white and blue arrowheads in (b) and (c) indicate the laser focal point and the bright fluorescence observed near the laser focal point after laser irradiation, respectively. The orange, green, and purple circles in (d) indicate the points of fluorescence intensity measurement for the evaluation of intracellular diffusion, approximately at 0 μm, 10 μm, and 20 μm from the laser focal point, respectively. The fluorescence intensity of whole cytoplasm area, surrounded by blue and brown dashed lines in (d) for the target cell and its adjacent cell, respectively, was measured for the evaluation of intercellular diffusion.

Fig. 3 .

 3 Fig. 3. Photoinjection efficiencies into (a) mannitol treated and (b) enzyme-and-mannitol treated tobacco BY-2 cells after irradiation with single 20 nJ pulse from fs laser amplifier under presence of objects with different particle sizes: FITC-20k (d = 3.2 nm), FITC-2M (d = 26.7 nm), GFNPs (d = 80 nm) and R-FNPs (d = 110 nm).

Fig. 4 .

 4 Fig. 4. Temporal evolution of fluorescence intensity in the cytoplasm for enzyme-and-mannitol treated target cells and their adjacent cells after single 20 nJ pulse fs laser irradiation. The graphs are averaged for three cells and normalized by the ΔI at 240 s of the target cell value. The error bar indicates standard deviation.

Acknowledgments This work was supported in part by the Grant-in-Aid for Scientific Research (C) (No. JP18042069), Scientific Research on Innovative Areas (No. JP18H05493) from the Japan Society for the Promotion of Science (JSPS). The collaboration between Japan and France was supported by JSPS and MEAE-MESRI under the Japan -France Integrated Action Program (PHC Sakura 40951TG).


Nicolò R Sgreva
email: nicolo.sgreva@u-psud.fr

Anne Davaille

Ichiro Kumagai

Kei Kurita

Interaction between a falling sphere and the structure of a non-Newtonian yield-stress fluid

Keywords: Yield stress, Falling sphere, SAP, Irregular motion, Intermittent motion

We present an experimental study using mixtures of aqueous superabsorbent polymers (SAP) where we systematically investigate the influence of the size of grains that make up the fluid structure on the effective rheology and its domain of validity. In water, these polymer powder grains swell up to 200 times and form gel grains whose size can be controlled by controlling the size of the initial powder. The rheology of this mixture (water and touching grains) combines viscous, elastic and plastic aspects and can be characterized using the free-fall of spheres of different diameters (between 3 and 25 mmdiameter) and densities (from 2200 to 15000 kg/m 3 ). As the typical size of the gel grains was varied between 1 and 8 mm, there is a range where it becomes comparable to the size of the falling spheres. We observe five different motion regimes. (1) A linear regime where the sphere has a rapid and linear fall and reaches a constant terminal velocity. (2) An irregular regime where the sphere's velocity varies around a constant value. (3) An intermittent regime where periods of no-motion and periods of irregular falls follow one another.

Introduction

Yield-stress fluids are a broad category of fluids which have been intensively studied in the last decades due to their key roles both in industry (e.g. drilling fluids, cement pastes, foams) and in natural phenomena (e.g. lava and mud flows, avalanches, landslides). Very often, the material consists typically of a suspension of particles in a liquid. When the particle volume fraction increases, they become in close contact with one another, and the material can face jamming. In this jammed state, the material can support stresses without flowing. A so-called yield-stress fluid does not flow if the applied stresses are not high enough to unjam the structure, but it does flow when the stresses become larger than a threshold value, the yield stress (σ Y ). The latter introduces a non-linearity into the material rheology, which strongly impacts on its dynamics [START_REF] Makse | Statistical mechanics of jammed matter[END_REF][START_REF] Bonn | Yield stress materials in soft condensed matter[END_REF][START_REF] Coussot | Rheometry of pastes, suspensions, and granular materials: applications in industry and environment[END_REF][START_REF] Barnes | A review of the slip (wall depletion) of polymer solutions, emulsions and particle suspensions in viscometers: its cause, character, and cure[END_REF]. An example of these materials are soft particle glasses. Above the yield stress, their flow results from the concurrence of particles elastic deformations and plastic events during which particles rearrange [START_REF] Seth | A micromechanical model to predict the flow of soft particle glasses[END_REF]. Plastic local events can influence each other and lead to non-local elastic relaxation. This collective rearrangement can be the reason of spatial variation of viscosity in an homogeneous stress field in emulsions [START_REF] Goyon | Spatial cooperativity in soft glassy flows[END_REF][START_REF] Goyon | How does a soft glassy material flow: finite size effects, non local rheology, and flow cooperativity[END_REF]. Other finite size effects to the bulk rheology of a jammed system appear due to concentration gradients of particles in non-Brownian particle suspensions [START_REF] De Cagny | Local rheology of suspensions and dry granular materials[END_REF] or due to structural heterogeneities in polymers [START_REF] De Kort | Yielding and flow of cellulose microfibril dispersions in the presence of a charged polymer[END_REF].

Several rheological models have been proposed to describe the rheology above the yield stress, such as the Bingham and the Herschel-Bulkley "HB" models. One issue is the domain of validity of these rheology models: how well do they describe the fluid behavior close to the yield stress and how much larger than the microstructure typical length scale must be the flow length scale? For example, if we take the simple (and well known in Newtonian fluids) problem of a solid sphere settling in a particulate fluid, when and how will the fluid microstructure interact with the sphere motion? There are two end-members depending on the size ratio between the fluid microstructure d g (e.g. the size of the particles in suspension in the liquid) and the falling sphere d s . When d g /d s >> 1, the sphere will fall through a porous medium, moving along the liquid channels between particles, and when d g /d s << 1, the sphere will see a continuum with an effective rheology (e.g. "HB" model).

But what happens for intermediate d g /d s ?

The d g /d s << 1 case has been well studied for viscoplastic fluids, such as Carbopol [START_REF] Bonn | Yield stress materials in soft condensed matter[END_REF]: the sphere reaches a steady-state downward motion only if a critical force is overcome, otherwise the fall is prevented [START_REF] Beris | Creeping motion of a sphere through a bingham plastic[END_REF][START_REF] Tabuteau | Drag force on a sphere in steady motion through a yield-stress fluid[END_REF][START_REF] Balmforth | Yielding to stress: recent developments in viscoplastic fluid mechanics[END_REF]. However, the accurate value of this critical force is still debated and it varies, for example, considering or not the elasticity of the material [START_REF] Fraggedakis | Yielding the yieldstress analysis: a study focused on the effects of elasticity on the settling of a single spherical particle in simple yield-stress fluids[END_REF]. Once the sphere is set in motion, crucial aspects are the determination of the drag coefficient [START_REF] Beris | Creeping motion of a sphere through a bingham plastic[END_REF][START_REF] Tabuteau | Drag force on a sphere in steady motion through a yield-stress fluid[END_REF][START_REF] Blackery | Creeping motion of a sphere in tubes filled with a bingham plastic material[END_REF][START_REF] Deglo De Besses | Sphere drag in a viscoplastic fluid[END_REF] and the determination of the shape and position of the yield surface that separates the solid-like region from the fluid-like one.

The yield surface and the drag force on the sphere have first been determined numerically and analytically for a Bingham fluid by Beris et al. [START_REF] Beris | Creeping motion of a sphere through a bingham plastic[END_REF] and subsequently for an HB fluid numerically [START_REF] Blackery | Creeping motion of a sphere in tubes filled with a bingham plastic material[END_REF][START_REF] Liu | Convergence of a regularization method for creeping flow of a bingham material about a rigid sphere[END_REF][START_REF] Deglo De Besses | Sphere drag in a viscoplastic fluid[END_REF] and experimentally by Tabuteau et al. [START_REF] Tabuteau | Drag force on a sphere in steady motion through a yield-stress fluid[END_REF] in Carbopol. Within the liquid-like region, the flow in Carbopol is confined nearby the sphere, within an area which size depends on the yield stress of the fluid [START_REF] Putz | Settling of an isolated spherical particle in a yield stress shear thinning fluid[END_REF], and can be either symmetric [START_REF] Dinkgreve | Carbopol: From a simple to a thixotropic yield stress fluid[END_REF] or asymmetric [START_REF] Putz | Settling of an isolated spherical particle in a yield stress shear thinning fluid[END_REF], depending on the sample preparation.

Similar results had been found for Laponite, a very thixotropic colloidal suspensions [START_REF] Gueslin | Sphere settling in an aging yield stress fluid: link between the induced flows and the rheological behavior[END_REF],

where the asymmetry of the flow pattern increases with the age of the fluid through the appearance of a negative wake (i.e. upward fluid motion in the sphere's wake [START_REF] Hassager | Negative wake behind bubbles in non-newtonian liquids[END_REF]) and in viscoelastic aqueous polyacrylamide solutions [START_REF] Arigo | An experimental investigation of negative wakes behind spheres settling in a shear-thinning viscoelastic fluid[END_REF]. Another peculiar characteristic of a falling sphere in yield stress fluids or in viscoelastic materials is that its velocity can depart from a constant value and show oscillations and irregularities. One of the first observations of an oscillating particle settling has been in hydroxyl propyl guar (HPG)

polymer gels [START_REF] Mollinger | An unexpected phenomenon observed in particle settling: oscillating falling spheres[END_REF], a yield-stress power-law family of fluids [START_REF] Jiang | The rheological characterization of hpg gels: Measurement of slip velocities in capillary tubes[END_REF]. In this work, the authors suggested that the irregularities in the descending motion were linked to the elasticity of the gel. Something similar is reported in entangled wormlike micellar fluids [START_REF] Jayaraman | Oscillations of a solid sphere falling through a wormlike micellar fluid[END_REF][START_REF] Kumar | Oscillatory settling in wormlike-micelle solutions: bursts and a long time scale[END_REF], in Laponite [START_REF] Fazilati | Quasi-periodic and irregular motion of a solid sphere falling through a thixotropic yieldstress fluid[END_REF], where the vertical velocity can oscillate in bursts and, again, in HPG [START_REF] Weidman | On the instability of spheres settling through a vertical pipe filled with hpg[END_REF]. In cornstarch suspensions [START_REF] Von Kann | Nonmonotonic settling of a sphere in a cornstarch suspension[END_REF][START_REF] Von Kann | Velocity oscillations and stop-go cycles: The trajectory of an object settling in a cornstarch suspension[END_REF], the sphere velocity never approaches a steady terminal velocity but instead it oscillates, decreases with time, to reach a series of stop-go cycles as the sphere becomes close to the experimental tank bottom. In wormlike micellar fluids, the cause of non-transient oscillations is to be found in the formation and successively breakage of flow-induced structures; in cornstarch suspensions, in the formation and dissolution of a jammed layer that, under sufficient stress, increases drag and slows the sphere down; while the irregular motion in laponite is due to the existence of flow instabilities and shear banding.

In this paper, we consider viscoplastic fluids that are constituted of an aqueous suspension of soft gel grains. To do so, we place ourselves in between the two end-members described above, using macroscopic super absorbent gels. Varying independently the size of the gels grains and the size of the falling spheres allows us to investigate the influence of the particle size on the effective rheology of the fluid and on the motion of spheres that are falling through it. The paper is organized as follow: in section 2 we describe the experimental set-up and the fluids used during the experiments. Section 3 presents the different regimes of motions that were observed. We close with section 4 and section 5 where results and their implications are discussed and summarized.

Experimental set-up and fluids

The fluids

According to the product company (Omiya Green Service Co., Ltd), the fluid used is a superabsorbent polymer (polyacrylamide) made by copolymerization of acrylic acid and acrylamide. It is part of the family of aqueous dispersions of superabsorbent polymers (SAPs), polymeric materials with large capacity of water absorption [START_REF] Zohuriaan-Mehr | Superabsorbent polymer materials: a review[END_REF][START_REF] Ramazani-Harandi | Rheological determination of the swollen gel strength of superabsorbent polymer hydrogels[END_REF] which are used in many different water absorbing applications in agriculture [START_REF] Guilherme | Superabsorbent hydrogels based on polysaccharides for application in agriculture as soil conditioner and nutrient carrier: A review[END_REF], in health care industry (e.g. in sanitary pads and baby diapers), in sealants and in air-fresheners [START_REF] Po | Water-absorbent polymers: a patent survey[END_REF][START_REF] Graham | Modern superabsorbent polymer technology[END_REF]. The dry material is a white and granular powder of particles with different irregular shapes. In water, the original small particles can swell up to 200 times and form gel grains whose size can be controlled by controlling the size of the initially dry grains. Beside the general behaviour of such materials that combines viscoelastic [START_REF] Ramazani-Harandi | Rheological determination of the swollen gel strength of superabsorbent polymer hydrogels[END_REF] and shear-thinning aspects [START_REF] Zhang | Using screening test results to predict the effective viscosity of swollen superabsorbent polymer particles extrusion through an open fracture[END_REF], the measurement of appropriate rheological properties at the swollen state is still quite complicated as typical commercial SAPs have large (up to few millimeters) and irregular swollen particles that render classical rheological methods difficult to use [START_REF] Ramazani-Harandi | Rheological determination of the swollen gel strength of superabsorbent polymer hydrogels[END_REF].

Starting from the provided raw material as the largest grains size endmember, we gradually ground the dry particles to form four thinner powders.

All the final gel samples were then obtained by letting 17.5 ± 0.1 g of dry product react with 3.5 l of distilled water for one night. The resulting swollen mixture have been gently stirred (100 rpm) for 2-3 days with an electric stirrer in order to remove air bubbles and to homogenize the mixture. Careful weight measurements of a given volume of the final aqueous suspension before and after it was completely drained from the free water allowed to estimate a particle volume fraction of 62±2%.

As it was for the starting dry powders, the five final fluids under scrutiny differ from each other only in the size of the grains that make up their structure. To obtain the grains size distribution of the swollen gel samples, we analyzed several high resolution images (4288 x 2848 pixels) of a laserilluminated vertical cross-section of the experimental tank (Fig. 1) by using Multispec, an image analysis software [START_REF] Biehl | Multispeca tool for multispectral-hyperspectral image data analysis[END_REF]. structure. Fig. 1b shows an example of how the fluid structure looks like.

The arrangement of the grains is easily recognizable and the sampling of each of them during the images analysis can be automatically or manually done directly from black and white images without needing more complex adjustments. Once the outline of the grains is defined from the image, the grain surface comes by counting the number of pixels of which it is composed.

We describe each grain by the equivalent diameter of a circle with the same area as the measured surface. The final grains size distributions are reported in Fig. 2 where data are normalized to a probability density function and in Table 1 where the characteristics of the distributions are listed.

Setup

All the spheres are released in a 50 cm deep (H) and 10 cm wide (D) cylindrical Plexiglas vessel, previously filled with approximately 3.5 liters of fluid (Fig. 1a). In order to reduce optical distortions during the recording of images illuminated laterally by a laser sheet, the vessel was placed in a rectangular tank filled with water. A vertical cross-section of the tank is then illuminated by a 532 nm laser sheet, and images are recorded using a video camera (Canon Legria HF S21 1080p) with frame rate of 25 frames/s.

The position of the center of the sphere on each frame was determined by plotting the spatio-temporal evolution of a vertical pixel line centered on the sphere, or by using the blob analysis method in Matlab [START_REF]Blob analysis in computer vision system tm matlab[END_REF] (fig. 4 and5a).

Both methods allow the derivation of the local vertical velocity of the sphere through time. And the blob analysis also provides the measurements of the horizontal position and velocity of the sphere.

We run experiments with spheres of different materials and sizes. Their properties (i.e. density and diameter) are summarized in Table 2. All the steel spheres (density of 7980 and 7970 kg/m 3 ) have been painted with black spray paint to avoid strong reflections, while the other spheres have been left

with their original surfaces. The use of spray paint could have increased the surface roughness that in turn can affect the velocity field in the fluid around the object as well as the shape of the yielded region [START_REF] Holenberg | Particle tracking velocimetry and particle image velocimetry study of the slow motion of rough and smooth solid spheres in a yield-stress fluid[END_REF]. However, we did not observe any discrepancies between the painted and the bare spheres.

Before each experiment, we stirred gently the fluid to remove any possible preferential path that might have been formed during previously runs. SAP, in fact, shows a quite strong hysteresis where the fluid structure does not recover completely its original state after it has been deformed. This aspect is observed when a second sphere is dropped from the same position as the previous one without stirring the fluid in between the two runs. For identical spheres, we have always measured higher velocities for those released later and falling through an already deformed path. Stirring the fluid between runs prevents this issue and allowed us to get reproducible velocities for identical spheres.

ρ s (kg/m 3 ) d s (mm) 14952

Experimental results

As first step of our investigation, it is useful to test the reproducibility and the accuracy of the experimental techniques and possible limits arising from the dimension of our setup. To do so, we measured the terminal ve-locities (v y ) of spheres settling in a Newtonian fluid at T c =15.2±0.

v Stokes = 2 9 R 2 (ρ s -ρ f luid ) η g, (1)
where R is the radius of the spherical object and g the acceleration due to gravity. In Fig. 3 we plot the velocities measured for spheres of same density and diameters in between 3 and 18 mm. Larger spheres are subjected to wall effects and in this case the sphere terminal velocity can be corrected by a coefficient that writes as

K -1 1 = 1 + α 1 d s D + α 2 d s D 3 + α 3 d s D 5 , (2)
where α 1 = -2.104443, α 2 = 2.08877, α 3 = -0.94813 and D is the tube diameter [START_REF] Happel | Low Reynolds number hydrodynamics: with special applications to particulate media[END_REF]. Experiments carried out in syrup show that the measured settling velocities are in good agreement with relative corrected Stokes velocities, that is v y /(v Stokes K -1 1 ) = 1, ensuring that our measurement system is accurate (Fig. 3).

We now move to the case where the fluid through which spheres fall is one SAP sample described in the previous section. In Fig. 4 we report a selection of falling profiles in gel C (grains mean diameter d g =3.3 mm) for four falling spheres with same density but different diameters. The same type of plot is also showed in Fig. 5a where this time we track spheres having same size and no-motion (downward-pointing triangles). Regimes (1), ( 2), ( 4) and [START_REF] Seth | A micromechanical model to predict the flow of soft particle glasses[END_REF] in the text, respectively. The insert shows the distance as function of time for regime (4) with a semi-log scale. The color line is the fit y(x)=a log(b x), where a=98.34 and b=0.66. The correlation coefficient is R=0.980.

but different densities settling in gel D (d g =4.8 mm). In both figures, the starting position along y-axis, h=0, does not refer to the top fluid surface but to the upper limit of the camera view at which movies were recorded (Fig. 1a). This is a consequence of the fact that all movies have been taken for a magnified portion of tube h=300 mm, starting from 70 mm under the surface. Camera position and therefore size of the view are kept constant for each run done in the same fluid. For each gel, depending on the diameter and density of the sphere, we distinguish the following five regimes (Fig. 4 and Fig. 5). (1) A linear regime: rapid and linear fall in which spheres reach a constant terminal velocity; [START_REF] Bonn | Yield stress materials in soft condensed matter[END_REF] an irregular regime superimposed to a linear one where spheres never stop during their way down (v y >0) but their local velocity fluctuates around a constant value; (3) an intermittent regime: another irregular regime but here periods of no-motion, where the sphere's vertical velocity goes to zero, and periods of irregular falls follow one another; (4) a logarithmic regime where the sphere's speed progressively decreases (inset in Fig. 4); and (5) a no-motion regime in which spheres do not move at all from their initial positions.

Regimes (1) and ( 5) are characteristic of a yield stress fluid [START_REF] Tabuteau | Drag force on a sphere in steady motion through a yield-stress fluid[END_REF][START_REF] Fazilati | Quasi-periodic and irregular motion of a solid sphere falling through a thixotropic yieldstress fluid[END_REF][START_REF] Ferroir | Motion of a solid object through a pasty (thixotropic) fluid[END_REF] and one can interpret them based on the state of the medium that surrounds the object. In the linear regime (1) the fluid around the sphere has a liquidlike behaviour that allows spheres to go down, whereas in regime (5) local stresses are below the yield stress (σ Y ) and the material around the sphere remains in the solid state, leading to the no motion of the object. Here the stress field is generated by the weight of the sphere minus the buoyancy and the no-motion regime is always observed for smallest diameters d s (Fig. 4) or density contrast ∆ρ (Fig. 5a, as in [START_REF] Tabuteau | Drag force on a sphere in steady motion through a yield-stress fluid[END_REF]). We have had no-motion whether the sphere was placed on the surface of the fluid column, or it was gently pushed within it and released deeper in the fluid. In both cases, we left the spheres in that position for days (up to a week for one run) and they never moved away from their starting positions.

Although the distinction between the two end-member cases (1) and ( 5) is easy to establish, to discern quantitatively between the first steady-state motion regime (1) and the more irregular ones is more delicate. We summarize how we define the different motion regimes in Table 3. From the falling profiles one can notice that, although sometimes it is still possible to identify an almost constant falling speed superimposed on fluctuations (circle symbols both in Fig. 4 and Fig. 5a), in some other cases, spheres seem to slow down logarithmically during their way to the bottom (square symbols in Fig. 4), making impossible the determination of a constant terminal velocity. We report this aspect in Fig. 6a by plotting the mean squared variations (MSE) obtained by fitting the falling profiles with a linear or a semilog fit as function of spheres buoyancy. The logarithmic regime ( 4) is defined for MSE linear > MSE log , that is when falls are better represented with a logarithmic fit. Hence, this regime describes a continuous decrease in sphere's speed which never reaches the steady-state. Although in the case of Carbopol [START_REF] Tabuteau | Drag force on a sphere in steady motion through a yield-stress fluid[END_REF] velocity decreases following a power law, in most of our cases it best fits a logarithmic curve. Although from Fig. 6a is evident that the logarithmic behaviour appears more often at lower buoyancies, the distinction between linear and non-linear regime is not straightforward. In a few cases, indeed, at same buoyancy the two behaviours coexist.

Once the separation between stoppage cases (logarithmic regime (4) and no-motion ( 5)) and the other cases is established, one can identify and recognize diversities between the steady-state motion (regime (1)) and the more irregular and chaotic regimes (2) and (3). Depending on the buoyancy of the object and hence on its velocity, the irregularities show up as fluctuations in the sphere local vertical velocity and in a progressive increase of its indicates a better linear fit whereas ♦ a better semilog fit (logarithmic regime (4), in the text). (b) mean falling velocities (v y ) for those falling profiles that have been linearly fitted. The error for v y is 0.25 mm/s, smaller than the symbols. We do not report the standard deviation for v y here (whereas it is in (c) and (d)) as being also linked to the physical fluctuation of the sphere, for small values of velocity it would make the plot unreadable. In (c) and (d) we separate those experiments that show a more irregular and chaotic motion ( ) from those that do not show any irregularities during the fall ( ) by plotting std(v y )/v y . The dashed line is for std(v y )/v y = 0.15, from the Newtonian case.

In (e) and (f) we separate the intermittent regime (3) (×) from the other irregular regime (2) ( ) based on the amount of time t s the sphere remains at v y err(v y ). t 0 is the total time the sphere takes to descend h. The dashed line indicates t s /t 0 =5. horizontal motion (Fig. 5b). Regimes (2) and (3) show some similarities: in both of them, spheres have a long term constant falling velocity on which is superimposed a more complex and irregular pattern of fluctuations. The latter are characterized by accelerations and subsequent decelerations in a short period of time and, in some cases, they can result in a succession of stops of the sphere. In the most drastic cases, we observe the complete arrest of the sphere after it has passed through most of the fluid (Fig. 5a).

To determine the boundaries between the irregular regimes ( 2) and [START_REF] Coussot | Rheometry of pastes, suspensions, and granular materials: applications in industry and environment[END_REF] and regime (1), one can look, as a first approximation, at the ratio between the standard deviation (std) of the local falling velocity and its mean value.

Variations of the local velocity for spheres in irregular regimes are much larger than the experimental uncertainty (Fig. 6c-d). The latter is due to the precision of the technique we used to detect the sphere position in time (hereafter err(v y )). It is estimated to be 0.25 mm/s from the syrup-case measurements and from local velocity measurements in the no motion regime. In the Newtonian case, std(v y )/v y always remains smaller than 0.15. We define the linear regime (1) using this value, that is for std(v y )/v y < 0.15 (dashed line in Fig. 6c-d). This limit represents therefore a boundary between the steady-state motion and what shows a more spread distribution of velocities.

Once the linear regime is bounded in this way, we distinguish between the irregular regime (2) and the intermittent regime (3) as follows. In the intermittent case the sphere remains during a certain cumulative amount of time, t s , at v y err(v y ). We then compare t s to the total time, t 0 , the sphere takes to descend 30 cm. The intermittent regime is defined as t s /t 0 ≥ 5%.

The remaining irregular regime (2) is therefore defined as t s /t 0 < 5% and Regimes (1) (2) (3) (4) ( 5)

v y = 0 ± err(v y ) x MSE lin >MSE log x MSE lin <MSE log & std(v y )/v y < 0.15 x MSE lin <MSE log & std(v y )/v y 0.15 x x MSE lin <MSE log & t s /t 0 ≥ 5%
x Table 3: Summary of the observed motion regimes and the way they are classified std(v y )/v y > 0.15 (Fig. 6e-f).

Fluctuations, accelerations and subsequent decelerations and the resulting general complex set of behaviours of the moving sphere, characterize the experiments both in the direction of fall and in the one perpendicular to it.

In fact, the irregular motion we have observed is not only present and visible in the vertical direction but it also affects the horizontal motion of spheres in a such a way that, along this direction, they collect a certain amount of traveled distance (L x ) with respect to their initial position, depending on their diameter (Fig. 7a), on the gel grains size (Fig. 7b), and on the vertical falling distance h (Fig. 8a). Fig. 8b shows that for each run, L x is roughly proportional to h. Fig. 7 shows that small spheres change frequently their horizontal position, oscillating continuously around the vertical axis, and traveling for larger distances L x than the large spheres. If we normalize L x by h and d s by the diameter of the gel grains (d g ), we achieve a good collapse for all our data along a line with slope -2. That is

L x /h = C exp (d s /d g ) -2 . ( 3
)
where C exp = 0.49±0.26. In this way the total travelled distance in horizontal direction turns to be a good candidate to underline the strong interaction of the sphere and the gel grains that make up the fluid structure when those two have similar sizes.

Discussion

Determination of the SAP effective rheology

The fall of a spherical object within certain boundary conditions is described by the combination of a proper set of equations of motion and the constitutive relation of the fluid through which the sphere is falling. For yield stress fluids a commonly used rheological model that relates stress and strain in a non-linear way is the Herschel-Bulkley model [START_REF] Bonn | Yield stress materials in soft condensed matter[END_REF][START_REF] Balmforth | Yielding to stress: recent developments in viscoplastic fluid mechanics[END_REF][START_REF] Coussot | Yield stress fluid flows: A review of experimental data[END_REF]:

     σ = σ Y + K v γn if σ > σ Y γ = 0 if σ ≤ σ Y ( 4
)
where σ is the stress, σ Y the yield stress, γ the shear rate, K v the consistency and n the shear-thinning index. Considering a fluid which flows following such model, for small Reynolds numbers (Re 1), the problem of a settling particle needs two parameters to be characterized [START_REF] Tabuteau | Drag force on a sphere in steady motion through a yield-stress fluid[END_REF][START_REF] Ansley | Motion of spherical particles in a bingham plastic[END_REF]. They are the

Bingham number

Bi = σ Y K v (v y /d s ) n , (5)
that compares yield stress with viscous stresses, and the yield number

Y = 3σ Y gd s ∆ρ (6)
that compares the yield stress with the buoyancy stress.

As described in section 1, most of the experimental work done to study this problem is carried out by using Carbopol. The main difference between experiments done with Carbopol and the present work regards the typical size of the fluid structure compared to the dimension of the object that is moving through it. Carbopol corresponds to the case where the fluid structure (grains) is much smaller than the spheres. In this case, only regimes

(1), ( 4) and (5) were reported (e.g. [START_REF] Tabuteau | Drag force on a sphere in steady motion through a yield-stress fluid[END_REF]) while neither the irregular regime (our regime (2)) nor the intermittent regime (3) were reported; and the fall of the sphere is always reported to be free of any sort of chaotic motion [START_REF] Tabuteau | Drag force on a sphere in steady motion through a yield-stress fluid[END_REF][START_REF] Putz | Settling of an isolated spherical particle in a yield stress shear thinning fluid[END_REF][START_REF] Holenberg | Particle tracking velocimetry and particle image velocimetry study of the slow motion of rough and smooth solid spheres in a yield-stress fluid[END_REF].

In regime [START_REF] Makse | Statistical mechanics of jammed matter[END_REF], that is when a constant terminal velocity can be determined, Tabuteau et al. [START_REF] Tabuteau | Drag force on a sphere in steady motion through a yield-stress fluid[END_REF] showed that their data are well-fitted by the following equation

1 Y = 7 + 8.52 Bi (7)
with n=0.5, and that there is a critical value of Y , Y c =0.145, above which there is no motion.

Taking into account the differences between SAP and Carbopol, one can use eq.( 7), defined for d s d g , and check whether it applies for a fluid with a much larger structure. Thus, assuming a framework in which eq.( 7) is valid also in our case, we can rewrite it as

gd s ∆ρ 3 = 7σ Y (ef f ) + 8.52K v(ef f ) v y d s n , (8)
where σ Y (ef f ) and K v(ef f ) represent the effective yield stress and consistency of the SAP through which the sphere is falling. In eq.( 8), v y is the average falling velocity defined for those experiments which do not show a decrease in speed, that is those in regime ( 1), ( 2) and [START_REF] Coussot | Rheometry of pastes, suspensions, and granular materials: applications in industry and environment[END_REF]. By plotting the buoyancy stress (gd s ∆ρ) as function of the effective strain rate, i.e. (v y /d s ) n , Fig. 9 shows for SAP the same linear relationship showed for Caropol [START_REF] Tabuteau | Drag force on a sphere in steady motion through a yield-stress fluid[END_REF].

Therefore the slope of the linear fit is related to an effective consistency through

K v(ef f ) = gd s ∆ρ 25.56 -0.86 σ Y (ef f ) v y d s -n . (9)
and, similarly, an effective yield stress is related to the intercept at zero effective strain rate through

σ Y (ef f ) = gd s ∆ρ 21 . (10)
In this way, by using the falling sphere system as reometer we can define an effective rheology, which has to have same form of eq.( 7), for each of our gel sample. Fits in Fig. 9 depend on the shear thinning exponent n used to calculate the effective shear rate. Minimum mean squared values for linear fit are found for n in between 0.5 and 0.7 (Fig. 11). Unless otherwise indicated, in what follows we will use n=0.5 since it is with this value that eq. ( 7) is defined. Interesting, experiments in both regime (1), ( 2) and ( 3) can be linear fitted by the same line. Moreover, effective yield stress and consistency increase with the gel grains size (Fig. 10).

For experiments involving large spheres, it is opportune to check whether wall corrections are needed or not. For sufficiently small sphere-to-tube diameter ratio, d s /D, we do expect no reductions of sphere's terminal velocity.

This means that the measured velocity, v y , is equal to the ideal terminal velocity of a sphere falling in an unbounded medium, v ∞ . On the other hand, if d s /D exceeds a critical value, (d s /D) crit , velocity decreases leading to v y /v ∞ <1. The critical diameter ratio for a yield stress fluid writes as [START_REF] Atapattu | Wall effect for spheres falling at small Reynolds number in a viscoplastic medium[END_REF] (

d s /D) crit = 0.055 + 1.114 Y if Y Y c . (11)
In our case, d s /D is smaller than (d s /D) crit as long as d s /D <0.1. The influence of walls causes a reduction in speed such that [START_REF] Atapattu | Wall effect for spheres falling at small Reynolds number in a viscoplastic medium[END_REF]

v y v ∞ = 1 -1.7 [(d s /D) -(d s /D) crit ] if d s /D > (d s /D) crit . (12)
For 0.1 d s /D 0.17 experiments are above critical conditions and, according to eq. ( 12), wall effects lead to 0.8< v y /v ∞ <1, that is a maximum decrease of v y /v ∞ of less than 20%.

Considering only the experiments with v y /v ∞ > 0.8 and ignoring the others, we recalculate the effective rheology in the same way described above Tabuteau et al. [START_REF] Tabuteau | Drag force on a sphere in steady motion through a yield-stress fluid[END_REF] approach, v y is the average falling velocity for those experiments where a steady decrease of the speed has not been observed (i.e. we do not consider regime (4)) and n=0.5 the power-law index. Colours represent the fluids (same as Fig. 7)

while symbols the regimes (circles for regime (1), squares for regime (2) and crosses for regime (3)). Empty symbols for steel sphere; symbol with a point inside for ρ s =7799 kg/m 3 ; with a plus sign for ρ s =14952 kg/m 3 and filled symbols for ρ s =2200 kg/m 3 . Solid lines are the linear data fits. (Fig. 12a). The new values of σ Y (ef f ) and K v(ef f ) do not differ much from the previous ones (Fig. 12b); which means that our set up is big enough.

Once we got the yield stress and the consistency of each sample, we plot Y -1 as function of Bi -1 in Fig. 13. For values of Bi -1 >0 we achieve a good collapse along eq.( 7) for experimental data that belong to regimes (1), ( 2) and (3). In this way, even considering the irregular regime (2) and the intermittent [START_REF] Coussot | Rheometry of pastes, suspensions, and granular materials: applications in industry and environment[END_REF] where the interaction between moving objects and fluid structure results in highly irregular falls of objects themselves, an effective rheology following the Hershel-Bulkley model with a classical shear-thinning index of 0.5 well represents the majority of our data. As we highlighted at the beginning of this section, equation ( 7) was derived for Carbopol when the size of the moving object is much larger than the size of the particles than make up the fluid structure. Inside this assumption, the critical condition is well predicted in our case too. In gel E, for example, we had no-motion for a large glass sphere (d s /d g >2.5 and ∆ρ=1200 kg/m 3 ) for Y =0.478 (Fig. 14), greater than Y c =0.145. It is only when the diameter of the sphere becomes comparable to the diameter of the gel particles that Y c moves to lower values and the stoppage regime (5) cannot be anticipated any more from eq.( 7). For the same gel E, a smaller steel sphere (d s /d g <2 and ∆ρ=6970 kg/m 3 ) is in stoppage for Y =0.13< Y c . We attribute this to the interaction between spheres and the gel structure: when d s →d g , the spheres "see" obstacles on their way and the effective rheology breaks down. open-symbols for ρ s =7970 kg/m 3 , point-filled for ρ s =7799 kg/m 3 , asterisk-filled for ρ s =3227 kg/m 3 and plus-filled for ρ s =14952 kg/m 3 . We show the good collapse of our data along eq.( 7) (solid black line). The inset shows the data zoomed around the critical conditions in order to better display the entrapment conditions. Fluid elasticity have been also linked to a peculiar characteristic of the fluid flow around the sphere, namely its asymmetry. The rise of a fore-arc asymmetry and the formation of a negative wake at the sphere's tail has been detected experimentally in Carbopol [e.g. 17, 38] and associated with fluid elasticity [START_REF] Fraggedakis | Yielding the yieldstress analysis: a study focused on the effects of elasticity on the settling of a single spherical particle in simple yield-stress fluids[END_REF][START_REF] Putz | Settling of an isolated spherical particle in a yield stress shear thinning fluid[END_REF] or the emergence of thixotropy [START_REF] Dinkgreve | Carbopol: From a simple to a thixotropic yield stress fluid[END_REF]. In SAP, we have observed a fore-arc asymmetry in the velocity field collected from PIV measurements. Indeed, Fig. 15 shows that a negative wake in the region at the sphere back can be recognized at least for the linear motion regime [START_REF] Makse | Statistical mechanics of jammed matter[END_REF].

Entrapment conditions

Although to discern the cause that breaks the fore-arc symmetry in SAP is not straightforward, its formation is more probably due to elasticity. In fact, between SAP grains there are not polymer chains that form system-spanning networks and the fluid, during a single experiment, does not undergo aging.

More details about will be given in a coming work where the fluid velocity field will be properly analyzed and explained.

Moreover, a further analysis on how the properties of a single elastic grain could influence σ Y (ef f ) and K v(ef f ) might be relevant to evaluate the role that its viscoelasticity plays for the bulk rheology. However, the extreme softness of the material makes it very difficult. After several attempts to measure forces and penetration depth by using a microindenter with a minimum normal force of 20 mN, we have not been able to get any reliable measurements since slip or rupture of the grain always occurred even for such small values of normal force.

Conclusion

In this work we have experimentally investigated the influence that the size of particles which build up the structure of an yield stress fluid has on the dynamics of a spherical intruder falling under gravity through it. As the typical size of the gel grains was varied between 1 and 6 mm and the sphere diameter was varied between 3 and 25 mm, there exists a range where they become comparable.

We have been able to expand from three (in Carbopol) to five (in SAP) the observed regimes of motion of the sphere: (1) a linear regime where the sphere has a steady-state motion; (2) an irregular regime superimposed to a linear one where spheres never stop during their way down but their velocity varies; (3) an intermittent regime in which periods of no-motion and periods of irregular falls follow one another; (4) a logarithmic regime where the sphere velocity progressively decreases; (5) no-motion. Beside the classical steadystate motion (1) and stoppage regimes [START_REF] Seth | A micromechanical model to predict the flow of soft particle glasses[END_REF], typical of viscoplastic fluids, the interaction between moving objects and the fluid structure results then in two additional regimes ((2) and ( 3)) where the motion becomes more chaotic.

The same strong interaction between the sphere and the gel grains is also observed from the sphere total travelled distance in the horizontal direction. Moreover, we have shown that the effective rheology of SAP follows the Herschel-Bulkley model with a power index n=0.50, making this simple setup a good rheometer for this type of fluids.

We also found that both yield stress and consistency of the gel increase with the gel grain size, suggesting that the effective rheology, due to the material heterogeneous properties, depends on the size of the particles which constitute the fluid structure. However, it is remarkable that the overall sphere motion predicted by the Herschel-Bulkley model begins to "see" the texture only for sphere diameters less than twice the SAP grain size.

Figure 1 :

 1 Figure 1: (a) Experimental set up. The cylindrical vessel containing the fluid is 50 cm high (H) and 10 cm wide (D). Videos are recorded in the tube section between the two black dotted lines (h). (b) Black and white image of the fluid structure illuminated by a laser sheet. The gel grains are in dark grey and the thin water film around them is in light grey.

Figure 2 :

 2 Figure 2: Grains size distributions for the five gels used in the experiments. Histogram representation of grains equivalent diameters with a normal distribution curve on top for gel A (a), gel B (b), gel C (c), gel D (d) and gel E (e). (f) All gels together.

Figure 3 :

 3 Figure 3: Measured terminal velocities, v y , in syrup (circles) and Stokes velocities, v Stokes , (crosses and dashed line) as function of spheres squared radius, R 2 s . Uncertainties for Stokes velocities come from the uncertainties in determining the syrup viscosity whereas the error for v y is 0.25 mm/s, smaller than the symbols. Star symbols are the ideal velocities considering the wall correction, eq. (2).

Figure 4 :

 4 Figure 4: Depth of spheres settling through gel C as function of time. Spheres have constant density (∆ρ=6980±10 kg/m 3 ). Symbols refer to the diameter of the sphere: d s =10 mm ( ), d s =8 mm (•), d s =7 mm ( ), d s =4 mm ( ). Varying sphere's diameter, one can recognize the following motion regimes: linear regime (upward-pointing triangles), irregular motion superimposed to a linear trend (circles), logarithmic regime (squares)

Figure 5 :

 5 Figure 5: (a) Depth of spheres settling through gel D as function of time. Here spheres have constant diameter (d s =6 mm) but different densities. As in Fig.4, we report the case of linear regime (upward-pointing triangles), no-motion (downward-pointing triangles) and, in addition, the intermittent regime (circles). (b) Sequence of snapshots every 100 frames showing the irregular motion in gel D for a sphere with d s =7 mm and ρ s =7970 kg/m 3 .

Figure 6 :

 6 Figure 6: Determination of the different regimes of motion. In (a), (c) and (e) data refer to experiments carried out in gel E. In (b), (d) and (f) data are for all fluids and colours indicate the grain size as in Fig.2. (a) best fit (M SE linear /M SE log ) of the falling profiles.

Figure 7 :

 7 Figure 7: Total distance traveled in the horizontal direction (L x ). (a) L x as function of the sphere diameter d s . Dashed line indicates the uncertainty in position. (b) Same plot with y-axis normalized by the vertical falling distance h and x-axis by the gel grain size d g . Colours indicate the different fluids as indicated in the legend. Symbols refer to the different motion regimes:for the linear regime (1); for an irregular regime superimposed to a linear one (2); × for the intermittent regime (3) and ♦ for the logarithmic regime[START_REF] Barnes | A review of the slip (wall depletion) of polymer solutions, emulsions and particle suspensions in viscometers: its cause, character, and cure[END_REF]. Empty symbols are for steel spheres while the symbol with a point inside refers to ρ s =7799 kg/m 3 and those filled with a plus sign are for ρ s =14952 kg/m 3 .

Figure 8 :

 8 Figure 8: Distance traveled in the horizontal direction in gel C. (a) Distance traveled in different 15-cm-vertical-sections of h, L x [Sect], as function of spheres buoyancy. Sections are Sect A =[h=0 cm, h=15 cm]; Sect B =[h=15 cm, h=30 cm]; Sect C =[h=7.5 cm, h=22.5 cm]. (b) Ratio between the total distance traveled in 30 cm, L x , and the distance travelled in 15 cm, L x [Sect].

Figure 9 :

 9 Figure 9: Sphere buoyancy as function of effective shear rate, (v y /d s ) n . Following

Figure 10 :

 10 Figure 10: Effective yield stress (σ Y ) and effective consistency (K v ) as function of gel grains size.

Figure 11 :

 11 Figure 11: Power law index (n) that best linearly fits data in Fig.9 for all gel samples. yaxis represents the mean squared errors (MSE) of the linear fit normalized by the maximum MSE found. The best range of n is in between 0.50 and 0.70.

Figure 12 :

 12 Figure 12: Influence of the velocity reduction due to wall effects on the effective rheology.

n=0. 5 .

 5 (a) The new fit that considers only experiments with v y /v ∞ > 0.8 (solid lines) is compared with the previous fit from Fig.9 (dashed lines). Colours and symbols here are those described in Fig.9. (b) Effective yield stress and consistency as function of gel grains size from the new fit (blue symbols). Red symbols are those in Fig.10.

Fig. 13 also

 13 Fig.13 also shows that for Bi -1 =0 spheres get entrapped also for Y -1 > Y -1 c . From the experimental work of Tabuteau et al. [11], the critical value of the yield number (Y c ) above which spheres do not move in a simple yield stress fluid is 0.145. This value is very close to 0.143, predicted in previous works of Beris et al.[10] and Blackery and Mitsoulis[START_REF] Blackery | Creeping motion of a sphere in tubes filled with a bingham plastic material[END_REF]. On the other hand, a more recent study[START_REF] Fraggedakis | Yielding the yieldstress analysis: a study focused on the effects of elasticity on the settling of a single spherical particle in simple yield-stress fluids[END_REF] points out how the critical value may increase as function of the Deborah number if elastic effects are taken into account. In our case, however, we observe sphere arrest also for Y < Y c . For example, we observed arrest even for Y = 0.042 ≈ Y c /3 (4 mm steel sphere).

Figure 13 :

 13 Figure 13: 1/Y as function of 1/Bi. Colours are for the different fluids. The dashed line indicates Y -1 c =0.145 -1 [11]. Symbols are filled in a different way to distinguish between spheres of different materials: full-filled symbols are for light spheres (ρ s =2200 kg/m 3 ),

Figure 14 :

 14 Figure 14: Range of validity of the effective rheology. Colours are for the different fluids. Symbols are those used in Fig.13. The dashed line indicates Y -1 c =0.145 -1 [11]. For d s /d g approaching the value of 2, we have no motions (triangles), even for yield numbers lower than the critical one commonly used.

Figure 15 :

 15 Figure 15: Fluid flow from PIV measurements. We display the vertical component of the right half part of instant velocity field (v z ) for a sphere (d s =19 mm) falling in the linear regime (1) through SAP.

Table 1 :

 1 Results from the imaging analysis of our five gel samples. Here we report the

	The sample was stirred between

Table 2 :

 2 Spheres density ρ s and diameter d s .

	6.00

  3 • C. The fluid is glucose syrup (Glucor 60/80) with ρ=1394.66 kg/m 3 and η=12.3±0.7

	Pa.s, density and viscosity at T c , respectively. Spherical objects falling in a
	Newtonian fluid at low Reynolds number reach a terminal velocity given by
	the Stokes velocity

Acknowledgments N.R. Sgreva is supported by the Initial Training Network No 642029-ITN CREEP, an Horizon 2020 -Marie Sklodowska-Curie Action.

Grains rearrangement and elasticity

The emergence of a yield stress in SAP can be related to the jamming of soft grains. Gel grains act as athermal elastic particles that interact by gliding and compressing each other. The collective behaviour of a fluid made up by these elastic grains or particles can shows elasto-viscoplastic responses [START_REF] Seth | A micromechanical model to predict the flow of soft particle glasses[END_REF][START_REF] Seth | Influence of short-range forces on wall-slip in microgel pastes[END_REF].

Therefore, one can expect that also the observed irregular motion regimes can be linked to the jamming nature of the fluid. Within this framework, going from the linear to the intermittent regime by decreasing the sphere size makes the elasticity and the yield stress seen by the sphere no longer the collective one of the bulk material, but rather the one of few gel grains.

Once this occurs, to allow the continuous fall of the sphere, grains in front of it must undergo rearrangement. The mechanics of the flow appears then similar to the one of soft glassy materials. In this case the flow occurs through a sequence of elastic (reversible) deformations and local plastic (irreversible) rearrangement events which induce a long ranged relaxation of stress [START_REF] Goyon | Spatial cooperativity in soft glassy flows[END_REF][START_REF] Mansard | Local and non local rheology of concentrated particles[END_REF][START_REF] Barrat | Heterogeneities in amorphous systems under shear, Dynamical heterogeneities in glasses[END_REF]. Grains rearrangements will be all the more common the higher the local shear stresses are (e.g. in the case of a larger sphere). Since the number of plastic events for unit of time controls fluid flow and relaxation [START_REF] Goyon | Spatial cooperativity in soft glassy flows[END_REF], a sufficient high number of them leads the system to be more fluid. The resulting sphere motion would be in the linear regime (1); whereas irregular and intermittent regimes might be the result of a less fluid system where the local rheology seen by the sphere becomes more relevant. In this way, the observed deformation induced by the sphere falling in the linear regime is mainly the elastic interaction between the grains, as it is for flow of soft particle glasses [START_REF] Seth | A micromechanical model to predict the flow of soft particle glasses[END_REF].