Planning in-hand dexterous micro-manipulation using 3-D rotations decomposition

Pardeep KUMAR, Redwan DAHMOUCHE, Michaël GAUTHIER

Abstract-This paper aims to contribute to the improvement of dexterity in contact micro-manipulation by performing inhand dexterous micro-manipulation planning. Previous experimental works on planar micro-manipulation showed that such an approach allows for large rotations of arbitrary shaped objects. Moving from planar to 3-D manipulation significantly increases the complexity of the manipulation planning, especially when considering the rolling of the fingers on the object during the manipulation. We propose in this paper a dexterous manipulation planning algorithm that leverages the complexity of 3-D manipulation planning by decomposing the desired 3-D rotations into three successive rotations within two different planes. Optimal paths of the manipulating fingers are thus obtained in the planar spaces and then combined to form the trajectories in the 3-D space. Besides the relevance of the approach, the simulation results show that exploiting adhesion forces improves the robustness of the manipulation and extends the manipulation capabilities, but at the expense of the computation time.

Index Terms-Micro-manipulation, Manipulation Planning, Dexterity, Grasping, In-hand manipulation.

I. INTRODUCTION

Dexterous manipulation at the macro-scale has been an active research topic in the last three decades [START_REF] Ozawa | Grasp and dexterous manipulation of multi-fingered robotic hands: a review from a control view point[END_REF]- [START_REF] Okamura | An overview of dexterous manipulation[END_REF]. Different operations like rolling [START_REF] Bicchi | Dexterous manipulation through rolling[END_REF], sliding [START_REF] Brock | Enhancing the dexterity of a robot hand using controlled slip[END_REF], and finger gaiting [START_REF] Han | Dextrous manipulation by rolling and finger gaiting[END_REF] have been studied to perform in-hand manipulation tasks [START_REF] Bullock | A hand-centric classification of human and robot dexterous manipulation[END_REF]. The demonstrated capabilities range from simple pick-andplace to more complex tasks, such as in-hand rotation of a cube [START_REF] Andrychowicz | Learning dexterous in-hand manipulation[END_REF], picking up a coin from the floor [START_REF] Bullock | A hand-centric classification of human and robot dexterous manipulation[END_REF], opening a bottle cap [START_REF] Ott | A Humanoid Two-Arm System for Dexterous Manipulation[END_REF], and many others.

At the micro-scale, manipulation methods are classified into two groups: (i) non-contact manipulation, in which the object trajectory is controlled with a physical field (e.g. electric, magnetic, etc.) and; (ii) contact-based manipulation, where the manipulator is in physical contact with manipulated object, just like the way humans manipulate common objects. This paper focuses on this second category, which is currently limited to basic tasks [START_REF] Xu | Design, fabrication, and testing of an mems microgripper with dual-axis force sensor[END_REF]- [START_REF] Sun | A novel piezo-driven microgripper with a large jaw displacement[END_REF].

In contact based micro-manipulations, two main types of robot architectures are proposed. The first one consists of a basic tweezers mounted on a precision manipulation robot. Yet, we face hurdles in achieving the desired rotations using this architecture, while achieving translations is usually more trivial. Indeed, the uncertainties and variations of the robot FEMTO-ST Institute, Univ. Bourgogne Franche-Comté CNRS, Besanc ¸on, France.(firstname.lastname@femto-st.fr) kinematic parameters and the existence of mechanical defects (backlash, out of roundness, etc.) affect the position of the rotation center of the tweezers and introduce undesired motion of the manipulated object [START_REF] Hériban | Improving rotation behaviour of robotic structures for micro-assembly[END_REF]. In addition, the entire robot has to rotate in order to manipulate the micro-object, which requires a large space. Unfortunately, such spaces are not available in many applications such as in Scanning Electron Microscope's chambers or in several medical applications. Such constraints limit the feasibility of this approach.

The second approach in contact micro-manipulation is to use two or more decoupled probes or fingers that have several Degrees of Freedom (DoF). In this approach, the trajectory planning of the fingers is more complex but the manipulation of the objects requires less space since it is done in the manipulating "micro-hand" [START_REF] Sjövall | Automatic dextrous microhandling based on a 6-dof microgripper[END_REF], [START_REF] Inoue | Micro hand with two rotational fingers and manipulation of small objects by teleoperation[END_REF]. However, most of the works that use this approach are limited to simple object shapes like spheres [START_REF] Chu | Automated dual-arm micromanipulation with path planning for micro-object handling[END_REF]- [START_REF] Liu | A novel dual-probe-based micrograsping system allowing dexterous 3-d orientation adjustment[END_REF]. Other works consider more complex shapes but the dexterity of the proposed approaches is limited. For instance, two fingers are used in [START_REF] Wason | Dextrous manipulation of a micropart with multiple compliant probes through visual force feedback[END_REF] to rotate a Lego-like object but the rotation angles are small since no re-grasp operation is performed. In [START_REF] Venkatesan | Path Planning and Micromanipulation Using a Learned Model[END_REF], the manipulated object lays on a substrate (2-D motion) and is pushed by four probes to get to the desired object's pose. The manipulation is thus limited to a horizontal plane.

Contact micro-manipulation is also made difficult because of surface forces (capillary, van der Waals and electrostatic). Those forces are often dominant over gravitational and inertial forces which make the manipulated objects stick to substrate and the tweezers [START_REF] Fearing | Survey of sticking effects for micro parts handling[END_REF]. For a long time, adhesion forces were seen as disturbing effects that we have to mitigate [START_REF] Dejeu | Reducing the adhesion between surfaces using surface structuring with ps latex particle[END_REF]. However, it was recently shown that such forces can be exploited to improve dexterous micromanipulation capabilities. Indeed, the adhesion forces allow the fingers to pull the object which enhances the stability of the grasps [START_REF] Seon | Planning trajectories for dexterous in-hand micro-manipulation using adhesion forces[END_REF]. In [START_REF]On the contribution of adhesion and friction in planning dexterous in-hand micromanipulation[END_REF], Seon et al. demonstrated that in-hand planar micro-manipulation of arbitrary shaped planar objects can be achieved in presence of adhesion forces. Dexterous manipulations with rotations above 220°for a wide variety of planar objects at the micro-scale were experimentally achieved.

The aim of this work is to enhance the capabilities of manipulating micro-objects with complex shapes in the 3-D space. The main challenge in performing 3-D dexterous manipulation in general, and micro-manipulation in particular, is the complexity of the manipulation planning algorithm that might increase exponentially with the dimension of the configuration space of the object and that of the hand.

In this article, we develop a dexterous 3-D manipulation algorithm with a moderate complexity. The main idea is to perform 3-D dexterous manipulation by using 2-D planners. The proposed approach consists in: i) decomposing 3-D rotations into three individual 2-D rotations; ii) planning the manipulation within the sub-spaces and; iii) combining them to form 3-D manipulation trajectories. The increase in complexity from 2-D manipulation to 3-D is thus moderate (three 2-D paths). The adhesion forces are exploited to stabilize the grasps and the A * search algorithm is used to generate optimal paths in the considered sub-spaces.

II. MODELING AND BACKGROUND

A. Grasping Forces

To grasp an object with F fingers, the fingers must make contact and apply some grasping force on the object's surface. As proposed in [START_REF] León | From robot to human grasping simulation[END_REF], we assume in micro-scale that the contact forces can be modeled using a combination of the Coulomb's Law in which the friction magnitude is independent of the velocity and contact area, and the pull-off force ( f po ) representing the force required to detach the finger from the object as:

f 2 x + f 2 y ≤ µ( f n + f po ) , (1)
where f x and f y are the tangential components of the force on Xaxis and Yaxis respectively, f n is the normal component of force, and µ is the friction coefficient.

Due to the presence of the pull-off force f po , it is possible to apply a negative force (pulling the object) in micro-scale since the applied force lies in the modified friction cone, whereas in macro-scale only the positive grasping forces (pushing the object) are possible.

