domain: nlin sdu sdu.astr

text:
Jean-Baptiste Durrive

Madhurjya Changmai

Rony Keppens

Pierre Lesaffre

Daniela Maci

Georgios Momferatos

Swift generator for three-dimensional magnetohydrodynamic turbulence

3 ) solenoidal vector fields in minutes to hours on desktop computers. Our model is inspired by recent developments in 3D incompressible fluid turbulence theory, where a Gaussian white noise vector subjected to a nonlinear transformation results in an intermittent, multifractal random field. Our B × C model has only few parameters that have clear geometric interpretations. We directly compare a (costly) DNS with a swiftly B × C-generated realization, in terms of its (1) characteristic sheetlike structures of current density, (2) volume-filling aspects across current intensity, (3) power-spectral behavior, (4) probability distribution functions of increments for magnetic field and current density, structure functions, and spectra of exponents, and (5) partial variance of increments. The model even allows to mimic time-evolving magnetic and current density distributions and can be used for synthetic observations on 3D turbulent data cubes.

I. INTRODUCTION

Fluids and magnetic fields are usually turbulent, and researchers often need to model and analyze turbulent data. Since fully nonlinear, turbulent, analytic solutions to the Navier-Stokes (hydro) or the magnetohydrodynamic (MHD) equations are unavailable, the most common tool to construct realistic models is by means of direct numerical simulations (DNS), which are-unfortunately-extremely expensive resource-wise [1]. This led to the creation of online turbulence databases (e.g., [2] for the Johns Hopkins Turbulence Database) where selected snapshots of isotropic hydro turbulent fields up to 8192 3 size, or 1024 3 incompressible MHD states, are stored for web-based access. To date, DNS models of increasingly larger size provide the only means to verify theoretical scaling laws, which for MHD in particular, are still subject of lively contemporary debate (e.g., see [3]). MHD turbulence, especially in three-dimensional (3D) incompressible settings, is discussed in many modern textbooks (see, e.g., [START_REF] Goedbloed | Magnetohydrodynamics of Laboratory and Astrophysical Plasmas[END_REF][START_REF] Biskamp | Magnetohydrodynamic Turbulence[END_REF][START_REF] Galtier | Introduction to Modern Magnetohydrodynamics[END_REF]), and these invariably emphasize its scaling and shape in power spectra, and its typical current-sheet dominated visual appearance.

In the quest for finding "exact" solutions to the incompressible Navier-Stokes equations, [START_REF] Chevillard | [END_REF] suggested an explicit, concise, and yet efficient, analytical expression for a random field which shares many properties of experimental and * rony.keppens@kuleuven.be numerical incompressible hydrodynamical turbulence (see also [8][9][10][START_REF] Chevillard | Habilitation à diriger des recherches[END_REF][START_REF] Pereira | [END_REF][13][14][15]). Intermittency (i.e., non-Gaussianity) in this model stems from the fact that the random field is constructed as products, i.e., a nonlinear transformation of Gaussian white noises. For this reason this approach belongs to the mathematical field called "Gaussian multiplicative chaos," first formalized by [16]. In order to build similarly parametric models for astrophysical environments [17], recently suggested an extension of the aforementioned model to magnetized fluids, mimicking MHD turbulence. In these constructed random fields, their statistics are controlled by a couple of free, physically motivated, parameters.

The approaches above have a threefold ambition: The random fields must (1) resemble real data as much as possible, (2) be physically motivated, and (3) be as numerically efficient as possible, to be worthwhile compared to DNS. They are useful in many ways, e.g., to quickly generate synthetic data (effective, "surrogate," models), to characterize turbulent data with few parameters for observers or experimentalists, and for constructing nontrivial (i.e., with at least self-similar and small-scale structuring) initial conditions for DNS. In the currently latest HD [START_REF] Chevillard | [END_REF] or MHD [17] efforts of this kind, objectives (2) and (3) are satisfyingly fulfilled, as the models are constructed from the physics of vortex stretching and flux tube shearing, and numerically they are several hundreds of times less resource-consuming than DNS. As for objective (1) to resemble real turbulent data, in the hydrodynamical case all efforts have focused on the statistics of the fields but not on the shape of the structures. Hence, while many statistical properties of the random incompressible velocity fields are fairly realistic, their 3D visualizations are far less convincing.

We here present a path to solve this problem, i.e., to build very efficiently [objective (3)] random fields that visually resemble DNS results [objective (1)]. We do this here directly for the MHD case, where the challenge is to get both current and magnetic field vector quantities behave in DNS-like fashion. We name our model B × C, standing for "magnetic fields from multiplicative chaos." Our reasoning is purely geometric, in the sense that we motivate our parametrized transformation mostly from getting visual correspondence with 3D turbulent magnetic vector fields. In practice, these parameters also relate, in a yet-to-be-quantified fashion, to the physical processes of vortex stretching and shearing (as we will "deform" spiral patterns based on gradient fields). The geometric parameters also are inspired by, and impact on, the statistical properties of the 3D turbulent states, and we provide various quantitative comparisons further on, notably in terms of energy spectra.

Incidently, it is straightforward to also adapt our model to the two-dimensional (2D) case, by starting with the wellknown 2D Biot-Savart's law and keeping the eddy modeling 2D as in Sec. II B. The interested reader may have a look at for example [14] who do work with fractional Gaussian fields (see definition below) in two dimensions. However, the strength of our model lies on its 3D nature, since 2D DNS are fairly cheap to run and B × C is an interesting complimentary tool to DNS only in the 3D case.

The paper is organized as follows. In the first part of the paper, we detail the construction of our model. After giving some background, we construct a formula mimicking an isolated eddy in two dimensions, as a set of constant-curvature spirals swirling around a single point. Then, in an efficient single mathematical step, we extend this formula to 3D sheets, with nonuniform curvature, randomly distributed throughout space. We also expose how to straightforwardly emulate a time evolution of our turbulent magnetic field. In the second part of the paper, we show an example of a 3D vector magnetic field and its current density built with our model, and compare them to a modern DNS result. The comparison is performed in multiple ways, inspecting several visual aspects and by means of quantitative statistical tools.

II. MAGNETIC FIELD CONSTRUCTION

A. Preliminaries

Biot-Savart's law expresses a magnetic field B in terms of its current density j as the convolution

B = N B R 3 j × r r 3 dV, ( 1
)
where N B ≡ μ 0 /4π , with μ 0 the vacuum permeability. Inside all integrals we use the usual short-hand notations r ≡ xy and r ≡ | r|, not to be confused later with the 2D (r, θ ) polar r coordinate.

The basic structure of the models in [START_REF] Chevillard | [END_REF] and [17] is the modified version of Biot-Savart's law

B = N B r L c × r (r 2 + η 2 ) h dV. ( 2
)
Compared to (1), the integration region is restricted to a ball of radius L, the kernel's fixed r -3 power-law behavior is set to vary freely with a power h, and the kernel's singularity at r = 0 is regularized with η, so that parameters L, h, and η, respectively, serve to control the large-scale cutoff, the slope, and the small-scale cutoff of the power spectrum of B ≡ | B|.

The normalizing constant N B will be used to control the total energy of the field (moving vertically the power spectrum; cf. Sec. III D). Finally, we write c instead of j, because due to the above modifications, c in (2) is not exactly the current density j anymore. The strategy is to first construct c, then deduce B through (2), and only then deduce j ≡ ∇ × B from B. Most importantly, the form (2) guarantees B to be divergence-free, for any c (so ∇ • c = 0 is not required) as long as η is large enough for the field to be smooth on small scales such that gradients are well approximated (e.g., [START_REF] Pereira | [END_REF]). At the same time we better take η to be small, to have a large inertial range, and we take as a trade-off between these two constraints η = 3/N at a resolution N.

The core of this model is to choose a relevant c. Hereafter we call s a Gaussian white noise vector, the tilde symbol reminding its random nature and "s" standing for "seed." The three components of s are Gaussian white noises, independent of one another, zero-averaged, and with unit variance. The simplest idea takes c equal to s and (2) reduces to

R ≡ N R r L R s × r r 2 + η 2 R h R dV, (3)
a field referred to as a fractional Gaussian field (fGf) [14]. We renamed it to R for reasons that will become clear further on and added subscripts R to the parameters in (3) as they will have different numerical values than in (2). Now, magnetic fields in nature are non-Gaussian (intermittent), while R is still a Gaussian field because the Biot-Savart operation (3) is a linear transformation (namely, a convolution) on a Gaussian field (namely, s). To build an intermittent model, we must find a nonlinear transformation on s, which is the purpose of the two following sections.

B. Designing an isolated eddy

