M Belhache

S Guillou

P Grangeret

A Santa-Cruz

F Bellanger

Fluid Structure Interaction of a loaded Darrieus Marine Current Turbine

Keywords: Vertical Axis Marine Turbine, Fluid Structure Interaction, Power Conversion System, Numerical Simulation, Wind Tunnel experiment

The present study aims to develop a model to simulate in 2 dimensions, the motion of the fluid and of the loaded Darrieus turbine without imposing the rotational speed. A resistive torque has been added to account for the load of a generator. Experimental and numerical studies were conducted and presented here. The experiment was set-up in the air first, on a wind turbine for reason of convenience. The simulations were made with Fluent (v14) for which a User Defined Function was developed. Then the experimental conditions were reproduced numerically for three cases. The study focuses on the determination of a factor (f) giving the relationship between the experimental and the 2 dimensions simulations resistive torque. The results are quite satisfying in this preliminary study as it demonstrate that the model developed works well. Some improvement has to be made to refine the value of f .

Introduction

A. Context

Owing to the exhaustion of classic energy resources (fossil fuels) exploited since the industrial revolution as well as the potential threat of climate change, the issue of energy supply is the challenge of the XXI century. Since the Earth summit in Rio de Janeiro in 1992, there is a consensus in the need to develop renewable energies. Renewable energies exit in various forms and some expressions can be found in the ocean. The most mature marine energy conversion system is offshore wind turbines and the next systems on the way towards commercialization are marine turbines. Harvesting tidal energy has the advantage to be predictable by opposition to the intermittency of wind energy.

Regarding marine turbine, no standard has emerged yet. If the horizontal axis has been favored for the most advanced prototypes as it could benefit from the experience and years of development invested in wind energy, the vertical axis type, patented in 1926 in France, undergo renewed interest. Indeed the vertical axis turbine presents some characteristics that are interesting for naturel flow. First, the turbine can accept flow from any direction which is convenient on account of flow reversal, besides high potential site may present rotary current [START_REF] Blunden | Initial evaluation of tidal stream energy resources at Portland Bill, UK[END_REF]. Moreover, in the marine environment flows are turbulent by nature and recent studies tend to demonstrate that the performance of vertical axis turbine are not altered by turbulent flow but on the contrary can be increased [START_REF] Simão | The near wake of the VAWT, 2D and 3D views of the VAWT aerodynamics[END_REF]. The present study focus on the numerical prediction of the behavior of a vertical axis marine turbine connected to a generator.

B. Darrieus characterisation : numerical and experimental modelling strategies

Numerical studies of marine and wind turbines always follow the same pattern: for a given fluid velocity, several rotational speed of the engine are tested in order to build the non-dimensional power coefficient curve (Cp-λ curve) that pictures the power coefficient (Cp) against the tip-speed ratio (λ= ̇/U 0 ) for a specific fluid velocity (U 0 ). The rotational speed ( ̇ is then forced and remains constant over the whole revolution of the turbine ([3]- [START_REF] Lain | Simulation and Evaluation of a straight-bladed Darrieus-type cross flow flow marine turbine[END_REF]). The curve, enable the identification of the point of maximum performance. The working condition of the turbine for the fluid velocity considered is then supposed to be somewhere on the curve depending on the load applied.

Experimental works ([6]- [START_REF] Maitre | Modeling of the flow in a Darrieus water turbine: Wall grid refinement analysis and comparison with experiments[END_REF]) apply the same procedure. In all cases a generator is used to drive the turbine at a desired rotation speed. Some questions arise from this method: does forcing numerically and experimentally the rotational speed are the right approaches? Are the force and moment (intensity and oscillations) well predicted?

Is the rotational speed constant during a rotation? Thus the objective of this study is to simulate the overall system behavior with a real interaction between the fluid, the turbine, and the load of a generator in order to better predict:

-The forces acting on the rotor for structural analysis -The turbine behavior -The power extracted and power quality In this preliminary work, numerical and experimental results will be correlated, first, on a wind turbine in "free conditions". Note that the final goal is to apply the method to a marine turbine. A vertical axis marine current turbine prototype has already been realized, it is optimized to have it driven by the fluid without any external help (Figure 1).

Experimental Set-Up

With a concern of speed and efficiency some preliminary tests were performed in an open wind tunnel on a threebladed Darrieus wind turbine (fig. 2). The experimental work has been favored in air first as it is easier to have the turbine working in free conditions with higher Reynolds number. Experiment in the fluid water brings more constraints; owing to water density the maximum velocity in the channel is reduced. As a matter of fact the size of the prototype is limited. Thus the prototype operates at lower Reynolds numbers and the side effects make it more difficult to have the turbine forced by the fluid.

A. Turbine characteristics

The turbine is a three-bladed vertical axis turbine of 1 meter diameter (D). Blades are straight, 1 meter high and show a curved profile with a chord length of 0.21 m (c). The core of the blade is made of acetal resin covered with an external envelope of Carbon (fig 3). The rotating shaft extends to the hub situated at rotor mid-height. Each blade is attached to the hub with two aluminium rigid arms that forms a trapezoid. The generator installed is a reversible permanent magnet synchronous motor for electrically assisted cycles

B. Synchronous generator characterization

The efficiency of a generator is not constant but varies with the level of torque transmitted to the system and with the rotational velocity owing to iron and mechanical losses. So the first step was to characterize the generator according to the mechanical power developed and the angular velocity of the turbine. The generator was evaluated on a test rig where it was driven by a motor shaft at constant speed. The motor rotated at a constant velocity and supplies a mechanical power P m , and the rotation of the generator produced an electrical power P e .

The efficiency was thus determined as the ratio of P e to P m . The performance curves of the generator according to the mechanical power for various rotational speeds are given on fig 4.

C. Wind Tunnel Experiment

The experimental test took place in an open wind tunnel (fig. 2) of 1m x 1m in section. The tunnel is 12 m long with honeycomb in order to laminerize the flow. The wind turbine was positioned on a base so that its axis was 0.6m from the tunnel outlet. The objective of the tests was to have different couple U 0 /P m corresponding to different generator efficiencies () with the following target values for the mechanical power: 100, 300 and 500W. For each test, a specific air velocity at the outlet of the tunnel is asked to the command system. The rotational velocity ( ̇) and the electrical power produced (Pe) were measured. The mechanical power (P m ) developed by the turbine is deduced from the two latter values and the generator characterization curves. A dimmer was then used to accelerate or slow down the turbine to have the duet ̇/Pe matching with the target mechanical power. Table I summarise the experimental test conditions. Flow velocity at the outlet of the wind tunnel was not recorded and the level of turbulence is unknown.

Modelling parameters

A. Domain and Grid Setting

The domain is 2 dimensional equivalent to a plan situated mid-height of the turbine. It extends 4D upstream and 14D downstream. The lateral sides are 2.8D distant from the rotation axis, which correspond to a blockage ratio of 0.18. The computational domain was divided in two parts, a stationary sub-domain and a rotating sub-domain separated by a sliding interface to enable the rotation of the turbine. The rotating sub-domain is a 3D circle including the turbine of which the position is updating at each time step.

The whole domain was discretized with 503318 elements composed of 13570 elements inside the stationary subdomain and 489750 in the rotating domain. The vertical axis turbine hydrodynamic is very complex, especially in stall regime: attached and separated flow as well as reattachment may occur simultaneously. Moreover the blades can interact with their own wake and the one created by the other elements of the rotor. On account of this unsteady nature, a fine discretization inside the rotor is essential, particularly around the blade. To capture the flow feature inside the boundary layer, a boundary layer cmesh was built around each blade (fig. 5). It develops half of a chord length normally to the blade and is composed of 56113 cells shared out in 120 layers.

B. Calculation parameters

An Unsteady Reynolds Average Navier-Stokes (URANS) approach was chosen for the turbulent modeling. For turbulence closure, the k-ω SST Model was used owing to its capability to deal with the specific flow features related to the operation of vertical axis turbine (strong adverse pressure gradient and flow reversal) [START_REF] Vermeer | Wind turbine wake aerodynamics[END_REF]. The velocity-pressure coupling corresponds to the SIMPLEC method and the spatial discretization is of second order. A symmetry condition was applied on the wall domain, the inlet is defined as a constant velocity inlet with 5% turbulence intensity and constant pressure condition was imposed on the outlet. Blades were defined as wall with no slip condition. The time step defined for each test is given in table II, it corresponds to a rotation angle of 0.75°. The relative convergence criterion is set to 0.05 with a maximum of 100 iterations per time step.

For each tests, only the velocity inlet and resistive torque were imposed, all other values were calculated. The velocity ordered to the wind tunnel control console was corrected with 3-d anemometer measurements. Each simulation was initialized with an imposed constant angular velocity case. The angular velocity was chosen close to the target one in order to reduce the computational time. The initialization cases were run on 5 revolutions until the forces and moment reach a periodic behavior. The Reynolds number and tip speed ratio corresponding to each case are given in table I.

C. Modeling strategy

In order to simulate the fluid structure interaction of a three-bladed vertical axis marine turbine connected to a generator with a defined resistive torque; a User Defined Function (UDF) has been developed in Ansys Fluent (v.14). At each time step the function gets the calculated force and moment applied on the blade and resolves equation 1 to determine the turbine rotational velocity ( ̇). Once the rotational speed is computed, it is passed on fluent for an update of the mesh position before the next time step. 6). The resistive torque of the generator (Cr) is defined as the ratio between the mechanical power and the rotational velocity derived from the experimental test multiply by a factor f.

The computation is 2 dimensions which means that there will be some differences in the forces computed with experiment as some elements of the turbine geometry are not represented (axe, support arms…) and 3 dimensions phenomenon are not taken into account. For each tests a first run was performed with the resistive torque from the experiment (f=1). If the velocity computed has not stabilized after 20 revolutions, f was increased or reduced depending on the slope. This process was repeated further until the velocity reaches a value showing an acceptable difference with the experimental one.

Results and Discussion