B. Grasping Equilibrium

The manipulated object must be in equilibrium during the whole manipulation process. The equilibrium condition for rigid body is that the sum of all the wrenches should be equal to zero. Whereas, a wrench vector is composed of forces and torques at a contact point as provided by:

w = f τ ( 2
)
where w is a wrench, f and τ represent the force and torque respectively. Grasp Equilibrium has been discussed in detail in [START_REF] Seon | Planning trajectories for dexterous in-hand micro-manipulation using adhesion forces[END_REF], [START_REF]On the contribution of adhesion and friction in planning dexterous in-hand micromanipulation[END_REF], [START_REF] Seon | Enhance in-hand dexterous micromanipulation by exploiting adhesion forces[END_REF]. A grasp using F fingers is stable if:

F ∑ i=1 w i + w ext = 0 , (3)
where w i is the grasping wrench applied by i th finger, and w ext is the external wrench applied to the object.

C. Impact of Pull-off forces on Finger Gaiting

Pull-off force has an important role in re-grasping and finger gaiting (i.e. fingers replacement). When a finger (e.g. the F th finger) is being detached, it will pull the object with a force corresponding to the pull-off force ( f po ). This pulloff force may disturb the grasping equilibrium. Thus, it is necessary for the remaining F -1 fingers to compensate for this f po to maintain the object's stability. Thus, removing a finger is possible only if the following equation is satisfied:

F-1 ∑ i=1 w i + w ext + w po = 0 , (4)
where w po is the pull-off wrench caused by releasing the F th finger.

III. PROBLEM FORMALIZATION

In order to perform fingers path planning, we assume that the object's shape is known through its Computer-Aided Design model (CAD model). Considering the small size of the manipulated object, we also assume that the fingers movement range is significantly larger than the object's dimensions. We also consider only two and three fingers grasps. The admissible contact points on the object that can be used to grasp it are obtained by sampling the object's surface. We then compute the set of all stable grasps, called Maps (M), as well as all the admissible finger gaiting configurations (3 to 2 fingers grasps that can resist a finger detachment). Therefore, there are three different sets of maps: i) initial grasp for which the object can be detached from the substrate (M IG ); ii) object translation and rotation while the fingers roll on the object without sliding (M RG ); iii) and reconfiguration where the number of fingers used for a grasp changes (M RC ).

To generate finger trajectories, the maps are converted to Graph (G), where each stable grasp is connected to its subsequent stable grasp(s). Each element of graph is considered as a Node (n) which are connected by Edges (e). Finding an optimal finger trajectory to rotate an object consists in finding a path in the Graph (G). This step can be performed using a graph search algorithm. Since, a node is an outcome of links of stable grasp configuration of manipulated object, which is based on contact points, fingers, and object orientation; thus, it is possible to represented a node by four parameters i.e. n = [i j k l ] where i, j, and k define the indexes of contact points on object of the first, second and third fingers respectively, and l defines the index of angular position of the object.

For finger path planning, we are using A * algorithm whose time complexity is O(b d ) where d is the depth (number of nodes to reach the goal) and b is the branching factor (average number of branches per node). There are two cases for branching factor, one is when two fingers are in use, and the second is when three fingers are used. For first case when two fingers are already in contact with object, the branching factor is sum of operations: clockwise rotation, counter clockwise rotation and addition of third finger over the remaining c -2 contact points. Whereas for the second case when three fingers are in use, the branching factor is sum of operations: clockwise rotation, counter clockwise rotation, and detaching one of the finger.

A. Reducing the complexity using Euler's angle

Beside the complexity of search algorithm, the major issue to deal with is the number of grasping possibilities (gp), which is the result of the number of contact points (c) of an object, the numbers of fingers (F) being used to grasp the object, and the possible angular positions/orientations (l θ ) of object (i.e. one turn of rotation (360°) divided by constant rotational step ∆θ).

In planar manipulation, where we consider to manipulate an object in a single plane, the number of grasping possibilities is computed as:

gp 2D = c F • l θ , (5)
Whereas, to manipulate an object in spatial case (S ), we imagine P numbers of planes over which the contact points are generated. In this case, the number of grasping possibilities is computed as:

gp SD = (P • c) F • l θ = P F • gp 2D (6)
This increase in the number of grasping possibilities, exponentially increases the complexity of search algorithm. One of the challenges in 3-D manipulation is to reduce the number of grasping possibilities and search algorithm complexity. According to Euler's Angles, it is possible to decompose any 3-D rotation into three individual 2-D rotations along any two orthogonal axes [START_REF] Euler | Problema algebraicum ob affectiones prorsus singulares memorabile[END_REF]. Thus keeping that in mind, we propose to only retain two orthogonal planes P 1 at XYaxis and P 2 at XZaxis intersecting the object over lines L 1 and L 2 as illustrated in Fig. 1. Since, the three individual 2-D rotations will be carried out in two different planes, we consider the two intersecting points I and J of lines L 1 and L 2 as a common link for these three individual rotations. The three individual rotations will be carried out as: R(z, θ 1 ) over P 1 , R(y, θ 2 ) over P 2 , and R(z, θ 3 ) over P 1 respectively as represented in Fig. 2, which will reduce this number of grasping possibilities to sum of grasping possibilities of each rotation carried out over these two planes as:

gp 3D = 2 • gp 2D L 1 + gp 2D L 2 . ( 7
)
The reduction of complexity introduced by two planes also reduces the genericity of our approach. Indeed the considered objects should have a geometry in which both planes P 1 and P 2 can be defined.

B. Sampling strategy

As we consider that the object's CAD model is known, we sample the manipulation lines L 1 and L 2 to generate contact points (c) for grasping. We consider the curvilinear abscissa as the coordinates for each contact point on the object. Whereas, for sampling on the object we considering that the rolling between two successive sample points corresponds to a constant rotation ∆θ of the object.

IV. IN-HAND MANIPULATION STRATEGY

Given the number of fingers, and object geometry, the finger trajectory generation is accomplished using two steps. The first is to compute the stable grasps and generate the Graph (G), while the second step is to define a path by traversing the graph(s) to achieve the desired configuration.

Since, the object modeling is restricted only to two planes P 1 and P 2 , and intersecting contact points I and J being the only common link induces some constraints for successive rotations. These constraints are:

• The first rotation should end at the intersecting points I and J. • The second rotation should start from these intersecting points I and J, and end at the same intersecting points. • The third rotation should start from these intersecting points I and J.

In previous work, Seon et al. [START_REF] Seon | Planning trajectories for dexterous in-hand micro-manipulation using adhesion forces[END_REF], [START_REF]On the contribution of adhesion and friction in planning dexterous in-hand micromanipulation[END_REF], [START_REF] Seon | Enhance in-hand dexterous micromanipulation by exploiting adhesion forces[END_REF] achieved only desired pose of object for planar manipulation, while to manipulate the object in 3-D we need to comply with these constraints. Thus, we propose an original method to provide finger trajectories using A * algorithm.

The A * algorithm uses a heuristic for traversing the graph(s), while ensuring that it computes a path with minimum cost through the nodes n. The algorithm optimizes the function f (n), which consists of the cost function g(n) and heuristic function h(n) as:

f (n) = g(n) + h(n). (8)
A. Cost Function g(n)

Cost function is used to characterize the distance between two nodes. There are three cases for cost function g(n); i) initial node cost g(n o ) i.e. when fingers grasp the object for first time, ii) when current node n is reached using finger rolling on object g(n) rolling , iii) when current node n is reached using finger gaiting g(n) gaiting .

Initial node cost g(n o ), is the distance between fingers' initial position, and first contact points on object for grasp. When the node n is reached using finger rolling (i.e. without reconfiguration), the cost can be defined as the sum of rolling distance d roll F i which is the curvilinear abscissa between two contact points covered by each finger F i , and previous cost g(n previous ) as:

