Fabrice Barbe
email: fabrice.barbe@insa-rouen.fr

Ivano Benedetti
email: ivano.benedetti@unipa.it

Vincenzo Gulizzi
email: vgulizzi@lbl.govd

Mathieu Calvat
email: mathieu.calvat@insa-rouen.fr

Clément Keller
email: clement.keller@insa-rouen.fr

Elucidating the effect of bimodal grain size distribution on plasticity and fracture behavior of polycrystalline materials

. The analysis of the computational results, compared to the experimentally characterised tensile properties of a bimodal 316L stainless steel, suggests that the elasto-plastic interactions taking place prior to micro-cracking may play an important role in the mechanics of fracture of this steel.

INTRODUCTION

As characterized in the seminal works of Hall and Petch [START_REF] Hall | The deformation and ageing of mild steel: III discussion of results[END_REF][START_REF] Petch | The Cleavage Strength of Polycrystals[END_REF], the yield stress and strength of a metal alloy is generally inversely proportional to the square root of the mean Grain Size (GS). This observation can be made for most metal alloys and for grain sizes larger than 1µm. Reducing the grain size however also results in the degradation of the ductility as could be characterised in [START_REF] Flipon | Tensile properties of spark plasma sintered AISI 316L stainless steel with unimodal and bimodal grain size distributions[END_REF][START_REF] Meyers | Mechanical properties of nanocrystalline materials[END_REF][START_REF] Pineau | Failure of metals III: Fracture and fatigue of nanostructured metallic materials[END_REF]. This has motivated the development of different strategies to overpass this strength-ductility dilemna by introducing in the material some microstructural characteristics acting as local stress relaxers. Tailoring the grain size distribution is one of these strategies [START_REF] Ramtani | A revisited generalized self-consistent polycrystal model following an incremental small strain formulation and including grain-size distribution effect[END_REF][START_REF] Berbenni | Impact of the grain size distribution on the yield stress of heterogeneous materials[END_REF][START_REF] Zhu | A statistical model for predicting the mechanical properties of nanostructured metals with bimodal grain size distribution[END_REF] and one remarkable case corresponds to the bimodal GS distribution, where Coarse Grains (CG) are present in a matrix of Fine Grains (FG).

One of the most efficient means to elaborate these kind of architected materials is the sintering of metal powders [START_REF] Flipon | Tensile properties of spark plasma sintered AISI 316L stainless steel with unimodal and bimodal grain size distributions[END_REF][START_REF] Srinivasarao | Synthesis of high-strength bimodally grained iron by mechanical alloying and spark plasma sintering[END_REF][START_REF] Dirras | Microstructure and mechanical characteristics of bulk polycrystalline Ni consolidated from blends of powders with different particle size[END_REF][START_REF] Zhang | Improvement of mechanical properties in SUS304L steel through the control of bimodal microstructure characteristics[END_REF][START_REF] Tingaud | Investigation of deformation micro-mechanisms in nickel consolidated from a bimodal powder by spark plasma sintering[END_REF]. It indeed enables to control the main metrics of a bimodal microstructure, viz. the mean GS of the CG population, the mean GS of the FG population, and the volume fraction of CG. In the experimental work [START_REF] Flipon | Tensile properties of spark plasma sintered AISI 316L stainless steel with unimodal and bimodal grain size distributions[END_REF][START_REF] Flipon | Elaboration of austenitic stainless steel with bimodal grain size distribution and investigation of their mechanical behavior[END_REF], different associations of these three metrics in a 316L austenitic stainless steel have been considered from the spark plasma sintering to the tensile fracture properties. Whereas unimodal cases confirmed the reduction of the elongation at break as GS was decreased, bimodal cases have shown that a good combination of strength and ductility could be obtained for a GS ratio larger than 10 and a FG size of the order of few µm.

Different modelling approaches have been proposed in the literature to predict the mechanical properties according to the GS distribution and to study the interactions taking place between FG and CG. Some of these proposals rely on mean field hypotheses of self-consistent homogenization approaches [START_REF] Pipard | A new mean field micromechanical approach to capture grain size effects[END_REF][START_REF] Ramtani | A bimodal bulk ultra-fine-grained nickel: Experimental and micromechanical investigations[END_REF][START_REF] Wang | A 3d crystal plasticity model of monotonic and cyclic simple shear deformation for commercial-purity polycrystalline ti with a harmonic structure[END_REF] or on the idealization of a GS population as a homogeneous phase [START_REF] Magee | Representation of a microstructure with bimodal grain size distribution through crystal plasticity and cohesive interface modeling[END_REF], or are restricted to the early stages of plastic deformation [START_REF] Lavergne | Effects of grain size distribution and stress heterogeneity on yield stress of polycrystals: A numerical approach[END_REF]. In the study [START_REF] Flipon | A full-field crystal-plasticity analysis of bimodal polycrystals[END_REF], a comprehensive presentation of a full field crystal plasticity modelling applied to bimodal polycrystals has been provided. Yet, although this modelling approach provides a very fine description of the elastoplastic mechanisms originally taking place in a bimodal polycrystal, its applicability is restricted to crystal plasticity. The present study aims at constructing the bases of a full field modelling approach accounting for both crystal plasticity and microcracking. For such preliminary steps, the micro-mechanics of a model bimodal polycrystal is studied considering first its crystal plasticity, by the finite element method, and then its intergranular fracture behaviour, by means of a boundary element framework [START_REF] Sfantos | Multi-scale boundary element modelling of material degradation and fracture[END_REF][START_REF] Benedetti | A three-dimensional cohesive-frictional grain-boundary micromechanical model for intergranular degradation and failure in polycrystalline materials[END_REF].

MODELLING

Material and Microstructure

The material under concern is a 316L stainless steel elaborated by spark plasma sintering. Its microstructure and tensile properties have been presented in [START_REF] Flipon | Tensile properties of spark plasma sintered AISI 316L stainless steel with unimodal and bimodal grain size distributions[END_REF] for different combinations of grain sizes. As illustrated by EBSD scans presented therein, CG are commonly observed to be embedded in the matrix of FG as isolated inclusions. Previous numerical analyses [START_REF] Flipon | A full-field crystal-plasticity analysis of bimodal polycrystals[END_REF] have shown that in such a configuration of isolated CG, a domain containing at least 9 CG should be considered to obtain quantitatively representative macroscopic properties; yet a 2D domain containing a single CG could provide qualitativaly representative properties, especially for the analyses of local CG-FG interactions. This is the retained configuration of study in this paper, also motivated by the reduction in computational resources and the comfort of analysis it provides. Three model polycrystals are considered (fig. 1): a bimodal polycrystal with GS of 5µm and 0.5µm resp. for CG and FG families and a unimodal microstructure containing 256 grains and with a mean GS corresponding to either CG or FG. They result from Laguerre-Voronoi tessellations with specific settings on the nucleation and growth processes [START_REF] Flipon | A full-field crystal-plasticity analysis of bimodal polycrystals[END_REF][START_REF] Flipon | Modelling of polycrystals using well-controlled voronoi-type tessellations and its applications to micromechanical analyses[END_REF]. The Neper software was used to generate these microstructure and to obtain the finite elements meshes [START_REF] Neper | Neper: polycrystal generation and meshing[END_REF][START_REF] Quey | Large-scale 3-D random polycrystals for the finite element method: Generation, meshing and remeshing[END_REF].

Crystal plasticity

The Schmid-based crystal plasticity constitutive model of Méric and Cailletaud is used [START_REF] Méric | Single Crystal Modeling for Structural Calculations. Part 2: Finite Element Implementation[END_REF]. The plastic slip rate γs of the slip system s is deduced from the resolved shear stress τ s and the isotropic hardening r s through the Norton power law 1, where K and n control the rate sensitivity:

γs = |τ s | -r s K n signe(τ s ) (1)
The hardening component r s evolves according to the saturation law 2 expressed in terms of p the cumulated plastic strain:

r s = R 0 + Q N ∑ i=1 h si 1 -e -bp i with ṗi = | γi | (2)
The matrix h si describes the interactions between slip systems, Q and b affect the rate of saturation and R 0 is the critical resolved shear stress, above which slip is activated. Hall-Petch-type grain size dependence of the macroscopic yield stress is introduced directly on the scale of the slip system with the definition of R 0 given in eq. 3:

R 0 = r 0 + k √ d (3)
d is the grain diameter, k controls the dependence to the grain size and r 0 is the critical resolved shear stress in a grain having large dimensions.

All the material parameters have been identified on the basis of experimental tensile stress-strain curves of the considered 316L steel elaborated by spark plasma sintering, for different mean GS of unimodal polycrystals. It is worth noting that the Hall-Petch dependence introduced in the slip system scale results also in a Hall-Petch dependence of the macroscopic scale of a unimodal polycrystal [START_REF] Flipon | A full-field crystal-plasticity analysis of bimodal polycrystals[END_REF].

Elastoplastic behaviour