15 simulations, were performed to assess the optimum value of f in each test. The latter represents the difference between the simulated and measured mechanical power. Figure 7 illustrates the behavior of the turbine in terms of rotational velocity to the changing values of f. For small values of f (curves A to C), the mean mechanical torque developed by the turbine is large in relation with the resistive torque. In consequence the turbine can accelerate and the mechanical torque decreases revolution after revolution to find an equilibrium with Cr in accordance with equation (1). On the contrary, for the greatest values of f, the resistive torque imposed remains higher to the mechanical so that the rotational velocity never stabilizes. The more important is the difference between the resistive and the simulated mechanical torque, the more significant is the acceleration or deceleration of the turbine. By proceeding in the manner described in section 3, the optimal value of f was determined for tests 1 to 3 (table III). Figure 8A shows the stabilized curves. The turbine reached a steady state between the 4 th and 6 th revolutions.

The rotational speeds calculated are around 0.1% higher than the experimental one which is very close. This discrepancy is in the same order of magnitude than the uncertainties of the measurements realized during the experiment.

Once the resistive torque balances with the mechanical torque, the latter does not evolve further. It regulates to a mean constant value, the rotational speed follows the same behavior. However a slight deviation of the curve is observed.

The optimum values of f are fairly close for the 3 tests. More tests should be investigated to obtain a clear tendency for f. It is not possible to say if it is possible to predict the complete behavior of the turbine in "free conditions with a fixed value of f or if it should depend on the operating point. The tests correspond to different tip speed ratio (λ) that may have an influence. Moreover, the form of the Cp curve for the various inlet velocities is not known; the position of the operating point on this curve may be a parameter to consider as it has been demonstrated for the 2 dimensions simulation of the power coefficient [START_REF] Maitre | Modeling of the flow in a Darrieus water turbine: Wall grid refinement analysis and comparison with experiments[END_REF].

These results highlight the pulsating nature of the rotation velocity (Figure 8B). This is consistent with the theory on Darrieus rotor. As the blades rotate around the rotational axis, the angle of attack changes resulting in variation of the hydrodynamic forces applied on the blades during a revolution. Two breaking zones and two acceleration zones can be identify inside the rotor. These oscillations have the same frequency than those of the mechanical torque, 3 oscillations per revolution which is in relation with the number of blades. The amplitude of these oscillations are very low (table III) and cannot modify critically the mechanical power in comparison to a constant rotational speed case. One can however wonder if the influence of these variations is still negligible in the fluid water with the density difference. To improve the evaluation of f, the experimental set-up has to be review to perform the tests in more controlled conditions. Indeed a certain number of uncertainties exist mainly related to the flow pattern experience by the turbine. Indeed, the wind tunnel is open on the outside so the turbine may be influenced by external wind; The second fan situated on the right side of the main duct could have an effect. Recording the level of turbulence could dissipate some uncertainties. This may explain the difference in the optimal f value found for the 3 tests.

The Cp-curve of the turbine should be determined numerically and experimentally for the fluid velocity considered. Different working conditions on these curves should be studied to investigate the sensibility of f to λ and to the position on the Cp curve.

Conclusion and Perspectives

In this preliminary study, experimental tests were conducted on a vertical axis wind turbine loaded with a generator of which characteristics were determined. Numerical simulation were performed in Fluent (v14) in the same conditions with the means of a User Defined Function developed to account for the fluid structure interaction of the turbine as well as the resistive torque of a generator. The comparison of the two set of data gave promising results and interesting clue to enhanced the method and assess f in a more precise way in order to predict the behavior of a Darrieus turbine connected to a generator in 2 dimensions. The validity of the UDF function developed for this study has been demonstrated as the turbine reacts as expected (evolution of the rotational velocity and mechanical torque) to the various resistive torques imposed. Some improvement will be done in particular to correct the slight deviation of the rotational speed when the mechanical torque is stabilized. First, it concerns the order of the scheme used. The sensibility of the factor f used in the 2 dimensions simulation to various parameters will also be investigated.

Fig 1 .

 1 Fig 1. Vertical axis marine current turbine prototype for water channel test.

Fig 2 .

 2 Fig 2. Experimental wind tunnel

Fig 3 .

 3 Fig 3. Example of the cross section of a NACA0012 blade builds with the same material.

Fig 4 .

 4 Fig 4. Generator performance behavior in relation to the mechanical power for different rotational speed. The points represent experimental data and the curves are fitted laws.

Fig 5 .

 5 Fig 5. Boundary layer C-mesh around a blade (left) and rotating sub-domain mesh (right)

Fig 6 .

 6 Fig 6. Parameterization of the turbine motionBlade (I 1 ) and generator inertia are taken into account for the computation, but blades arms inertia is assumed negligible. Generator inertia is split in two parts, I 2 and I 3 , representing respectively stator and rotor inertia. The coefficient 0.23 is the speed-reduction gear ratio between the rotor and the stator. Ӫ denotes the acceleration of the turbine, m and r are the mass of a single blade and the radius of the turbine. ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ is the moment resulting from the hydrodynamic forces acting on each blade related to point O, origin of the coordinate system (fig.6). The resistive torque of the generator (Cr) is defined as the ratio between the mechanical power and the rotational velocity derived from the experimental test multiply by a factor f.

Fig 7 .

 7 Fig 7. Evolution of the rotational velocity for decreasing values of the resistive torque (A-D). The straight dash line (E) represents the target velocity corresponding to the experimental data.

Fig 8 .

 8 Fig 8. Rotational velocity ratio ̇ ̇ ⁄ . The time lap corresponds to one revolution for test 1. The arrows identify one revolution for test 2 and 3.

Table I .

 I -Experimental test conditions

	Test	U 0 (m/s)	̇ (rad/s)	Pe		P m
	1	7.13	26.3	54	0.54	100
	2	9.9	44.6	210	0.7	300
	3	10.8	55.6	210	0.7	300

Table II

 II

			. -Numerical parameters
		U 0 (m/s)	 m (rad/s)
	Test	imposed on the	used for the initialization	∆t (s)	λ	Re c
		inlet	case
	1	7.13	26.3	5.10 -4	1.84 2.7.10 6
	2	9.9	44.6	3.10 -4	2.25 4.7.10 6
	3	10.8	55.6	2.5.10 -4 2.57 5.38.10 6

Table III .

 III -Main parameters extract from the numerical results

	Test	f op	̇ (rad.s -1 )	̇ (rad.s -1 )	oscillation
	1	0.9	26.33	26,3	±0,029 %
	2	0.5	44,66	44.6	±0,0174 %
	3	0.5	55,63	55.6	±0,0135 %

Acknowledgement

This study was founded in association by the Low Normandy region and the society WIBEE (Wind Building Engineering). We would like to thank the Syndicat Mixte du Cotentin for financing the computational resources.

A B


Aurélien Gourrier
email: aurelien.gourrier@cnrs.fr

Chenghao Li

Stefan Siegel

Oskar Paris

Paul Roschger

Klaus Klaushofer

Peter Fratzl

Stehphan Siegel

the mineral nanoparticles in fluorotic bone using a stack of cards model

Keywords: Bone, Nanoscale, Ultrastructure, Mineral particles, Synchrotron radiation, Microbeam, Small-angle scattering, SAXS, Scanning X-ray scattering imaging, Modelling

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The structure of bone tissues at the nano-scale is currently receiving considerable attention. This partly results from recent findings showing that the clinical indicators measured at higher hierarchical levels and currently used to characterize various pathologies are, in some cases, insufficient to establish a reliable diagnosis. More specifically, there is growing evidence supporting the view that the micro-architecture and bone mineral density (BMD) alone are not sufficient to determine the risk of fracture in elderly people (Watts, 2002). Therefore, bone material quality is increasingly recognized as a major determinant of bone fragility (Felsenberg & Boonen, 2005;Seeman & Delmas, 2006). A better understanding of the fundamental structure of the tissue is thus required in order to improve medical treatments and therapies.

At this level, the structure of bone is quite heterogeneous due to complex physiological remodelling activities occurring throughout the lifetime of the tissue. Mineralized collagen microfibrils of ~ 100 nm in diameter act as building blocks, packing into complex arrangements of larger scale structures in the form of bone packets or osteons with characteristic dimensions of 100 µm in diameter (see, e.g., [START_REF] Currey | Bones: Structure and Mechanics[END_REF]Fratzl & Weinkamer, 2007). There is a general consensus over the fact that the mineral phase within and around the fibrils is found in the form of thin calcium phosphate platelets of nanometer dimensions (Weiner & Traub, 1986, 1989;[START_REF] Landis | [END_REF]Fratzl et al., 2004). The size of these nanoparticles in the longitudinal and transverse directions is generally estimated to ~ 50 nm and ~ 25 nm respectively in mature human bone (Robinson, 1952;Weiner & Traub, 1992). The causes for the dispersion in the measured values most likely reside in the constrained mineralization inside the collagen matrix [START_REF] Hodge | Aspects of Protein Structure[END_REF][START_REF] Glimcher | Treatise on Collagen[END_REF] as well as in the difficulty in sample preparation and analytical interpretation (Ziv & Weiner, 1994). However, most studies also indicate that the distribution in thickness is relatively narrow, usually falling in the range 2 -6 nm (Weiner et al., 1991;Rubin et al., 2003). Remarkably, this characteristic length scale has also been found in a wide range of species and organs such as tooth dentin or rein-deer antler [START_REF] Currey | Bones: Structure and Mechanics[END_REF]. Furthermore, the size, structure and distribution of mineral particles has been shown to play a major role in the mechanical behaviour of the bone material (Fratzl et al. 2004;Roschger et al. 2008). As a result, the thickness of the mineral platelets, their relative orientation and organization across the tissue are important parameters of interest for fundamental and medical studies of bone.