g(n) rolling = 3 ∑ i=1 d roll F i + g(n previous ) . (9)
Whereas when the node n is reached using a finger reconfiguration (i.e. finger attachment or detachment), the cost g(n) gaiting will be:

g(n) gaiting = g r + g(n previous ) , ( 10
)
where g r is a constant and estimated as the minimal distance applied by the actuator to guarantee that the finger is detached. In our case, we choose the value of g r as three times the finger radius.

B. Heuristic Function h(n)

In search algorithms, heuristic function h(n) estimates the minimum cost (in terms of distance for our case) from current node n to the goal node n g . As it is important to develop a good heuristic function that guarantees the shortest path; for which the heuristic must underestimates the actual cost. Since, a node is dependent on contact points, fingers being used, and angular position of object, thus we take these parameters into account to formalize our heuristic as:

h(n) = 3 ∑ i=1 d roll F i + h r (p), (11)
where d roll F i is the estimated rolling distance covered by finger F i from current node n to goal node n g , and h r is the heuristic for reconfiguration, function of parameter p:

p = (i + l -i g -l g ) 2 + ( j + l -j g -l g ) 2 + (k + l -k g -l g ) 2 (12)
where current node is n = [i j k l ] and target node is n g = [i g j g k g l g ].

Parameter p checks whether a path between current node n and goal node n g without any reconfiguration exists or not. If such path exists then p will always be zero which will result in h r (p) being zero as well. Otherwise, to reach the goal node n g we will need at least one reconfiguration, and the value of h r (p) in such case will be three times the finger radius.

V. RESULTS

The proposed methodology presented in the previous sections, has been simulated and implemented to generate the finger trajectories for three objects with different curvatures in 3-D space i.e. Ellipsoid, Convex shaped object, and Concave shaped object provided in Fig. 3, using three spherical fingers having diameter of 10µm each. For the simulations, we have considered the physical properties of silicon for all fingers and objects, a pull-off force of 1.5 µN [START_REF] Gauthier | Analysis and specificities of adhesive forces between microscale and nanoscale[END_REF], and a friction coefficient of 0.3.

As described in section IV, we propose to decompose the movement in 3 successive rotations as: R(z, θ 1 ) over XYplane P 1 , R(y, θ 2 ) over XZplane P 2 , and R(z, θ 3 ) over XYplane P 1 respectively.

A. Characteristics of the graph G

The Table I represents the total number of nodes for finger path generation against number of contact points over manipulation lines L 1 and L 2 when adhesion forces are not present. The difference in number of nodes generated is 2 -4 times higher when adhesion forces are available in comparison to when there is no adhesion force. Indeed, when there is adhesion, the object is more stable and a higher number of grasping is stable.

B. Example of finger trajectory

Fig. 456represent an example of the optimal finger path to perform three successive optimal rotations (respectively 40°, 70°, 20°) of Concave Shaped Object. The object is expended on Z-Axis, which is the case for most of the objects at micro-scale built using clean-rooms micro-fabrication methods. We show that the proposed method enables to determine a finger path planning to rotate the micro-object in 3-D space.

C. Performance of the path planning

In order to illustrate the performance of our proposed methods, we determined several optimal finger paths considering the objects presented in Fig. 3. In all the cases, we have considered two possibilities: with adhesion and without adhesion to evaluate the impact of adhesion on fingers paths. Various large rotations: for ellipsoid object (40°, 30° †, 50°), (40°, 140°, 70°), (130° †, 150°, 90°), (90°, 120°, 140°), for convex shaped object (40°, 70° †, 30°), (40°, 140°, 80°), (90°, 160°, 30°), (80°, 90°, 30°), and for concave shaped object (50°, 140°, 30°), (40°, 70° †, 20°), (40°, 90°, 20°), (60°, 150°, 30°) have been carried out. The rotations with ( †) represent particular cases where the desired rotation cannot be reached due to the instability of object in the absence of adhesive forces. Indeed, some of the rotations are not accessible without adhesion. It comes from the fact that the graph G is smaller without adhesion.

a) b) c) d) e) f) g) h) i) j) k) l) m) n)
The technical data of three simulations are described in Table II. With adhesion forces, the reconfigurations are required only for the 2 nd rotations, while in the case of without adhesion forces (60°, 150°, 30°for concave shaped object), the other rotations (1 st or 3 rd ) may require reconfiguration(s). The average time to generate: the finger path of 1 st rotation for all various object is around 50 seconds (with adhesion) and 1 second (without adhesion), the finger path of 2 nd rotation is 3 seconds (with adhesion) and 4 seconds (without adhesion), and the finger path of 3 rd rotation for all objects is 1 second in both cases with and without adhesion. The time taken to generate the finger path of 1 st rotations(s) is higher due to the initial grasp possibilities from fingers' reference position (off-contact) to initial nodes on object (on-contact).

VI. CONCLUSION

We proposed in this paper a new method to perform finger path planning for 3-D dexterous manipulation of microobjects. To leverage the complexity of 3-D dexterous micromanipulation, we proposed to perform three planar rotations about two perpendicular axes. To ensure the continuity of the manipulation process over the three rotations, we proposed an algorithm that is able to rotate the manipulated object over three successive angles while starting and/or ending the rotations with predefined fingers positions on the object. The simulation results showed that most trajectories are generated within few seconds. Exploiting adhesion forces enables more feasible trajectories with a lower number of reconfigurations but at the expense of a higher number of nodes (stable grasps) and a longer calculation time. The focus for our future work will be on the fingers collision avoidance and the experimental validation of the proposed method.

Fig. 1 :Fig. 2 :

 12 Fig. 1: 3-D object model of ellipsoid and two intersecting orthogonal planes.

Fig. 3 :

 3 Fig. 3: Different objects a) Ellipsoid, b) Convex Shaped Object, and c) Concave Shaped Object, each with isometric, and 2-D top and side views; used for in-hand dexterous manipulation in 3-D space.

Fig. 4 :

 4 Fig. 4: Concave Shaped Object (2-D top view of XY plane): Sequence of first rotation of 40°along XY plane/z-axis; (a)initial grasp of object (b)-(e) rotation of 40°with rotational step of 10°.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Concave Shaped Object (2-D side view XZ plane): Sequence of second rotation of 70°along XZ plane/y-axis; (a) position of fingers on the intersection points (I, J) (sec. III), (b) rotation of rotation of 10°, (c)-(f) reconfiguration of finger composing of finger detachment and finger addition, (g)-(h) rotation of 20°, (i)-(j) finger detachment to go back to its reference position, (k)-(n) rotation of 40°.

TABLE I :

 I Simulated numerical data for both (i) considering adhesion with a pull-off force of 1.5 µN and (ii) without adhesion. Finger diameter is 10µm, friction coefficient (µ) is 0.3, and rotational step (∆θ) is 10°.

		Number	Number of nodes in the generated
	Object	of contact points c	graph (in millions) With Adhesion Without Adhesion
		(L 1 , L 2 )	Nodes n (L 1 , L 2 )	Nodes n (L 1 , L 2 )
	Ellipsoid	362, 362	337.8, 337.8	110.7, 110.7
	Convex	288, 302	530.1, 595	268.8, 228.5
	Concave	440, 302	1817.4, 595	416.2, 228.5

TABLE II :

 II Comparison of statistical data for an object for both cases with and without the impact of adhesive forces, where W.A., Wo.A.,O.L., C.L., stands for with adhesion, without adhesion, open list (explored by the A * algorithm) and closed list (describing the optimal path) respectively, and θ i represents i th rotation to be carried out for θ°.

	Object	Rotations	Number of nodes in: W.A. Wo.A.	Number of required reconfiguration(s) Time to generate path (sec)
		θ i	θ°O.L.	C.L.	O.L.	C.L. W.A.	Wo.A.	W.A.	Wo.A.
		θ 1	40	147,612	5	6,522	5	0	0	45	< 1
	Ellipsoid	θ 2	140	2,536	19	1,689	19	4	4	5	4
		θ 3	70	946	6	946	6	0	0	< 1	< 1
		θ 1	80	143,824	9	8,696	8	0	0	48	1
	Convex	θ 2	90	650	14	930	16	4	6	2	3
		θ 3	30	682	4	682	4	0	0	< 1	< 1
		θ 1	40	224,458	5	6,958	5	0	0	63	< 1
	Concave	θ 2	90	428	12	930	16	4	6	2	3
		θ 3	20	693	3	693	3	0	0	< 1	< 1