The polycrystals are tensile strained to 1%. The related axial stress and strain responses are presented on two scales: macroscopic behaviour in fig. 2 and field contours in fig. 3. The presence of the CG in the FG phase, although it occupies only 27% in volume, significantly decreases the macroscopic stress response. This can be explained from the role that the CG plays in the bimodal system: the CG is not only the base of highest strains localisation, it is responsible for a significant decrease of stress in the FG at its upper and lower neighbourhood: in this relaxation channel, the stresses in the FG population is much lower than in the unimodal FG polycrystal.

The effect of the presence of the CG within the FG matrix is further analysed by means of the mean responses per grain in fig. 4: the case of the bimodal distribution is compared to the case of the unimodal FG distribution. According to the grains stress-strain curves, a large dispersion of the tensile responses is observed in both cases, and FG grains do not seem to reach higher stresses in the bimodal case than in the unimodal case. The statistical distributions in the inserts show that the presence of the CG significantly broadens the dispersion of stresses, creating a bimodal distribution in the range of 60% to 160% from the macroscopic bimodal polycrystal stress. It also increases the number of FG grains at both very small strains (∼ 25% the mean polycrystal strain) and large strain (∼ 150% the mean polycrystal strain).

To quantify the mean effect of the CG on the FG behavior, we consider in fig 5 the mean behavior of each grain population in the bimodal polycrystal (resp. the CG and the FG population) divided by its corresponding mean properties in its counterpart unimodal polycrystal. This is done for each considered variable: axial strain and axial stress. One can then observe the particular behavior of the CG as it interplays with the FG grains: (i) a strain which FIGURE 3: Spatial distributions of strain (top) and stress (bottom) at 1% imposed strain: from left to right, Unimodal CG, Bimodal, Unimodal FG reaches 180% the mean strain reached in the unimodal CG polycrystal, but which however remains in a range that many grains reach in the unimodal CG polycrystal (cf. fig. 4), and (ii) a stress which abruptly increases to 120% and then slowly decreases to 110% the mean unimodal CG polycrystal. The coincident effect of this CG behavior on the FG population is a decrease to 62% and 78% the mean strain and stress (resp.) of the unimodal FG polycrystal. These increases or decreases of the mean behavior of a grain population in the bimodal system with respect to what it would be in a unimodal configuration results in the intermediate macroscopic stress-strain response given in fig. 2, as could be expected in any kind of "soft inclusion"-"hard matrix" composite. A peculiarity is however observed concerning the early stages of plasticity: three deformation regimes corresponding to (1) elasticity, (2) a first elastoplastic hardening regime and (3) a second elastoplastic hardening regime. These three stages can be distinguished on the curves of stress, strain and plasticity evolutions in each grain population (fig. 5; "plasticity" corresponds to the sum of the slips on the slip systems and thus quantifies a mean plastic slip activity): plasticity in the CG of the bimodal polycrystal starts at higher imposed strain than in the unimodal CG polycrystal and this moment of the load corresponds both to the CG stress relaxation and to the transition from the first to the second deformation regime; the second regime lasts until the FG stress saturates to 78% the mean stress in the FG unimodal polycrystal. It is further noticed that during the second regime, FG of the bimodal polycrystal display much more plasticity than in the unimodal configuration. It actually starts halfway during the second regime. FIGURE 4: Stress-strain curves of all the grains in the Unimodal FG polycrystal (left) and in the Bimodal polycrystal. The inserts present the distribution of the grains stresses and strains at the end of the load. Each grain response (subscript "GR") has been normalized by its corresponding mean polycrystal response (subscript "PX").

FIGURE 5:

The ratio of the mean CG (resp. FG) response in the bimodal polycrystal to the mean response of its unimodal counterpart is plotted against the imposed axial strain. As an example, the mean CG strain in the bimodal polycrystal is the same as in the CG unimodal polycrystal for a strain less than 0.0015 and becomes 1.8 times larger at the end of the test.

Intergranular fracture propagation

To help elucidate the behaviour of polycrystalline aggregates with a bimodal grain size distribution, the analyses are completed by an investigation into their fracture behaviour, employing a recently developed computational framework based on the use of a boundary integral representation of the polycrystalline micro-mechanics in conjunction with cohesive-frictional interfaces [START_REF] Sfantos | Multi-scale boundary element modelling of material degradation and fracture[END_REF][START_REF] Benedetti | A three-dimensional cohesive-frictional grain-boundary micromechanical model for intergranular degradation and failure in polycrystalline materials[END_REF]. The method, which is built on a Laguerre-Voronoi representation of the microstructure, allows to study the progressive degradation and micro-cracking of aggregates subjected to quasi-static loading. It employs the boundary element method [START_REF] Aliabadi | The boundary element method[END_REF] for discretising the so-called displacements boundary integral equations associated with, or collocated at, the mesh nodes of the discretised individual grains. The process, whose details are described in Refs. [START_REF] Benedetti | A three-dimensional cohesive-frictional grain-boundary micromechanical model for intergranular degradation and failure in polycrystalline materials[END_REF][START_REF] Gulizzi | An enhanced grain-boundary framework for computational homogenization and micro-cracking simulations of polycrystalline materials[END_REF], leads to an aggregate system of equations of the form

A • X I (X, d * ) = B • Y(λ ) 0 ( 4
)
where the matrix blocks A and B collect constant entries stemming from the boundary element numerical integration, the vector X collects unknown components of displacements and tractions associated with points lying on the surfaces of the individual grains, and it is then suitable for representing intergranular degradation and micro-cracking processes, the vector Y collects known components of displacements or tractions associated with points lying on the external surface of the aggregate, thus stemming from the boundary conditions, and the non-linear matrix block I (X, d * ) implements the intergranular cohesive-frictional laws, accounting also for the loading history through an irreversible damage field whose nodal values associated with the interface pairs are collected within in the vector d * . The methodology allows to study the microstructural evolution capturing the transition of individual interfaces from a pristine to a damaged and eventually cracked status. It has been successfully employed to model quasi-static micro-cracking [START_REF] Gulizzi | An enhanced grain-boundary framework for computational homogenization and micro-cracking simulations of polycrystalline materials[END_REF][START_REF] Gulizzi | Micro-cracking of brittle polycrystalline materials with initial damage[END_REF][START_REF] Gulizzi | Modelling intergranular and transgranular micro-cracking in polycrystalline materials[END_REF], small-strains crystal plasticity [START_REF] Benedetti | A grain boundary formulation for crystal plasticity[END_REF], hydrogen embrittlement [START_REF] Benedetti | Grain-boundary modelling of hydrogen assisted intergranular stress corrosion cracking[END_REF], high-cycle fatigue [START_REF] Benedetti | A grain-scale model for high-cycle fatigue degradation in polycrystalline materials[END_REF], micro-cracking of piezoelectric aggregates [START_REF] Benedetti | A microstructural model for homogenisation and cracking of piezoelectric polycrystals[END_REF] and multi-scale modelling of polycrystalline components [START_REF] Benedetti | Multiscale modeling of polycrystalline materials: A boundary element approach to material degradation and fracture[END_REF].

In this work, the framework is employed to contribute to better understanding how the bimodal grain size distribution may affect the micro-mechanics of this class of aggregates. At this stage, no crystal plasticity is considered in the boundary element model and the individual grains behave linearly, with the non-linear processes confined at the intergranular interfaces only: this is done to obtain a preliminary assessment of purely morphological effects on the micro-mechanics of cracking, while more general effects will be included extending the study presented in Ref. [START_REF] Benedetti | A grain boundary formulation for crystal plasticity[END_REF].

Fig. 6 shows the comparison between the micro-cracking patterns obtained for the unimodal FG polycrystal and the bimodal polycrystal. The FG populations have the same mean grain size. The two polycrystals are subjected to a tensile load along the vertical direction, enforced in displacement control conditions. In the unimodal case, if the same cohesive properties are assigned to all the intergranular interfaces, the crack initiates at the grain boundaries lying on planes orthogonal with respect to the loading direction. The interplay between grains anisotropy and orientation mismatch may render some interfaces more susceptible to crack initiation, although the interface orientation remain the dominating initiation factor (in the mentioned modelling assumptions). In the bimodal case, the micro-cracking starts at the interfaces between the larger grain and the smaller ones lying at the top or bottom of the large grain: this is consistent with the unimodal analyses as, in the considered morphology, such interfaces lye on planes normal to the loading direction. After initiation, the micro-cracks follow the border of the large inclusion up to locations at which the external loads starts aligning with the interface, thus reducing the effective traction on the interface. In the bimodal case, also multiple cracking is observed. Fig. 7 shows the macro stress-strain curves obtained from the two analysed cells. In the adopted assumptions, a light effect of the bimodal morphology is detectable only on the softening branch of the curves. It is apparent that the absence of plastic effects has an important effect on such curves, as all the dissipative processes (damage) are confined over the intergranular interfaces, instead than within all the volume (plasticity). More comprehensive investigations will include the effect of crystal plasticity in the boundary element framework.