Small-angle X-ray scattering (SAXS) has been used for a long time to characterize size and habit of bone mineral crystals [START_REF] Fratzl | [END_REF](Fratzl et al., , 1992(Fratzl et al., , 1994(Fratzl et al., , 1996(Fratzl et al., , 1997(Fratzl et al., , 2004)). In most of this work, parameters were determined based on invariants of the scattering which do not depend on specific assumptions about the mineral particles, such as their general shape or distribution. The advantage of this approach is that it is fairly generic, while the main drawback is that it only provides a very rough description of the mineral crystals. In particular, the size of the mineral particles has been characterised by the T-parameter which describes the smallest dimension of a particle regardless of its shape. Recently, more refined models were proposed where one includes additional information to describe the SAXS data (Fratzl et al., 2005;Burger et al., 2008;Bünger et al., 2010). Based on electron microscopic evidence (Weiner & Traub, 1992;Landis et al., 1993), these models assume plate-shaped mineral crystals with a regular local parallel arrangement within the collagen microfibrils while a lesser degree of correlation can be expected between the fibrils. The models differ by their assumptions about the size distribution of mineral plates and their spacing but they have the advantage of giving some indication on the distance distribution between particles, in addition to their thickness T.

Recent development of synchrotron x-ray sources with X-ray beams of micrometer size now allows performing scanning-SAXS (sSAXS) experiments to characterize the nano-structure with a scanning resolution comparable to the size of typical heterogeneities in the tissue (Fratzl et al., 1997;Riekel, 2000;[START_REF] Gourrier | [END_REF][START_REF] Paris | Biointerphases 3[END_REF]. This is a clear advantage over TEM or AFM for medically oriented studies, for which the measure of an average particle size, spacing and relative orientation over large sample areas is usually preferable to a precise characterization of the particle organization of a very localized part of the tissue.

The application of such techniques is proving extremely useful, for example in studying the effect of osteoporosis treatments on bone material quality (Li et al., 2010), or analysing the structure forming around metallic implants [START_REF] Bünger | [END_REF].

This study reports the use of one of the stack of cards models (Fratzl et al., 2005) for the semi-automatic analysis of sSAXS from pathological bone treated with NaF. Previous reports (Fratzl et al., 1996(Fratzl et al., , 1997) ) on such samples qualitatively indicate that important structural modifications occur in the bone deposited during the treatment in the remodeling process. Here we present results that allow to discriminate between healthy and abnormal tissue and quantitatively characterize the changes in the packing of the mineral particles.

Materials and methods

The sample used in this study was prepared from a transiliac bone biopsy of an 82 yrs old woman with postmenopausal osteoporosis treated for 4 years with NaF (25mg/d). The specimen was fixed and dehydrated in a graded series of ethanol and embedded in methyl methacryalate. The surface of the sectioned bone area was polished using a grinding / polishing device (PM5 Logitech, Glasgow, Scotland) for investigations by quantitative backscattered electron imaging (qBEI), from which the bone mineralization density distribution (BMDD) was determined (Roschger et al. 1995;Roschger et al., 1998). Subsequently, a thin section of 20 µm in thickness, containing the same surface, was prepared for sSAXS measurements using a slow-speed diamond saw (Buehler Isomet, Lake Bluff, USA) and sand paper for grinding.

The scattering experiments were conducted at the microfocus beamline (μSpot) at BESSY II, Berlin, Germany (Paris et al., 2007). The X-ray beam was monochromatized to a wavelength λ = 0.984 Å (E = 12.6 keV) with a W/Si multilayer monochromator. Focussing was achieved by means of a toroidal mirror combined with a beamdefining pinhole to a diameter of 15 µm in diameter at sample position. The parasitic scattering was reduced with a guard pinhole of 20 µm in diameter placed close to the sample. The specimens were mounted on a high precision translation stage and visualized through an on-axis microscope. Optical calibration of the beam position was achieved with an accuracy of 2 μm by placing on the sample stage a high-resolution CdWO4 crystal fluorescing in the visible range. The microscope could be remotely interchanged with a 16-bit charge couple device (CCD) area detector (MarMosaic 225, Rayonix, Evanston, USA) with a 225 x 225 mm 2 square converter screen. The images were recorded with a readout time of 5 s at full resolution, i.e. 3048 x 3048 pixels with a pixel size of 73.24 x 73.24 µm 2 . The sample-to-detector distance, the detector tilt and the beam centre were calibrated using a silver behenate standard (Blanton et al., 1995).

The 2D data were azimuthally and radially integrated using the FIT2D software package (Hammersly, 1997). The 1D profiles were corrected for the dark current of the CCD and for the parasitic scattering. This procedure and the subsequent analysis were carried out using a dedicated SAXS analysis library written in Python by one of the authors (A. Gourrier). This code essentially relies on the SciPy [START_REF] Jones | SciPy: Open Source Scientific Tools for Python[END_REF] and Matplotlib [START_REF] Barret | Matplotlib -A Portable Python Plotting Package[END_REF] libraries for the analysis and imaging parts. The relative degree of orientation (so-called ρ parameter) as well as the average direction of the mineral particles were calculated using the azimuthal profiles as described in previous studies (Fratzl et al., 1996). This procedure involves fitting the azimuthal profile using Gaussian functions and a linear base line to account respectively for the fractions of oriented and unoriented particles. The relative degree of orientation is then expressed as the ratio of the areas of the peak to the intensity below the base line. Finally, the average orientation is calculated as the average value of the Gaussian peak positions taking into the account the symmetry of the SAXS signal. The analysis of the radial profiles is more complex and will therefore be described in detail in the following section. The programme MPFIT was used to fit our data. This procedure involves a least-square minimization of the statistically weighted total  2 value based on a Levenberg-Marquardt algorithm [START_REF] Markwardt | Non-Linear Least Squares Fitting in IDL with MPFIT[END_REF]. The goodness of the fit was assessed by calculating the coefficient of determination (R 2 ).

Theoretical description of the organization of the nanoparticles in bone

In dealing with the SAXS analysis of bone, we recall the basic assumption usually made on the tissue consisting of a two-phased media with nanometer size mineral inclusions embedded in a collagen matrix. This organization is described schematically in Fig. 1. The particles are oriented in preferred directions, which lead to an elliptical SAXS pattern on a 2D detector [START_REF] Fratzl | [END_REF](Fratzl et al., , 1992)). Experimentally, the signal acquisition is limited at small angles by the size and position of the beamstop and at large angles by the outer diameter and location of the CCD camera. Typically, the resolution in reciprocal space is in the order of q = qmin -qmax = 0.2 -3 nm -1 corresponding roughly to d = 1 -30 nm in direct space. The length of the scattering vector q is defined as q=4sin where 2θ is the scattering angle.

Estimation of the particle width using integral parameters

The average thickness of the particles in the volume probed by the beam can be expressed as the ratio of the volume fraction of mineral  to the total surface per unit volume of tissue σ:

(1)

Assuming a two-phased system with sharp interfaces and a mineral volume fraction  of 50 %, this quantity can be calculated from the radial SAXS profile using the following expression [START_REF] Fratzl | [END_REF](Fratzl et al., , 1992(Fratzl et al., , 1996):

(2), where (3) and (4)

This so-called T-parameter was first introduced by Porod (Porod, 1951(Porod, , 1952) ) and has the advantage of being independent of the shape and arrangement of the particles.

The first limitation in this procedure comes from the determination of the Porod constant P which can only be calculated using a restricted part of the scattering curve at high q values known as the Porod region. Futhermore, the scattered intensity I(q) also contains a constant contribution from the incoherent and inelastic scattering.

Hence, Porod's law can be expressed as I(q) = P/q 4 + L, where L is a constant. Therefore, P and L are respectively determined by linear regression of the curve I(q)q 4 vs q 4 at high q values. Although the choice of the lower limit of the Porod region is known to be critical for the determination of P, this limitating step can be partly overcome by fitting the log(I(q)) vs log(q) curve in the corresponding region where a slope of -4 is expected.

An additional limitation lies in the experimental restrictions of the accessible q-range. This implies that the only part of the integral parameter J, defined in equation ( 4), which can be measured in this range is:

. At high q, the experimental curve can be extrapolated using Porod's law (equation ( 3)). In this case, . However, there is no rigorous way to account for the missing data at low q. Two approximations are generally used assuming either a linear decrease or a constant extrapolation of the intensity at q = 0, where gives respectively , or

. From this, it follows that: J = Jmin + Jexp + Jporod.

In order to overcome these limitations, an alternative aproach has been developed, making use of models which will be described in the following sections.

Analysis of the mineral arrangement using a stack of cards model

Assuming that the T-parameter has been calculated, all the information about the shape and arrangement of the nanoparticles is contained in a rescaled function of the dimensionless parameter x = qT (Fratzl et al., 1996), defined as:

(5), where (6).

Systematic changes in G(x) curves have been found with differing species, bone age or volume fraction of mineral (Fratzl et al., 1996). Such differences seem to be very consistent throughout the existing literature and indicate some systematic differences between tissue structures. These observations entailed the need for the development of structural models to describe the G(x) curves and thus the organization of the mineral phase.

Since the particles are mostly found in the form of parallel platelets, a model describing the arrangement of particles as stacks of cards has been developed to fit the G(x) curve and thus provide sets of parameters describing the organization at this structural level (Fratzl et al., 2005) and later modified with different assumptions by Burger et al. (2008). Here we use the first of the two models which has fewer parameters and simpler analytic expressions, but, yet, still provides good fits of the data and provides quantitative indicators for the nanoparticle ordering. This model assumes that the mineral particles are positioned along the direction normal to the platelets (indicated by an arrow in Fig. 1). The function G(x) is related to the normalized (onedimensional) correlation function g(ζ) by a Fourier transform (Fratzl et al., 2005) (7)

The simplest assumption, which we are considering here, is to describe the correlation function as a damped oscillation:

(8)

The period of the oscillation describes the typical spacing between particles (in units of T) and the exponential decay (parameter 1/α) describes the distance over which the periodicity in particle spacing is damped. The normalization of G(x) given by eqs.( 5) and ( 6), which can equivalently be written as (Fratzl et al., 2005) and , immediately defines the other two constant as and

. This leads to simple model functions for and for in the form:

(9) Previously, it had been shown that the absence of positional correlation between the plates would result in (Fratzl, 1994): ( 11) and ( 12), which corresponds to the special case of equations ( 9) and ( 10) where α = 2 and β = 0. In cases where the particles are strongly disorganized, the G(x) curve can,therefore, be expected to have a Lorentzian shape.

Model based refinement of the particle size and order parameters

From a practical point of view, fitting the stack of cards model to the G(x) curves obtained from the experimental data I(q) using equation ( 10), implies that the parameters α and β will be biased by the approximation at low q in the calculation of Jmin. In order to avoid this error propagation, it is therefore preferable to fit the unbiased experimental data using equation ( 5) rewritten as: (13) using the expression of G(x) given in equation ( 10) with x = qT. P is then a fixed parameter, since it is calculated independently from the model, while T is allowed to vary to account for the unknown value of Jmin. It can thus be expected that this procedure will yield a better estimate of T than under the initial assumptions.

From the previous considerations, it also follows that the parameter β/α defined previously is insufficient to characterize the degree of ordering since it does not take into account the uncertainty in Jmin. We propose to use a new parameter, κ, defined as: (14) where Tref is a reference value of T that can be obtained either from a reference sample (e.g. control), from a region of interest of the scan considered to be representative of the healthy tissue or even from the mean T value of the scan. A previous study using this model to fit SAXS data from healthy bone lead to values of β/α ~ 1 (Fratzl et al., 2005), thus implying that κ ~ 1, since T ~ Tref in this case. Inversely, a strong disorganization should yield κ = 0 since, as seen in section 3.2, β = 0 in such a case.

Results

A quantitative back-scattered electron microscopy image (qBEI) of the sample used for this study is shown in Fig. 2a. The grey levels are proportional to the local mineral content of the bone matrix, brighter pixels indicating a higher mineral density. The trabecular feature displays a central part (roughly running vertically through the sample) which strongly differs from the adjacent bone packets. The compact and homogeneously mineralized bone matrix in the central part has the aspect of normal, healthy bone. The peripheral tissue, on the other hand, exhibits a higher mineral density albeit with a greater amount of mineralization defects. Altogether, the presence, throughout the whole biopsy, of an important fraction of more highly mineralized bone packets results in a significant shift in the bone mineralization density distribution (BMDD) curve towards higher mineral content as shown in Fig. 2d (Roschger et al., 2008). (Roschger et al., 2008) is also indicated as a dashed line and the grey contour represents a 95% prediction interval (P.I.) for healthy bone.

To analyse the structural differences across the sample, a region of interest of 600 (h) x 105 (v) µm 2 indicated by the rectangle in Fig. 2a was selected using an online microscope and scanned in steps of 15 (h) x 15 (v) µm 2 . Since the size of the scan steps corresponds to the beam diameter at sample position, this indicates a full resolution map of the area. The SAXS patterns obtained in the central and surface parts, shown respectively in Fig. 2b and Fig. 2c are quite different.

The first is typical for strongly oriented bone showing an overall elliptical aspect, while the circular shape of the second tends to indicate a lesser degree of orientation in the volume probed by the beam. This is a good indication for a stronger texture in the centre and a much less organized structure at the edges of the sample. This observation is confirmed by the examination of the corresponding Kratky plots, q 2 I(q) vs. q, shown in Fig. 3. The data Fig. 3a (corresponding to position 1 in Fig. 2a) exhibits a clear maximum and seems to decrease towards low q-values while the curve in Fig. 3c (position 3 in Fig. 2a) has a pronounced Lorentzian shape with a maximum at q = 0. Similar observations were made in previous studies, pointing towards a disorganization of the ultrastructure of the tissue in the second case (Fratzl et al., 1994(Fratzl et al., , 1996;;Roschger et al., 1997). The data in Fig. 3b (position 2 in Fig. 2a.) are somewhat intermediate since no clear maximum is found but there is an inflection point indicating a possible decrease below the resolution limit. Representative set of experimental Kratky data, I(q)q2 vs. q (crosses), and result of the fitting procedure (solid curves) using the solution of the 'stack of cards' model described by equations ( 13) and (10).

The difference profile calculated in the data range is shown below each plot. Additionally, the usual approximations made at low q in the calculation of Jmin are are indicated by dotted lines. The data in (a), (b) and (c) were respectively acquired at positions 1, 2 and 3 in Fig. 2a.

In order to derive quantitative information on the particle arrangement, the Kratky curves q 2 I(q) vs. q were fitted using equation ( 13) within a range of q = 0.2-3 nm -1 . In order to test the robustness of the fits, this procedure was repeated with a large range of initial parameters αi and βi. This showed that equation ( 13) was sufficient to provide a good fit for all the data, despite substantial differences in the shapes of the Kratky curves. This is illustrated in Fig. 3a, 3b, 3c where the fitted parameters were respectively α = 1.44, β = 1.63 (R 2 = 0.962), α = 1.37, β = 0.77 (R 2 = 0.985) and α = 2.00, β = 0 (R 2 = 0.940). It is worth noting that the later case corresponds to the conditions leading to equation ( 12) indicating a total loss of positional correlation of the platelets.

A summary of the results is shown in Fig. 4 in the form of a plot of 2π/ α vs. d = T.2π/β. This implies that the fits leading to a solution where β = 0 were discarded in these graphs. The first striking observation is that the data obtained in the central part of the sample (open circles) can clearly be distinguished from those obtained in the rest of the scan (closed circles). Interestingly, the first group of data are centred within a narrow range about a value of d ~ 17 ± 2 nm while the second are much more scattered towards higher values of d = 25-50 nm. This tends to show that there is a very sharp distribution of inter-particle distances in the central part while a progressive loss of correlation is observed elsewhere in the sample. Note that, as previously stated, these two groups yield different mineral volume fractions. Secondly, most of the data fall in a range of 2π/ α = 3.5 ± 1.5 which implies that, even for the central part of the sample where the platelets have a high degree of spatial coherence, there are some fluctuations in the extent of the short-range order. To put these results in perspective, a map of the T-parameter calculated from the model is shown in Fig. 5a. The values in the central part of the sample are close to T ~ 4 nm, a typical value for healthy bone (Fratzl et al., 2004;Roschger et al., 2008). In the rest of the sample, 'patches' can be observed where T ~ 4.5-5.5 nm indicating larger crystal sizes. These patches are separated by regions where T >> 10 nm appearing in bright on the map due to scaling. Such values are quite unusual and most likely result from an abnormal mineralization process.

Based on these observations, the degree of ordering, κ, was calculated from equation ( 14) with Tref = 4 nm (Fig. 5b). A clear distinction can be observed between the central region where κ ~ 0.9-1.0 and the rest of the scanned area where κ ~ 0 for > 90 % of the data, indicating a total disorganization of the mineral platelets. The remaining values fall in a range 0 < κ < 1. A complementary analysis of the orientation of the particles was undertaken as described in section 2. This is illustrated by the vector map in Fig. 5c. The particles are arranged in a very parallel fashion with a high degree of alignment in the central part of the tissue while a close to random orientation is found near to the sample surfaces, on the extreme left and right hand side of the image. In volumes closer to the centre, there are still regular patterns in the orientation of the particles, although gradually loosing the original vertical direction towards the edges of the sample. It is worth noting that the spatial resolution of the scan, given by the size of the X-ray beam at sample position, is sufficiently high in this study (15 µm) to obtain fine maps of the bone packets. The results shown in Fig. 5 can therefore be interpreted at the level of the individual bone packets which are the areas replaced during a single remodelling event.

Discussion

We have investigated the biopsy of an osteoporotic patient treated for 4 years with NaF using microbeam smallangle X-ray scattering. The results presented in the previous sections reveal that the ultrastructure of the vertical central part of the sample is close to that of healthy bone, while the rest of the tissue is much more disorganized.

It is therefore most likely that the former part was produced prior to the medical treatment, and the latter was formed during the ingestion of NaF by the patient. Interestingly, this implies that the central part was essentially unaffected by the treatment, as could have been postulated, e.g. by diffusion of important quantities of Fluoride ions into the hydroxyapatite particles. This is in agreement with previous work conducted in the field (Grynpas, 1990;Bang & Baud, 1990). Only those parts of the bone which were remodelled during the treatment are affected and show a strongly modified nanocomposite structure with larger mineral crystals (Fratzl et al., 1994).

The model proposed in this study turns out very successful in discriminating between healthy and pathological bone. Fig. 3 shows that the Kratky plot of the data (and hence G(x)) is significantly different in fluorotic bone. Nevertheless, the fitting procedure was found to be very stable, even in the case of strongly disorganized bone, leading to a rapidly converging result. Furthermore, the fact that this model provides a relatively simple analytical solution made it possible to process a large quantity of data in batch mode, thus allowing us to produce maps of the degree of ordering, κ, on a scale comparable to histological images.

The use of this model has yet another important advantage in providing a more realistic estimate of the Tparameter, and thus of the integral value J below the resolution limit of the instrument, than would be obtained under the usual assumptions described in section 3.1. This is schematically illustrated in Fig. 3 where the dotted lines represent the two approximations generally used in the calculation of Jmin. : the horizontal line is a constant extrapolation of qmin 2 I(qmin) at q = 0 (Jmin2) and the diagonal line assumes a linear intensity decrease to q = 0 (Jmin1). In the case of healthy bone, Jmin1 and Jmin2 will provide lower and upper bounds of the model value: Jmin1 < Jmodel < Jmin2 (Fig. 3a). However, if severe disorganization of the tissue occurs, then the two extrapolations will clearly underestimate the model value: Jmin1 < Jmin2 << Jmodel (Fig. 3c). More recently, Bünger et al. (2010) proposed a new model to calculate the T-parameter in which neither J nor the Porod constant P need to be derived. This method could provide an advantage in dealing with data with poor statistics at high q values as, for example, measurements performed with laboratory sources. However, with the generalization of the use of synchrotron radiation, we believe that the estimation of the Porod constant can be performed with sufficient statistical accuracy. A difference map between the T-parameter calculated with the stack of cards model and with the classical approach is shown in Fig. 6. Overall, the values obtained using the model are always higher. However, the difference in the central part of the sample is very small, with ΔT < 0.5 nm, i.e. < 10 %, which means that the classical approach provides reasonable estimates for healthy bone. On the other hand, the values obtained in the rest of the sample are generally higher by ΔT = 0.7-1. nm, i.e. 18-25 % and up to ΔT > 1.5 nm, i.e. > 38 % in the brightest regions. The use of our model therefore seems clearly beneficial when analysing bone with severe structural pathological alteration.

Conclusion

The procedure described in this paper represents an improvement in the quantitative description of the organization of bone tissues at the nanometer scale. A more reliable measure of the average mineral particle size was obtained using a model describing the arrangement of the platelets as a regular stack of cards. This model is very close to the description of the tissue's ultrastructure found in the literature and is relatively simple, yielding an analytical solution for the scattering curves. In particular, new parameters were derived to characterize the fluctuations in positional correlation of the particles. Using X-ray microbeams from a second generation synchrotron source, we were thus able to discriminate between healthy and pathological bone with a resolution comparable to histological observations. This allowed to evidence severe structural changes in the bone packets deposited during a anti-osteoporotic treatment based on the administration of NaF.

Figure 1 .

 1 Figure 1. Schematic illustration of the fibril structure showing the internal arrangement of mineral platelets separated by layers of collagen molecules in a 'stack of cards' fashion. The short-range order results in a strong positional correlation of the mineral platelets along the direction normal to the faces indicated by the arrow. The fluctuations in position and thickness of the platelets and organic layers along this axis can adequately be described by a one-dimensional Poisson distribution. In the other two directions, there is no experimental evidence for such short-range order.

  construction of the model, the parameter β describes the degree of spatial correlation between successive plates, separated by an average distance d = T.2π/β. The other parameter, α, controls the damping of the oscillations and, therefore, provides an indication on the relative extent of the ordering. Smaller values of 2π/α point to a reduction in the short-range ordering of the platelets (greater disorder). As a result, at constant T values, β/α can be seen as a reasonable (scalar) indicator of the degree of ordering since the degree of spatial correlation is then modulated by the relative extent of the short-range order.

Figure 2 .

 2 Figure 2. (a) qBEI image of a region of interest of a sample section prepared from a transiliac bone biopsy of an 82 years old osteoporotic woman treated with NaF. The rectangle indicates the area scanned with an X-ray beam of 15µm in diameter. The size of the region probed by the microbeam is indicated by the square dots enlarged by a factor of 2 for clarity. The scale bar indicated at the bottom is 200µm. Examples of the 2D SAXS pattern obtained in the centre of the sample and on the left edge are shown respectively in (b) and (c). The BMDD deduced from the grey level histogram of the qBEI image of the whole sectioned compartment of the biopsy is shown as a solid line in (d). A reference curve(Roschger et al., 2008) is also indicated as a dashed line and the grey contour represents a 95% prediction interval (P.I.) for healthy bone.

Figure 3 .

 3 Figure 3.Representative set of experimental Kratky data, I(q)q2 vs. q (crosses), and result of the fitting procedure (solid curves) using the solution of the 'stack of cards' model described by equations (13) and (10).The difference profile calculated in the data range is shown below each plot. Additionally, the usual approximations made at low q in the calculation of Jmin are are indicated by dotted lines. The data in (a), (b) and (c) were respectively acquired at positions 1, 2 and 3 in Fig.2a.

Figure 4 .

 4 Figure 4. Plot of 2π/α vs. d = T 2π/β. The open and closed circles represent respectively the results of the fit of the data collected in the central part of the sample and in the rest of the scanned area.

Figure 5 .

 5 Figure 5. A map of the T-parameter in the scanned area is shown in (a) and corresponding parameter κ in (b) with appropriate scaling. Similarly, the average orientation of the mineral particles at each position of the scan is shown in (c) in the form of vectors. The direction of the vectors indicate the average orientation of the mineral particles in the volume probed by the beam while their length and colour correspond to the relative degree of orientation ρ. The scale bar indicated at the bottom is 200 µm.

Figure 6 .

 6 Figure 6. Difference map between the T-parameter calculated with the stack-of cards model and the values obtained with the integral parameters using equation (2). The scale bar at the bottom is 200 µm.


Florian Wolf
email: florian.wolf@cea.fr

Franc ¸ois Dehmas Vincent Berg
email: vincent.berg@cea.fr

Valérian Mannoni
email: valerian.mannoni@cea.fr

Sébastien De Rivaz
email: sebastien.derivaz@cea.fr

Multi-Frequency Phase Difference of Arrival for Precise Localization in Narrowband LPWA Networks

Keywords: Localization, LPWA network, MF-PDoA, frequency hopping, narrowband localization, phase synchronization

In Low Power Wide Area (LPWA) networks, radio localization based on Time Difference of Arrival measurements has major benefits: energy consumption of the node and spectrum usage are contained. Temporal resolution and hence positioning accuracy is however limited by the bandwidth of the LPWA narrowband signals. Multi-channel ranging, that relies on multiple narrowband signals that are sequentially transmitted on different and discontinuous channels, has recently been proposed to improve precision. However, the concept has not been adapted yet to differential uplink localization techniques. This paper studies and proposes adaptations to multi-channel ranging called Multi-Frequency Phase Difference of Arrival (MF-PDoA) for this purpose. It analyzes the impact of imperfect synchronization on its performance and concludes that MF-PDoA is robust and well adapted to the LPWA scenario.

I. INTRODUCTION

Low Power Wide Area (LPWA) networks are emerging communication technologies that enable wireless connectivity on a large variety of objects in the context of the Internet of Things (IoT). Long-range communication of both, proprietary (e.g. Long Range from Semtech (LoRA)) and standardized solutions (Narrow-Band IoT (NB-IoT)), is achieved thanks to low levels of receiver sensitivity obtained by low data rates and narrowband modulation schemes.

Wearable health monitoring is one example of a promising application of LPWA networks [START_REF] Lehne | D3.1: Healthcare Vertical Trial Requirements Definition and Execution[END_REF]. It allows to remotely monitor a patient and in case of an emergency (e.g. heart stroke) alert services to quickly react. LPWA radios are suited for these applications thanks to their small size, low-power and low-cost, while providing indoor and outdoor coverage. Localization of the patient is essential to facilitate the work of emergency services. Global Navigation Satellite System (GNSS) solutions are unsuitable because of extra energy consumption and unavailability for indoor operation. An energy efficient approach should derive position information from the same radio signal as used to transmit health-monitoring data. However, current approaches are not accurate enough to meet the requirements of the application [START_REF] Lehne | D3.1: Healthcare Vertical Trial Requirements Definition and Execution[END_REF].

Time of Arrival (ToA) based localization is bound by the signal bandwidth used to estimate ranging information [START_REF] Skolnik | Radar Handbook[END_REF]. This makes precise LPWA localization challenging and difficult to resolve in a multipath environment. Practical implementations and field trial measurements in LoRA networks typically achieve mean positioning errors of 200 m in an urban environment [START_REF] Podevijn | Experimental TDoA Localisation in Real Public LoRa Networks[END_REF]. A new technique, called coherent multi-channel ranging, has therefore been proposed to significantly improve precision. It relies on multiple narrowband signals that are sequentially transmitted on different and discontinuous channels to virtually increase the bandwidth [START_REF] Wolf | Coherent Multi-Channel Ranging for Narrowband LPWAN: Simulation and Experimentation Results[END_REF] [START_REF] Wolf | Accurate Narrowband LPWA Ranging: Principles, Performance in AWGN and Multipath Channels[END_REF]. This technique, based on multiple phase measurements or Phase of Flight (PoF), significantly improves temporal resolution and ranging precision, while preserving the narrowband modulation necessary for long-range communication. The technique has been so far evaluated in two-way ranging scenarios, where elements (nodes or base-stations) sequentially perform a two-way packet exchange in order to remove absolute time synchronization constraints [START_REF] Wolf | Coherent Multi-Channel Ranging for Narrowband LPWAN: Simulation and Experimentation Results[END_REF] [START_REF] Wolf | Accurate Narrowband LPWA Ranging: Principles, Performance in AWGN and Multipath Channels[END_REF]. When LPWA base stations are strictly synchronized using e.g. GNSS time reference, Time Difference of Arrival (TDoA) is preferred as it does not require a bi-directional exchange and therefore preserves the energy of the IoT node and limits spectrum occupancy.

Application of TDoA positioning techniques to multi-frequency PoF ranging or Multi-Frequency Phase Difference of Arrival (MF-PDoA) is however not straightforward. MF-PDoA in the context of Radio-frequency identification (RFID) passive tags has already been considered and proved efficient [START_REF] Qiu | Multifrequency Phase Difference of Arrival Range Measurement: Principle, Implementation, and Evaluation[END_REF]. However, Phase Difference of Arrival (PDoA) in this context considers perfect phase, frequency and time synchronization between distributed antennas leading to major simplifications because frequency and phase is generated from a single Radio Frequency (RF) source (i.e.: the RFID reader). Furthermore, only very short ranges are possible. Promising implementation of MF-PDoA in a short-range indoor context is considered in [START_REF] Pichler | Multi-Channel Distance Measurement With IEEE 802.15.4 (ZigBee) Devices[END_REF] using 2.4 GHz IEEE 802.15.4 PHY-layers, however, synchronization between access points is performed using a reference device transmitting from a known location. This transmission signal is time interleaved with the transmission of the devices to be located to avoid large synchronization drifts. This approach is not recommended for LPWA networks where frequency resources are scarce. Besides, 80 MHz bandwidth and ranges up to a few meters are considered in [START_REF] Pichler | Multi-Channel Distance Measurement With IEEE 802.15.4 (ZigBee) Devices[END_REF].

This paper evaluates and derives performance of coherent

Analog Numerical s 0 f w [T] c , φ R [T] c PA, g R [T] c s R [T] c h [T,R i ] (τ) r R [T,R i ] c LNA, g R [R i ] c f [R i ] w , φ [R i ] R r ADC [T,R i ] c A D r IF [T,R i ] c f c r [T,R i ] c
LPWA transmi er Base station receiver multichannel ranging using a new MF-PDoA approach adapted to the LPWA scenario. A new phase synchronization mechanism imposes some architecture constraints on the base stations while significantly relaxing requirements on the nodes. This paper is organized as follows. The system and signal model is introduced in Section II. The MF-PDoA estimator is presented in Section III. Section IV evaluates and compares performance of MF-PDoA with TDoA. Section V concludes the paper and outlines perspectives.

II. SYSTEM AND SIGNAL MODEL

Range difference estimates between a transmitter T and multiple base stations R i , ∆d [Ri,R l ] , are used by hyperbolic multilateration to derive the position of the transmitter. Position accuracy can be approximated as a product of range difference estimate standard deviation and geometric dilution of precision [START_REF] Spirito | On the Accuracy of Cellular Mobile Station Location Estimation[END_REF]. Hence, precise and accurate range difference estimates are key to achieve precise localization.

Multi-channel ranging relies on sequentially transmitting a narrowband signal, s 0 (t), over different frequencies

f wc = f w + f c , c ∈ [0, C -1].
The frequency, f w , is called the central frequency and f c , the channel frequency and is usually equal to c∆f . C is an integer equal to the number of channel frequencies considered by the algorithm.

In order to evaluate the quality of the localization metrics, a signal model with impairments is considered. The transmitter, channel and receiver architecture used to derive the signal model is given in Fig. 1. An arbitrary LPWA RF transmitter architecture is allowed on the node side, while on the base station, the reception stage is done through an Intermediate Frequency (IF) followed by a digital conversion and a second conversion stage performed digitally to keep coherence of phase between frequencies. It should be noted that this architecture model is compatible with cellular Internet of Things (IoT) base stations where receiver bandwidth is often much larger than a single channel.

Since transmitter and base stations are generally not synchronized with respect to the absolute reference time t, their local time may be expressed as:

t [X] = 1 + δ [X] f t + t [X] 0 , (1)
where δ

[X] f is the normalized relative frequency offset of element X ∈ {T, R i }. In the following and without loss of generalization, δ 0 are set to zero. As the localization scenario only considers one transmitter and multiple receivers, performance estimation of the ranging metrics may be done considering that the transmitter time is the reference time and imperfect synchronization of the receivers is relative to the transmitter time.

With these assumptions, when transmitted over carrier frequency f wc , the signal at the antenna of the IoT node becomes:

s R [T ] c (t) = g R [T ] c s 0 t -t [T ] d e j(2πfw c t+φR [T ] c ) , (2)
where g R [T ]
c , is a complex number value that models RF hardware gains, t c , an arbitrary phase term function of carrier f wc . The RF signal is then transmitted over a wireless RF channel so that at base station

R i , r R [T,Ri] c (t) = s R [T ] c (t) * h [T,Ri] (τ ) + ν(t), (3)
where h [T,Ri] (τ ) is the Channel Impulse Response (CIR) between transmitter T and base station R i and is assumed stationary during the communication, * , the convolution product and ν(t), Additive White Gaussian Noise (AWGN). Lets now assume that the coherence bandwidth of the channel, h [T,Ri] , is significantly larger than the bandwidth of the baseband signal s 0 (t). This assumption allows to consider flat fading conditions for the transmitted signal s R

[T ] c (t), c ∈ [0, C -1]
and the convolution product in (3) can be replaced by a simple complex product.

r R [T,Ri] c (t) = α [T,Ri] c e jϕc [T ,R i ] s R [T ] c t -τ [T,Ri] 0 + ν(t), (4)
where α

[T,Ri] c , ϕ c [T,Ri] and τ [T,Ri] 0
are real numbers respectively equal to the amplitude, phase and delay of the channel at frequency f wc .

At the base station receiver, R i , the received signal is amplified by a low noise amplifier, then mixed by f

[Ri] w = (1 + δ [Ri] f )f w , to be downconverted to IF. r ADC [T,Ri] c (t) = g R [Ri] c r R [T,Ri] c (t)e -j 2πf [R i ] w t+φ [R i ] R + ν Rc [T,Ri] (t), (5)
r [T,Ri] c [k] = g R [T ] c g R [Ri] c α [T,Ri] c • s 0 δ [T,Ri] T kT S -t A [T,Ri] c • e j 2πδ [T ,R i ] f (fw+fc)kTS+φA [T ,R i ] c + n [T,Ri] c [k] (8)
where g R

[Ri] c

is the complex gain of the receiver, φ

[Ri] R , the unknown initial phase value of the receiver, and ν

[T,Ri] Rc (t), noise ν(t) at the receiver IF. Signal r ADC [T,Ri] c (t) is then sampled every kT [Ri] S
by the local clock of receiver R i , where k is an integer:

r IF [T,Ri] c [k] = r ADC [T,Ri] c kT S -t [Ri] 0 1 + δ [Ri] f . ( 6
)
The samples are finally numerically down-converted to baseband with f c such that

f c + f w = f wc . r [T,Ri] c [k] = r IF [T,Ri] c [k] e -j2πfckTS + n [T,Ri] c [k], (7)
where n

[T,Ri] c [k] is the equivalent AWGN noise at the baseband of receiver R i .
The baseband received signal in [START_REF] Pichler | Multi-Channel Distance Measurement With IEEE 802.15.4 (ZigBee) Devices[END_REF] can finally be expressed as function of the baseband waveform s 0 (t) in ( 8) by combining (1) to [START_REF] Qiu | Multifrequency Phase Difference of Arrival Range Measurement: Principle, Implementation, and Evaluation[END_REF]. The following parameters are introduced:

1) Time dilation: Time dilation is a number close to 1 and accounts for relative sampling time errors between transmitter and receiver.

δ [T,Ri] T = 1 1 + δ [Ri] f (9)
2) Time of Arrival:

the transmitted signal, s R [T ] c , is received at time, t A [T,Ri] c , at base station R i : t A [T,Ri] c = 1 + δ [Ri] f τ [T,Ri] 0 + t [Ri] 0 + 1 + δ [Ri] f t d [T ] . (10)
3) Relative frequency offset: The received baseband signal at R i is modulated by the remaining relative frequency offset, δ

[T,Ri] f (f w + f c ), where: δ [T,Ri] f = -δ [Ri] f 1 + δ [Ri] f (11)
4) Phase of Arrival (PoA): The phase of arrival is a term dominated by the rotation of the carrier frequency over the duration of the time of arrival. Impairments on this term include relative frequency offset of the transmitter, the channel phase, ϕ

[T,Ri] c
, phase offsets due to time offsets between transmitter and receiver, as well as initial phase.

φ A [T,Ri] c = -2π (f w + f c ) τ [T,Ri] 0 + ϕ [T,Ri] c + 2π δ [Ri] f f w -f c 1 + δ [Ri] f t [Ri] 0 + φ R [T ] c -φ [Ri] R . (12)
The latter term is particularly critical as it is the main limiting factor of quality of the MF-PDoA estimation process. The phase signal, φ

[T,Ri] c

[k], is then introduced as the phase of the received signal r

[T,Ri] c [k]: φ [T,Ri] c [k] = 2πδ [T,Ri] f (f w + f c ) kT S + φ A [T,Ri] c . (13)

III. RANGE DIFFERENCE ESTIMATION

A. Multi-Frequency Phase Difference of Arrival (MF-PDoA)

The principle of MF-PDoA is based on the quasi linear relationship between phase φ Ac and delay τ 0 in an homogeneous transmission medium (e.g. free space). Measuring the phase delay over multiple frequencies provides two main benefits: it improves range resolution and increases the maximum unambiguous range [START_REF] Wolf | Coherent Multi-Channel Ranging for Narrowband LPWAN: Simulation and Experimentation Results[END_REF]. Range resolution is inversely proportional to the overall signal bandwidth while range ambiguity is inversely proportional to the frequency resolution. Fig. 2 gives an example of MF-PDoA spectrum occupation for an LPWA transmission. Packets of narrowband instantaneous bandwidth B are sequentially transmitted on C channels with channel spacing ∆f . An overall virtual bandwidth B virt = ∆f (C -1) is therefore used.

MF-PDoA relies on a pair of base stations (R i , R l ) to estimate the difference of PoAs from [START_REF] Etsi | Evolved Universal Terrestrial Radio Access (E-UTRA); User Equipment (UE) radio transmission and reception[END_REF]. As for TDoA, MF-PDoA constant but unknown terms in PoA are eliminated by the difference.

The PDoA measurement between the transmitter T and the pair of base stations (R i , R l ) for channel c can be derived by subtracting (13) for both base stations:

∆φ [Ri,R l ] c [k] = φ [T,Ri] c [k] -φ [T,R l ] c [k] (14)
Eq. ( 14) can then be extended in details and analyzed.

∆φ [Ri,R l ] c [k] = -2π (f w + f c ) τ [T,Ri] 0 -τ [T,R l ] 0 (15a) + ϕ [T,Ri] c -ϕ [T,R l ] c (15b) -φ [Ri] R -φ [R l ] R (15c) -2π t [Ri] 0 1 + δ [Ri] f - t [R l ] 0 1 + δ [R l ] f f c . (15d) + 2π δ [Ri] f t [Ri] 0 1 + δ [Ri] f - δ [R l ] f t [R l ] 0 1 + δ [R l ] f f w (15e) -2π δ [Ri] f 1 + δ [Ri] f - δ [R l ] f 1 + δ [R l ] f (f w + f c ) kT s (15f)
1) Perfect base station synchronization: When the base stations are perfectly synchronized in both time and frequency, i.e δ can be estimated and (15c) has almost no impact on the range difference estimation.