Acknowledgment: This work was funded by the EUR EIPHI program (Contract No. ANR-17-EURE-0002), and by HEC Pakistan's HRDI-UESTP program in collaboration with Campus France.


Emna Dellali

François Lanzetta
email: francois.lanzetta@univ-fcomte.fr

Sylvie Begot

Eric Gavignet

Jean-Yves Rauch

Data reduction of friction factor, permeability and inertial coefficient for a compressible gas flow through a milli-regenerator

A regenerator of a Stirling machine alternately absorbs and releases heat from and to the working fluid which allows to recycle rejected heat during theoretical isochoric processes. This work focuses on a milli-regenerator fabricated with a multiple jet molding process. The regenerator is a porous medium with a total length L = 60 mm and an internal diameter D = 5 mm filled with a dense pillar matrix. The pillars have a geometrical lens shape with a form factor FF = 0.5 (aspect ratio width/length) and a height h = 2.5mm. Two metallic layers (chromium and copper) with a total thickness of 800 nm are deposited on the polymer pillars to increase heat transfer inside the regenerator. We performed experiments on different milli-regenerators corresponding to three porosities ( = 0.80, 0.85 and 0.90) under nitrogen steady flow and oscillating flows (oscillating Reynolds number in the range 0 < Re < 60 and Reynolds number based on the hydraulic diameter ReDh,max<6000) for different temperature gradients (T < 100°C). Temperature, velocity and pressure experimental measurements are performed with microthermocouples (type K with 7,6 µm diameter), hot-wires and miniature pressure sensors, respectively, to estimate the friction factor, the permeability and the inertial coefficient for steady and oscillating compressible flows of nitrogen.

Introduction

The regenerator plays a key role on the performances of Stirling machines (engines, coolers, heat pumps). It stores and releases heat during the isochoric displacements. In order to provide high efficiency of the machine, this regenerative heat must be exchanged in a regenerator with a high thermal efficiency. In an ideal Stirling cycle, this efficiency is 100 % and the thermal efficiency of the machine corresponds to the Carnot efficiency. In this case, all thermodynamic processes are thermodynamically reversible, there is no pressure drop in the regenerator, not heat conduction loss from the warm end to the cold end, there is an infinite rate of heat transfer between the working gas and solid matrix and the solid matrix presents an infinite heat capacity. In reality, the regenerator efficiency is less than 100 % and the fluid flow across the regenerator creates a pressure drop and the heat transfer between the gas and the solid matrix are not reversible and infinite. A regenerator is a solid porous matrix with a high surface ratio in order to increase the heat transfer with the working gas. To increase the heat transfer, the speed of the fluid is also increased. But a high velocity flow increases the pressure losses and finally decreases the performance of the regenerator. The regenerator design must be optimized for each gas (Nitrogen, Helium, Hydrogen) considering different geometrical parameters (length, porosity, hydraulic diameter) and thermal parameters (temperature gradient between the two ends, thermophysical properties of the solid matrix).

Many experimental and numerical work has been performed in order to analyse the performances of regenerators [1][START_REF] Ibrahim | Stirling Convertor Regenerators[END_REF][START_REF] Gheith | Fluids Engineering Division Summer Meeting[END_REF][START_REF] Sadrameli | [END_REF][5]. In Computational Fluid Dynamics analysis, the regenerator can be modelled using commercial codes (ANSYS FLUENT, Star CCD, Open Foam) upon different numerical methods (SIMPLE, PISO) [1,6]. The regenerator is divided into several sections (from three to five) and the results are focused on spaces and time variations of the temperature gradient along the regenerator and heat transfer between the gas and the solid matrix. Experimental works concern flow resistances and heat transfer characteristics of regenerators with different geometrical forms and materials (wire mesh, metal foams, microchannels, involute foils, sponge metals). Authors performed experiments to determine the friction factors the Darcy coefficients, the inertial coefficients). Tanaka et al. obtained the efficiency and heat transfer coefficient by measuring the pressure drops and the periodic temperature variation at both ends of the regenerator (wires, sponge metal) [7]. They determined the friction factor empirical relationship based on hydraulic diameter considered as a representative length. Ishii et al. determined flow characteristics and pressure drop for matrix composed of sintered metal fiber with high porosity [START_REF] Ishii | Proc. 17th ISEC[END_REF]. They obtained an empirical equation for the friction factor function of the hydraulic diameter and the porosity of the matrix. They formulated an empirical equation for the Nusselt number function of the Reynolds number and they estimated the reheat loss and thermal efficiency. Kumar et al. proposed a correlation-based methodology to quantify the Darcy permeability and Forchheimer's inertial coefficient of porous structures for Stirling cryocoolers [START_REF] Kishor Kumar | [END_REF]. They performed a CFD simulation model (2D axis-symmetric) with Ansys Fluent and determined the pressure drops at different mass-flow rates at both ends of a regenerator. They identified and compared a correlation to reported experimental pressure drop. Landrum et al. calculated Darcy permeability and inertial coefficients for different mesh screens from CFD simulations compared to experimental study [10].

In this study, we will determine the permeability and inertia coefficients for unidirectional and oscillating compressible gas flows. We will identify semi-empirical correlations from experimental measurements performed in milli-regenerators of different porosity.

Materials and methods

Experimental setup

Fabrication of the milli-regenerator

An additive manufacturing process, a Multiple Jet Molding, associated to a 3D printer ProJet™ HD 3500 is used to produce the different geometries of regenerators. The laminating pitch is 0.2 mm. The thickness of the layers deposited is 32 μm for a resolution of 375 x 375 x 790 DPI (X Y Z). The precision of the machine is of the order of 0.1 -0.2% of the dimensions of the printed object. The material is a UV-curable acrylate polymer of the type "Visijet 39 Crystal and their thermophysical properties are presented in the Table 1. Each regenerator consists of a cylinder with a rectangular channel with a matrix of pillars runs through. Two of three regenerators with different porosities ( = 0.85 and 0.90) are presented in Figure 1. The regenerator is a porous medium with a total length L = 60 mm and an internal diameter D = 5 mm filled with a dense pillar matrix. The pillars have a geometrical lens shape with a form factor FF = 0.5 (aspect ratio width/length) and a height h = 2.5mm. Two metallic layers (chromium and copper) with a total thickness of 800 nm are deposited on the polymer pillars to increase heat transfer inside the regenerator.

Test bench

The experimental setup is shown in Fig. 2. An electric motor (1) drives two pistons (2 and 3) by means of a connecting rod/crank mechanism to generate the alternating flow (Fig. 2a). In our experiments, the angle between the two pistons is set to 180° in order to generate only an alternating transfer flow between the two volumes of the chambers. The rotational speed of the electric motor varies between 10 to 600 rpm. The stroke of each piston varies between 5 and 25 mm and the inner diameter of the cylinder of each chamber is 40 mm. The reference crank angle is set to 0° when the cold side piston is at its Bottom Dead Centre (BDC) and the hot side piston at its Top Dead Centre (TDC) (Fig. 2b). The milli-regenerator ( 6) is in the middle of the setup between two heat exchangers (4 and 5). These heat exchangers can be cooled (4) or heated (5) by a cold or hot water flow respectively in order to impose a gas temperature gradient T between the two ends of the milli-regenerator and 0 < T < 70°C. The pressure (P1, P2), velocity (V1, V2) and temperature (T1, T2) of the airflow are measured at both ends of the milli-regenerator. The diameter of the connecting pipes between the heat exchangers and the regenerator is 5 mm.