A characteristic feature of fully developed MHD turbulence is the omnipresence of current density sheets swirling around throughout space. This is-in a sense-the MHD counterpart of the vortex tubes from hydrodynamics, where eddies of varying sizes communicate nonlinearly. In the simplest nontrivial model to mimic a swirling current sheet, we are led to the Archimedean spiral (such spiral shapes may represent the outcome of some underlying physical process, such as a Kelvin-Helmholtz roll-up pattern in a shear flow (e.g., [START_REF] Goedbloed | Magnetohydrodynamics of Laboratory and Astrophysical Plasmas[END_REF]; this Archimedean spiral also returns as the Parker spiral of interplanetary magnetic field [18] in an important historical model for solar wind magnetic fields). The latter is described in the 2D plane by the polar equation

r(θ ) = c 0 + d θ, ( 4
)
where r and θ are the usual polar coordinates. The parameter c 0 moves the center point of the spiral outward from the origin, while d controls the distance between the spiral arms. Actually, in the top-left panel of Fig. 1, we rather show r p ≡ 1 -√

x 2 + y 2 instead of r (for pedagogical reasons only, to focus on a clump rather than a void region), together with θ ≡ 1 π atan2(y, x) where the atan2 function generalizes the relation θ = arctan(y/x) which holds only for x > 0. The 1/π factor simply keeps the field in the normalized range [-1, 1] for convenience.

To construct an actual spiral-shaped scalar field in the plane, we consider λ ≡ r pc 0d θ , a local length that measures how far a given point is from the spiral (4). This local length is then given as argument to a suitable filter, for example the smoothed top hat

T (s) ≡ 1 2 tanh s + w/2 -tanh s -w/2 , ( 5
)
which is a function such that T (s) equals 1 in a region of width w near the origin s = 0, and equals 0 elsewhere with a smooth transition from 1 to 0 of thickness controlled by the length . The field T (λ) is a field with a spiral shape, because T selects the regions of space where λ is close to 0, up to a certain width w. Most importantly, so far r p is a 2D field, but in the next section we will replace it by a 3D field related to R from (3), turning T (λ) into a 3D scalar field with spiral-shaped sheetlike structures. Anticipating this, we refer to T (λ) as a sheet. At this point, we obtain further guidance from the current density field j as obtained in actual DNS studies, where it appears relevant to distinguish two types of sheets in the modeling of turbulent magnetized flows. Indeed, they suggest clearly a bimodality in (1) intense (i.e., high | j| regions), which are thin, and relatively rare sheets (i.e., intermittent) and (2) more diffuse weaker | j| regions distributed in thicker, and more abundant sheets (i.e., more volume-filling), which surround the intense sheets.

Therefore, we define an intense filter T i and a diffuse filter T d , which are identical to T in (5) with differing numerical values for the parameters (w i , i ) and (w d , d ) respectively: w i and i are smaller than w d and d , to mimic the fact that intense sheets are thinner and less blurry than diffuse sheets. Second, as seen in the top-middle panel of Fig. 1, the atan2 function introduces a discontinuity where θ = ±1. As a simple work-around to avoid jumps in our magnetic field model, we impose a spatial dependence to the width w i of intense sheets through the prescription (recall that θ ∈ [-1, 1])

w i = w max i cos(πθ/2), (6)
where w max i is a constant. As a result, wherever θ is discontinuous, intense sheets become infinitely thin, and therefore vanish. Intense sheets are then also less volume-filling, and hence more intermittent, as it appears in DNS simulations. It turns out to be unnecessary to do the same for the width w d of diffuse sheets, because these sheets have weak amplitudes, so their discontinuities are smoothed out when taking the Biot-Savart law (2) in the last step of our construction. Finally, as illustrated in the top-right panel of Fig. 1, to model the fact that diffuse sheets are numerous and surround intense sheets, we generalize our filtering to

S(λ) ≡ T i (λ) + T d [cos(k d λ)]. ( 7
)
The first term corresponds to an intense sheet, and the second to several diffuse sheets. Indeed, instead of T d (λ) we consider

T d [cos(k d λ)],
which gives rise to as many sheets as there are zeros in cos(k d λ), i.e., k d controls the number of diffuse sheets. Moreover, in regions where λ ∼ 0 this cosine does not vanish so that diffuse sheets are absent, which adequately gives room to the intense sheet T i (λ) sitting there. Lastly, the free parameter , assumed to be small, makes diffuse sheets more diffuse than intense sheets by controlling their relative amplitude.

To sum up, for our 2D field r p and angle θ as in the first two panels of Fig. 1, S given by ( 7) is a field of nested, 2D spirals where the central one is intense, as shown in the top-right panel of that same figure. This constitutes the basic structure of an (isolated) eddy in our model. The key point of the next section is that we will insert in (7) a 3D (random) scalar field instead, such that S will indeed be a field of 3D sheets with artificially constructed spiraling behavior. Note that thus far, our spirals have constant curvature, to be remedied in what follows as well.

C. Randomly distributing eddies

We now present an efficient way (i.e., a simple single step) to simultaneously (1) extend from two to three dimensions the above considerations, (2) introduce nontrivial spatial variations of the curvature of the sheets, and (3) distribute eddies in the whole domain, with the properties of the sheets (size and wiggliness) controlled by a few parameters.

As mentioned in the preliminaries, our fractional Gaussian field R given by ( 3) is a poor stochastic model for a turbulent magnetic field. In the bottom-left panel of Fig. 1 we show a 2D cut of a realization of its norm, R ≡ | R|. The 3D scalar field R consists of an ensemble of nearly spherical clumps of various sizes, randomly distributed throughout space. The fact that this field does not resemble actual turbulent structures is related to the (visual appearance) shortcoming we alluded to in our introduction of present multiplicative chaos models for hydro turbulence. The clumpiness of R and the typical size of its largest clumps are readily controlled by the Hurst parameter h R and the cutoff L R in (3), respectively. Having noticed this, we will now use this clumpy field to build spiral-shaped structures swirling around intense clumps. Hence, we are not going to use R as a magnetic field vector B model, but as our foundation to build a current c, to plug in the formula (2) for B.

We now have a natural "radius field" R, but in analogy with the construction of 2D spirals, it remains to find a relevant angle θ . This is indeed possible noticing that we may also write θ = atan2(∂ y r, ∂ x r)/π , a relation that becomes clearer after checking that it does reduce to the standard arctan(y/x) for x > 0. With this viewpoint, it is now natural to define, for the 3D case,

θ R ≡ 1 π atan2(∂ y R, ∂ x R). ( 8
)
Finally, we redefine the length λ as

λ R ≡ R -c 0 -d θ R . ( 9
)
Our motivation for these peculiar definitions is purely geometrical, in the sense that we introduce them independently of the dynamical equations. However, an expression such as (8) should not be surprising, since dot products between fields and gradients (and therefore angles) are omnipresent in (magneto-)fluid dynamics, notably with the advection operator v • ∇. Note that, in this 3D case, we could likewise consider a second angle, inspired from the φ angle of spherical coordinates, but we deliberately keep our model as elementary as possible.

All in all, our magnetic field model B is the modified Biot-Savart law (2) with the "current" vector field in it taken as

c ≡ S R, (10)
i.e., c starts from the fractional Gaussian field R given by (3), scaled by a sheetlike field with a spiral structure S given by [START_REF] Chevillard | [END_REF], where the top-hat functions T i and T d are given by ( 5), the angle θ R by (8) and the length λ R by (9). We name our model B × C, which stands for "magnetic fields from multiplicative chaos" in reference to notably [START_REF] Chevillard | Habilitation à diriger des recherches[END_REF]16,17].

We can motivate our construction as follows. Evidently, the core of turbulence studies is to understand and be able to model the intricate interactions between scales in turbulent fields. A classical paradigm is to consider as total field a split into a sum of fields of different nature, e.g., constituted as an ordered (strong background) plus a turbulent field, or an equilibrium plus a perturbed field. An archetypical example is the mean-field dynamo theory where the magnetic and velocity fields are split into large-scale, mean-field parts and small-scale, fluctuating parts [19]. In this paper, we introduce another procedure when we use the fGf R field. We effectively introduce a scale splitting linked to the correlation length scale of R: inside each "blob" of R (cf. bottom-left panel of Fig. 1) a spiral-shaped eddy forms, while on larger scales, beyond R's correlation length, the eddies decorrelate. Since we expect the statistics of our field to become Gaussian on large scales (see also PDFs of increments further shown in Fig. 9 below), it seems appropriate to use a Gaussian field, such as a fGf.

D. Mimicking a time evolution

A particular feature of the present type of modeling is that it consists in applying a deterministic transformation to a given white noise. Being deterministic, once a realization of the white noise is chosen, we can transform the magnetic field smoothly by varying continuously the parameters (L, h, η, , . . . ). This can be used to emulate a (artificial) time evolution: to each parameter p we give a simple time dependence p = p + σ p sin(ω p t + φ p ), i.e., the value of p oscillates around a mean value p, with an amplitude σ p , at a frequency ω p , and a phase shift φ p . It is paramount to choose different phase shifts for the various parameters. The oscillations will then be out-of-phase, which avoids spurious periodicities. In other words, we thus move continuously in a rather chaotic way into the parameter space. An animation exemplifying this can be found at [20] for a 512 3 resolution.

III. COMPARING OUR MODEL TO A DNS

In this section we analyze a realization of a magnetic field built with our B × C model, as well as its corresponding current density field, and we compare them to a realization of a magnetic field and current density generated using a DNS, to assess the realism of our model.