[Ri] f = δ [R l ] f and t [Ri] 0 = t [R l ] 0 ,
2) Imperfect base station synchronization: In the more realistic case of imperfect synchronization between the base stations, e.g. when base stations are synchronized through GNSS signals [START_REF] Loraalliance | Geolocation Whitepaper[END_REF], the time and frequency offset terms of (15) have to be considered. The inter-base station time synchronization offset ∆t 0 ≈ t

[Ri] 0 -t [R l ] 0
in (15d) directly impacts range difference measurement in (15), as it depends on the channel frequency f c and is therefore generally not compensated for.

The term (15e), although directly impacted by the relative time and frequency errors t

[Ri] 0 , t [R l ] 0 , δ [Ri] f and δ [R l ]
f , is independent from f c and hence remains constant during the reception process. It thus has very little impact on the range difference estimation.

Finally, term (15f) is a phase rotation that is function of each channel c, the carrier frequency f w and the channel frequency f c . It is due to the relative Carrier Frequency Offset (CFO) between the base stations R i and R l . Although estimation of the CFOs, δ

[Ri] f and δ [R l ]
f , is possible, the level of precision required is such that it would lie below the theoretical bound. This approach is therefore excluded. To compensate for this last term a symmetric frequency hopping scheme is thus introduced as illustrated in Fig. 2.