Uncertainty analysis

We followed the method presented in the Guide to the expression of Uncertainty in Measurement (GUM) to determine the experimental uncertainties [12]. The pressure and velocity measurement are the average of 40 successive cycles. The maximal type A uncertainties are observed at the lowest frequency (2 Hz), they were of 0.022 m.s -1 or 1.45% for the maximal velocity at this frequency and 8.75x10 -5 bar or 1.48% for the maximal pressure drop at this frequency. This worst-case uncertainty is used in the calculation of the total uncertainty. For the type B uncertainties, we use the supplier characteristics of the calibration apparatus for the velocity and pressure (resp. +/-2% and +/-2.5x10 -4 bar). The total relative velocity uncertainty varies from 3.7 % at 2 Hz to 2.4% at 10 Hz. The relative pressure uncertainty varies from 5.7 % at 2 Hz to 3.9 % at 10 Hz. The temperature sensors are calibrated with a reference probe and the uncertainty is +/-0.1°C. We consider that the uncertainty in the thermophysical properties due to the temperature measurement uncertainty is small compared to the other uncertainties. The dimensions uncertainties are 0.2% due to the resolution of the 3D printer. The porosity uncertainty is assumed to be 1%. The errors are supposed uncorrelated. Therefore, applying the law of propagation of uncertainties leads to relative uncertainties for the Reynolds number, pressure drop and friction factor. The relative uncertainties depend on the frequency. They vary from 3.7% to 2.6% for the Reynolds number, from 5.7 to 0.4% for the pressure drop, and for the friction factor, from 7.8% to 3.6%.

Results and discussions

We performed experiments on different milli-regenerators with three porosities ( = 0.80, 0.85 and 0.90) under nitrogen steady and oscillating flows (oscillating Reynolds number in the range 0 < Re < 60 and Reynolds number based on the hydraulic diameter ReDh < 6000) for different temperature gradients (T < 100°C).

Steady flow

A continuous controlled mass flow of nitrogen passes through the regenerator. The maximal pressure drop P is measured as a function of the Reynold number ReDh based on the hydraulic diameter Dh with ReDh < 6000 (Fig. 3). Steady flow can be considered as the zero-frequency limit of oscillating flow at maximum. Darcy's law requires correction when the particle-diameter-based Reynolds number exceeds unity. Darcy's law means that the drag is linearly proportional to the velocity. In the fluid momentum equation written for compressible flow, the permeability K and the inertial coefficient F need to be evaluated for each kind of regenerator porous medium. Vafai and Tien [13] and Pamuk and Ozdemir [14] presented semi-empiric correlation equations for fluid flow through packed bed of spherical balls based on the model of Ergun [15]. Even if the milli-regenerator is filled with pillars (not spherical balls), we consider a quadratic relation for pressure drop such as the Darcy-Forchheimer quasi steady-flow form of the 1-D fluid momentum equation:

P VF V LK K    =+ 2 (1)
This equation represents the balance between the two drag forces and the pressure force needed to overcome them, to maintain average velocity in the porous medium. Considering maximal pressure drop P, Reynolds number ReDh, permeability K and inertial coefficient

F, equation (1) becomes:

h h h Dh D D F D P Re L V K K   =+ 22 ( 2
)
This equation is written in terms of the Darcy friction factor Cf for steady flows in microchannels [14 ; 16-20] as follows :

h h f Dh h D FD P K C Re V L K D   = = = + 2 2 2 2 2 (3)
Friction factor Cf correlations for various type of regenerator materials (random fiber, wire screen, microchannels, involute foils) are of the general form we based on hydraulic diameter Dh [7,16]:

c std d f ,st th h d sD D a b Re Re C =+ (4)
Equation ( 4) is a three-term composite correlation equation that best-fit the experimental data in steady flow. We define the friction factor Cf,std, the permeability Kstd and the inertial coefficient Fstd :

h std std h f ,std Dh std D K F D C Re K = = + 2 2 2 (5)
The coefficients Kstd and Fstd are obtained from Equations (2) to (5) and: The three constants a, b and c are identified from experimental data [21]. Figures 4a,b,c,d show the friction factor after constants identification (Tab. 4).

Oscillating flow

In order to compare results with steady flows, we conducted experiments in oscillating flows with oscillating Reynolds number Re= Dh 2 / the range 0 < Re < 60 for different temperature gradients (T < 100°C), based on the hydraulic diameter ReDh < 6000. The exchangers are temperature-controlled and keep a temperature gradient between the two ends of the regenerator of ∆T = 30°C, 40°C and 50°C. The strokes of the pistons are S = 24 mm and S =30 mm. The different tests were made for five frequencies: 2, 4, 6, 8 and 10 Hz. Figure 5 shows the principle of measurement acquisition (Fig. 2). From 0 to 180°, the fluid leaves the cold exchanger (CHX) and enters the regenerator: this is the "Cold Blow" phase. Conversely, from 360° to 180°, the fluid leaves the hot exchanger (HEX) and enters the regenerator: this is the "Hot Blow" phase. Normally we should have negative velocity oscillations from 0 to 180° and positive from 360 to 180°. On figure 5, these oscillations are positive because the velocity sensor measures only positive quantities. There is a phase shift of 18° between the piston and the fluid. This phenomenon is caused by the inertia of the compressible fluid at each change of direction. We have observed that this phase shift varies between 18° ( = 0.90) and 22° ( = 0.90) depending on the porosity. Figures 5 and6 show that the maximum velocity occurs when the pressure drop is maximum and the piston is zero. Figure 6 shows the influence of temperature gradients on the pressure drop. It can be seen that the temperature gradient only affects the Cold Blow phase when the fluid reaches its highest negative pressure change (at angle 120°). This slight expansion changes the thermophysical properties of the fluid (density and viscosity) enough to create a pressure change. During the Hot Blow phase, the fluid is sufficiently hot and the slight compression (at angle 290°) does not generate a sufficiently large overpressure variation. Usually, the temperature gradient between 30°C and 70°C does not affect the pressure drop. We can now continue the study from the measurements made in isothermal regime for T = 0°C.

The rotation frequency changes the pressure drop within the regenerator (Fig. 7). The pressure drop increases with the square of the fluid speed and therefore with the frequency. The inertial effects coupled with the compressibility of the fluid create a phase shift of about 22° between the maximum pressure drops at 2 and 10 Hz (for the Hob Blow and the Cold Blow periods). The viscous frictions decrease with the porosity of the regenerator (Fig. 8). Indeed, the greater the porosity, the smaller the contact surface area with the fluid. A slight phase shift of around 10° is observed between the extreme values of pressure gradient. As the porosity decreases, the effects of compressibility and friction increase and tend to delay the appearance of the pressure peak.

In oscillating flows, velocity, pressure and temperature vary with time and it is obvious that friction factor varies with time too. In order to compare our results with literature, we consider velocity and Reynolds number ReDh,max corresponding to maximum pressure gradient value Pmax (Fig. 9). According to Eq. ( 4), the friction factor Cf,osc will be based on the maximal value of the velocity measured on a period of oscillation, i.e. the maximal value of the Reynolds number ReDh,max: Based on equations ( 4) and ( 6), the equivalent oscillating permeability Kosc and inertial coefficient Fosc at maximal Reynolds number ReDh,max become: Permeabilities in steady and oscillating flows, Kstd and Kosc respectively, are of the same magnitude order (from 2.10 -10 to 5.10 -10 m 2 ) and 0.66 < Kosc/ Kstd < 1.37. In oscillating flows, the inertial coefficients Fosc are larger than in steady flows and Fosc/ Fstd vary between 57 at ReDh,max = 100 and 72 at ReDh,max = 5000. For high values of Reynolds number, oscillating flows provide the highest inertial effects because of the higher velocity through the regenerator.

Conclusion