CONCLUSION

The mechanical problem of a single Coarse Grain (CG) embedded in a matrix of Fine Grains (FG) has been considered. It is inspired from the experimental realisations of 316L stainless steel polycrystals with bimodal grain size distribution, which have led to an enhancement of the combination of strength and deformability by comparison to unimodal polycrystals. In this set configuration of study, crystal plasticity multiscale analyses show that interactions between the CG and the FG could not be returned to a simple mixture rule of the properties: with the presence of the CG, a large stress relaxation channel appears along the tensile direction, high strains localize within the CG, stress distribution is broadened, plasticity activity starts at lower imposed strain and elastoplastic hardening evolves according to two distinct regimes.

This new situation exercised by the FG can have different consequences on the fracture properties, depending on the mechanisms controlling damage and fracture. The investigations on the intergranular fracture are based on the assumptions of an elastic linear behaviour assigned to the grains and a dissipative process related to damage at grain boundaries. They show that the presence of the CG does not significantly modify the macroscopic inception of softening as compared to the unimodal case. This can be explained by the fact that, in the adopted framework, the orientation of a grain boundary with respect to the tensile direction is a dominant factor in the micro-cracking process, and it would remain such both in the uni-modal or bi-modal case, as almost all orientations are represented in both cases.

Yet, as a plane perpendicular to the tensile direction is the most favourable location for crack initiation, microcracking naturally occurs first at the CG-FG interface in the bimodal case, on the top or the bottom of the CG. This also corresponds to the place where crystal plasticity predicts the smallest states of stress within the FG matrix. As experiments on 316L stainless steel have shown the positive effect of the presence of CGs on the elongation to fracture, these patterns of stress and strain interactions between FG and CG as provided by crystal plasticity may be considered as key factors in the micro-cracking process. Fracture analyses of the considered 316L steel should thus include not only the morphological characteristics of the grain boundaries but also the heterogeneous nature of stresses and strains resulting from the polycrystalline plasticity. This will correspond to future developments.

FIGURE 1 :

 1 FIGURE 1: The three considered polycrystalline model microstructures. (From left to right) Unimodal Coarse Grain (grain size 5µm); Bimodal, with a single CG (5µm) embedded in a matrix of FG (0.5µm); Unimodal Fine Grain (grain size 0.5µm).

FIGURE 2 :

 2 FIGURE 2: Tensile stress-strain curves of the three polycrystals.

FIGURE 6 :

 6 FIGURE 6: Micro-cracking patterns computed with the boundary element framework for polycrystalline micro-mechanics in Refs.[START_REF] Benedetti | A three-dimensional cohesive-frictional grain-boundary micromechanical model for intergranular degradation and failure in polycrystalline materials[END_REF][START_REF] Gulizzi | An enhanced grain-boundary framework for computational homogenization and micro-cracking simulations of polycrystalline materials[END_REF] for (a) the unimodal and (b) the bimodal polycrystal. The present simulations do not consider crystal plasticity effects and provide a preliminary assessment of morphological effects of the bimodal distribution on the micro-cracking mechanics in case of linear elastic grains with damaging interfaces.

FIGURE 7 :

 7 FIGURE 7: In the assumption of linear grains with damaging intergranular interfaces, the bimodal morphology has a light effect on the softening branch of the macro stress-strain curves.

ACKNOWLEDGMENTS

Ivano Benedetti acknowledges the support of INSA Rouen Normandie, where he has been invited as Foreign Visiting Professor between November and December 2019.


Thomas Michalek
email: tomas.michalek@fel.cvut.cz.

Aude Bolopion
email: aude.bolopion@femto-st.fr.

Hurak Zdenek

Michaël Gauthier

T Michálek

Z Hurák

Control-oriented

Control-oriented model of dielectrophoresis and electrorotation for arbitrarily shaped objects Tomáš Michálek, Aude Bolopion, Zdeněk Hurák, Senior Member, IEEE, and Michaël Gauthier Abstract-The most popular modeling approach for dielectrophoresis (DEP) is the effective multipole (EM) method. Yet nowadays, in an unrestricted source electric field, a general multipolar description of the potential caused by their charge distribution is only available for spherical objects. We present a method for computing such multipolar representation of arbitrarily shaped objects polarized by arbitrary electric fields. The method is based on numerical solutions of Laplace equation for the potential and the orthonormality of spherical harmonics. We show how the superposition principle holding for the potential can be utilized to construct a basis of solutions, from which the total potentical can be composed in real time. It is then possible to calculate the DEP forces and torques for arbitrarily shaped objects under the influence of arbitary electric fields in fractions of a second, which is needed in model-based control applications. This contrasts with the Maxwell stress tensor-based numerical method. We validate the proposed method against reference numerical solutions and analyse the importance of the higher-order multipolar moments using a sample case of a Tetrisshaped micro-object placed inside a quadrupolar micro-electrode array and exposed to electrorotation. The implementation of the model in Matlab and Comsol are offered for free download.

Index Terms-dielectrophoresis, non-spherical objects, arbitrarily-shaped objects, control-oriented model, electrorotation

I. INTRODUCTION

E LECTROKINETIC effects, when exploited at micro- scale, constitute a fundamental principle of numerous non-contact micro-manipulation devices. They use the electric field to impart forces and torques either directly on the objects of interest or indirectly on the liquid medium surrounding them, which then induces a fluid flow carying the objects to the desired locations.

A widespread representative of the former is the dielectrophoresis (DEP). First described by Pohl [START_REF] Pohl | Dielectrophoresis: The Behavior of Neutral Matter in Nonuniform Electric Fields[END_REF], it is a physical phenomenon enabling actuation of electrically neutral objects. By placing them in an nonhomogenous external electric field, they polarize and the newly emerged charge distribution inside them interacts via Coulomb forces with the source field. Alternating electric fields are commonly used to eliminate the (unwanted) phenomena of electrolysis and electrophoresis (interaction between the possible intrinsic charge in the object and the external field). Two major modes of dielectrophoresis are distinguished: conventional DEP (cDEP) and travelingwave DEP (twDEP). The former arises from a spatially varying magnitude of the field while the latter is due to a spatially varying phase of the field. The closely related concepts of electroorientation and electrorotation, both causing a torque acting on the object, are usually treated separately in the literature. A shared name-generalized DEP (gDEP), coined in [START_REF] Wang | A unified theory of dielectrophoresis and travelling wave dielectrophoresis[END_REF], can be used to jointly describe all the above mentioned polarization-related phenomena. The review of gDEP may be found, for example, in [START_REF] Lei | Review of the theory of generalised dielectrophoresis[END_REF], where the authors consider not only the quasi-static but also arbitrary time-varying fields. The latter is denoted as the transient gDEP or the polarization history (or crossing trajectory) effect. Since the period of the AC electric field is usually much higher than the time scale related to the object's movement, the quasi-static theory is usually sufficient.

A mathematical model of gDEP enables us not only to perform various simulations and analyses helping us to understand the described physics, but it is also necessary when it comes to applying gDEP to a precise position and orientation control of the micro-objects. In such cases, apart from the model accuracy, also its computational time becomes important as it is detailed in [START_REF] Zemánek | Feedback control for noiseaided parallel micromanipulation of several particles using dielectrophoresis[END_REF], where a device for independent position control of several micro-spheres is described. Following this motivation, we aim at developing a fast computational scheme for gDEP force and torque applicable not only to spheres, but to objects of arbitrary shapes.

State of the art

There are two basic approaches for modeling DEP: effective moments method (EM) [START_REF] Jones | Multipolar dielectrophoretic and electrorotation theory[END_REF]- [START_REF] Morgan | AC Electrokinetic: Colloids and Nanoparticles[END_REF] and Maxwell stress tensor method (MST) [START_REF] Wang | General expressions for dielectrophoretic force and electrorotational torque derived using the Maxwell stress tensor method[END_REF], [START_REF] Rosales | Numerical comparison between Maxwell stress method and equivalent multipole approach for calculation of the dielectrophoretic force in single-cell traps[END_REF]- [START_REF] Michálek | Dipole and multipole models of dielectrophoresis for a non-negligible particle size: Simulations and experiments[END_REF]. The former approximates the polarized object by electric multipoles (the higher the order of the multipole, the higher the accuracy of the approximation) and provides an analytical formulas for the force and torque computation. Its use is, however, limited to spherical objects. The latter is applicable to completly general shapes, but imposes huge time and memory requirements due to the underlying finite element (FEM) simulations. For control purposes, we need the speed of EM and the possibility to treat also non-spherical objects provided by MST. Partial solutions exist for ellipsoidal objects [START_REF] Jones | Electromechanics of Particles[END_REF], [START_REF] Stratton | Electromagnetic Theory[END_REF], where analytical formulas for induced dipole usable in EM are known, and also for arbitrary cylindrically symmetrical objects [START_REF] Green | Numerical determination of the effective moments of non-spherical particles[END_REF], [START_REF] Ogbi | Numerical Identification of Effective Multipole Moments of Polarizable Particles[END_REF], where it is shown how multipolar moments up to the 9th order can be obtained from numerical simulations and then used in EM. In both cases, the symmetry requirement stems from the fact that only the linear multipoles are used for the description of the polarized object when all the charges' locations are constrained to a single line. Another related strong limitation of both approaches is that they are valid only if the external electric field is rotationally symmeteric.