A frequency hopping scheme repeated with time reversal symmetry is thus considered (see Fig. 2). The time reversal symmetry imposes that the signal transmitted at time t = kT s is the same as the one at time (Nk)T s . It is assumed that the duration of the narrowband signal s 0 is equal to N c T s . Every N c T s seconds, the carrier frequency of the signal is therefore changed. The overall C channels are covered in CN c T s seconds and the overall signal duration is equal to N T s seconds, with N = 2CN c . Besides, if at time t 1 = kT s , the channel frequency of the narrowband signal is f c , then it is also the channel frequency of the channel at time t 2 = N T st 1 = (Nk)T s since channels are assumed stationary.

∆φ [Ri,R l ] c [k] and ∆φ [Ri,R l ] c [N -k] are combined into a new variable Φ c [k]
assuming the property of time reversal in the transmitted sequence.

Φ c [k] = ∆φ [Ri,R l ] c [k] + ∆φ [Ri,R l ] c [N -k] = -4π (f w + f c ) τ [T,Ri] 0 -τ [T,R l ] 0 (16a) + 2 ϕ [T,Ri] c -ϕ [T,R l ] c (16b) -2 φ [Ri] R -φ [R l ] R (16c) -4π t [Ri] 0 1 + δ [Ri] f - t [R l ] 0 1 + δ [R l ] f f c (16d) + 4π δ [Ri] f t [Ri] 0 1 + δ [Ri] f - δ [R l ] f t [R l ] 0 1 + δ [R l ] f f w (16e) -2π δ [Ri] f 1 + δ [Ri] f - δ [R l ] f 1 + δ [R l ] f (f w + f c ) N T s . ( 16f
)
Only term (16d) of ( 16) remains dependent of time as it changes every N s symbols (i.e. every time the channel frequency is modified). This term is dependent on the synchronization between both receivers.