Gaseous flow characteristics of nitrogen in milli-regenerators were studied. Each regenerator, fabricated with a multiple jet molding process, is a cylinder with a rectangular channel with a matrix of pillars with three porosities ( = 0.80, 0.85 and 0.90). The regenerator is a porous medium with a total length L = 60 mm and an internal diameter D = 5 mm filled with a dense pillar matrix. The pillars have a geometrical lens shape with a form factor FF = 0.5 (aspect ratio width/length) and a height h = 2.5mm. Two metallic layers (chromium and copper) with a total thickness of 800 nm are deposited on the polymer pillars to increase heat transfer inside the regenerator. The regenerators were tested under nitrogen steady flows and oscillating flows (oscillating Reynolds number in the range 0 < Re < 60 and Reynolds number based on the hydraulic diameter ReDh < 6000) for different temperature gradients (T < 100°C). The major findings can be drawn as follows.

In steady flow: -We identified a three-term composite correlation equation for the friction factor based on a Darcy-Forchheimer flow model that best-fit the experimental data: In oscillating flow: -The compressibility and inertia of the gas affect the characteristics of the oscillating flow. A phase shift occurs between the velocity of the piston and the velocity of the gas. The maximum gas velocity is delayed by an angle of about 18° with respect to the maximum piston velocity. Similarly, as the frequency of rotation increases, the pressure maxima are shifted by an angle of about 22°.

-The temperature of the gas (between 30°C and 70°C) do not affect the pressure oscillations.

-We have adopted the same kind of friction coefficient correlation as for the steady flow. The Reynolds number ReDh,max corresponds to the maximal velocity obtained in an oscillating cycle, at the maximal pressure of the flowing gas.

Fig. 1 .

 1 Fig. 1. Milli-regenerator structures made by Multiple Jet Molding (Length L = 60 mm, internal diameter D = 5 mm, porosities  = 0.90 and  = 0.85.

Fig. 2 . 1 :

 21 Fig. 2. Experimental setup -1: Electric motor; 2 and 3: Pistons; 4: Cold Heat Exchanger (CHX); 5: Hot Heat Exchanger (HHX); 6: Milli-regenerator; 7 and 8: Data acquisition system Gas pressures are measured with ultraminiature pressure transducers (Kulite XCQ-055, 1.7 BARA-8068 and bandwidth is 210 kHz). They were calibrated using a Drück PV621 Pressure Station within the pressure range [Patm-3 bar]. The velocity measurements are achieved with hot wire anemometers (TSI IFA300 range 0.15-200 m.s -1 , 600 kHz bandwidth). Each probe measures the axial velocity component and was calibrated in the laboratory. The fluid temperatures are measured with home-made 12.7 μm diameter type K microthermocouples (accuracy ± 0.1 °C and cut-off frequency 30 Hz). The static response is not presented in this work but the probes were placed in a regulated oven within the range [30-300°C] and the measured temperatures were compared with a precision thermometer (Pt100 platinum reference resistance, accuracy ± 0.005 °C). Each thermocouple presents a Seebeck coefficient of 40 µV/°C at 20°C. The experimental dynamic characterization of the sensors consists in the determination of their temporal and frequency characteristics (response time, cut-off frequency and phase shift of signals) for different types of convective excitations in non-stationary regimes [11]. The measured signals are registered with the data acquisition system (7 and 8) operating at 200 kHz maximum sample rate.

Fig. 3 .

 3 Fig. 3. Maximal pressure drop P for a steady unidirectional flow versus Reynolds number ReDh

Fig. 4 .

 4 Fig. 4. Friction factor Cf,std for the three regenerator porosities ( = 0.90_std ; 0.85_std and 0.80_std respectively) in steady unidirectional flow versus Reynolds number Results indicate that the inertial coefficient Fstd of the compressible gas in steady flow decreases when the porosity increases and, for each porosity, increases with the Reynolds number ReDh, (Tab. 2).

Fig. 5 .Fig. 6 .

 56 Fig. 5. Friction gas velocity at the hot outlet of the regenerator ( = 0.80, S = 30 mm, isothermal flow at 22°C) Fig. 6. Influence of the temperature gradient on the pressure drop ( = 0.80, S = 30 mm, T=0 (isothermal flow), 30°C, 40°C and 50°C respectively).

Fig. 7 .Fig. 8 .

 78 Fig. 7. Pressure drop of the regenerator fluid flow at different frequencies ( = 0.80, S = 30 mm, isothermal flow T = 0°C)

Fig. 9 .

 9 Fig. 9. Maximal pressure drops Pmax for an oscillating flow versus maximal Reynolds number ReDh,max

Fig. 10 .

 10 Fig. 10. Friction factor Cf,osc for the three regenerator porosities ( = 0.90_osc ; 0.85_osc and 0.80_osc respectively) in oscillating flow versus maximal Reynolds number ReDh,max

2 -

 2 coefficients Kstd and inertial coefficient Fstd are: The inertial coefficient Fstd of the compressible gas decreases when the porosity increases and increases with the Reynolds number ReDh for each porosity. The steady flow show permeability coefficient of the same magnitude order.

  =+where the permeability coefficients Kosc and inertial coefficient Fosc are:

2 -

 2 The inertial coefficients Fosc are larger than in steady flows and increase with the Reynolds number with an amplitude much larger.

Table 1 .

 1 Thermophysical properties of the regenerator material.

Table 2 .

 2 Friction factor coefficients, permeability K and inertial coefficient F for the regenerators at steady flow.

	Porosity	Hydraulic diameter Dh (µm)	a	b	c	Permeability Kstd (m 2 )	Inertial coefficient Fstd ReDh = 100 ReDh = 5000
	 = 0.80	268	392.2	0.8181	0. 0614	3.6626 x 10 -10 1.384 x 10 -3 1.759 x 10 -3
	 = 0.88	445	980.7	0.0547	0.3532	4.0384 x 10 -10 1.418 x 10 -4 5.648 x 10 -4
	 = 0.90	622	2589.3	0.000027	1.1448	2.9980 x 10 -10 1.034 x 10 -6 91.14 x 10 -6

Table 3 .

 3 Friction factor coefficients, permeability Kosc and inertial coefficient Fosc for the regenerators at oscillating flow.

	Porosity	Hydraulic diameter Dh (µm)	aosc	bosc	cosc	Permeability Kosc (m 2 )	Inertial coefficient Fosc ReDh,max ReDh,max
							= 100	= 5000
	 = 0.80	268	574.5	8.428	-0.2503	2.5004 x 10 -10 78.52 x 10 -3 29.49 x 10 -3
	 = 0.88	445	721.6	13.707	-0.314	5.4885 x 10 -10 849.7 x 10 -4 248.8 x 10 -4
	 = 0.90	622	3832	0.000011	1.275	2.0192 x 10 -10	44.5 x 10 -6	6537 x 10 -6

  This work has been supported by the EIPHI Graduate School (contract ANR-17-EURE-0002) and the Region Bourgogne-Franche-Comté

	C	Factor	Subscripts
	CHX	Cold Heat eXchanger	atm	atmospheric
	D	diameter (m)	Dh	hydraulic diameter
	E	effusivity (J.K -1 .m -2 .s -1/2 )	f	friction
	f	frequency (Hz)	max	maximal value
	F	inertial coefficient	osc	oscillating
	FF	form factor	std	steady
	GUM	Guide to the expression of
		Uncertainty in Measurement
	h	height (m)	Greek Symbols
	HHX	Hot Heat eXchanger		difference
	K	permeability (m 2 )		porosity
	L	length (m)		thermal conductivity (W.m -1 .K -1 )
	P	pressure (bar)		cinematic viscosity (m 2 .s -1 )
	Re	Reynolds number		density (kg.m -3 )
	S	stroke (m)		angular frequency (rad.s -1 )
	Nomenclature
	a, b, c	constants	T	temperature (°C)
	BCD	Bottom Dead Center	TDC	Top Dead Center
	cp	specific heat (J.kg -1 .K -1 )	V	velocity (m.s -1 )


Remi Mahmoud

Pierre Casadebaig

Nadine Hilgert

Noemie Gaudio

Intercropping for sustainability: Research developments and their application

Keywords: Complementarity effect, dominance effect, biodiversity effect, N fertilisation, intercropping

proposed a framework to estimate the net biodiversity effect (over-or under-yielding) on the productivity of plant communities and to further decompose this effect into dominance and complementarity effects. Here we applied this method to analyse the performance of eight cereal-legume intercropping experiments grown in diverse environmental conditions in France, with three combinations of species mixtures: i) durum wheat/faba bean ii) soft wheat/pea and iii) durum wheat/pea. For each species mixture, we tested whether N fertilisation influenced dominance or complementarity effects using a Bayesian approach. Our results showed that the Biodiversity Effect (BE) was positive in the three intercrops under unfertilised conditions, mainly thanks to complementarity effect. The relative importance of BE and both of its components depends on the species intercropped. Finally, we showed that the impact of N fertilisation on BE of intercrops strongly depends on the competitive ability of the species intercropped.