Contribution

We propose a method of computing the multipolar moments of the polarized objects utilizing the numerical solution of Laplace equation for potential and the orthogonality property of the spherical harmonics. No limitation regarding the shape of the object, inhomogenity or direction of the source electric field are considered. We furhter use this method to construct a hybrid (consisting of some precomputed parts as well as expressions evaluated in real-time) control-oriented model, which can be evaluated in fractions of a second compared to units of hours in the case of a MST simulation. It serves as a solid foundation for future design of control algorithms with the forseen goal of achieving simultaneous DEP based position and orientation control of arbitrarily shaped micro-objects in fluidic mediums.

II. MATERIALS AND METHODS

A. Summary of the effective multipole method

To compute the gDEP force and torque acting on any polarized object placed in an inhomogeneous external electric field, we can simply use the formulas derived by Jones and Washizu in [START_REF] Jones | Multipolar dielectrophoretic and electrorotation theory[END_REF]

F = N n=1 F (n) (1)
T = N n=1 T (n) , (2)
where N is the order of approximation and

F (n) = 1 2 Re p(n) [•] n (∇) n Ẽ * n! (3)
T (n) = 1 2 Re 1 (n -1)! p(n) [•] n-1 (∇) n-1 × Ẽ * . ( 4
)
These expressions represent the EM method mentioned in the introduction. The multipolar moments represent the electric potential arising from the charge distribution inside the polarized object being a consequence of the induced charge displacement. From a sufficiently far distance, this potential field seems to be generated just by a single point charge with a magnitude given by an integral over the whole charge distribution. If the object was electrically neutral, all the posive and negative charge introduced by the polarization principle should mutually cancel out and the mentioned integral should be equal zero. However, as we approach closer to the charge structure, it will be possible to distinguish some positive and negative areas, and an electric dipole could be a good approximation of such charge distribution. Even closer, more asymmetries in potential would be apparent, and thus quadrupoles or even higher order multipoles will have to be used to describe the considered charge distribution (and its potential field) accurately enough.

F (n) , resp. T (n) is
Formally, this can be shown by taking an arbitrary charge distribution ρ(r) locally confined to an interior of some sphere S of radius R located, without the loss of generality, at the origin of the coordinate system. Outside this sphere, there is no charge. This situation may, for example, represent the alone polarized micro-object. The potential Φ(r) being created by such charge distribution (by the polarized object) may then be written in the following form

Φ (r) = 1 4π 0 S(R) ρ(r) |r -r | d 3 r , (5)
where 0 is the permittivity of a vacuum and r is the location of the test particle. By expanding the term 1 |r-r | in (5) in a Taylor series about the origin, we get

Φ (r) = 1 4π 0 •    p (0) r + p (1) • r r 3 + 1 2 3 i,j=1 r i r j r 5 p (2) ij + • • •    , (6)
where r = r 2 , r i is the ith element of vector r, p (n) are the so-called Cartesian multipolar moments and p

(2) ij s are the individual elements of quadrupolar moment p (2) . Since the sum is infinite, all the higher order moments are omitted from the listing. Notice that, in agreement with our following discussion, the higher the order of multipole is the faster it decays with the distance r from the object's position. If the potential being described was complex and represented by a phasor Φ (r), also the multipolar moments would be complex and denoted by phasors p(n) as for example in Eq. (3) and Eq. (4) .

B. Overview of the proposed modeling scheme

The problem of using the EM method for arbitrarily shaped objects was that up to now there was not any way for fast computation of their multipolar moments. Here, we show, how to overcome this limitation by employing the two key procedures: numerical computation of the multipolar moments and construction of a solution basis enabling us to do all the time demanding numerical simulations in advance and make the rest of the model evaluatable in real time. In Fig. 1, there is s slightly oversimplified overview of the whole model highlighting its main working principles. We choose a set of orthogonal source electric fields (field_1, field_2, . . . , field_n) and use them as a non-complete basis for an approximate description of the external electric field applied to the object of interest. We then use the coefficinets α 1 ,α 2 ,. . . α n of the corresponding linear combination to compute the multipolar moments according to the diagram by using a basis table. This table is obtained in advance by a numerical method (described later on) and contains the multipolar moments of the object when polarized individually by fields field_1, field_2, . . . , field_n. In the following sections, we will detail the individual blocks of the presented diagram.

C. Extraction of multipolar moments from numerical simulations

The multipolar moments are usually used in their Cartesian form as they appear for example in expressions Eq. (3) and Eq. ( 4) and throughout the previous subsection. In the following development we will, however, need to utilize some advantageous mathematical properties of the so-called spherical harmonics and therefore we will use a so-called spherical form of multipolar moments fitting better in the given framework.

These are obtained by performing the same expansion as in Eq. ( 6), but in spherical coordinates. Doing so, we get

Φ(r) = 1 4π 0 m ∞ l=0 l m=-l 4πq * l,m 2l + 1 • Y l,m (θ, φ) r l+1 , (7)
where q * l,m = r l ρ(r, θ, φ)Y l,m (θ, φ)dV are the components of the spherical form of multipolar moments. l is the order of the moment, and m = -l . . . l denotes its elements. Y l,m (θ, φ) = (2l+1)(l-m)! 4π(l+m)! P l,m (cos θ)e imφ are the spherical harmonics defined using Legendre polynomials P l,m (cos θ). i = √ -1 denotes, as usual, the imaginary unit. Since q l,m is already a complex quantity, the extension of the above formula for complex potential Φ (r) is not so trivial as it was in the Cartesian case. Special care will have to be taken later on when complex potential will be needed.

Although the spherical form of multipolar moments varies from the Cartesian one, they both represent the same quantity. We can find expressions for mutual conversion between these two description by comparing the two expansions in Eq. ( 6) and Eq. ( 7) (as it is done for example in [START_REF] Lea | Spherical multipole moments[END_REF]). In the following text we will, for the reasons stated above, mostly utilize the spherical form and we will stress when the Cartesian one will be used.

For further development the property of orthonormality formally written as

Y l,m (θ, φ)Y * l ,m (θ, φ)dΩ = δ ll δ mm , (8)
is of particular importance. In the expression above, δ ij is the Kronecker delta and dΩ = sin θdθdφ, since the integration is over a unit sphere. Other two useful properties are the following conjugation rules for spherical harmonics and spherical multipolar moments

q * l,m = (-1) m q l,-m and Y l,-m = (-1) m Y * l,m . (9)
Applying Eq. ( 8) and Eq. ( 9) to Eq. ( 7) enables us to express the multipolar moments.

At first, let's apply the first one to the potential expression in Eq. [START_REF] Jones | Basic theory of dielectrophoresis and electrorotation[END_REF]. We get

Φ(r) = 1 0 m l=0,∞ m=-l,l (-1) m q l,-m 2l + 1 Y l,m (θ, φ) r l+1 . ( 10
)
Next, we change the order of summation and apply the second conjugation rule from Eq. ( 9). Then, we multiply both sides of the equation by Y a,b (θ, φ) and integrate them over a sphere S(R). We obtain

S(R) Φ(r)Y a,b (θ, φ)dS = = S(R) 1 0 m ∞ l=0 l m=-l q l,m 2l + 1 Y * l,m (θ, φ) r l+1 Y a,b (θ, φ)dS. (11)
Most of the terms on the right-hand side are independent on the integration variable and can be factored out in front of the integral. Thanks to the orthonormality of the spherical harmonics expressed in Eq. ( 8) and modified for a case of integration over a nonunit sphere to S(R) Y * a,b Y a,b dS = R 2 , the only summation term that remains is the one for l = a, m = b. This enables to express the multipolar moments as

q a,b = (2a + 1) 0 m R 1-a S(R) Φ(r)Y a,b (θ, φ)dS (12)
just based on the potential Φ(r) of the polarized object, which is easier to obtain than the unknown charge distribution that generates it. To get the potential, we simply numerically solve the corresponding Laplace equation for the two cases-one in a setting without the presence of the object of interest (just to get the source potential field) and the other including the object. By subtracting the source field from the second solution, we get the potential field Φ(r) corresponding to the polarized object itself. The whole procedure is illustrated in Fig. 2 for an example of a Tetris "S/Z"-shaped object placed in between two planar electrodes. We use this example to demonstrate and prove the validity of the above described method. The object of interest is alongside with its dimensions in detail shown in Fig. 3. We used the method at first to compute its multipolar representation and then we plugged them into the Eq. ( 7) and compared the results against the solution obtained numerically using Comsol Multiphysics 5.1. External electric field with intensity of 1.1238 • 10 20 V/m pointing along the x-axis was used here. Fig. 4 showes the results for various orders of approximation.