Φ c [k] of ( 16) is then used to construct, Hc , a compound channel transfer function of the MF-PDoA system:

Hc = A c e jΦc = H [T,Ri] c 2 H [T,R l ] c 2 * , (17)
where

A c = α [T,Ri] c α [T,Ri] 2C-1-c α [T,R l ] c α [T,R l ]
2C-1-c . The time domain equivalent of Hc is the Compound CIR (CCIR), and may be derived by Inverse Discrete Fourier Transform (IDFT):

h(τ ) = h [T,Ri] * h [T,Ri] ⊗ h [T,R l ] * h [T,R l ] (τ ), ( 18
)
were ⊗ is the cross-correlation product. The range difference ∆d [Ri,R l ] is extracted by estimating the delay of the maximum amplitude in the CCIR of (18). The time reversal symmetry introduced in the transmitter to mitigate frequency error results in the auto-convolution of the CIR h [T,Ri] (τ ) and h [T,R l ] (τ ), while the cross-correlation product is introduced by the phase difference between transmitter and both base stations. This leads to unwanted path components in the CCIR and degrades the performance of estimation of the differential range.

B. Time Difference of Arrival (TDoA)

In comparison, the TDoA method estimates time difference of arrival from ToA estimates of [START_REF] Wolf | Improved Multi-Channel Ranging Precision Bound for Narrowband LPWAN in Multipath Scenarios[END_REF]. When base stations, (R i , R l ) are synchronized, i.e. t