Note from the outset that we expect our model to be primarily useful (1) to generate extremely high-resolution fields (including a mock time evolution) that are out of reach of DNS and (2) to reduce drastically the time needed to create nontrivial initial conditions for DNS. We therefore will assess whether our model can reproduce with much reduced resources various aspects also present in a given DNS. It is to this end that we ran a full DNS. The latter will constitute some reference data, considered as "realistic," and in this part of the paper we show by means of a series of side-by-side comparisons, that our model shares many properties of this DNS, both qualitatively (notably sheetlike structures with appealing visual aspects) and quantitatively (notably providing evidence of intermittency, and the expected shape for power spectra, namely, a well-defined power-law-behaving inertial range between clear large-and small-scale cutoffs). Naturally, since our model is a fast parametric model, future work could easily extend it with an automated systematic parameter survey, such as Monte Carlo Markov chain analyses.

This part of the paper is organized as follows. We give details of how we implemented numerically our DNS and our model magnetic fields. We carry on by comparing the DNS and B × C fields in five ways. First we compare the resources required to generate them, then we inspect their visual aspects (2D slices as well as 3D appearance, with both scalar and vector visualizations), after which we provide several quantitative comparisons using the standard statistical tools of turbulence studies, namely power spectra, PDFs of increments for B and j, structure functions and spectra of exponents, supplemented with a Partial Variance of Increments analysis.

A. Numerical implementation

Throughout the paper, unless otherwise stated, the fields have a resolution corresponding to N 3 = 1024 3 collocation points.

The DNS data set considered is a snapshot at the temporal peak of total dissipation from a pseudospectral simulation of decaying 3D isotropic MHD turbulence that was performed with the ALIAKMON code [START_REF] Momferatos | The sites of extreme turbulent dissipation in the diffuse interstellar medium: Structure & properties[END_REF]. The nonlinear terms in the equation were de-aliased using the standard two-thirds rule, while advancement in time was performed by a fourthorder Runge-Kutta method. The product of the maximum wave number that was represented in the simulation with the Kolmogorov microscale was at all times kept above 2. At the temporal peak of total dissipation, the Taylor microscale Reynolds number is approximately equal to 270, while the Reynolds number based on the integral length scale is approximately equal to 2000. The initial condition used is a superposition of a large-scale Arnol'd-Beltrami-Childress (ABC) flow at |k| = 2 and a Gaussian random field with an exponentially decaying energy spectrum.

For B × C, we detail the reasoning that lead us to the chosen numerical values of the parameters. Note that in our code the box size is taken equal to unity, so the numerical values of the lengths below should be read as percentages of the box size.

First we chose the parameters for the fGf R, because the fGf directly controls the typical size of the large eddies, as illustrated by the bottom panels of Fig. 1. In order to obtain about ten large eddies along each direction of our data cube, we took L R equal to about a tenth of the box size, specifically L R = 0.075. Then, for the inertial range to be as large as possible, we needed to choose η R as small as possible, but as previously mentioned, at the same time this parameter should be large enough for the gradients of this random field to be well approximated. A usual trade-off in such models (e.g., [START_REF] Pereira | [END_REF]17]) is to take η R = 3dx, where dx = 1/N is the pixel size on our grid of size N = 1024. Finally, the Hurst parameter h R of the fGf controls how smooth R is, and consequently, given the construction, it controls how wiggly the sheets are (cf. bottom panels of Fig. 1). Since in our reference DNS data the sheets are particularly smooth, sometimes even almost flat, we were led to choose a very small numerical value for this Hurst parameter, and took h R = 0.05.

For the Biot-Savart law (2) we took L = 0.3 to integrate on sufficiently large regions for the magnetic field to span on large scales, as in our DNS. The choice h = 2 was based on enabling the magnetic field to have enough power at small scales, since Hurst parameters control the slope of the power spectrum, and the standard value η = 3/N was chosen as for η R above.

For the properties of the sheets, we focused on the parameters controlling the spiral shapes. Given the properties of Archimedean spirals, we chose c 0 = 0.3 to offset the centerpoint of the spirals from the origin to avoid having sheets converging artificially at the same points, and d = 0.2 for the sheets to be well separated, as in the DNS. Otherwise, in order to make our intense sheets very thin, as in our reference DNS, we chose a width several orders of magnitude smaller than the box size, namely w max i = 3 × 10 -5 and i = 5 × 10 -3 . We then constructed the diffuse sheets relatively to the intense ones: In our reference data diffuse sheets appear typically an order of magnitude wider, hence w d = 0.05, and being "diffuse" translates into d = 0.2 to be an order of magnitude larger than in intense sheets (the filter thus being far less steep). In addition, from [START_REF] Chevillard | [END_REF] it is clear that k d controls the number of sheets swirling inside a given eddy, measured in multiples of 2π . The choice k d = 6π leads to a few diffuse sheets and an appropriate volume-filling aspect. To make sheets diffuse and hence less intense, it is natural to weigh their amplitude relative to the intense sheets by a number of the order of a percent, hence = 5 × 10 -3 .

Finally, in Biot-Savart's law we smoothed the truncation of the integration region to a ball of radius L, by multiplying the kernel by 0.5[1tanh (r -Lc )/ c ] with c = 0.23L. Second, noticing that in places where | ∇R| ∼ 0, such as in the origin of the plane in the top row of Fig. 1, many sheets converge in a relatively artificial manner, we multiplied S by 1exp(-| ∇R| 2 /10), and found that this improved slightly the results.

B. Comparison 1: Required resources

The important difference between B × C and our DNS run is the resources used: the reference DNS required about 50 000 core hours (on an HPC system with eight-core Intel E5-2670 Xeon processors running at 2.60 Ghz) while a magnetic field with our code is generated in less than half an hour on a desktop with 40 logical cores, as detailed in Fig. 2. Our model is very fast because (2) and (3) are nothing but convolution products, i.e., simple products in Fourier space. In contrast to the original HD and MHD models of this kind, where much more intricate nonlinearities were used to mimic turbulence statistics, this aspect makes our geometric, parametrized construction scalable to extreme resolutions, beyond those achievable by DNS on modern supercomputers, and only bound by local memory requirements. Order 500 3 realizations FIG. 6. Statistical comparison of the DNS and B × C fields. In the top panel we plot the power spectrum P(k) of the norm of the magnetic field, from three orthogonal 2D slices passing through the center of the data cube. Continuous lines correspond to B × C and dotted lines to the DNS, where in red the data used are from the slice with fixed x = 512 (the resolution being N = 1024), in green with fixed y = 512 and in blue with fixed z = 512. Red, green, and blue curves of a given data set match because the fields are statistically isotropic. The bottom panel is the same with the current density field. The important point is that the power spectra of the B × C fields have the characteristic shape of turbulent fields, with a clear power-law inertial range delimited by a large-scale cutoff at small k and a smallscale cutoff at large k.

are feasible on any laptop, while modern desktops can easily generate far larger fields.

C. Comparison 2: Visual aspects

In Figs. 345, the left columns correspond to the DNS and the right columns to our B × C model. The first row of Fig. 3 shows the magnetic fields, while all the other figures correspond to the current density fields, which B × C aims at reproducing. We insist that j here is computed, as it should, by taking the curl of the magnetic field (2): it does not simply correspond to c given by (10), because ( 2) is a modified Biot-Savart formula.

In Fig. 3 we start by exhibiting 2D slices of the norms of B and j, in the top and middle rows respectively. It appears that the B × C fields have a fluid aspect in the sense that FIG. 7. Same as Fig. 6, but where the spectra have been compensated.

smooth variations alternate sudden concentrated structures. The overall size distribution of larger and smaller patches, in both the magnetic field magnitude variation and in the current intensity, is fairly similar between DNS and our model. An advantage of this construction is that the properties of the sheets are easily controlled by a few parameters: c 0 and d in λ given by ( 9) control the (deterministic) shape of individual spirals constituting the sheets, while L R and h R in the fGf R given by (3) control, respectively, the typical size of the swirling regions and how wiggly the sheets are. For example, in the B × C realization shown here, we deliberately chose a very small numerical value for the Hurst parameter h R . In this way R is very smooth (cf. bottom-left panel of Fig. 1) so that the sheets are not very wiggly, as we observed in our reference DNS data.

Then, in the bottom row of Fig. 3, to provide some vector information, we show a 2D slice (the same as for the two rows above) of j x , the x component of j. In the blue regions j x is positive, while it is negative in the red regions. Comparing the DNS and our analytic reproduction (the left and right columns) we conclude that B × C reproduces, qualitatively and statistically speaking, the orientation of the sheets. This is an extremely important finding, since we targeted this 3D turbulent vector correspondence from the outset, not just a scalar reproduction. In Fig. 4 we show isocontours of j at 60%, 30%, and 10% of its maximal value, in the top, middle, and bottom rows respectively. These 3D visualizations confirm that the B × C current density field is indeed composed of 3D sheets with nontrivial shapes (nonuniform curvature and wiggly edges). The distribution (i.e., the positions, the orientations and the volume-filling aspect) of the sheets is rather realistic, in the sense that intense j regions are not volume-filling, which is one known facet of MHD intermittency.

In Fig. 5 we show yet more vector information, complementing the bottom panel of Fig. 3. The top row is a zoom into a 3D vector visualization, while the bottom row is a zoom on a 2D vector visualization, both displaying regions with many sheets. We again conclude that the look and feel of B × C is convincing, and it should be noted that we have not yet attempted to optimize the free parameters involved in any way. This can probably be done in follow-up work, but it is to be stressed that we can easily generate many realisations within hours on desktop resources, which in principle are equally likely, just by changing our starting Gaussian noise model.