In this section we have shown how a multipolar representation can be obtained for arbitrarily shaped and oriented objects located in arbitrary electric fields. After performing conversion from the spherical to Cartesian form, these multipolar moments could be used for DEP force and torque computation according to Eq. ( 3) and Eq. ( 4). Although this procedure may still be faster than using the Maxwell Stress Tensor q 00 = -1.81e-05 q 10 = 1.74e-08 q 01 = 0.011 -0.001i q 20 = -6.33e-13 q 21 = -7.35e-14 -8.29e-14i q 00 = -7.32e-06 q 10 = 3.97e-08 q 01 = -0.001 + 0.013i q 20 = -5.06e-12 q 21 = -2.11e-14 -4.28e-13i q 00 = -2.53e-05 q 10 = -2.23e-08 q 01 = 0.01 + 0.012i q 20 = -5.70e-12 q 21 = -9.75e-14 method, it does not eliminate the necessity for running a couple of numerical simulations every time the object changes its position or orientation, or when there is a change in the external electric field (for example the voltage signal applied to any of the electrodes changes). We will focus on this problem in the following section.

D. Construction of the solution basis

For the sake of simplicity, we will at first consider the orientation of the object to be fixed. In order to avoid repeated numerical computations, it is possible to utilize the principle of superposition holding for an electric potential (and as a consequence of Eq. ( 12) also for the multipolar moments). This means that knowing the multipolar descriptions of the object polarized by external fields E 1 , resp. E 2 , we can get the multipolar moments of the object polarized by αE 1 +βE 2 , α, β ∈ R by the same linear combination of the corresponding known multipoles. See Fig. 5 for a simple illustrational example. It is depicted here how the multipolar moments for the field aligned along the diagonal direction can be obtained from the already known multipolar moments for the two separate cases of fields oriented along the horizontal and vertical directions.

# q S 0,0 q S 1,0 q S 1,1 q S 2,0 q S 2,1 q S 2,2 • • • 1 0 1 0 0 0 0 • • • 2 0 0 1 0 0 0 • • • 3 0 0 i 0 0 0 • • • 4 0 0 0 1 0 0 • • • 5 0 0 0 0 1 0 • • • 6 0 0 0 0 i 0 • • • 7 0 0 0 0 0 1 • • • 8 0 0 0 0 0 i • • • . . . . . . . . . . . . . . . . . . . . . . . .
The basic idea is therefore to run the above described computational scheme for a suitable, rich enough set of source external electric fields and use the obtained results to form a basis-a table containing couples of the source electric field and the corresponding multipolar representation of the objects polarization. This should then enable us to perfom all the subsequent computations in real time without the need of any further numerical computations just as it is shown in Fig. 5. The only remaining problem is that in order to compute the multipolar representation of the object polarized by an arbitrary electric field (not only uniform one), the basis would need to have an infinite number of elements. Since we anyway have to chose and use some finite order of approximation n in EM based force and torque computations, it gives no sense to consider source electric fields with nonzero spatial derivatives of an order higher than n -1 in the process of basis construction. That is because polarization caused by such electric field could not be anyway represented by the considered multipolar moments.

Instead of representing the source electric field by its value and the value of its various spatial derivatives in the point of interest a more compact, and from the viewpoint of the numerical solutions to be performed, also more natural way of representing the source electric field is through the boundary potential conditions that generate it. For the purpose of the basis creation, we encapsulate the object by a large enough sphere S(R ) (so that the boundary effects would be negligible) centered at its center of gravity (or some other point we choose to compute the DEP force in) and by setting the potential on its interior surface we establish the source electric field. To get an orthonormal basis, we once again utilize the spherical harmonics and their orthonormality property. We use the expression Eq. ( 7) evaluated on the spherical boundary of the domain. By setting the qs, according to Table I and evaluating the above described scheme, we generate a set of mutually orthogonal basis elements. Note that the monopole term q 0,0 is always equal to zero since we assume that the object is electrically neutral. Further note that there are no negative indices related to qs. The reason is to save both computational time and storage space, as these can be determined by the above stated conjugation rules (Eq. ( 9)). Since spherical multipolar moments are in general complex quantities, we have to also take into account their imaginary parts as can be seen for example in the 3rd, 6th or 8th row of Table I. An exception here are the elements q a,0 , a ∈ Z, where the imaginary part is always equal to zero (this can be seen again from the conjugation rules in Eq. ( 9)).

The potential, ΦS , resulting from such boundary conditions can alone generally be also complex even though the boundary potential condition applied to the sphere surface is purely real. The shift in phase leading to complex values is caused by the complex permittivity of the polarized object. As a consequence, two sets of multipolar moments have to be extracted from every simulation-one based on the real part of the potential field of the polarized object and the other based on its imaginary part. In contrast to the Cartesian representation of the multipolar moments, we can not use here the phasor notation (since the spherical multipoles are already complex quantitites) and we have to keep the two mentioned parts separately during storage as well as during subsequent computations making the notation slightly complicated. To summarize, from each simulation we will obtain and store 2

• (1 + 2 + • • • + n) = n(n + 1
) complex values being the elements of our basis.

Considering all of the above, we will use for a basis description a notation q B,cde,f a,b , which denotes the element b of the spherical multipolar moment of order a resulting from a source electric field caused by a potential on S(R ). This potential is given by Eq. ( 7) with q l,m = 1 when l = c, m = d and e = r or by q l,m = i when l = c, m = d and e = i (q l,m = 0 for l = c or m = d) and based on whether f = r or f = i, only its real or imaginary part is considered. For example q B,10i,r 2,1 is the second element of the second order spherical multipolar moment extracted from a simulation when the source potential on S(R ) was given by real part of Eq. ( 7) evaluated with q 1,0 = i and all other qs being equal to zero.

In this section, we showed how the principle of superposition enables us to construct in advance a basis based on which multipolar moments resulting from any external field can be then computed in real time.

E. Decomposition of the source electric field

Having the basis and some arbitrary source electric field, for which we want to know the multipolar description of the polarized object, we now need to get the coefficients of the linear combination of the basis elements. As it was already stated above, these coefficients are exactly the same as are used in the linear combination of the corresponding source electric fields (see Fig. 5).

The source electric field is typically generated by some set of micro-electrodes. Knowing the voltage signals applied to the electrodes, we can compute the electric field at the object's location. In some specific cases of electrode designs an analytical solution may be known (for example for the so-called inter-digitated electrode arrays [START_REF] Sun | Analytical solutions of ac electrokinetics in interdigitated electrode arrays: Electric field, dielectrophoretic and traveling-wave dielectrophoretic forces[END_REF], [START_REF] Chang | Closed-form solutions in the electrical field analysis for dielectrophoretic and travelling wave inter-digitated electrode arrays[END_REF]) in others one has to settle for the numerical solution. We used a partially analytical method based on Green's functions [START_REF] Gurtner | Green's function-based control-oriented modeling of electric field for dielectrophoresis[END_REF], which allowed us to compute high enough (in our case up to the 5th order) spatial derivatives of the electric field. In contrast to [START_REF] Gurtner | Green's function-based control-oriented modeling of electric field for dielectrophoresis[END_REF], we used a slightly different discretization scheme. Instead of using semiinfinite rectangles suitable for the electrode array presented there, we divided the electrode plane into small enough square bins enabling us to approximate arbitrary electrode shapes.

For the purpose of basis generation we represented the source electric fields by the spherical moments q S determining the potential boundary conditions for electric field generation. Now, we need to represent an arbitrary source field in the same way. In other words, we need to find such values of moments, which after plugging into expression Eq. ( 7) and evaluating it on the sphere S(R ) would give us a boundary condition corresponding to the source field.

In order to simplify things a little bit, let's for now assume that the source electric field is just real (the extension to the general complex case will follow in the subsequent section), which is the case when all the harmonic voltage signals applied to the electrodes have exactly the same phase. The multipolar description can then be found by at first using the well-known analytical expressions known for the Cartesian multipolar moments of a sphere [START_REF] Jones | Multipolar dielectrophoretic and electrorotation theory[END_REF] (expression 22):

p E,(n) = -4π 0 m R (2n+1) (2n -1)!! ∇ (n-1) E, (13)
and then just converting it to its spherical form. We denote the result as q E a0 , q E a1 , . . . , q E aa . These are then the desired coefficients of the linear combination of our precomputed basis elements. The resulting moments describing the real, respectively the imaginary part of the potential caused by the polarized object are then given by

q lr a,b = n c=1 c d=0
Re q E c,d q B,cdl,r a,b

, where l ∈ {r, i},

q li a,b = n c=1 c d=0 Im q E c,d q B,cdl,i a,b (14)
, where l ∈ {r, i}.