[Ri] 0 = t [R l ] 0 and δ [Ri] f = δ [R l ]
f , the unknown and arbitrary transmitter time offset t

[T ] 0 cancels and the resulting range difference is equal to

∆d [Ri,R l ] = d [T,Ri] -d [T,R l ] = c 0 t [T,Ri] A -t [T,R l ] A . (19)
Due to narrowband signals, ToA estimates t A will not be as precise as MF-PDoA, but will serve as performance comparison to the proposed MF-PDoA technique.

IV. PERFORMANCE COMPARISON BETWEEN TDOA AND MF-PDOA FOR LPWA IOT SCENARIOS

Range difference derivations of Section III are used to evaluate and compare performance between TDoA and MF-PDoA for LPWA IoT scenarios. Two different LPWA scenarios have been considered for the evaluation: LoRA and NB-IoT. The conducted numerical MATLAB simulations assume AWGN channels. In the following, the Signal to Noise Ratio (SNR) level of the received signal has been considered identical for both receivers R i and R l . This is practically rarely the case, but constitutes a worst case scenario.

A. Cramer Rao Lower Bound (CRLB) for TDoA and MF-PDoA

The CRLB for MF-PDoA in presence of AWGN channel can be derived from the coherent multi-channel bound given in [START_REF] Wolf | Improved Multi-Channel Ranging Precision Bound for Narrowband LPWAN in Multipath Scenarios[END_REF] where the CRLB is derived for PoA. MF-PDoA is a difference of two PoA metrics, φ

[T,Ri] c and φ [T,R l ] c
, as defined in (14). Therefore the variance of the delay estimation of the difference of both phases is equal to the sum of the variances since φ

[T,Ri] c and φ [T,R l ] c
are uncorrelated. We can thus write:

σ M F -P DOA ∆τ,CRLB = 1 2π 2 Es N0 BW 2 RMS ( 20
)
with root-mean-squared bandwidth

BW RMS = +∞ -∞ f 2 |S 0 (f )| 2 df +∞ -∞ |S 0 (f )| 2 df (21)
and Es N0 , the symbol to noise power density over the transmitted sequence. Hence, under AWGN, the derivation of the CRLB of MF-PDoA is similar to the derivation of the CRLB of TDoA but with different values of Es N0 and BW RMS . CRLB however only considers performance with perfect synchronization. Numerical simulations are thus considered to compare performance of non perfectly synchronized receivers with the derived CRLB.

B. Application of MF-PDoA to LoRA numerology

For this application, one transmitter and two base stations are considered. The transmitted signal consists on a Gold code sequence of length N c = 32 coded onto a Binary Phase Shift Keying (BPSK) modulation at the rate of 125kHz (T s = 8µs). The transmitted sequence is repeated on C = 16 different frequencies with a channel spacing ∆f equal to 200kHz and a time gap of 4 symbol periods (or 32µs). Time reversal symmetry is provided as shown in Fig. 2. The overall sequence duration is therefore less than 9.2ms. Simulations have considered arbitrary and independent carrier frequency errors of the receivers within ±10 ppm relative to the transmitter. Carrier frequency error is compensated at the receiver by the Mengali algorithm [START_REF] Mengali | Data-Aided Frequency Estimation for Burst Digital Transmission[END_REF].

Fig. 3 gives the standard deviation of the range difference error as a function of the SNR for the various simulated cases. The theoretical CRLB is also given for both TDoA and MF-PDoA. When perfect time synchronization between the base stations is considered, as expected, simulations converge towards CRLB for higher SNR levels and PDoA significantly outperforms TDoA. Although performance is limited by the carrier frequency estimator algorithm, the CRLB is reached. For a SNR of 6dB, the standard deviation error of MF-PDoA is equal to approximately 1m while the performance of TDoA is more than 30-times larger. Then time synchronization between receivers is relaxed. For each receiver, time synchronization errors are introduced using a normal distribution of increasing standard deviation, respectively 10ns, 30ns and 50ns on each receiver. For MF-PDoA, a performance floor is observed. This performance floor exactly matches the accuracy of the base station synchronization, as the performance floor converges towards √ 2σ

t [R i,l ] 0 c 0 , where σ t [R i,l ] 0
is the standard deviation of the time synchronization error between base stations. When the synchronization error is set to 10ns (respectively 50ns), the RMS range error floor is equal to 4m (respectively 21m). Since TDoA performance is not as close to these limits, it only gets affected by synchronization errors at higher SNR values and only for the 50ns. The MF-PDoA behavior remains however very robust in presence of synchronization impairments.

C. Application of MF-PDoA to NB-IoT numerology

A second simulation scenario is geared towards application of MF-PDoA to NB-IoT. For this case, we assumed that signaling is transmitted over a Narrowband Physical Uplink Shared Channel (NPUSCH) modulated using a single carrier . modulation (π/2-BPSK) at a rate equal to 15kHz. The transmitted signal consists of a Barker sequence of N c = 7 bits in order to fit the coded sequence within the duration of a NB-IoT frame slot. The transmitted sequence is repeated on C = 16 different frequencies with a channel spacing ∆f = 180kHz. The complete localization signal could therefore be transmitted over two 15kHz single carrier NB-IoT frames (or 16ms). As for the LoRA simulations, the CRLB is estimated for this scenario. The carrier frequency error is set to ±1 ppm as NB-IoT requires transmitter to synchronize within ±0.2 ppm [START_REF] Etsi | Evolved Universal Terrestrial Radio Access (E-UTRA); User Equipment (UE) radio transmission and reception[END_REF] and therefore receivers are not expected to compensate for large frequency errors. The carrier frequency error however requires compensation at the receiver and the Mengali algorithm is also used [START_REF] Mengali | Data-Aided Frequency Estimation for Burst Digital Transmission[END_REF]. The time synchronization hypotheses are the same as for the previous scenario. It should be noted that although the numerology of NB-IoT has been carefully considered, frequency hopping is currently not supported by the standard and would require some amendments.