D. Comparison 3: Power spectra

Finally, in Fig. 6 we computed detailed statistical information to be more quantitative. In the DNS and B × C code simulations, we consider the power spectrum as a tool to quantify and compare the statistics of the scale dependence of the fluctuations. The power spectrum (P) is defined as the change in kinetic energy (E ) as a function of wave number (k), P(k) = dE/dk. From the isotropic incompressible 3D data of the field, we generate a one-dimensional (1D) radially averaged power spectrum [START_REF] Miville-Deschênes | [END_REF] from 2D slices along coordinate directions. The 2D field f (x, y) is Fourier transformed, yielding the 2D power spectrum from the amplitude defined as P(k x , k y ) = | f (k x , k y )| where f denotes the Fourier transform of the field. The collapsed 1D radial average of P(k x , k y ) between k and k + dk, where k = k 2

x + k 2 y , yields the power spectrum P(k)dk. This is shown in Fig. 6. The fact that the red, green, and blue curves of a given data set-which differ in their slice orientation-overlap each other [23], stems from the statistical isotropy of the fields. This behavior would obviously change if we were to combine a B × C prescription with background guide fields, or given spatiotemporally varying, smooth background magnetic field models. This figure shows that the B × C fields, in particular the magnitudes B and j, have similar power spectra than that of the DNS. Indeed, they provide a clear proof of concept, and produce fields with power spectra that have the characteristic shape of turbulent fields, namely a large-scale cutoff at small k corresponding to the injection scale, a power-law inertial range at intermediate k which, physically speaking would correspond to the energy cascade, and a clear small-scale cutoff at large k which mimics the effects of dissipation.

Note that we can always scale the magnetic field strength in the B × C to match the DNS power spectrum at a specific length scale, notably with N B . The comparison between power spectra in Fig. 6 is repeated in compensated form in Fig. 7, showing a very acceptable level of agreement, given that no parameter optimization has been performed. The numerical values of the powers in k for the compensations were chosen such that the regression lines of the inertial ranges fit become horizontal. In Fig. 8 we provide three examples of how the power spectra vary when varying the values of the parameters of our model, and how sensitive they are to such variations. Specifically, in the top, middle, and bottom panels we vary, respectively, h, η, and L R while keeping all the other parameters to their values of the reference run. We indicate with black arrows how varying these parameters may help tweak a given power spectrum: The Hurst parameter h is a convenient degree of freedom to modify the slope of the spectrum, while η and L R enable refining the cutoffs at the small and large scales, respectively. Figure 8 simply illustrates there are enough degrees of freedom in our model to fit DNS data rather precisely, but this possible optimization is out of the scope of this paper.

E. Comparison 4: PDFs of increments, structure functions, and spectrum of exponents

As turbulent fields are in general not Gaussian fields, power spectra cannot fully characterize a turbulent state. Hence, we now supplement our analysis with the most common tools of diagnosis in turbulence studies which reveal the existence of intermittent corrections to the scaling of the increments of the fields and their moments (structure functions and spectrum of exponents) with respect to length scale. Specifically, let us define the increment over a lag of the norm B of the magnetic field as the quantity δ B( x) ≡ B( x + ) -B( x). [START_REF] Chevillard | Habilitation à diriger des recherches[END_REF] In the following we will also consider the norm j of the current density field, and consider the same expression, replacing B by j.

A first traditional way to reveal intermittency is to compare the probability density functions (PDFs) of the increments of the considered field to those of a Gaussian field. Indeed, the PDFs of increments in intermittent fields undergo a continuous deformation as the norm of the lag is decreased, the PDF having an almost Gaussian shape at large lags but large tails at small lags. This behavior is a typical signature of intermittency, and the large tails are often called "non-Gaussian wings." Now, as we saw in the previous section, our fields (both DNS and B × C data) are statistically isotropic since the power spectra of the three slices (cuts along x, y, and z directions) are very close to one another. Therefore only the norm of the lags matters, and here we compute the PDFs of the increments for = 4, 7, 10, 13, 17, 22, 30. For larger s the PDFs are near Gaussians. In addition, we use this isotropy to improve our statistics as follows. In practice we compute the PDFs for each direction considering the x, y, and z slices as independent realizations of a single process, and we show in Fig. 9 the median PDF at each lag, with the gray areas indicating the standard deviation from this median. In this figure, the top row corresponds to PDFs (normalized to unit variance) of B and the bottom row to PDFs of j, the left and right columns corresponding to the DNS and B × C data, respectively. In all those plots the departure from Gaussianity is evident, with the aforementioned characteristic continuous deformation when varying the lag. Comparing the two columns, it is manifest that the B × C data do reproduce well the statistics of the reference DNS.

A second usual way to identify intermittency in isotropic turbulence studies is to analyze the so-called structure functions, and reveal their power-law behavior with respect to scale. Let us define the nth order structure function as the nth moment of the absolute value of magnetic field increments, namely,

S n ( ) ≡ |δ B| n , ( 12
)
where brackets • denote the expectation value [START_REF] Frisch | Turbulence. The Legacy of A. N. Kolmogorov[END_REF]. A similar expression stands when using the current density j instead of B. In practice we compute the structure functions using the above PDFs of magnetic field increments, and we show them in the top row of Fig. 10. In this figure it appears that in the range of lags considered here the seven first structure functions do behave as power laws. We find that we do not need to invoke extended self-similarity, which consists in considering structure functions as functions of the third order structure function S 3 , rather than of the lag in order to widen the power-law-behaving range. Hence, we have S n ∝ ζ n , where ζ n is called the spectrum of exponents (in the following we will normalize it with the third exponent, i.e., we will consider ζ n /ζ 3 ). The dependence on n of ζ n quantifies the intermittency, as the field is intermittent if and only if ζ n depends nonlinearly on n. The spectra of exponents deduced from our DNS and B × C magnetic fields are shown in the bottom row of Fig. 10. As they should, they strongly deviate from Kolmogorov's linear scaling. From these plots it appears once more that both magnetic fields are non-Gaussian and that they have similar statistical properties.

F. Comparison 5: Partial Variance of Increments: Correlating intermittent current sheets with discontinuities in magnetic fields

Figure 9 quantified increments in norms of magnetic field and current density, further used in Fig. 10 for structure function analysis. Now we will use similar incremental magnetic field changes along a parametrised path (using path parameter s) written as B(s, s) = B(s + s) -B(s) and current density j(s, s) = j(s + s)j(s) to produce a statistical analysis to identify intermittent turbulent structures (i.e., current sheets) by analyzing the discontinuities present in magnetic fields. This time, we measure the normalized partial variance of increments (PVI) where • = (1/l ) l •ds denotes a spatial average over the entire length l of the path considered (concatenated paths across the domain), and s is the spatial lag. The square of the above quantity I s, s is referred to as PVI as given in [START_REF] Greco | [END_REF]26]. We follow the idea presented in [27] to detect discontinuities along a tangential 1D path traced within 2D simulation cuts along the X , Y , and Z directions of the 3D simulation data. As shown in the top of Fig. 11, we sample the simulation along this 1D path, which is 14 • with respect to X , Y , and Z directions for Z, X , and Y cuts, respectively. The periodicity of the data enables us to consider the entire 1D path along the domain, where the path reenters the opposite periodic side. We have adopted this 14 • angle such that the offset distance between the path reentry is greater than the integral scale of the data. Along this 1D path, we measure the PVI. This shows the correlation between current structures formed due to the turbulence and intermittent PVI events along each cut for the For lags s = 1, 10, 50, the PVI series calculated for the DNS and B × C are plotted for X , Y , and Z cuts in Fig. 12. The PVI series can easily measure the presence of intermittent events relating to current sheets or magnetic reconnection. In a turbulent flow, the non-Gaussian events fill up the space in addition to these very rare intermittent events, whose values lie above the standard deviation of the sample. By applying a threshold method to the PVI analysis of numerical simulations [27,28], found a direct correlation between PVI events satisfying the threshold parameters to the non-Gaussian and intermittent events of a flow. In our analysis, the threshold parameter θ is set to 3σ , where σ is the standard deviation calculated across the PVI series. The increment of the threshold parameter leads to separating even higher intermittent events from the sample. We find the distinct regions of intermittent (rarely occurring) and non-Gaussian events (frequently occurring) in the turbulent flow above and below this threshold, respectively. In Fig. 12 the PVI signal for DNS (B × C) data is plotted on the left (right) column for values of different lag, s = 1, 10, 50. Considering a threshold of PVI > θ, the smallest lag of s = 1 captures the highest intermittent events compared to higher lags for each cuts of the two cases. PVI > θ captures both intermittent and non-Gaussian events as we increase the lag and as such the information gets saturated with lower intermittent events which we see in all the plots. As shown in Fig. 12, we clearly expect the B × C to provide similar information about the discontinuities present in the magnetic field compared to the DNS data.