These spherical moments are then converted back to their Cartesian form q lm 0,n , q lm 1,n , . . . q lm n,n → p (n),lm where l, m ∈ {r, i} and assembled the following way to just one complex Cartesian moment

p (n) =(p (n),rr + p (n),ri ) + i(p (n),ir + p (n),ii ). (16)
In this form they can be finally used to compute the dielectrophoretic force and torque.

F. Considering complex potentials

In a case when source potential itself can be complex (e.g. when there are phase-shifted signals applied to the electrodes), we will have instead of just one set of linear combination coefficients q E a0 , q E a1 , . . . , q E aa , two such sets: q E,r a0 , q E,r a1 , . . . , q E,r aa and q E,i a0 , q E,i a1 , . . . , q E,i aa representing the real respectively the imaginary part of the external electric field. The linear combination of the basis elements is then simply extended to

q klr a,b = n c=1 c d=0 Re q E,k c,d q B,cdl,r a,b
, where k, l ∈ {r, i}, ( Im q E,k c,d q B,cdl,i a,b

, where k, l ∈ {r, i}, [START_REF] Chang | Closed-form solutions in the electrical field analysis for dielectrophoretic and travelling wave inter-digitated electrode arrays[END_REF] and the resulting spherical moments are then again converted to their Cartesian form and combined the following way

p (n) = (p (n),rrr + p (n),rri ) + i(p (n),rir + p (n),rii ) + i (p (n),irr + p (n),iri ) + i(p (n),iir + p (n),iii ) (19)
to get the one resulting Cartesian moment suitable for DEP force and torque computations along Eq. ( 3) and Eq. ( 4).

G. Considering orientation of the object

So far, we have assumed that the orientation of the object is fixed. In practice, however, the manipulated object can revolve freely in the space and thus change its orientation. This slight complication can be simply solved by expressing the electric field and its necessary spatial derivatives in a new rotated frame of reference always attached to the object. After performing the force and torque computation, the results have to be again translated into the original global coordinate system. This can be accomplished by simply multiplying the resulting vectors by a corresponding rotational matrix.

III. RESULTS AND DISCUSSION

To prove the validity of the model, we compared its outcomes against the MST reference numerical simulations done in Comsol Multiphysics 5.1. A Tetris "S/Z" piece depicted in Fig. 3 was used as an example of an object, for which the realtime force and torque computations were so far unavailable. The dimensions of the object are included in the picture. As a test scenario for validation of the results, we used the so-called quadrupolar electrode array, which is usually used for torque generation in electrorotation experiments, but it can equally easily generate a DEP force. It is shown together with all the relevant dimensions in Fig. 6. As the driving signals we used four sinusoidal waveforms having amplitudes 10, 10, 50 and 50 V, respectively, and phase shifts 0 • , 90 • , 180 • and 270 • , respectively. The object's orientation was set such that from its initial orientation (as depicted in Fig. 6) it was rotated by 45 • and 25 • subsequently about the x and y axes. Its geometric center was placed on the z axis. Several simulations for levitation heights varying between 60 µm and 150 µm were performed. To keep the errors negligible, we used moments up to the 5th order (32-poles). Note that to accurately represent polarization of non-spherical objects the higher order multipolar moments will be necessary even when the external electric field is uniform as could be already seen in the potential comparison shown in Fig. 4.

In Figures 7 and8, we can see that results from our model resemble well the reference solution when high enough number of multipoles is used. As the object approaches the edges of the electrodes and the electric field around object becomes more inhomogeneous, the accuracy of the solution gets worse, since it would be necessary to use even more multipolar moments in such case to represent accurately the object's polarization.

During the simulations made above, we measured the times required for computation of DEP forces and torques by both the MST method and our proposed model. A conventional PC (Intel Core i5, 3.30 GHz, 8 GB RAM, 64-bit, Win 7) was used. The results for different orders of approximation are stated in Table II. The evaluation of our model is clearly much faster than performing the numerical simulation for MST method. Nevertheless, with the increasing number of multipoles, the .37n , which means that the highest possible multipole computable under the time limit of 1 s would in the current implementation be the 256-pole (n = 8). Since all the necessary numerical simulations can be done in advance (and just once for every new shape of the object), the real-time evaluation of the model takes only fractions of a second in comparison with hours of MST computation. Our Matlab implementation of the described mathematical model can be found at c1 www.example.com.

IV. CONCLUSION

We introduced a simulational scheme based on EM method capable of real-time computation of the DEP forces and torque acting on objects of arbitrary shapes under arbitrary orientations and located in arbitrary electric fields. The enabling ingredience is the usage of numerical solutions of the Laplace equation for getting the multipolar description of objects under such general scenarios. A specific way of how to construct a solution basis (from precomputed solutions) is then described so that all further computations can be done in real time. As such, the described control-oriented model could be used in future for simultaneous control of position and orientation of non-spherical objects with the envisioned biology and microassembly applications.

Fig. 1 .

 1 Fig. 1. Simplified diagram of the whole model.

FEMFig. 2 .Fig. 3 .

 23 Fig.2. Illustration of the procedure for extracting the spherical multipolar moments q a,b of an arbitrary order for objects having an arbitrary shape.

Fig. 4 .

 4 Fig. 4. Comparison of potentials along the circle laying in the xy-plane and enclosing the object as depicted in Fig. 3.

Fig. 5 .

 5 Fig. 5. Illustration of the superposition principle utilization.

Fig. 6 .

 6 Fig. 6. Dimensions of the electrode array and Tetris-shaped micro-object used in the validation simulations.

Fig. 7 .

 7 Fig. 7. Comparison of the force prediction obtained by the described model and the reference MST model.

4 Fig. 8 .

 48 Fig. 8. Comparison of the torque prediction obtained by the described model and the reference MST model.

TABLE I VALUES

 I OF qS USED IN COMPUTATION OF THE BASIS ELEMENTS simulation

TABLE II COMPARISON

 II OF THE COMPUTATIONAL TIMES OF DIFFERENT MODELS computational time t (in seconds) rises exponentially approximatelly accoring to t = 9.768 • 10 -6 e 1

			Model	Computational time
			n = 1 (dipole)	0.1 ms
	EM	method	n = 2 (quadrupole) n = 3 (octupole) n = 4 (hexadecapole)	0.2 ms 0.6 ms 2.1 ms
			n = 5 (32-pole)	9.3 ms
			MST method	1 h 40 m 57 s

ACKNOWLEDGMENT

This work has been supported by the Labex ACTION project (contract "ANR-11-LABX-0001-01"), by the French RENATECH network and its FEMTO-ST technological facility, by the ANR project CoDiCell (grant No ANR-17-CE33-0009), by the MiMedi project funded by BPI France (grant No. DOS0060162/00) and the European Union through the European Regional Development Fund of the Region Bourgogne-Franche-Comte (grant No. FC0013440) and by the Czech Science Foundation within the project P206/12/G014 (Center for advanced bioanalytical technology, http://www.biocentex.cz).


Ronald Terrazas
email: ronald.terrazas@ulb.ac.be

Adrien De Maeijer

Aude Bolopion
email: aude.bolopion@femto-st.fr

Michaël Gauthier
email: michael.gauthier@femto-st.fr

Michel Kinnaert
email: michel.kinnaert@ulb.ac.be

Pierre Lambert
email: pierre.lambert@ulb.ac.be

Capillary Dipoles: Towards Thermocapillary Micromanipulation of Multiple Particles Floating at the Free Surface

Keywords: Thermocapillary micromanipulation, capillary forces, optofluidcs, feedforward, capillary dipole, Cheerios effect

This paper presents a contribution to thermocapillary micromanipulation of floating particles at the free surfaces. Based on a physical model for the so-called capillary dipole interaction, the force between two floating particles is predicted with a Bessel function, whose magnitude can be calculated from the measured distance between particles. This force is then considered as a disturbance force and rejected in the control scheme thanks to a feedforward strategy. We provide theoretical results for this strategy as well as an experimental proof of concept.

I. INTRODUCTION

Within the framework of non-contact micromanipulation, this paper presents a new strategy towards handling multiple particles floating at the free surface thanks to the so-called thermocapillary micromanipulation technique.

Non-contact micromanipulation uses remote force fields to control the position of the objects, such as for instance magnetophoresis [START_REF] Oulmas | 3 d closed-loop motion control of swimmer with flexible flagella at low reynolds numbers[END_REF], dielectrophoresis [START_REF] Gauthier | Analytical formulation of the electric field induced by electrode arrays: Towards automated dielectrophoretic cell sorting[END_REF], acoustophoresis [START_REF] Zhou | Controlling the motion of multiple objects on a Chladni plate[END_REF] or optical effects [START_REF] Yin | Highbandwidth 3d force feedback optical tweezers for interactive biomanipulation[END_REF]. However, the manipulation in parallel of several particles independently is a challenge since all the objects are affected by the remotely induced force fields. Different approaches have been proposed, such as anchoring the objects that should not move [START_REF] Pawashe | Multiple magnetic microrobot control using electrostatic anchoring[END_REF], controlling the force using objects with slightly different physical properties [START_REF] Diller | Control of Multiple Heterogeneous Magnetic Microrobots in Two Dimensions on Nonspecialized Surfaces[END_REF], using the noise inherent to the micromanipulation systems [START_REF] Zemnek | Feedback control for noise-aided parallel micromanipulation of several particles using dielectrophoresis[END_REF] or shaping the force field [START_REF] Wong | Independent control of identical magnetic robots in a plane[END_REF], [START_REF] Chowdhury | Automated manipulation of biological cells using gripper formations controlled by optical tweezers[END_REF].

More recently, some authors have proposed to use the socalled thermocapillary effect. This consists in generating a temperature gradient along the interface between two fluids, leading to a surface stress thanks to the so-called Marangoni effect (surface tension depends on temperature). This stress is then used to trigger flows able to displace particles. In [START_REF] Vela | Non-contact mesoscale manipulation using laser induced convection flows[END_REF]- [START_REF] Munoz | Optimizing the speed of single infrared-laser-induced thermocapillary flows micromanipulation by using design of experiments[END_REF], 1 -100 µm particles are moved in the bottom of a water reservoir by heating the top free surface with a laser. In [START_REF] Hu | Interactive actuation of multiple opto-thermocapillary flow-addressed bubble microrobots[END_REF], the generation of bubbles in oil is controlled, which in turn creates flows in the oil, able to move particles. A similar approach was applied to bubbles in water [START_REF]An opto-thermocapillary cell micromanipulator[END_REF].

Comparing all these principles, [START_REF] Mallea | Thermocapillary micromanipulation: laser induced convective flows towards controlled handling of particles at the free surface[END_REF] showed the micromanipulation of particles with the laser-induced thermocapillary convection to be an alternative method with less or different drawbacks w.r.t. other non-contact micromanipulation techniques. Especially, four key points can be emphasized: (1) using laser heating to generate thermocapillary convection allows to generate fast localized flows which allows to exert significant forces on the particles and also opens the possibility to perform manipulation in parallel; (2) particles manipulated at the interface (i.e. the free water-air surface) move faster than particles manipulated in the fluid bulk. However, this imposes wettability conditions for the particles to float and remain at the fluid interface; (3) using a closed-loop controller allows to control the unstable system defined when manipulating particles at the interface and also improves the performance of the system; (4) the method is based on a local effect, which opens the way towards multiple particles handling.

However, multiple particles handling remains challenging because of the capillary interactions between them (known as the Cheerios effect). This paper proposes a strategy to reject the capillary perturbation induced by a neighboring particle thanks to a feedforward control scheme. It is structured as follows. Section II presents the physics of lateral capillary forces, Section III recalls the feedback controller presented in [START_REF] Mallea | Thermocapillary micromanipulation: laser induced convective flows towards controlled handling of particles at the free surface[END_REF] and presents the new feedforward controller of this study. Simulation and experimental results are reported in Section IV and discussed in Section V.

II. LATERAL CAPILLARY FORCES

As shown in Fig. 1, lateral capillary forces acting on floating particles can be attractive or repulsive, according to the sign of the capillary dipoles Q i = r i sin φ i . Fig. 1: Lateral capillary forces between particles [START_REF] Kralchevsky | Chapter 7 -Lateral Capillary Forces between Partially Immersed Bodies[END_REF]. a) Attracted particles with same wettability (r 1 sin φ 1 • r 2 sin φ 2 > 0) (b) Repelled particles with opposite wettability (r 1 sin φ 1 • r 2 sin φ 2 < 0).

Kralchevsky [START_REF] Kralchevsky | Chapter 7 -Lateral Capillary Forces between Partially Immersed Bodies[END_REF] proposed a model for the lateral capillary force F LatCap , which can be used in an appropriate controller:

F LatCap = 2 π σ Q 1 Q 2 L -1 cap K 1 L -1 cap ||X part 1 -X part 2 || (1)
where σ is the surface tension, L cap is the capillary length σ ρfluid g (m) and K 1 is the modified Bessel function of the second kind, first order. The operator ||X|| refers to the norm of vector X. One of the assumptions here is that particles should not be very close to each other (typically not closer than L cap ).

Equation ( 1) was experimentally validated using the experimental setup shown in Fig. 2. A spherical AISI 304 steel sphere with 500 µm diameter (1) was glued on a Femto-Tools FT-S1000-LAT (50 nN-1000 µN) force sensor, and approached from a IPL-780 resin 2.1 mm diameter perturbing cylinder [START_REF] Zhou | Controlling the motion of multiple objects on a Chladni plate[END_REF]. Since the perturbation arises from the free surface deformation by the perturbing particle (characterized by the capillary charge), equivalence can be achieved by considering equal capillary charge Q i for both situations. The cylinder diameter was selected such that the resulting lateral capillary force F LatCap could be in the range of the force sensor.

Fig. 3 compares equation ( 1) with the experimental lateral capillary force F LatCap as a function of the obstacle-particle distance r obst-part . As foreseen by the model assumptions, the theoretical model provides a good estimate for obstacleparticle distances larger than 2 mm, and is indeed inaccurate for smaller distances. Since the particle is unlikely to be closer from disturbing neighbors than these 2 mm, the following Fig. 2: Force measurement setup: a spherical steel sphere with 500 µm diameter (1) was glued on a force sensor (2), and approached from a perturbing cylinder (3) Fig. 3: Lateral capillary forces between particles. Each measured point represents the average of 50 continuous signals coming from the force sensor. The fitted model has been introduced to better estimate the force at small obstacleparticle distance, for which equation ( 1) is expected not to be accurate. empirical fitted model was used to design the feedforward controller:

F LatCap = 9.9 • 10 -6 NK 1 800.56 m -1 ||X obst-part || (2)
where X obst-part is the center-to-center obstacle-particle vector [m] and the resulting force F LatCap is given in N.

Nevertheless, the thermocapillary actuation force generated by the laser was shown to be about 10 nN [START_REF] Mallea | Thermocapillary micromanipulation: laser induced convective flows towards controlled handling of particles at the free surface[END_REF], which is much lower than the disturbing lateral capillary force F LatCap at small obstacle-particle distances (at 2 mm, the force is still about 2.37 µN). Therefore in the following, the particle will only be controlled at distances larger than 10 mm. According to the model, two 500 µm diameter steel spheres could be manipulated, provided the minimum distance between them should be larger than 4.7 mm, distance at which the F LatCap ≈ 9.95 nN. This limitation could be removed by increasing the laser power.

III. CONTROL STRATEGY WITH AN OBSTACLE

A. Feedback Controller w/o Perturbation Disturbance

The thermocapillary flows actually push particles. Therefore, particle velocity is likely to deviate from the current target direction (Fig. 4). Based on the measured particle position Fig. 4: Based on the measured particle position and velocity, a new laser position is calculated to realign the particle velocity towards the target. X part and particle velocity Ẋpart , the closed loop controller shown in Fig. 5 computes the velocity correction Ẋc part to realign the particle velocity towards the target. As shown in [START_REF] Mallea | Laser-induced thermocapillary convective flows: A new approach for non-contact actuation at microscale at the fluid/gas interface[END_REF], Ẋc part depends on the laser-particle distance, so that the inverse model of this empirical characteristics provides a new laser-particle distance, and therefore the targeted laser position X lasD . The Controlled Mirror System can then be commanded for the targeted laser position to be set. As shown in [START_REF]Closed-loop particle motion control using laser-induced thermo-capillary convective flows at the fluid/gas interface at micrometric scale[END_REF], the transfer function of the Controlled Mirror system is considered to be equal to identity so only the control of the thermocapillary system has to be recalled here.

As represented in Fig. 7, the particle velocity magnitude and direction are actually decoupled, leading to the so-called steady-state particle velocity magnitude rpart QSS and direction θ corr . The targeted laser position X lasD is obtained from the inverse model of the empirical characteristics [START_REF] Mallea | Laser-induced thermocapillary convective flows: A new approach for non-contact actuation at microscale at the fluid/gas interface[END_REF] (r las-part in mm, rpart QSS in mm s -1 ) :

r las-part = -0.732 rpart QSS + 4.177 (3)
The orientation θ las-part of the laser-particle vector X las-part is then computed based on the correction angle θ corr θ las-part = θ part-targ -θ corr [START_REF] Yin | Highbandwidth 3d force feedback optical tweezers for interactive biomanipulation[END_REF] where θ part-targ is the desired orientation, the orientation of the particle-target vector X part-targ . From the laser-particle distance r las-part and laser-particle orientation θ las-part , the laser-particle relative position X las-part is computed in the "Vector Comp" block (see Fig. 7). The corresponding desired laser position X lasD = X part -X las-part is then computed.

B. Feedforward Controller w/ Perturbation

A feedforward controller (green box in Fig. 7) has been added to the system, to reject the disturbance of the lateral capillary force. As sketched in Fig. 6, the disturbing capillary force is calculated with (2) from the obstacle-particle distance. This disturbing force F LatCap is then converted into a disturbing velocity using the so-called Stokes' Law (commonly used to compute a viscous force from a particle velocity, according to [START_REF] Flynn | Parallel transport of biological cells using individually addressable {VCSEL} arrays as optical tweezers[END_REF]- [START_REF] Zhong | Trapping red blood cells in living animals using optical tweezers[END_REF]). In other words, an equivalent velocity rpart QSS FF is computed, which is proportional to the disturbing force:

rpart QSS FF = F LatCap 6πµR part ( 5
)
where µ is the liquid dynamic viscosity and R part the particle radius.

Fig. 6: Based on the measured particle and obstacle positions, a disturbing force F LatCap is computed from ( 2) and converted into a disturbing velocity from ( 5), which can be compensated with an appropriate laser position to realign the particle velocity towards the target.

Beside the disturbance magnitude, a correction angle θ corr FF is also computed in the "FF Direction action" block of Fig. 7 using the following relation:

θ corr FF = θ part-targ -θ obst-part (6)
This equation represents the direction opposite on which the lateral capillary force is acting w.r.t. the particle-target direction θ part-targ . Consequently, we see that this feedforward controller computes an estimated particle velocity at quasi steady-state rpart QSS FF and a correction angle θ corr FF (FF denotes the values of the inputs computed by the feedforward controller) based on the measured obstacle-particle position X obst-part . The system inputs computed by the feedforward and feebdack controllers are then added, as shown in Fig. 7. Using these resulting values, a desired laser position X lasD is computed using the linearization procedure recalled in the previous subsection.

To define the feedforward controller, an inversion based approach was used that enabled to define rpart QSS FF which is computed in the "FF Mag action" block (green box in Fig. 7). This block has as input the obstacle-particle distance r obst-part which is used to compute the lateral capillary force F LatCap acting on the particle using eq. 2. The inputs computed by both controllers are added in the "Input combinator" block as if they were defining a vector with polar coordinates: rpart QSS xx and θ corr xx , according to equations (7-8) on page 4.

The resulting rpart QSS (once limited to 4.5 mm/s in the saturation block) and θ corr are then used to compute a desired laser position X lasD . The latter one allows reaching the goals of the feedback and feedforward controllers. Indeed, for particles close to the obstacle (Fig. 8a), the capillary force F LatCap Fig. 7: Block diagram representation of the entire closed-loop system with the feedforward controller. This controller has two outputs: an estimated particle velocity magnitude at quasi steady-state rpart QSS FF and a correction angle θ corr FF . These outputs are combined with the outputs of the feedback controller which for the sake of distinction are referred as rpart QSS FB and θ corr FB . The combination is performed in the block "Input combinator" which computes the estimated particle velocity magnitude at quasi steady-state rpart QSS and the correction angle θ corr FB which are used to compute the desired laser position X lasD .

rpart QSS = r2 part QSS FB + r2 part QSS FF + 2 rpart QSS FB rpart QSS FF cos(θ corr FF -θ corr FB ) (7)
θ corr = θ corr FB + arctan rpart QSS FF sin(θ corr FF -θ corr FB )

rpart QSS FB + rpart QSS FF cos(θ corr FF -θ corr FB ) (8)
is very strong and the feedforward controller dominates the feedback controller (whose goal is to take the particle to the target location X targ ). Therefore, the desired laser position X lasD generates a total radial force F dr σ r counteracting the large lateral capillary force F LatCap . In Fig. 8a, the laser-particle vector X las-part is almost aligned with the lateral capillary force F LatCap . Conversely, for particles far away from the obstacle (Fig. 8b), the capillary force F LatCap is small and the feedback controller dominates the feedforward controller. Therefore, the desired laser position X lasD leads to a total radial force F dr σ r displacing the particle towards the target location X targ . That is why in Fig. 8b, the laser-particle vector X las-part is almost aligned with the particle-target vector X part-targ . Note well that we considered here the case of repelling particles only.

Fig. 8: (a) For small particle-obstacle distance, the FF action dominates and leads to a force F dr σ r mainly counteracting the large lateral capillary force F LatCap ; (b) For larger particleobstacle distance, the FB action dominates and leads to a resulting force F dr σ r pushing the particle towards the target.

IV. RESULTS

A. Simulation Results

The control strategy was first tested in simulation using the Matlab/Simulink software. The starting and target positions of the sphere are respectively (0, -15) mm and (0, 20) mm (fig. 9). The disturbing cylinder is located in (7, 0) mm.

The model for the particle velocity magnitude ṙpart can be seen as a cascade of a non-linear steady-state characteristics and a linearly time-variant dynamics. The first one is given by : rpart QSS = 26.086 exp(-1.433 r las-part )

The second one links the steady-state output rpart QSS to the particle velocity ṙpart , according to the following second-order transfer function:

VelDyn(s) ≡ ṙpart rpart QSS = 52.027 s 2 + 7.6458s + 52.0274 [START_REF] Vela | Non-contact mesoscale manipulation using laser induced convection flows[END_REF] The feedforward controller was first turned off. We can see in Fig. 9a that the laser primarily aims at pushing the particle towards the target location. Consequently, the X-axis deviation can reach up to 1.22 mm in the worst case. The feedforward controller was then turned on to counteract the effect of the lateral capillary force F LatCap while also displacing the particle towards the target location. As a result (Fig. 9b), the particle is only deviated by a distance of 29 µm in the Xaxis. The lateral capillary force F LatCap could however not be counteracted completely, because of saturation of the actuator.

B. Experimental Results

Experimental results were obtained with a AISI 304 steel 500 µm diameter particle floating on distilled water in a plastic container placed above a white LED panel improving the particle position acquisition in real time with a camera. The laser is moved by a mirror using piezoelectric actuators [START_REF] Mallea | Laser-induced thermocapillary convective flows: A new approach for non-contact actuation at microscale at the fluid/gas interface[END_REF]. As already devised, the perturbation is produced with a fixed IPL-780 resin cylinder (2.1 mm diameter, produced by Nanoscribe Photonics GT). The start and target points have the same Y coordinate approximately only because it is not possible to control exactly the initial position of the floating particle. On its trajectory to the target, the closest point from the cylinder is at 12.8 mm. As shown in Fig. 10, a maximum deviation of 3.05 mm is observed without the feedforward while it has been reduced to 1.27 mm with the feedforward activated.

V. DISCUSSION AND CONCLUSIONS This work opens perspectives towards parallel manipulation of multiple floating particles as long as they are not too close to each other, despite the intrinsic unstable physical principle. Indeed, due to the significant interaction capillary force, the actuator is saturated, which is accounted for by introducing a modle of this saturation in the controller. A more powerful laser would probably create a stronger surface tension gradient, and therefore allow the system to be faster. Of course, the main limitation of this work is that the (single) disturbing object is fixed. A next step would be to test the controller having two free moving particles. The capillary dipole interaction is indeed a 2D analogue of the Coulomb law in electrostatics [START_REF] Kralchevsky | Chapter 7 -Lateral Capillary Forces between Partially Immersed Bodies[END_REF]. Therefore, different forces arising from different particles could be added to estimate the force acting on a given particle. If these disturbing particles are fixes (=obstacles), the control scheme is probably still valid. If they are however free to move as well, the action of the laser beam must be distributed among the different particles to control, and this probably requires a modification of the control scheme. Another perspective would be to mitigate the strong capillary force with another force field. For instance, magnetic repulsion could reduce the capillary attraction. To show the effect of the controllers, the laser and particle positions are plotted for the time instants t = 2, 3, 4, 5, 6, 7 s. Note: the particle is only displaced towards a target location, no trajectory tracking is performed. Step response of the particle in the presence of an obstacle. The particle starts at the right of the cylinder and is controlled towards a target position (magenta cross) with a Y coordinate similar to its origin. The laser and particle positions are plotted for the time instants t = 0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30 s.

X-axis [mm]

Fig. 5 :

 5 Fig. 5: Simplified schematic representation of the control scheme.

ACKNOWLEDGMENT Supported by 7 /

 7 38 MicroMAST IAP (Belspo, Belgium).

Fig. 9 :

 9 Fig.9: Simulation results of a controlled manipulation of a 500 µm diameter steel sphere in presence of a cylindrical obstacle. To show the effect of the controllers, the laser and particle positions are plotted for the time instants t = 2, 3, 4, 5, 6, 7 s. Note: the particle is only displaced towards a target location, no trajectory tracking is performed.

  Fig.10:Step response of the particle in the presence of an obstacle. The particle starts at the right of the cylinder and is controlled towards a target position (magenta cross) with a Y coordinate similar to its origin. The laser and particle positions are plotted for the time instants t = 0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30 s.