, σ t0 [R i,l ] = 0ns TDoA, σ t0 [R i,l ] = 50ns CRLB PDoA PDoA, σ t0 [R i,l ] = 0ns PDoA, σ t0 [R i,l ] = 10ns PDoA, σ t0 [R i,l ] = 30ns PDoA, σ t0 [R i,l ] = 50ns
Fig. 4 gives the performance results for this scenario. Since the duration of the sequence is much shorter than in the previous case, performance of the TDoA and MF-PDoA algorithms converges for larger SNR values. For a SNR of 10dB, the standard deviation error of MF-PDoA is equal to approximately 2m while the standard deviation error of TDoA is equal to 150m. Floor performance behavior is the same as in the LoRA scenario when synchronization between transmitter and receivers is relaxed: when synchronization error is set to 10ns (respectively 50ns), the RMS range error floor is equal to 4m (respectively 21m). Both TDoA and MF-PDoA are robust to synchronization timing errors between receivers. However, since the signaling sequence is much shorter than in the previous application, MF-PDoA remains significantly better than TDoA. The TDoA performance is much worst (i.e. larger than 120m) than the error due to time synchronization (i.e. 21m).

V. CONCLUSION This paper proposes and investigates an important evolution to coherent multi-channel ranging designed to limit energy consumption based on MF-PDoA for LPWA. The resulting metric provides difference in distances between the IoT node and two reference receivers (base stations) and provides localization with hyperbolic trilateration. As MF-PDoA relies on phase, time and frequency synchronization between base stations, the impact of imperfect synchronization between receivers is first analyzed and then simulated and compared against CRLB performance assuming perfect synchronization. Simulations performed for both LoRA and NB-IoT scenarios demonstrate the benefits of MF-PDoA in this context. Further work should consider field trial implementations of MF-PDoA with experimental hardware to confirm the hypothesis used for this simulation work.

Fig. 1 .

 1 Fig. 1. Architecture model of the LPWA transmitter, channel model and base station receiver. Numerical stages of the base station are assumed coherent in phase.

Fig. 2 .[

 2 Fig. 2. Time versus frequency spectral representation of the transmitted MF-PDoA signal.

Fig. 3 .

 3 Fig. 3. LoRA scenario. Range difference error as a function of the SNR at the receiver and for different levels of synchronization errors σ t [R i ,l] 0 .

Fig. 4 .

 4 Fig. 4. NB-IoT scenario. Range difference error as a function of the SNR at the receiver and for different levels of synchronization errors σ t [R i ,l] 0

ACKNOWLEDGMENT

This work was supported by the European Union's Horizon 2020 research and innovation programme under 5G-HEART project (grant agreement No 857034).


domain: sde shs

text:
Forêt de Mormal

Grand Bois l'Evêque

Petit Bois l'Evêque

Bois de Ribeaucourt

  adjudicataires de vente à l'unité au Grand Bois l'Evêque (1600-1650)


text:
Petit Bois l'Evêque

Bois de Ribeaucourt

  adjudicataires de vente à l'unité au Grand Bois l'Evêque


text:
Terre d'Avesnes Forêt de Mormal

Prévôté du Quesnoy

Prévôté de Valenciennes Te rr e d e C a m b ré s is

Bois de Ribeaucourt

L ' E s c a u t


Christian Ersche

Ercüment Çelik

Wiebke Keim

Veronika Wöhrer

Editorial to SPECIAL ISSUE: Circulating Knowledge. Empirical Examples from the Social Sciences

Introduction

This special issue of Transcience contributes to the current debates about the circulation of knowledge in the social sciences. The circulation of practically applicable knowledge, as well as abstract theoretical knowledge; of worldviews, ideologies, and epistemologies, is not at all a new phenomenon. In recent decades several critical studies of such circulation processes have focused especially on the historical as well as the contemporary asymmetry of knowledge transfers. As social sciences were institutionalized during the historically short period of European world domination, it was easy for European scholars to declare non-European forms of knowledge irrelevant. Observing global trends in the social sciences of today, we can therefore find an interestingly ambivalent phenomenon: On the one hand, the European research area and its achievements still enjoy a high standing outside of Europe. On the other hand, the worldwide influence of the European theoretical tradition with its inherent claim to universality is increasingly being perceived as overbearing and presumptuous.

The papers collected in this issue add empirical studies and theoretical reflection to this topic from sociological, historical, and anthropological perspectives. They attempt to assess the acceptance potential of knowledge beyond its relative context of origin and reflect on the preconditions of acceptance or rejection. A closer look at specific empirical examples, however, shows that there are not only different forms of transfer and circulation of knowledge. More generally, this collection shows that there are different types and forms of knowledge which can and should bear relevance for the social sciences themselves. This issue brings together five articles from international scholars who conducted their research as fellows or team members of the research project "Universality and the Acceptance Potential of Social Science Knowledge: On the Circulation of Knowledge between Europe and the Global South", which is funded by the BMBF (German Federal Ministry of Education and Research) and affiliated with the Institute of Sociology at the University of Freiburg, Germany. The idea for this special issue emerged at a workshop entitled "For a Generational Shift in Analysing Circulating Knowledge", which was organized in the framework of this project in Freiburg in 2012. 1 The aim of the workshop and, in consequence, of this special issue is to add new perspectives to the study of processes of production and circulation of knowledge.

In his contribution, Cristiano Lanzano problematizes the popular concept of 'indigenous knowledge', an issue that connects with academic debates as well as practical development cooperation. In an ethnographic case study on 'traditional' environmental practices and sacred sites, which he conducted in the area of the Comoé-Léraba reserve in Western Burkina Faso, he shows the complexity of these practices and places as social and cultural phenomena. Although he acknowledges that the concept of 'indigenous knowledge' has been useful for an appreciation of the knowledge of local populations in the global south which has been devalued since colonial times, he also points to the ambiguity that lies in this concept's tendency to abstract from broader social and cultural aspects and reduce this knowledge to its effectiveness.

Ercüment Çelik explores the circulation of knowledge between academic and extra-academic fields in his paper by focusing on the intellectual engagements with the labour movement in South Africa in the 1970s and 1980s. In a case study on the Institute for Industrial Education (IEE) and the South African Labour Bulletin (SALB), Çelik shows the multiple functions and identities of intellectuals and the characteristics of their engagement with the labour movement. After reviewing several theoretical approaches to the role of intellectuals in the labour movement he argues for a complementary understanding of the engagements of different groups of intellectuals with this movement. Instead of focusing on the social and cultural distance between white intellectuals and black working class, he presents a view that looks for areas where this distance could be bridged.

In her field study with Muslim students of the Malabar Christian College in Kozhikode (Calicut) and their families, Barbara Riedel outlines a specific and locally rooted cosmopolitanism in a region that has been the stage for economic, cultural, and religious interaction for centuries. She shows the crucial role the Mappila Muslims of Kerala played in the cosmopolitan past of this region and the new kind of cosmopolitanism emerging among them. She describes "cosmopolitanism at work" as a process occurring not only in global cities but also in peripheral places and in small scale everyday situations. It is a process that extends far into history. It is neither exclusively modern, nor does it develop along a single timeline, nor is it irreversible.

In his contribution on intellectuals and knowledge production in the 18th Century Salonika, a province of Ottoman Empire, Irfan Kokdas demonstrates how knowledge production goes beyond a simple centerperiphery relation and becomes multi-centered. He underlines the importance of rising mobility of intellectuals and multiplication of political channels in the circulation of knowledge in a broad geography. This knowledge production related to the new political debates on the new political administration not only in the Ottoman lands, but also in various parts of Europe where these intellectuals traveled in this century. Kokdas also explains how all of these intellectuals and their production of knowledge developed alongside religious, ethnic and language ties, patronage networks, and land and money ownership. This century also experienced a democratization of knowledge production through the opening up of high culture to a large number of people. Anika Meckesheimer's paper deals with two levels of the circulation of knowledge: on the one hand she elaborates on the circulation of knowledge produced and used by activists and academics and on the other hand on the circulation between Mexico and Germany. In her exploration of de-colonial research, she argues for a view beyond academia and beyond Europe for new models: She describes feminist debates in Mexico, where indigenous women refuse to let their words be "taken away" by researchers. In her opinion, it is mostly due to structural obstacles that de-colonial research is so hard to achieve. She concludes that universities and research institutions should rethink their requirements for research (including theses) and offer the structural frameworks that allow for participatory research.

All of these papers draw our attention to knowledge production that is happening beyond today's supposed academic centers and classical canons in the social sciences. They stress the importance of locality, translatability and social embeddedness of knowledge production. The contributions of activist groups, ethnic or religious minorities, small towns and non-Western merchants or politicians broaden our perspective on what social science knowledge can mean. Furthermore, they make us question the social sciences' monopoly on the interpretation of society. Especially the connection between academic and non-academic actors in producing knowledge through dialogue and engagements is to be highlighted here. The perception of a variety of actors and locations as 'valid' persons and places of knowledge production and the acknowledgment of the various and often entwined paths of circulation of knowledge is on the agenda of all the presented authors. This is their major contribution to our understanding of what sociology on a global perspective might look like.

This workshop was part of a broader conference on knowledge circulation from which another publication with a related but slightly different focus emanated: 'Global Knowledge Production in the Social Sciences: Made in Circulation'(Keim et al. eds.:

forthcoming).