I s, s = | B(s, s)| 2 √ | B| , ( 13)
As a follow-up study, we present our analysis to correlate the intermittent events found by analyzing the increments of the magnetic field to the presence of magnetic reconnection or current sheet events. According to [26], the more prominent peaks of current density correspond statistically to more significant peaks of PVI. It is because of this that the PVI method can describe and identify the strong magnetic gradients. We compare the spatial signals of PVI 2 (in red) to J 2 / < J 2 > (in dashed green) for Z cut in Fig. 13. We analyze both the signals for a lag of s = 1. The reference DNS data show distinct statistical peaks to be in phase of the PVI signal and the current density. As so, the data demonstrate that the two quantities have a positive cross-correlation. The same can be interpreted for the B × C data, which presents similar statistical results in Fig. 13, demonstrating that the PVI method in this case is capable of successfully relating the magnetic field discontinuities to estimate the intermittency in current density for the B × C. A further statistical study to show the relation of PVI and current is shown by the joint PDF in Fig. 14. We plot the kernel density estimate of the joint PDF for PVI value compared to the J/J rms at the smallest spatial separation of s = 1 for the Z cut of DNS [Fig. 14 [Fig. 14(b)]. For both cases, a positive correlation is seen with the extreme values of PVI corresponding to the extreme values of current density and the bulk of the PVI population at lower PVI values corresponds to the lower current density values. The Pearson correlation coefficient between the variables is 0.62 for the DNS data and 0.71 for the B × C data. Thus, it shows how PVI helps in identifying these extreme events from magnetic discontinuities and in doing so relates them to the sharp gradients in current density effectively for B × C as it should be for the DNS data.

IV. OUTLOOK

We introduced a geometrically controlled, parametrized way to generate mock turbulent MHD fields, emphasizing the magnetic field and current density variations in 3D space as typically encountered in high-resolution DNS data. Our B × C model was shown to correspond visually, as well as statistically, with typical isotropic turbulent magnetic fields. In contrast to DNS models, our tool is not computationally intensive, and has direct parametric control on the spectral properties embedded in these turbulent fields. By generalizing this proof-of-concept to cases with also background organized fields, our model may become a direct tool for testing rivaling MHD (anisotropic) turbulence theories, and for inspecting their visual appearance. Potential applications of this tool are numerous [29], with the distinct advantage that laptop resources suffice. This can then quickly generate turbulent magnetic data cubes, to study, e.g., polarized light propagation through astrophysical turbulent media (Faraday effect); or for fitting our geometric parameters to match actual 3D DNS fields, that can then be artificially "upscaled" to ever larger sizes N 3 . Note that we can likewise generate pure hydro fields, where vorticity-velocity vectors behave like our current-magnetic vector fields, and hence produce data cubes for both incompressible flow and magnetic field vectors, for input to full MHD simulations with particular turbulent properties. Future work can try to generate a suitable generalization of this model for isotropic MHD turbulent fields, to those encountered in situations with a clear organized guide field, where differences in behavior parallel versus perpendicular to the guide field can be explored.

FIG. 1 .

 1 FIG. 1. From deterministic spirals in two dimensions to random sheets in three dimensions. (Top row) 2D setup: Using the deterministic r p field in the left panel and the θ field in the middle panel, we construct with (7) the spiral-shaped field S on the right. This S could mimic an isolated eddy. (Bottom row) 3D setup, generalizing the top row: Using the random R field [norm of (3)] in the left panel and the θ field (8) in the middle panel, we construct similarly the field S with swirling sheets on the right. This S is used in (2) to mimic a distribution of eddies.

FIG. 2 .

 2 FIG. 2. Resources required to generate a magnetic field realization with our B × C code (written in Python). Computing time as a function of resolution is plotted in blue (left y axis), and the required RAM memory in red (right y axis), performed with a 40-logicalcores desktop. Continuous lines are for double precision (float64) calculations, and dashed lines for single precision (float32). Hence, 1024 3 data are generated in about 10 to 25 min, depending on the precision needed.

FIG. 3 .

 3 FIG. 3. Figures 3-5 are visual comparisons of fields generated with a DNS (left column) to fields generated with our B × C model (right column). In the present figure, (a) in the top row are slices of B, the norm of the magnetic field, (b) in the middle row are slices of j in logarithmic scale, the norm of the current density, and (c) in the bottom row are slices of j x , the x component of the current density, which shows some vector information (orientation of j). From Figs. 3-5 we conclude that, while B × C fields are generated using several orders of magnitude fewer resources, they have a similar visual aspect to the DNS.

FIG. 4 .

 4 FIG. 4. Continuation of Fig. 3. Isocontours of j, the norm of the current density, are shown for values of 60% (top row), 30% (middle row), and 10% (bottom row) of the maximal value. The volume filling and the shape of the structures of the B × C field at different amplitudes of j matches qualitatively that of the DNS.

FIG. 8 .

 8 FIG.8. Three illustrations of how the power spectrum of the norm of a B × C-generated magnetic field changes when varying some parameters of the model (namely, h, η, and L R from top to bottom panels), while keeping the other parameters to their reference values. These examples were made at a resolution N = 512, and, as in the rest of the paper, the parameter L R is measured in box-size units and dx = 1/N. The black arrows suggest how tweaking these parameters may help fitting a given power spectrum.

FIG. 9 .

 9 FIG. 9. (a) Top row: PDFs of increments of the norm of the magnetic field generated with our DNS (left column) and that generated with our B × C model (right column), at lags = 4, 7, 10, 13, 17, 22, 30. The dotted black curves correspond to unit-variance Gaussian PDFs. As the lag decreases, the curves deviate from Gaussianity, which is characteristic of intermittency. (b) Bottom row: Same plots using the norm of the current density field instead of B. The fact that the plots on the left and right columns look like each other indicates that B × C generated fields have rather realistic statistical properties.

FIG. 10 .

 10 FIG. 10. (a) Top row: The first seven structure functions S n ( ) of the norm of the magnetic field generated with our DNS (left column) and that generated with our B × C model (right column), obtained using the PDFs from Fig. 9. (b) Bottom row: Spectra of exponents resulting from fitting the power-law behaviors of the structure functions shown in the top row. The dashed blue lines correspond to the spectrum of exponents of a nonintermittent field (Kolmogorov scaling), while the dashed black curves correspond to the DNS and B × C data in the left and right panels, respectively. The pronounced departures of the black curves from the linear laws manifest intermittency, and the similarity between the plots of both columns highlights B × C's ability to mimic DNS data.

FIG. 11 .

 11 FIG. 11. (a) An example of our 1D data path in a X -Y cut used in a PVI quantification, and (b) the locally sampled data series below as function of path parameter s.

FIG. 12 .

 12 FIG.12. PVIs computed for s = 1, 10, 50 (top, middle, bottom). The dashed horizontal lines indicated for each cut (marked in the same color) represent the threshold of PVI θ = 3σ , calculated separately for each PVI series. data sets. The PVI events for each separate tangential path on the sampled DNS data set are shown in the bottom of Fig.11.For lags s = 1, 10, 50, the PVI series calculated for the DNS and B × C are plotted for X , Y , and Z cuts in Fig.12. The PVI series can easily measure the presence of intermittent events relating to current sheets or magnetic reconnection. In a turbulent flow, the non-Gaussian events fill up the space in addition to these very rare intermittent events, whose values lie above the standard deviation of the sample. By applying a threshold method to the PVI analysis of numerical

FIG. 13 .

 13 FIG.13. PVI 2 and the square of the J normalized to its mean value is plotted as a spatial signal with lag s = 1 for the Z cut in (a) the DNS data and (b) the B × C data.

  FIG. 14. Kernel density estimate of the joint PDF of the magnitude of the current J/J rms and the PVI values for the detected PVI events. The PVI signal has been computed on a spatial separation s = 1 and the Pearson correlation coefficient is 0.62 for DNS (a) and 0.71 for B × C (b).

ACKNOWLEDGMENTS

We thank the referees for their constructive comments. R.K. and J.B.D. are supported by Internal Funds KU Leuven, Project No. C14/19/089 TRACESpace. R.K. received further funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant No. 833251 PROMINENT ERC-ADG 2018) and a FWO Project No. G0B4521N. P.L. acknowledges support from the European Research Council, under the European Community's Seventh Framework Programme, through the Advanced Grant MIST (FP7/2017-2022, No. 742719).


Vincenzo Di Paola
email: vincenzo.dipaola@edu.

Stéphane Caro

Matteo Zoppi

DETC2023-XXXX SENSITIVITY ANALYSIS OF TENSION DISTRIBUTION ALGORITHMS FOR CABLE-DRIVEN PARALLEL ROBOTS

Keywords: Sensitivity, Tension Distribution, Cable-Driven Robot NOMENCLATURE u 𝑖 Cable direction Equality constraints g Inequality constraints 𝝁, 𝝀 Lagrangian multipliers ℒ Lagrangian 𝜎 Sensitivity index

Cable-Driven Parallel Robots (CDPRs) constitute a class of robots which displace a load by means of cables. In particular, the use of cables to guide the end-effector continues to be a challenging task. Indeed, several Tension Distribution Algorithms (TDAs) have been developed so far to properly steer the load. Although comparisons between them have been made, no performance index is established to evaluate numerically such methods. Therefore, this paper aims at defining a performance index establishing the sensitivity of a TDA by looking at the maximum tension variation that stems from a change in the external wrench. This would help both comparing and choosing the proper TDA among the existing ones. The definition of the index pops out naturally from the linearization of the equations defining the optimization problems. Commonly used TDAs are compared using the defined index. In order to show an application of it, a planar CDPR is considered and its Wrench-Feasible Workspace (WFW) is characterized by the index defined.