Introduction

Plant diversity promotes primary production through complementarity and dominance effects, both involved in the net biodiversity effect (BE; [START_REF] Loreau | Partitioning selection and complementarity in biodiversity experiments[END_REF]. Dominance effect measures the part of the BE due to the excessive competition of one given species, whereas the complementarity effect measures the part of the BE due to niche complementarity and/or facilitation. In intercrops, these effects were already demonstrated, highlighting that productivity is higher in cereal-legume intercrops, mainly due to complementarity effect (e.g. [START_REF] Pelzer | Pea -wheat intercrops in lowinput conditions combine high economic performances and low environmental impacts[END_REF] on durum wheat/pea intercrops). However, the influence of the choice of the species intercropped or of the agronomic practices on these effects is not well understood.

In this study, we aimed to quantifying the BE and its components (i.e. complementarity and dominance effects) in three cereal-grain legume intercrops in unfertilised conditions (durum wheat/ pea, soft wheat/pea and durum wheat/faba bean), and to test whether differences exist depending on the species intercropped. For this purpose, we used a database containing eight field experiments located at two sites (Angers and Auzeville) in France. We then studied how nitrogen (N) fertilisation influenced these effects.

Materials and Methods

Field experiments

In order to explore the net BE in intercrops, we gathered eight factorial experiments in two locations in France (Angers, Auzeville), covering a wide range of climate conditions and agricultural managements. The data were formatted [START_REF] Wickham | Tidy Data[END_REF], homogenised and pooled in the same database.

Environmental conditions

Climate conditions of each experiment were characterised using variables retrieved from the NASA POWER API, i.e. the sum of precipitation (mm) and the mean temperature (°C) during the crop cycle (from sowing to harvest date). The experiments involved winter crops, with precipitation ranging from 286.4-712.5 mm and mean temperature ranging from 7.3-11.3°C during the crop cycle.

Agricultural management

All the experiments included cereal-grain legume intercrops of two annual crop species and their corresponding sole crops for which grain yield was measured (t ha -1 ). Cereals were represented by two species: durum wheat (Triticum turgidum L.) and soft wheat (Triticum aestivum L.); legumes were represented by two species: faba bean (Vicia faba L.) and pea (Pisum sativum L.). Finally, we had three intercrop combinations: durum wheat/faba bean, soft wheat/pea and durum wheat/pea. Intercropped species were sown and harvested at the same time. Within an experiment, the different treatments differed in i) the number of genotypes cropped by species (from one to five genotypes), ii) the relative sowing densities of each species, sole crops being characterised by densities of 1 or 0.5 relatively to the recommended density, and intercrops being characterised by densities ranging from 0.33-1 and iii) the nitrogen (N) fertilisation, with unfertilised and N-fertilised treatments (from 30-180 kg N ha -1 ). 89% of the intercrops were grown in a substitutive design (i.e. when the relative sowing densities of the two species intercropped sum to 1) and 11% of the intercrops were grown in an additive design (i.e. when the relative sowing densities of the two species intercropped sum higher than 1). Considering all these variations, the database contained 70 intercrop experimental units (site × year × mix of genotypes × relative densities × N-treatment), among which eight were in additive design and 62 in substitutive design, and 69 sole crop experimental units. The dataset was unbalanced because the experiments gathered were made for different purposes (i.e. different number of observations between groups) regarding many factors (N-treatment, mixture design, etc.). This heterogeneity constrained the statistical analyses we made.

Components of intercrop performance

For each experimental unit, we computed BE as the difference between the observed and expected yield in intercrops (Eq. 1):

BE = YO C + YO L -(YE C + YE L ) (Eq. 1)
where YO C (resp. YO L ) is the observed yield of the cereal (resp. legume) grown in intercrop, and YE C (resp. YO L ) is the expected yield of the cereal (resp. legume) grown in intercrop.

Expected yield was assessed from the yield of the species in sole crop weighted by the scaled relative density in intercrop (Eq. 2):

YEC = MC

and YEL = ML (Eq. 2)

where Mc (resp. M L ) is the yield of the cereal (resp. legume) grown in sole crop and RD C (resp. RD L ) is the relative density (i.e. the density of each species in the intercrop relatively to the sowing density in the reference sole crop) of the cereal (resp. legume) in intercrop.

In the case of substitutive designs, the expected yield is simply for cereal (resp. for legume):

BE can be partitioned into dominance effect (DE, Eq. 4) and complementarity effect (CE, Eq. 5) [START_REF] Loreau | Partitioning selection and complementarity in biodiversity experiments[END_REF][START_REF] Li | Yield gain, complementarity and competitive dominance in intercropping in China: a meta-analysis of drivers of yield gain using additive partitioning[END_REF], where:

A positive DE means that the species with the greater yield in sole crop benefits the most from intercropping. On the opposite, a negative DE means that the species with the lowest yield in sole crop has a higher relative yield gain. In the definition of CE, we introduced the classical Land Equivalent Ratio [START_REF] Willey | A Competitive Ratio for Quantifying Competition Between Intercrops[END_REF]. Hence, CE is equal to the LER minus 1, multiplied by the average yields in sole crops. Thus, the CE measures the part of the BE due to niche complementarity and or facilitation, without being able to disentangle these processes [START_REF] Loreau | Partitioning selection and complementarity in biodiversity experiments[END_REF] .

Data processing and analysis

All analyses were done under the Bayesian approach. Bayesian inference is based on the reallocation of credible values for a parameter (posterior distribution), given the prior knowledge (prior distribution) and the adequacy of the data to the model (likelihood). Under the Bayesian approach, it is possible to get information about the probability of a hypothesis being true given the data (P (hypothesis/data)). Bayesian estimation for the difference of group means [START_REF] Kruschke | Bayesian estimation supersedes the t test[END_REF] is an alternative to the classical Student's t test to compare the means of two groups. This method allows to get a posterior distribution for the values of mean differences between the two groups, and to derive a 95% highest density interval (HDI), described as the 95% most credible values of the parameter [START_REF] Kruschke | Rejecting or accepting parameter values in Bayesian estimation[END_REF].

For each species mixtures, we carried out Bayesian estimation for the difference of group means between the values of the components of the BE in fertilised and unfertilised treatments. We defined the null hypothesis (H0) as the hypothesis of equality of the mean BE components between fertilised and unfertilised intercrops. We applied the following decision rule according to the position of the null value (0) relatively to the position of the 95% HDI: reject H0 if 0 is completely outside the 95% HDI vs do not reject H0 if 0 falls inside the 95% HDI.

All indicator computations and statistical analyses were carried out with R 4.0.0 software (R Core Team, 2020). The data was ordered, reshaped and homogenised using the collection of R packages tidyverse [START_REF] Wickham | Welcome to the tidyverse[END_REF]. Bayesian statistical analyses were performed using the R package BEST [START_REF] Kruschke | BEST: Bayesian Estimation Supersedes the t-Test[END_REF].

Results

Species choice influenced BE and its components in unfertilised conditions