INTRODUCTION

Nowadays, Cable-Driven Parallel Robots (CDPRs) constitute a widely known class of robots where rigid links are replaced by flexible cables [START_REF] Bruckmann | Cable-Driven Parallel Robots[END_REF]. This design allows reducing costs and inertia of the system thus guaranteeing the ability to displace a load quickly and within a large workspace [START_REF] Albus | The NIST ROBOCRANE[END_REF][START_REF] Cone | Skycam: An Aerial Robotic Camera System[END_REF][START_REF] Bosscher | A Concept for Rapidly-Deployable Cable Robot Search and Rescue Systems[END_REF][5]. These peculiarities make them suitable for many applications ranging from common pick-and-place [6] to search and rescue in natural disasters [START_REF] Bosscher | A Concept for Rapidly-Deployable Cable Robot Search and Rescue Systems[END_REF], rehabilitation and collaborative purposes [7].

In particular, the load is steered through the use of cables and winches driven by electric motors; a simplified scheme of these robots is proposed in Fig. 1. The tensions exerted along the cables must be positive and bounded above to avoid controlling issues while accomplishing a task. Thus, choosing tension becomes a problem of crucial importance for addressing a given task. Several Tension Distribution Algorithms (TDAs) have been developed so far to solve mentioned problems [8][9][START_REF] Pott | Closed-form Force Distribution for Parallel Wire Robots[END_REF][START_REF] Taghirad | An Analytic-Iterative Redundancy Resolution Scheme for Cable-Driven Redundant Parallel Manipulators[END_REF][START_REF] Gosselin | On the determination of the force distribution in overconstrained cabledriven parallel mechanisms[END_REF][START_REF] Paola | Analytic Center Based Tension Distribution for Cable-Driven Platforms (CDPs)[END_REF]. Most of the time, because of different requirements dictated by the task, the solution to this problem relies on formulating an optimization problem. It is then clear the relevance of this topic for this class of robots and therefore why big efforts have been made hitherto in this direction. The main peculiarities of existing methods were summarized and compared in [START_REF] Pott | Closed-form Force Distribution for Parallel Wire Robots[END_REF][START_REF] Paola | Analytic Center Based Tension Distribution for Cable-Driven Platforms (CDPs)[END_REF]. However, none of the previous works investigated the sensitivity of these approaches to external disturbances or uncertainties. This question arises naturally since it is of interest assessing how much changes in the parameters of an optimization problem affect the optimal solution [START_REF] Castillo | Sensitivity analysis in optimization and reliability problems[END_REF].

In general, sensitivity analysis was already employed for CDPRs. The purposes were somehow different w.r.t. what this paper discusses. Indeed, most of the presented works, regard the evaluation of kinematic performances [START_REF] Patel | Manipulator performance measures-a comprehensive literature survey[END_REF]. Generally speaking, the first indices defined were the manipulability index [START_REF] Yoshikawa | Analysys and control of robot manipulators with redundancy[END_REF] and the dexterity index [START_REF] Salisbury | Articulated hands: Force control and kinematic issues[END_REF] well-known in kinematics. Today, many other works have been done in the kinematic context [5,[START_REF] Olds | Global indices for kinematic and force transmission performance in parallel robots[END_REF][START_REF] Cardou | Kinematic-sensitivity indices for dimensionally nonhomogeneous jacobian matrices[END_REF]. In particular, understanding how much cable errors (or joint position errors) affect the pose of the end-effector was the key aspect of mentioned contributions. In fact, the definition of different indices makes possible the comparison of various architectures under different perspectives. Thus, these indices help engineers during the design process of the robots.

Other works focus, for instance, on the effects of the assembly errors and dimensional tolerances on the performance of the robots [START_REF] Aref | Dynamic and sensitivity analysis of KNTU CDRPM: A cable driven redundant parallel manipulator[END_REF].

A recent work that aims to evaluate how a tension error in one cable affects the overall distribution of tensions in the other cables has been presented in [START_REF] Mattioni | Forcedistribution sensitivity to cable-tension errors in overconstrained cable-driven parallel robots[END_REF]. However, this methodology is developed regardless of the existing TDAs.

Hence, a preliminary analysis of the TDAs' sensitivity is hereby proposed. Indeed, the aim of this work consists in introducing an index measuring the maximum cable tension variation within the Wrench Feasible Workspace (WFW) [START_REF] Bouchard | On the ability of a cable-driven robot to generate a prescribed set of wrenches[END_REF] when the external wrench changes. In other words, this index helps estimate the maximum tension variation due to external disturbances. Moreover, it enables comparisons between TDAs. Therefore, as with kinematic indices, this study can guide the choice of the TDA to be used to fulfil a given task.

This paper is structured as follows. Section 2 recalls the main equations governing the CDPRs. Section 3 briefly describes the TDAs used for this work. In Section 4 the methodology used to define the index is described. Section 5 employs a planar CDPR to investigate the sensitivity of the TDAs. Finally, the results are summarized in Section 6.

SYSTEM MODELING

Let's consider CDPRs with 𝑚 cables and a load with 𝑛 Degree of Freedom (𝐷𝑜𝐹). Moreover, assume 𝑚 > 𝑛 thus, redundant CDPRs are considered. Generally, the Degree of Redundancy (𝐷𝑜𝑅) is defined to quickly characterize these classes of CDPRs as 𝐷𝑜𝑅 = 𝑚-𝑛 ≥ 1. Thus, 𝐷𝑜𝑅 returns the number of redundant cables w.r.t. the 𝐷𝑜𝐹. Now, the equilibrium of the load is given by

W𝝉 + w 𝑒 = 0, (1)
where w 𝑒 is the external wrench, 𝝉 is the cable tension vector and W is the wrench matrix1 which is defined as

W = ( u 1 . . . u 𝑚 ) , (2)
here u 𝑖 represents the 𝑖 𝑡 ℎ cable direction unitary vector. Thus, if 𝐷𝑜𝑅 ≥ 1, there exist infinite solutions of Eq.( 1) grouped in the following set

Σ = { 𝝉 | W𝝉 + w 𝑒 = 0 } . (3)
However, to maintain the equilibrium of the platform the cable tension limits have to be taken into account. Hence, the 𝑚-dimensional convex hypercube Π that defines the domain of the feasible tensions is

Π = { 𝝉 | 0 < 𝝉 ≤ 𝝉 ≤ 𝝉 } , (4)
where 𝝉, 𝝉 ∈ R 𝑚,+ are positive tension vectors limits containing the 𝑖 𝑡 ℎ lower and upper cable tension limits. Consequently, the set of feasible solutions Γ satisfying both Eq.(1) and Eq.( 4), is

Γ = Σ ∩ Π.
(5)

TENSION DISTRIBUTION ALGORITHMS

As anticipated, the problem of choosing a set of cable tensions among infinite solutions was extensively studied. Consequently, several TDAs exist, each with its own peculiarities. As expressed in [START_REF] Pott | An improved force distribution algorithm for over-constrained cable-driven parallel robots[END_REF], comparing them is not straightforward as there are many aspects, regarding the algorithms, to be considered. Generally, existing approaches are compared based on [START_REF] Pott | An improved force distribution algorithm for over-constrained cable-driven parallel robots[END_REF]: realtime capability, workspace coverage, continuity of the solution, degree of robustness2 and generality of the approach based on the achievable 𝐷𝑜𝑅.

However, many other aspects should be considered and this motivates the necessity to characterise each method succinctly and effectively. Although not all the features can be compared at once, some of them can be evaluated by means of an index; introduced in the following Section 4.

The aim of this section is to recall a few methods that will be investigated in terms of sensitivity. Among the existing approaches to compute cable tension solving an optimization problem, the choice lies in the Linear Programming [8], Quadratic Programming [START_REF] Taghirad | An Analytic-Iterative Redundancy Resolution Scheme for Cable-Driven Redundant Parallel Manipulators[END_REF], 4-norm Programming [START_REF] Gosselin | On the determination of the force distribution in overconstrained cabledriven parallel mechanisms[END_REF] and the Analytic Centre method [START_REF] Paola | Analytic Center Based Tension Distribution for Cable-Driven Platforms (CDPs)[END_REF]. These methods have distinct features in terms of robustness and continuity of the solutions that make them interesting to examine.

Linear Programming (LP)

The Linear Programming method is employed to ensure a robust solution and because of its computational speed in solving the tension distribution problem with generic (any 𝐷𝑜𝑅) CDPRs. In other words, the aim is to find a tension vector inside Γ as far as possible from the boundary of the hypercube Π. To achieve this, the objective function is built such that the distance 𝑑 (𝝉) of the tension vector 𝝉, from the facets of Π, is maximized. Formally, the optimization problem reads as

argmax 𝝉 ∈ Π 𝑑 (𝝉), W𝝉 + w 𝑒 = 0 (6)
where the objective function is

𝑑 (𝝉) = min{𝜏 -𝜏 1 , 𝜏 1 -𝜏, . . . , 𝜏 -𝜏 𝑚 , 𝜏 𝑚 -𝜏}. (7)
This formulation ensures reaching the maximum, among the existing techniques, robustness of the tension profiles. However, tension profiles are prone to discontinuities when a small change in the end-effector position occurs.

Quadratic Programming (QP)

A widespread approach solving the problem of tension discontinuity [8] while allowing reducing the energy consumption amounts to minimize the 2-norm of the cable tension vector argmin

𝝉 ∈ Π ||𝝉|| 2 2 . W𝝉 + w 𝑒 = 0 (8)
Despite the energy consumption is reduced, the system inherits low stiffness and therefore, the accuracy of the task can be affected by vibrations and oscillations of the platform. Moreover, the robustness of this method is the lowest possible as the solution always lies on the boundary of Π.

4-Norm Programming

An optimisation problem that turns out to be a suitable candidate for clarifying some of the characteristics of the index that will be defined in what follows, consists in minimising the 4-norm instead of the 2-norm (QP). The mathematical formulation takes the following form argmin 𝝉 ∈ Π ||𝝉|| 4 4 ,

W𝝉 + w 𝑒 = 0 (9)
which is quite similar to the previous one (QP). In other words, one could generically formulate all the above problems by using the 𝑝-norm 1 ≤ 𝑝 ≤ ∞. However, for sake of clarity, the methods were synthetically, although explicitly, reported.

Analytic Centre Solution

Guaranteeing robust, continuous and differentiable tension profiles for a system with generic 𝐷𝑜𝑅 is, generally, the purpose intended to be achieved by the various existing TDAs. However, they do not succeed in satisfying them all contemporary. This motivated the introduction of the Analytic Centre. Indeed, it was introduced as an attempt to take the strengths of the various methods presented in the literature. Essentially, to achieve mentioned goal, a new objective function was chosen, it exploits the properties of logarithms (or barrier functions), as explained in [START_REF] Paola | Analytic Center Based Tension Distribution for Cable-Driven Platforms (CDPs)[END_REF]. Formally, the objective function has the form

𝜙(𝝉) = 𝑚 ∑ 𝑖=1 -log(𝝉 -𝜏 𝑖 ) -log(𝜏 𝑖 -𝝉), (10)
which is one particular case of the general one while the problem can be formulated as argmin 𝜙(𝝉), W𝝉 + w 𝑒 = 0 [START_REF] Taghirad | An Analytic-Iterative Redundancy Resolution Scheme for Cable-Driven Redundant Parallel Manipulators[END_REF] Observe that, apart from ensuring continuity, differentiability and robustness of the tension profiles, it eliminates tension inequalities simplifying the problem formulation and numerical resolution.

The main differences among the reported methods can be visualized in Fig. 2 where the level-sets generated by the different objective functions are explicitly depicted. Their shape is fundamental to understand how the tension solution is selected when intersected with Eq. (1).

PERFOMANCE INDEX

The introduction of an index measuring the sensitivity of the solution provided by the TDAs is a necessary condition to make comparisons among them.

To reach this aim, it becomes of primary importance to estimate the variation of the tension vector 𝛿𝝉 w.r.t. a variation of the external wrench 𝛿w 𝑒 . Therefore, taking advantage of the Karush-Kuhn-Tucker (KKT) theorem [START_REF] Boyd | Convex optimization[END_REF], it becomes possible to transform a generic optimization problem into a system of equations of the form

𝑓 (𝝉, W, w 𝑒 ) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝜕ℒ 𝜕𝝉 = 0 h(𝝉) = 0 g(𝝉) ≤ 0, (12)
where ℎ(𝝉) and 𝑔(𝝉) represent the equality and inequality constraints, respectively while ℒ is the Lagrangian function, namely

ℒ = 𝜓(𝝉) + 𝝁 𝑇 g(𝝉) + 𝝀 𝑇 h(𝝉), (13)
with 𝜓(𝝉) : R 𝑚 → R objective function or cost function to be minimized (or maximized). It is generally related to the optimization problem considered and can assume different forms, e.g. expression of 𝜙(𝝉) Eq. ( 10) used in Eq. [START_REF] Taghirad | An Analytic-Iterative Redundancy Resolution Scheme for Cable-Driven Redundant Parallel Manipulators[END_REF]. Remaining symbols 𝝁 and 𝝀 are known as the KKT multipliers. Now, linearizing Eq.( 12) gives

𝜕 𝑓 𝜕𝝉 𝛿𝝉 + 𝜕 𝑓 𝜕w 𝑒 𝛿w 𝑒 = 0. (14)
Rearranging this equation, the link between 𝛿𝝉 and 𝛿w 𝑒 is found, in compact form

𝛿𝝉 = - ( 𝜕 𝑓 𝜕𝝉 ) † 𝜕 𝑓 𝜕w 𝑒 𝛿w 𝑒 , (15)
then, calling

S = ( 𝜕 𝑓 𝜕𝝉 ) † 𝜕 𝑓 𝜕w 𝑒 (16)
one has

𝛿𝝉 = -S 𝛿w 𝑒 , (17)
where S is the sensitivity matrix, formally

S = ⎛ ⎜ ⎜ ⎜ ⎝ 𝜕𝜏 1 𝜕𝑤 𝑒,1 . . . 𝜕𝜏 1 𝜕𝑤 𝑒,𝑛 . . . . . . . . . 𝜕𝜏 𝑚 𝜕𝑤 𝑒,1 . . . 𝜕𝜏 𝑚 𝜕𝑤 𝑒,𝑛 ⎞ ⎟ ⎟ ⎟ ⎠ 𝑖 = 1, . . . , 𝑚 𝑗 = 1, . . . , 𝑛 (18)
where the indices 𝑖, 𝑗 identify the components of S as well as those of the vectors 𝝉 and w 𝑒 . At this point, all the elements necessary for the definition of the performance index have been collected. Consequently, borrowing the idea from [START_REF] Olds | Global indices for kinematic and force transmission performance in parallel robots[END_REF] and considering 𝛿𝝉 as if it were the variable of the joint space while 𝛿w 𝑒 represents the term in the Cartesian space, respectively, the sensitivity index can be introduced as

𝜎 = max | | 𝛿w 𝑒 | | 2 =1 ||𝛿𝝉|| ∞ . (19)
This, in complete analogy with the kinematic indices [START_REF] Olds | Global indices for kinematic and force transmission performance in parallel robots[END_REF][START_REF] Cardou | Kinematic-sensitivity indices for dimensionally nonhomogeneous jacobian matrices[END_REF][START_REF] Mottola | Effect of actuation errors on a purely-translational spatial cable-driven parallel robot[END_REF]. It was proven that optimization problem Eq.( 19) can be directly solved substituting Eq.( 17) into Eq.( 19) to obtain

𝜎 = max | | 𝛿w 𝑒 | | 2 =1 ||S 𝛿w 𝑒 || ∞ = ||S|| 2,∞ . ( 20
)
where the matrix norm at RHS is known as mixed Hölder norm [START_REF] Bernstein | Matrix mathematics[END_REF].

Practically, since the columns of the matrix S identify the variation in the tension of each cable w.r.t. a variation of the wrench in a specific direction, this index estimates the maximum (upper bound) tension variation induced by a change in the external wrench. This can help quantify how much can be the variation in the tension along the cables. Observe that measuring cable tension experimentally is a hard task. Therefore, the possibility to gain some information (upper bound on 𝛿𝝉) based on knowledge of some disturbances can be useful when measuring cable tension with sensors. [START_REF] Mottola | Effect of actuation errors on a purely-translational spatial cable-driven parallel robot[END_REF]. [START_REF] Salisbury | Articulated hands: Force control and kinematic issues[END_REF] pops out to be similar to the closed-form solution 𝛿𝝉 = -W † 𝛿w 𝑒 . However, the procedure here developed can be considered a generalization of it in the context of sensitivity. Indeed, taking the derivatives of Eq.( 12) appearing in the definition of matrix S Eq.( 16), one can find the explicit form

Remark 4.1 Taking a point-mass as end-effector only requires taking into account the resultant force acting on the load. Hence no dimensional unit issue arises. In case the load is a rigid body, then two separate indices can be used

Remark 4.2 The linearized Equation

S = ⎛ ⎜ ⎜ ⎜ ⎝ ∇ 2 𝝉 𝜓 W I -I ⎞ ⎟ ⎟ ⎟ ⎠ † ⎛ ⎜ ⎜ ⎜ ⎝ 0 I 0 0 ⎞ ⎟ ⎟ ⎟ ⎠ , (21)
where ∇ 2 𝝉 𝜓 is the second-order gradient3 of the objective function 𝜓(𝝉) w.r.t. the tension vector 𝝉 while I identifies the identity matrix. Hence, from Eq.( 21), one can observe that S reduces to W † when no optimization problem is considered, i.e. when 𝜓(𝝉) and 𝑔(𝝉) are not considered. Observing this analogy stresses the familiar relationship between the solution provided by the closed-form (pseudo-inverse) and the QP approach. Indeed, the solution with pseudo-inverse leads to minimize the distance from the origin 𝝉 = 0. The main difference with QP, then is the absence of tension limits 𝝉, 𝝉.

NUMERICAL EXAMPLE

In this section, some simulations are performed to compare the TDAs recalled in Section 3. First of all, to evaluate the index, Copyright © 2023 by ASME the workspace of the robot is discretized and only the points inside the WFW are considered. Subsequently, for every feasible point, the index is computed. The value assumed by the index 𝜎 is assigned to the point considered and a map is generated within the workspace. This process is repeated for each TDA reported above.

In the following discussion, a planar CDPR with four cables is considered. The load is a point mass (𝐷𝑜𝑅 = 2) of mass 𝑚 = 1𝑘𝑔. Cable tension limits are set to be 𝝉 = 10𝑁 and 𝝉 = 100𝑁, respectively. The architecture of the robot and its WFW are depicted in Figure 3 for sake of clarity. Results of the simulations and the characterization of the WFW by means of 𝜎 are depicted in Fig. 4.

First, from Fig. 4 one can understand that the AC method results to be more sensitive than all the other methods. This means that, around each point of the WFW, the LP, QP and 4norm generate a ||𝛿𝝉|| that is lower than the one of the AC for a given ||𝛿w 𝑒 || 2 = 1.

To fully grasp why this happens, one should focus on both how the TDAs work and what happens inside the sensitivity matrix, i.e., how the index works. With this in mind, let's pair the methods by affinity in order to shed light on the results and discuss the above-mentioned aspects. Hence, consider the QP and 4-norm together and LP and AC consequently 4 . Now, starting from QP and 4-norm and analysing their results, one sees that some of the cable tensions assume the lowest values possible. In particular, for this case study, two tensions 𝜏 1 , 𝜏 2 (attached at the bottom of the robot) take the lowest bound 𝝉 as they do not bear any wrench. Therefore, only the remaining components 𝜏 3 , 𝜏 4 variate while the load moves inside the workspace. However, 𝜏 3 , 𝜏 4 vary as less as possible (because of how the problem is defined), they just have to maintain the equilibrium of the load. This explains why the sensitivity 𝜎 assumes similar values (𝜎 < 0.38) inside the WFW. However, looking at how QP and 4-norm approaches work, is not enough to understand why maps Figs. 4 (b)-(c) are different. Indeed, this is due to the values assumed by ∇ 2 𝝉 𝜓 inside matrix S. Computing the 4 The affinity can be understood by looking at the level-set reported in Figure 2.

second-order gradient for both

∇ 2 𝝉 𝜓 𝑄𝑃 = 2I 𝑚 ∇ 2 𝝉 𝜓 4-𝑛𝑜𝑟 𝑚 = 12 diag(𝜏 𝑖 ) 2 , (22)
shows that ∇ 2 𝝉 𝜓 𝑄𝑃 is constant, i.e., it does not depend on the solution 𝝉 while ∇ 2 𝝉 𝜓 4-𝑛𝑜𝑟 𝑚 does. As a consequence, the 𝜎-map is symmetric for the QP and it does not change if, for example, the external wrench changes.

Analogously, the same reasoning can be done for LP and AC. However, this time the values of 𝜎 are significantly different: the sensitivity of AC is far more higher than LP (QP and 4-norm as well). The reason why lies again in the second-order gradient (or matrix S). Indeed, in both cases the solution 𝝉 is robust but

∇ 2 𝝉 𝜓 𝐿 𝑃 = 0 (23)
results to be null while

∇ 2 𝝉 𝜓 𝐴𝐶 = diag ( 1 (𝝉 -𝜏 𝑖 ) 2 - 1 (𝜏 𝑖 -𝝉) 2 ) (24)
depends on 𝝉. Hence, because of the shape of ∇ 2 𝝉 𝜓 𝐴𝐶 , its values take weight inside S generating the highest 𝜎 among LP, QP and 4-norm.

Other useful data can be acquired by computing the (numerical) gradient of 𝜎 namely, ∇𝜎 reported in Fig. 5. This information can be relevant when the task requires precision: no vibration or oscillation of the platform. In other words, it can be useful when a small variation of ∇𝜎 is required. Indeed, in practice, ∇𝜎 quantifies the possible tension variation between two points of the workspace. Not surprisingly, the AC method has smaller ∇𝜎 than the others (smallest and homogeneous arrows inside the WFW). In fact, the AC homogeneity, in terms of ∇𝜎, is probably linked to the property of smoothness of the tension profiles. Recall that smoothness of tension profiles means that there is a small variation in the cable tensions from one point of the workspace to another.

Collecting the information acquired from this case study, one can conclude 𝜎 allows for estimating an upper bound of 𝛿𝝉, assessing robustness as well as observing the continuity of the 5

Copyright © 2023 by ASME solution looking at ∇𝜎. Hence, many of the aspects listed in Sec.3 can be somehow visualized with this index. Although the definition of 𝜎 is analogous to the wellestablished performance index for kinematics, some information is lost due to the linearization process. Indeed, with this definition, for example, it becomes impossible to assess the sensitivity of the single cables since it is absorbed by the use of the norm ||𝛿𝝉||.

CONCLUSION

A first performance index estimating the sensitivity of TDAs was defined in this paper. The results of the simulations and the maps reported show that the index is able to furnish information about the variation of tension due to a change in the external wrench. Moreover, it tells that the homogeneity of the 𝜎-gradient is linked to the ability of the TDA in providing smooth tension profiles. Although it seems to be an effective index that can help in choosing the TDA, some information is lost. Therefore, it will be necessary to dedicate future efforts to improve or edit the definition in order to recover information about single cable sensitivity. Other tests are also necessary to better understand the peculiarities of this index for instance one should consider the orientation of the platform and investigate a spatial case as well. Moreover, experimental tests should be undertaken to validate the proposed method.

FIGURE 1 :

 1 FIGURE 1: PLANAR CDPR: GENERIC PLANAR ARCHITECTURE OF A CDPR WITH A POINT-MASS LOAD.

FIGURE 2 :

 2 FIGURE 2: LEVEL SET: THE CONTOURS SHOW THE LEVEL-SETS INSIDE THE 2D BOX Π (RED LINES). IN PARTICULAR, FIGURE (A) DEPICTS THE LP METHOD, (B) THE QP METHOD, (C) THE 4-NORM AND (D) THE ANALYTIC CENTRE. OBSERVE THAT THEIR SHAPE INFLUENCES THE SOLUTION OF THE OPTIMIZATION PROBLEM WHEN INTERSECTED WITH THE AFFINE SPACE EQ. (1).

FIGURE 3 :

 3 FIGURE 3: CASE STUDY: (A) PLANAR CDPR WITH FOUR CABLES AND A POINT-MASS LOAD WHILE (B) SHOWS THE WFW OF THE CDPR. OBSERVE THAT THE WFW IS FOUND BY USING THE CAPACITY MARGIN µ [24], THE PURPLE CURVE µ = 0 DELIMITS ITS BOUNDARY.

FIGURE 4 :

 4 FIGURE 4: SENSITIVITY: σ-MAP WITHIN THE WFW FOR THE LP METHOD (A), THE QP METHOD (B), THE 4-NORM (C) AND THE ANALYTIC CENTRE ONES (D).

Only planar CDPRs with point-mass load are considered in this work.

By definition, a TDA is said to be robust if the provided solution is far from the boundary of Π.

Observe that computing the derivative of the Lagrangian function w.r.t. 𝝉 resembles in computing the second-order gradient of the objective directly.

ACKNOWLEDGMENTS

This work has been partially supported by ROBOTEX 2.0 (Grants ROBOTEX ANR-10-EQPX-44-01 and TIRREX ANR-21-ESRE-0015).

FIGURE 5: GRADIENT: ∇σ-MAP WITHIN THE WFW FOR THE LP METHOD (A) THE QP METHOD (B) THE 4-NORM (C) AND THE ANALYTIC CENTRE ONES (D). OBSERVE THAT THE SIZE OF THE ARROWS IDENTIFIES

THE MAGNITUDE OF THE GRADIENT. and Information in Engineering Conference (2005)DOI 10.1115/DETC2005-84324. [5] Merlet, J. and Daney, D. "A portable, modular parallel wire crane for rescue operations." IEEE International Conference on Robotics and Automation (2010)DOI 10.1109/ROBOT.2010.5509299. [6] Kawamura, Sadao, Choe, Won, Tanaka, Satoshi and Kino, Hitoshi. "Development of an Ultrahigh Speed Robot FAL-CON Using Parallel Wire Drive Systems." Journal of the Robotics Society of Japan (1997). [7] Métillon, Marceau, Charron, Camilo, Subrin, Kévin and Caro, Stéphane. "Performance and interaction quality variations of a collaborative Cable-Driven Parallel Robot." Mechatronics (2022)DOI https://doi.org/10.1016/j.mechatronics.2022.102839. [8] Borgstrom, Per Henrik, Jordan, Brett L., Sukhatme, Gaurav S., Batalin, Maxim A. and Kaiser, William J. "Rapid Computation of Optimally Safe Tension Distributions for Parallel Cable-Driven Robots." IEEE Transactions on Robotics (2009)DOI 10.1109/TRO.2009.2032957. [9] Gouttefarde, Marc, Lamaury, Johann, Reichert, Christopher and Bruckmann, Tobias. "A Versatile Tension Distribution Algorithm for 𝑛-DOF Parallel Robots Driven by