The BE was positive in 100% of the cases (Table 1). The average BE of intercropping was equal to 0.67 ± 0.11 t ha -1 (mean ± standard error). This was mainly due to the complementarity

DE = � � �� � � � � � �� � �� � ��� � � � � � � � � � �� � �� � ��� � �� � �𝑀𝑀 𝐶𝐶 � 𝑀𝑀 𝐿𝐿 � (Eq. 4) CE = � � �� � � � � � � � � � �� � �� � ��� � � � � � � � �� � �� � ��� � � � 𝑀𝑀 � � �𝐿𝐿�� � 1� (Eq. 5) LER = YEC = MCRDC (resp. YEL=MLRDL) (Eq. 3) (𝑀𝑀 � )
effect which represented in average 77% of the BE. The remaining 23% of the BE was due to the dominance effect. The relative contribution of complementarity and dominance effects in BE differed between the species intercropped (Table 1) but complementarity effect was predominant in all unfertilised intercrops. N fertilisation influenced BE and its components depending on species intercropped Biodiversity effect and its components were not influenced by N fertilisation in the same way depending on crop mixture (Fig. 1).

Considering durum wheat/pea mixtures, 0 was outside the 95% HDI (0.19; 0.74) for the values of mean BE differences between the unfertilised and N-fertilised treatments. Thus, BE was higher in unfertilised treatments, resulting from both an effect of dominance and complementarity. For the dominance effect, 0 was outside the 95% HDI (-0.5; -0.21) for the values of mean differences between the unfertilised and N-fertilised treatments, meaning that N fertilisation increased the dominance effect. For the complementarity effect, 0 was outside the 95% HDI (0.57;1.07), meaning that N fertilisation decreased the complementarity effect. Hence, in unfertilised durum wheat/pea mixtures, the complementarity effect compared to N-fertilised crops was balanced by the higher complementarity effect, leading to a greater BE.

Conversely, considering soft wheat/pea mixtures, 0 was included in the 95% HDI for the difference of means of BE between unfertilised and N-fertilised treatments meaning that N fertilisation had no significant influence on BE. Moreover, N fertilisation had no significant influence neither on the complementarity effect nor on the dominance effect, the null value being inside the 95% HDI.

In the same way, considering durum wheat/faba bean mixtures, N fertilisation had no significant impact on the BE or complementarity effect. Null value was also a credible value for the difference of means dominance effect between N-fertilised and unfertilised treatments. A trend was distinguishable in the posterior distribution of the difference of means, 93% of the posterior values of μ N being greater than μ N0 , meaning that N fertilisation tended to increase the dominance effect.

Discussion and Conclusion

In this study, we showed that the overall biodiversity effect of unfertilized cereal-legume intercrops was positive (pooling together all experimental units), mainly due to complementarity effect (77%), meaning that intercrops were more efficient considering yield than their corresponding sole crops, which is in accordance with previous findings (e.g. [START_REF] Li | Yield gain, complementarity and competitive dominance in intercropping in China: a meta-analysis of drivers of yield gain using additive partitioning[END_REF][START_REF] Pelzer | Pea -wheat intercrops in lowinput conditions combine high economic performances and low environmental impacts[END_REF][START_REF] Yu | A meta-analysis of relative crop yields in cereal/legume mixtures suggests options for management[END_REF]. However, this effect strongly depends on species intercropped and on N fertilisation management. The results of Bayesian estimations of group means showed that N fertilisation increased dominance effect and decreased complementarity effect in durum wheat/pea mixtures, which is a behaviour typically highlighted in existing literature for cereal-legume intercrops (e.g. [START_REF] Corre-Hellou | Effect of root depth penetration on soil nitrogen competitive interactions and dry matter production in pea-barley intercrops given different soil nitrogen supplies[END_REF][START_REF] Naudin | The effect of various dynamics of N availability on winter pea-wheat intercrops: Crop growth, N partitioning and symbiotic N-2 fixation[END_REF]. In unfertilised conditions, the superior competitive ability of the cereal forces the legume to rely more on its atmospheric N-fixing ability [START_REF] Ghaley | Intercropping of Wheat and Pea as Influenced by Nitrogen Fertilization[END_REF]. The complementarity for the use of N increased as the symbiotic N 2 fixation of the legume is enhanced leaving more mineral N for the cereal [START_REF] Corre-Hellou | Effect of root depth penetration on soil nitrogen competitive interactions and dry matter production in pea-barley intercrops given different soil nitrogen supplies[END_REF]. In fertilised conditions, dominance is boosted because durum wheat takes a competitive advantage over the legume. Our results showed, however, that if we change one of the species in the intercrop and the fertilisation conditions, we can move away from this textbook case.

Replacing durum wheat by soft wheat in wheat/pea mixtures, N fertilisation did not influence the components of the BE of the intercrop anymore. This could be due to a lower competitive ability of soft wheat compared to durum wheat. However, such a conclusion has to be tempered as N fertilisation levels in soft wheat/pea mixtures were lower than those of durum wheat/pea mixtures (40 kg N ha -1 vs 60-140 kg N ha -1 ), and also globally lower to conventional N fertilisation levels in wheat crops (170 kg N ha -1 ; Mercier-Poirier, 2014)). Thus, the amount of N brought by fertilisation might be insufficient to upset complementarity and dominance occurring in unfertilised conditions. Finally, when replacing pea by faba bean in the mixtures, we showed a tendency to an increase of dominance effect without decreasing complementarity effect in N-fertilised intercrops (in contrast to durum wheat/pea mixtures). The absence of impact of N fertilisation on the complementarity effect could be due to a stronger competitive ability of the faba bean compared to the pea in unfertilised conditions and/or to a lack of N availability for the cereal, even with the presence of the symbiotic fixation of the faba bean. Therefore, our study highlighted that BE was positive in the three unfertilised cereal-legume intercrops considered, mainly thanks to complementarity effect. However, the relative importance of BE and both of its components depends on the species intercropped. Using a Bayesian alternative to the t test, we showed that the impact of N fertilisation on BE of intercrops strongly depends on the competitive ability of the species intercropped.

Fig. 1 .

 1 Fig. 1. Posterior distribution of the difference of mean biodiversity effect (left panels), dominance effect (center panels) and complementarity effect (right panels) between unfertilised (μN0) and fertilised (μN) intercrops of durum wheat (Triticum turgidum L.)/pea (Pisum sativum L.), soft wheat (Triticum aestivum L.)/pea and durum wheat/faba bean (Vicia faba L.). The dotted lines represent the null value of the posterior difference of means.

Table 1 .

 1 Biodiversity effect, complementarity effect and dominance effect in all intercropping conditions, and especially in intercrops of durum wheat (Triticum turgidum L.)/pea (Pisum

	sativum L.), soft wheat (Triticum aestivum L.)/pea and durum wheat/faba bean (Vicia faba
	L.), in unfertilised conditions (mean ± standard error). The relative contribution (%) of
	complementarity and dominance effects to the biodiversity effect was precised. For each
	experimental unit, we averaged yield values across all replicates of the treatment
		(two to four replicates by treatment)
		Biodiversity	Complementarity	Dominance	No. experimental
		effect (t ha -1 )	effect (t ha -1 )	effect (t ha -1 )	units
	All intercrops	0.67 ± 0.11	0.83 ± 0.11 (77%)	-0.17 ± 0.05 (23%)	28
	Durum wheat/pea	0.52 ± 0.08	0.75 ± 0.07 (77%)	-0.22 ± 0.06 (23%)	18
	Soft wheat/pea	1.21 ± 0.51	1.03 ± 0.61 (69%)	0.18 ± 0.17 (31%)	4
	Durum wheat/faba bean	0.76 ± 0.21	0.99 ± 0.29 (83%)	-0.23 ± 0.09 (17%)	6

Acknowledgements

This work was supported by the French National Research Agency under the Investments for the Future Program, referred as ANR-16-CONV-0004.