Simon Chesne Lamcos

Indirect boundary force measurements in beam-like structures using a derivative estimator

Keywords: Numerical differentiation, Boundary identification, Indirect Force measurement

Indirect boundary force measurements in beam-like structures using a derivative estimator

Introduction

Source characterization is a major issue in the analysis of structural vibrations. In many vibro-acoustic problems, structures are excited by their boundaries. The identification of boundary forces (such as shear forces or bending moments in beam-like structures) is not straightforward, because these physical quantities are not readily accessible and sensors can rarely be positioned at suitable locations. From a more general point of view, the identification of these forces/sources provides information related to the type of vibration generated. They make it possible to: analyse stress inside the structure, describe how vibration energy propagates, through the use of structural intensity calculations [START_REF] Pavic | Measurement of structure borne wave intensity, part i: Formulation of the methods[END_REF], [START_REF] Noiseux | Measurement of power flow in uniform beams and plates[END_REF], and localize external forces applied to the structure.

Force localisation and estimation is a wide subject. A lot approaches exist. These approaches include transfer-function-based methods [START_REF] Leclere | Indirect measurement of main bearing loads in an operating diesel engine[END_REF], numerical model-based methods [START_REF] Parloo | Force identification by means of in-operation modal models[END_REF], modal analysis [START_REF] Maes | Identification of axial forces in beam members by local vibration measurements[END_REF], methods using discretised differential motion equation [START_REF] Pezerat | Two inverse methods for localisation of external source exciting a beam[END_REF], [START_REF] Zhang | Measuring the structural intensity and force distribution in plates[END_REF], [START_REF] Pezerat | Force analysis technique: Reconstruction of force distribution on plates[END_REF], [START_REF] Leclère | Vibration source identification using corrected finite difference schemes[END_REF], or using an operator obtained with the Finite Element Method [START_REF] Renzi | Vibratory source identification by using the finite element model of a subdomain of a flexural beam[END_REF], etc... An interesting review of certain methods is proposed in [START_REF] Sanchez | Review of force reconstruction techniques[END_REF]. Even if many force identification methods exist, few allow to identify boundary forces. Indeed, methods using discretized differential motion equation, which are very successful in a lot of applications, cannot be applied at structure boundaries. This is due to the fact that it is difficult to approximate the derivatives of a quantity, when this quantity is defined on only one side of a point of interest. The techniques described above are therefore unsatisfactory. The use of the annihilator operator, as proposed in this work, avoids this drawback during the derivative estimation step, which is a necessary step for this kind of approach. Other approaches as transfer-function-based [START_REF] Leclere | Indirect measurement of main bearing loads in an operating diesel engine[END_REF] or transmissibility-function-based methods [START_REF] Lage | Force identification using the concept of displacement transmissibility[END_REF] work well if a measurement is located at the point where the force has to be estimated. This can be an important drawback. Indeed, the displacements measured at a structure's boundaries contain the highest levels of error, because the structure's absolute displacements are often very small at these locations, i.e. the signal-to-noise ratio is relatively low. Measurement errors in the vicinity of a structure's boundaries generally lead to large variations in the estimation of the boundary forces.

The goal of the technique proposed in the present study is to identify shear forces and/or bending moments at any point in a structure, in particular at its edges, where such forces can in some cases be associated with the sources. Few force identification methods focus on the boundaries. In the literature, boundary condition identification is a topic that is usually treated independently to sources identification.

Many methods for boundary conditions identification use Finite Element Approach. The general procedure is to apply idealized constraints or to make use of previous experience in modelling apparently similar boundaries. Boundary conditions may be determined by using conventional identification or updating methods [START_REF] Mottershead | Model updating in structural dynamics: A survey[END_REF]. In many cases, for this kind of inverse problems, measurement errors and ill-conditioning lead to unsatisfactory estimates of the spatial model by the equation error method [START_REF] Ahmadian | Regularisation methods for finite element model updating[END_REF], [START_REF] Friswell | Combining subset section and parameter constraints in model updating[END_REF], [START_REF] Majkut | Identification of beams boundary conditions in ill-posed problem from measured displacements[END_REF]. Indeed measurement errors at the boundaries produce large variations in estimates of boundary parameters. One of the most popular methods in the determination of boundary parameters is the sensitivity method [START_REF] Ghandhi | Structural optimisation with frequency constraint: a review[END_REF]; the difference between model predictions and test observations is defined using linearised first order sensitivities. Another approach [START_REF] Ahmadian | Boundary condition identification by solving characteristic equations[END_REF] considers a method for the identification of structural boundary conditions which relies on the measured natural frequencies of the structure. From an another point of view, identification of boundary condition can be linked to the estimation of the boundary spatial derivatives.

If the spatial derivatives of the structural displacements are known, they can be used to determine these shear forces and/or bending moments. One derivative approximation technique involves the use of finite difference schemes [START_REF] Pezerat | Two inverse methods for localisation of external source exciting a beam[END_REF], [START_REF] Pezerat | Force analysis technique: Reconstruction of force distribution on plates[END_REF], spatial Fourier transforms [START_REF] Zhang | Measuring the structural intensity and force distribution in plates[END_REF], and modal expansions [START_REF] Gavric | Finite element method for computation of structural intensity by normal mode approach[END_REF]. The major difficulty with this type of approach is that derivation amplifies the influence of noise in the recorded data, and the higher the order of the derivative, the greater the sensitivity of spatial derivatives to measurement errors. It is thus not possible to use such techniques, since they are, and must include, regularization developments, such as wave number filtering [START_REF] Pezerat | Force analysis technique: Reconstruction of force distribution on plates[END_REF], modal truncation [START_REF] Gavric | Finite element method for computation of structural intensity by normal mode approach[END_REF], regularization techniques [START_REF] Renzi | Vibratory source identification by using the finite element model of a subdomain of a flexural beam[END_REF] which require the adjustment of regularization parameters. To estimate derivatives, some approaches present a natural regularizing aspect using an integral formulation for boundary characterization [START_REF] Chesne | Identification of boundary forces in beams from measured displacements[END_REF][21] [START_REF] Chesne | Distributed piezoelectric sensors for boundary force measurements in euler bernoulli beams[END_REF][23], using polynomial approximation for damage detection [START_REF] Chochol | An original differentiation tool for identification on continuous structures[END_REF], via operational calculus allowing algebraic derivative [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF][25] [START_REF] Fliess | An introduction to algebraic discrete-time linear parametric identification with a concrete application[END_REF] or based on a least square approximation [27][28]. This kind of differentiation by integration methods allows an accurate and robust estimation of the derivative. The approach developed for temporal signals in [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF] is used and adapted in the method detailed in this paper. Since boundary forces are proportional to the spatial derivatives of the structure's displacements, these approaches can be helpful to estimate these forces. But their quantification remains highly sensitive to uncertainties, especially at the boundaries, when post-processing techniques are applied.

To summarise, as for the force identification methods, two major difficulties are encountered when determining these spatial derivatives: they are highly sensitive to measurement uncertainties, and are difficult to measure at the structure's boundaries. The boundaries are often the transmission routes for vibrations, and measurements made at these points are also a major issue in terms of source identification.

It is known that classical characterization methods which use finite difference techniques are not suitable at boundaries. In [START_REF] Chesne | Identification of boundary forces in beams from measured displacements[END_REF], the spatial derivatives of a beam were determined using a weak form of the equation of motion. The principle consists in multiplying the equation of motion by a specific test function, and then integrating the resulting product by means of several successive integrations by parts. This approach allows the shear force (third spatial derivative) and/or the bending moment (second spatial derivative) to be extracted at one boundary of the integration domain. The integration is carried out using discretization, in which the displacements are measured with classical point sensors such as accelerometers, laser vibrometers, etc. The main drawback of this technique is the need to measure a spatial domain, which contains at least one wavelength at the angular frequency under consideration. An extension of this technique to the analysis of a plate structure presenting the same characteristics has also been proposed [START_REF] Chesne | Identification of plate boundary forces from measured displacements[END_REF]. This approach was recently implemented directly, through the use of a distributed PVDF sensor. Since the output of a piezoelectric sensor corresponds to the weighted average of the surface strains in the region covered by the sensor electrodes, discrete integral calculations are not needed [START_REF] Chesne | Distributed piezoelectric sensors for boundary force measurements in euler bernoulli beams[END_REF].

The method presented in the present paper is based on a numerical differentiation technique proposed by Mboup [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF], which was initially developed to handle the case of temporal signals involving a very large number of data-points. This approach has been adapted to the case of structural vibrations, and modified to allow the derivatives to be estimated through the use of iterative integrals of the spatial displacement. The resulting estimation of a structure's spatial derivatives allows its boundary forces to be identified. This novel approach to the physics of the problem has an impact on the method's accuracy. Various limitations, such as the relationship between the size of the integration domain, the vibration wavelengths, sensitivity to noise, and the order of the Taylor expansion on which the approach is based, are identified. The first section of this paper describes the theory of the proposed method. The second section makes use of numerical simulations to evaluate its advantages and drawbacks. The last section describes the experimental validation of this technique, in the case of boundary shear force identification.

Theory

System equations

In this paper, the transverse vibrations of a beam are modelled using the standard Euler-Bernoulli theory for flexural motion. Harmonic oscillations are considered:

v(x, t) = v(x)e jωt (1)
where v(x) is the transverse displacement and ω is the angular frequency.

For the sake of clarity, the time dependence e jωt is simplified in the following. The equation of motion can thus be written as follows:

EI ∂ 4 v ∂x 4 (x) -ρSω 2 v(x) = F (x) ( 2
)
where E is the complex Young's modulus, I is the flexural inertia, ρ is the mass density, S is the cross-section of the beam and F (x) is the excitation. 2014

The boundary forces identified in the present study can be computed using the spatial derivatives of v(x):

T (x) = EI ∂ 3 v(x) ∂x 3 M (x) = EI ∂ 2 v(x) ∂x 2 (3)
where T (x) is the shear force and M (x) is the bending moment.

Estimation of the boundary forces

Exact displacements

It is known from beam theory [START_REF] Guyader | Vibration in Continuous Media[END_REF] that transverse displacements verify Eq.( 2) and can be written as the sum of two propagative and two evanescent bending waves on the left and right sides of the excitation point:

v 1 (x) = A 1 cos(kx) + B 1 sin(kx) + C 1 e -kx + D 1 e kx v 2 (x) = A 2 cos(kx) + B 2 sin(kx) + C 2 e -kx + D 2 e kx (4)
where

v(x) = v 1 (x) for x ∈]0, x f ], v(x) = v 2 (x) for x ∈ [x f , L[), k
is the natural wave number, satisfying the dispersion equation:

k 4 = ρS EI ω 2 (5)
and

A 1 , B 1 , C 1 , D 1 , A 2 , B 2 , C 2
, and D 2 are the wave amplitudes. These must be computed in order to satisfy boundary and continuity conditions at x = 0, x = x f and x = L.

Displacement approximation

The first step in this approach is to consider that in a domain Ω = [0,

x Ω ] (with x Ω < x f ), the exact solutions (4) can be approximated by a truncation of the Taylor series expansion of the displacements. Assuming v(x) to be an analytic function over Ω, it can be approximated by:

v(x) N i=0 ∂ i v ∂x i (0) x i i! ( 6
)
where N is the truncation number. Fig. 1 illustrates this notion and its limitations for the first flexural mode of a simply supported beam, with N=3, N=5 and N=7. For the purposes of simplicity, the modes are normalized such that the maximum amplitude is equal to 1. In the following examples, the proposed method is presented using these 3 truncation values, thereby illustrating the dependence of the computed approximation on this parameter. This figure contains various types of information. As could be expected, the greater the value of N , the more accurate the approximation. The approximation can nevertheless be accurate, whatever the truncation numbers is, if a small domain close to the boundary is considered. Although Ω includes the whole beam in this example, in practice this is not a good choice: it can be difficult to measure physical parameters over the whole structure, and at higher frequencies the Taylor approximation is clearly inaccurate (even with N = 7), as shown in Fig. 2 for the 2 nd vibration mode.

It has been shown in previous studies ( [START_REF] Chesne | Identification of boundary forces in beams from measured displacements[END_REF]) that with this type of integral technique, the key parameter is not simply the angular frequency or the size of the domain Ω, but rather the number of wavelengths contained in the domain Ω.

Since a smaller number of wavelengths is present in Ω when the size of this domain is reduced, the accuracy of the Taylor series approximation is improved. The most well-adapted representations use the number of wavelengths included in Ω, as shown Fig. 3, which do not depend on the size of Ω, nor on ω (the angular frequency of the excitation): under such circumstances, the accuracy of the identification is related directly to the number of wavelengths present in the domain Ω.

Fig. 4 illustrates the difference between the signal and its truncated Taylor expansion, in a domain limited to one half wavelength. The Taylor series approximation can be seen to deviate from the exact solution at λ/4 when N = 3, at 3λ/8 when N = 5, and at λ/2 when N = 7. It illustrates also the interest to use the number of wavelength in the considered domain to plot the results. Indeed the Figs. 3 or 4 can correspond to the whole beam at its first natural frequency (Fig. 1) or to half of the beam at its second natural frequency (Fig. 2).

Derivative estimation

It has been shown that, through the use of an algebraic approach, a Taylor series expansion (6) can be used to numerically differentiate a quantity of interest [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF]. This process can be summarised as follows:

In the Laplace domain, ( 6) can be written:

v N (s) = v(0)/s + v (1) (0)/s 2 + ... + v (N ) (0)/s N (7)
where v (i) (x) is the i th derivative of the beam's displacement at the coordinate x. In order to identify the boundary shear force, the third spatial derivative at x = 0, which is v (3) (0) in the Laplace domain, needs to be estimated. All of the terms v (p) (0) in the expansion, for which p = 3, should thus be rejected. This can be achieved by applying various derivations or multiplications with the Laplace variable s. This operator is referred to as the 'annihilator' [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF]. With standard finite element methods, an Euler-Bernoulli beam element is described by a third order polynomial function.

In order to be coherent with the represented physical phenomenon, and to be able to compute the function's third derivative, a truncation number N ≥ 3 has to be chosen. This method is illustrated in the following example, using N = 3. If both sides of expression [START_REF] Zhang | Measuring the structural intensity and force distribution in plates[END_REF] are multiplied by the annihilator Π, defined as:

Π = s 6 ∂ ∂s s 2 ∂ ∂s s 2 ∂ ∂s s, (8)
(7) then becomes:

6s -4 v(s) + 18s -3 v (1) (s) + 9s -2 v (2) (s) + s -1 v (3) (s) = -6s -8 v (3) (0) (9)
This can be re-written in the spatial domain as:

6 (4) Ω v(x)dx 4 -18 (3)
Ω xv(x)dx 3 + 9

(2)

Ω x 2 v(x)dx 2 (10)
-

Ω

x 3 v(x)dx = -6A 7 /(7!) ∂ 3 v ∂x 3 (0)
where A is the length of the integration domain, defined by Ω = [0; A], and

(α) is the iterated integral of order α.

The Cauchy formula can be used to transform this iterated integral into a simple one, giving:

(α) Ω (-1) β x β v(x)dx α = 1 (α -1)! A 0 (A -x) α-1 (-1) β x β v(x)dx (11)
Expression [START_REF] Renzi | Vibratory source identification by using the finite element model of a subdomain of a flexural beam[END_REF] can thus be written as:

∂ 3 v ∂x 3 (0) = -7! 6A 7 A 0 (A -x) 3 -9(A -x) 2 x + 9(A -x)x 2 -x 3 v(x)dx (12)
The integral expression [START_REF] Lage | Force identification using the concept of displacement transmissibility[END_REF] can be used to estimate the third derivative of the displacements, and has the substantial advantage of being based on integral calculations, thus making it potentially less sensitive to measurement noise.

As previously shown, when the truncation number N used for the Taylor approximation is increased, the accuracy of the approximation v N (x) is improved. The third derivative can by computed using a longer Taylor expansion, and ∂ 3 v N ∂x 3 (0) can be computed as:

∂ 3 v N ∂x 3 (0) = A 0 P N (x)v(x)dx (13)
where P N (x) are the polynomial functions given in Appendix A, for N =3, 5 and 7.

Analytical results

Three main sources of error are naturally present in this method:

• the influence of the truncation number N corresponding to the displacement approximations.

• the influence of discretization of the integral (13), when using discrete measurement points.

• the influence of the noise present in the input data, corresponding to measurement uncertainties. 2014 The aim of the present section is to identify the influence of the truncation number N on identification, in the absence of errors arising from discretization of the integral, or noise on the measured displacements.

F(x f ) Ω=[0,A]
The displacements are computed analytically using Maple software, according to expression (4) for the case of the cantilever beam shown in Fig. 5. The integral of expression ( 13) is also computed analytically. Fig. 6 shows the estimated shear force at the boundary x = 0, computed using three different Taylor series expansions (N = 3, 5 and 7), expressed as a function of the number of wavelengths in the integration domain. Fig. 7 shows (on a logarithmic scale) the corresponding relative identification error, given by the ratio between the identified and exact values.

It should be noted that the same representation is used, whatever the frequency (see [START_REF] Chesne | Identification of boundary forces in beams from measured displacements[END_REF]). By comparing Fig. 4 with Fig. 7, it can be seen that the identification error is directly related to the error in the truncated Taylor series expansion. The identification is correct for A < λ/4 when N = 3, for A < 3λ/8 when N = 5, and for A < λ/2 when N = 7 (A being the size of the integration domain Ω). These values are similar to those observed in Fig. 4 and appear to be appropriate limits for the selection of the order of the polynomial test function P N (x). Although a more accurate estimation is achieved for greater values of N , when the analytical expression is used in the absence of noise, it is shown in the following section that the estimation's robustness with respect to noise follows a different pattern.

It should be noted that neither the structural parameters of the beam, nor the size of the domain Ω, change the behaviour of these errors. The parameters selected for the simulation and used here as reference values are: Steel beam, length: 1m, width: 5cm, height: 1cm, A = 0.95m.

Numerical Simulations

In practice, displacements must be measured at distinct locations, and the integral (13) must be discretized. This is one of the approximations associated with the present method, even when theoretically exact displacements are used. When the displacements are measured, the measurement uncertainties introduce additional errors into the identification.

Numerical integration

Various numerical techniques can be used to discretize an integral. In [START_REF] Chesne | Identification of boundary forces in beams from measured displacements[END_REF] or [START_REF] Chesne | Identification of plate boundary forces from measured displacements[END_REF], the trapezoidal method was studied for similar problems, but was found to be highly sensitive to measurement noise. This technique thus requires too many measurement points to be of any practical interest in the case of the present study. The Gauss-Legendre method is known to be exact when integrating polynomial functions if the discretization uses a sufficiently large number of points. In the general case, the integral approximation can be written:

A 0 f (x)dx A/2 n i=1 W i f (x i ) (14)
where f (x) is the function to be integrated, n is the number of points used, and W i are the weights applied to each of these points x i . In practice, the Gauss-Legendre method has the advantage of producing a better approximation, but requires the implementation of an accurate, irregular mesh. Moreover, this method does not necessitate the knowledge of the function's value at the boundaries. This can be an advantage when it is not possible to make measurements at the boundaries of the beam. It is known that an n-point Gaussian quadrature integral yields an exact result for polynomials of degree 2n -1 or less. This approach appears to be naturally more suitable for the analysis of the present problem. In fact, N is the truncation parameter of the Taylor expansion, but is also the degree of the weight function P N (x). Consequently, Eq.( 13) can be thought of as the integral of a polynomial function of degree 2N . The Gauss-Legendre integration is used in the analyses presented in the following.

Influence of noise

The aim of this section is to demonstrate the influence of noisy displacements (which can be expected when using experimental data) on the reconstructed shear force, and to show that this method is not characterized by the usual instabilities of inverse problems. The case studied here is based on that used in the previous section (Fig. 5), but noise has been added to 2014 the input data. In order to take these errors into account, and to illustrate their influence on the method, the integral ( 13) is computed using noisy displacements. The displacement noise is of the form:

v noisy N (x i ) = v exact N (x i ).∆v m .e j∆φ ( 15
)
where ∆v m is a real random number (with a Gaussian distribution) having a mean value equal to unity and a standard deviation equal to 1/100 x the magnitude of the displacement. ∆φ is another Gaussian real random number having a zero mean value and a standard deviation equal to 1 deg. This type of noise is representative of the measurement accuracy. The simulations use only 8 Gaussian integration points. With this number of points, the Gauss-Legendre integration is theoretically exact for the integration of a 15 degree polynomial function. This number was chosen with respect to Eq.( 13). The highest degree polynomial test function presented in this study is N = 7, also implying the use of a 7 th degree Taylor expansion, and leading to a globally 14 th degree function to be integrated. As a consequence, the errors are not directly related to a problem of discretization. Fig. 8 shows the shear forces at the boundary x = 0, estimated using noisy displacements for three different Taylor series expansions (N = 3, 5 and 7), as a function of the number of wavelengths in the integration domain. Fig. 9 shows, on a logarithmic scale, the resulting relative error of the identification. These figures show that for a small number of wavelengths in the integration domain, the method becomes more sensitive for greater values of N . In the case of a small value of N (N = 3), the polynomial function P 3 (x) is smooth (see annex) and acts like a low pass filter, thereby reducing the influence of noise. When N increases, the polynomial's natural ability to filter noise appears to be shifted towards slightly higher wave numbers. As a consequence, the choice of the correct parameter N and its associated test function P N (x) results from a compromise between a good Taylor series approximation and an adequate filtering capability. In Fig. 9, the vertical dotted lines (set at 0.4 and 0.9 wavelengths) provide an approximate indication of the regions in which it is recommended to use each of the three polynomial functions, in order to obtain a global error below 3dB. These choices were applied to the reconstruction of shear forces in the cantilever beam, using each of the 3 polynomial functions over various ranges. As the limits are not accurately defined, basic weight functions were designed, as shown in Fig. 10. These weights were applied directly to the spectrum of Fig. 8, to reconstruct the new estimation shown in Fig. 11. The simulation is accurate over the observed range of frequencies. It should be noted that for a higher number of wavelengths, the estimation is necessarily inaccurate, since the method is based on the fact that the displacement field can be approximated by a limited Taylor series expansion.

Experimental validation 4.1 Experimental setup

Various experiments were carried out, in order to validate the proposed method with measured displacements, in the case of shear force identifications. A Gauss-Legendre integration was implemented, using 8 points in

Indeterminate boundary condition

Shaker

Measurement area 25cm

Clamped force sensor 12). All displacements were measured with a laser vibrometer and random band excitation was used to study the reconstruction over a wide range of frequencies. As the method is local, the quality of the boundary fixation and location of the shaker (outside the measured domain Ω) had no influence on the results. A force sensor was placed at the lower boundary, between the beam and a rigid support. This measurement, corresponding to the identified boundary shear force, was used as a reference. The beam was made of steel, and its dimensions were: width l = 5cm, height h = 0.5cm.

Results

Displacement field measurement and design of the weight function W

In order to estimate the forces over a large range of frequencies, the first step involves defining frequency limits and appropriate weight functions W. Simulations showed that these limits could be set at approximately 2014 0.4 NoWID and 0.9 NoWID (see Fig. 10). These limits can be found in terms of frequencies when looking at the shape of the displacement fields. At 190Hz and 1100Hz (see Figs. 13 and14) the displacement fields using 8 Gauss points correspond to these criteria. In practice, these limits are not exact, and can vary slightly. Due to the proximity of the boundary, the analytical wavelength (Eq.5) is not repre-sentative of the observed shapes, and this parameter had to be chosen by observing the displacement field as done previously. The cutoff frequencies of the weight functions W were simply computed by adding or removing 10% of the chosen frequency limits. These weight functions are equal to 1 in the frequency range of interest. At the cutoff frequency, they linearly decrease until zero, and are equals to 0.5 when crossing the chosen frequency limits identified in Fig. 13 and14. Fig. 15 shows these 3 weight functions. Figs. [START_REF] Majkut | Identification of beams boundary conditions in ill-posed problem from measured displacements[END_REF] and 17 show the amplitude and phase of the reconstructed shear forces, which are compared with the reference measurements. As predicted by the simulations, the amplitude and phase are very close to the values given by the reference sensor over the observed range of frequencies. The estimations given by independent test functions can be found in the annex. The small differences observed in Figs. [START_REF] Majkut | Identification of beams boundary conditions in ill-posed problem from measured displacements[END_REF] and 17 can be explained by the fact that the reference sensor does not only measure the shear force from the Euler Bernoulli theory. Indeed, due to the presence of this sensor, a point force is created at the location where the Saint Venants principle is not locally respected. Differences also appear at frequencies for which the test function is shifted. These are due mainly to a loss of phase information in the estimated complex frequency response function.

Conclusion and discussions

This paper shows that the identification of boundary forces can be understood in terms of the determination of boundary spatial derivatives of transverse displacements. It is proposed to approximate the displacements Published in Journal of Sound and Vibration, Vol. 333 [START_REF] Chochol | An original differentiation tool for identification on continuous structures[END_REF], 2014, p6438-6452 in the vicinity of the boundary by Taylor series expansions, thus allowing the spatial derivatives to be estimated using specific point-wise derivative estimators. This approach makes it possible to extract these quantities through the use of a weighted spatial integration of the displacement field. Numerical simulations using exact and noisy data have revealed the relationship between the size of the integration domain, wavelength, and noise sensitivity, allowing an appropriate truncation order to be selected. Experimental measurements have demonstrated the feasibility and accuracy of the proposed method.

In comparison with others force identification methods, the proposed approach presents some drawbacks. The use of Gauss points to discretize and to estimate integral (13) can appear as a difficult task. The location of each measurement point has to be as accurate as possible. Depending on the measurement set up, it can be time consuming. But, nowadays, the use of technologies like scanning vibrometer makes this part easy and fast. The point coordinates are automatically generated and the measurement points are very accurately located. Nevertheless the different points are not measured simultaneously and real time estimation of the boundary forces is not allowed with this technology. The main limitation of the proposed approach is the fact that frequency range is directly linked to the size of the measurement domain and the truncation order used for the Taylor expansion estimation. Compared to numerical model-based methods as [START_REF] Parloo | Force identification by means of in-operation modal models[END_REF], the main advantage of the proposed approach is its local aspect. No global measurement is needed, only a few number of point near the boundary is required. The only restriction is that the equation of motion (2) (leading to eqs (3)) has to be verified inside the measurement domain Ω. Another advantage is that there is no particular equation to solve or a regularization step to tune after the measurement step (contrary to [START_REF] Majkut | Identification of beams boundary conditions in ill-posed problem from measured displacements[END_REF] or [START_REF] Ahmadian | Boundary condition identification by solving characteristic equations[END_REF]). Some simple weights have to be applied at each measurement points to obtain directly an estimation of the boundary forces. We remind that these weights are associated to Gauss points or to the weight functions W linked to the truncation order used for the Taylor expansion estimation. From this standpoint the approach is simple to apply. using various test functions without weights. The measured behaviour can be seen to be very similar to that determined by simulation. At low frequencies, a high value of N leads to a greater sensitivity to noise. As the frequency increases, the term N at which the Taylor series expansion is truncated also needs to be increased, in order to maintain the accuracy of the estimation.

C Amplitude reconstruction using 15 Gauss integration points

In order to test the proposed method, an additional experiment was carried out using 15 Gauss integration points. With this number of points, improved reconstruction quality could be expected. Table 1: Locations and the weights of the Gauss points for N=8 on a 25 cm length domain should however be accompanied by measurements characterised by a very low level of noise, since a higher truncation number implies an increase in noise sensitivity.

D Gauss points locations and weight

The tables 1 and 2 give the locations and the weights of the Gauss points for N=8 and N=25, on a 25 cm length domain.

Figure 1 :

 1 Figure 1: Taylor series approximation for beam displacement, at x=0, for the first mode of a simply supported beam: ω = (πL) 2 EI/ρS

Figure 2 :Figure 3 :

 23 Figure 2: Taylor series approximation for beam displacements at x=0, for the second mode of a simply supported beam: ω = (2πL) 2 EI/ρS

Figure 4 :

 4 Figure 4: Error in the Taylor approximation for x = 0, Ω = [0; λ/2], ∀ω

Figure 5 :

 5 Figure 5: An Euler-Bernoulli cantilever beam, Ω represents the measurement domain. x f > A

Figure 6 :

 6 Figure 6: Boundary shear force estimation (analytical computations)

Figure 7 :

 7 Figure 7: Relative error in the boundary shear force estimation (analytical computations)

Figure 8 :

 8 Figure 8: Estimation of the boundary shear force, using various truncation numbers and their corresponding test functions (noisy data)

Figure 9 :

 9 Figure 9: Relative error (dB) on the boundary shear force estimation, using various truncation numbers and their corresponding test functions (noisy data).

Figure 10 : 7 0

 107 Figure 10: Weight functions: W 3, W 5 and W 7

Figure 11 :

 11 Figure 11: Boundary shear force estimation using various truncation numbers and their corresponding test functions (noisy data), following application of the weight functions.

Figure 12 :

 12 Figure 12: Experimental setup

Figure 13 :Figure 14 :

 1314 Figure 13: Measured displacements at 8 Gauss points in the domain Ω, at 190Hz

Figure 15 :

 15 Figure 15: Weight functions W 3(f ), W 5(f ) and W 7(f )

Figure 16 :

 16 Figure16: Amplitude of the reference force sensor and force reconstruction, using 8 measurement points and 3 polynomial test functions with their associated weights.

Figure 17 :

 17 Figure 17: Phase of the reference force sensor and force reconstruction, using 8 measurement points and 3 polynomial test functions with their associated weights.

Figure A. 1 :

 1 Figure A.1: Polynomial test functions P 3 (x), P 5 (x) and P 7 (x)

Figure B. 2 : 7 x 1 x 2 x 3 x 4 x 5 x 6 x 7 x 8

 278 Figure B.2: Amplitude spectra corresponding to the reference force sensor measurements, and the simulated force reconstruction, using 8 measurement points and the test functions P 3 , P 5 and P 7

Figure C. 3 :

 3 Figure C.3: Amplitude spectra corresponding to the reference force sensor and the force reconstruction simulations, made using 15 measurement points and 3 polynomial test functions, associated with weights applied to the high frequency range.

Table 2 :

 2 Locations and the weights of the Gauss points for N=15 on a 25 cm length domain

	Point	x 1	x 2	x 3	x 4	x 5	x 6	x 7	x 8
	Location(cm) 0.0015 0.0078 0.0190 0.0344 0.0536 0.0757 0.0999 0.1250
	Weight	0.0038 0.0088 0.0134 0.0174 0.0208 0.0233 0.0248 0.0253
	Point	x 9	x 10	x 11	x 12	x 13	x 14	x 15
	Location(cm) 0.1501 0.1743 0.1964 0.2156 0.2310 0.2422 0.2485
	Weight	0.0248 0.0233 0.0208 0.0174 0.0134 0.0088 0.0038

Published in Journal of Sound and Vibration, Vol. 333[START_REF] Chochol | An original differentiation tool for identification on continuous structures[END_REF], 2014, p6438-6452

Bibliography References

A Polynomial functions

B Shear force estimation with various test functions


Explaining the slow diffusion of new renewable energy in the Argentine electricity market : a wrong policy mix or an unfavourable context ?

Overview

The production and consumption of electricity in Latin America has grown strongly in recent decades (about 4% per year) with an increasing share coming from fossil fuels, which has led to an increase in the carbon intensity of the electricity production. Large hydro still represents a substantial part of the electricity mix in most Latin-American countries. However, the construction of new dams has slowed mainly due to their local environmental consequences. In the last decade, most of these countries showed a growing interest in developing renewable energy technologies (RETs) for power generation, especially wind, solar, biomass, geothermal and small hydroelectric dams. This interest is explained primarily by the need of diversifying the power mix and increase security of supply. Additionally, other policy objectives have been considered, such as the electrification of isolated rural areas, the decrease of energy imports, the creation of new jobs and the reduction of GHG emissions. The latter goal became especially important after the COP21 (Paris, 2015), in which most countries agreed to follow decarbonisation pathways for their economies which means, among other measures, an increased effort to develop green energy. In this context, governments have set relatively ambitious targets and implemented public policies to encourage investment in RETs and thus take advantage of the great potential available. Different policy instruments have been implemented: tax exemptions, feed-in tariffs, feed-in premium, auction systems, tradable certificates, etc. However, even though many years of government effort and public resources have been invested in order to speed up the development, diffusion and implementation of RETs, experiences in different countries show that this is a very slow process. The current share of RETs is still low (or extremely low depending on the country), especially when compared to the ambitions of policy objectives. Support policies for RETs in Argentina is an interesting case to analyse the effectiveness of incentive mechanisms in a context of high risk perception. Recent experience in the electricity sector shows that the application of several theoretically effectives instruments did not produce the expected results despite the large potential available. Could it be explained by a failure in the design and implementation of the main promotion tools or by one unfavorable economic and institutional context and the related barriers?

Theoretical framework and method

There are two main Economic theories to explain the slow diffusion of RETs. The first is the neo-classical economic paradigm that highlights market failures as the main obstacle. Getting the prices right (by means of taxes, subsidies or similar measures) and public subsidies to compensate for private under-investments in RETs R&D are often proposed solutions. The second group of explanations is based on the evolutionary approaches coupled with a systemic and historical perspective of technological change (Maréchal, 2007). This paradigm highlights the systemic character of innovations and diffusion of new technologies. The technological lock-in process makes it unlikely that traditional cost-efficient measures aimed at internalising external cost will be sufficient to bring about the required radical changes in the field of energy, because they fail to address structural barriers. "Non-economic barriers"institutional, technical, cultural-are thus an important part of the explanation and would require a wider range of policies This paper analyses the interactions between policy tools and barriers from an historic and evolutionary perspective. After a theoretical introduction, it first explains the organization and functioning of the Argentine electricity market. Then, it analyzes the main policy instruments, their implementation and the various barriers that may determine the success or failure of the promotion policy. To achieve this, several step were performed: a deep literature review, an analysis of promotion laws and interviews with key stakeholders (investors, policy makers, energy analysts).

EXPLAINING THE SLOW DIFFUSION OF NEW RENEWABLE ENERGY IN THE ARGENTINE ELECTRICITY MARKET: A WRONG POLICY MIX OR AN UNFAVOURABLE CONTEXT?

Bersalli Germán, University of Grenoble-Alpes, +33 4 56 52 85 31, german-ariel.bersalli@univ-grenoble-alpes.fr

Results

Three main policy instruments have been implemented in the Argentine electricity sector to promote RETs: -The first policy tools for wind and solar energy were approved in 1998 (Law 25.019). It was a feed-in premium system that guaranteed an additional payment above the market price for each kWh of energy produced and injected into the national interconnected grid. The law also provided some tax exemptions for wind and solar projects. Such a system failed due to the freeze on electricity tariffs that occurred after the deep devaluation in 2002. In fact, the imbalance generated by the frozen tariffs facing rapidly growing inflation misrepresented the basis of calculation of the grant. Thus, it provoked an insignificant payment and consequently the inapplicability of this incentive policy.

-In 2006 a new law on renewable energy was approved. It has set a target of 8% of electricity production from renewable sources (excluding large hydro) for 2016. To achieve this target, the new law provides incentives including tax exemptions and the payment of premium remuneration. In contrast to the regulatory approach adopted by Chile, the target of 8% in the Argentine law is only indicative, as "obligation" will fall on any particular actor. Moreover, regulation and law-enforcement only came in 2009, without real results.

-In 2010 the GENREN auction program was implemented. From this program, competitive auctions to incorporate 1000 MW of renewable power to the national electricity system were organised. In the first tender, projects representing 1,436 MW of total power were proposed, that is to say, an amount 40% higher than requested, which shows the initial great interest of investors in the program. However, most of these projects have yet to be realised.

The main obstacle encountered by companies responsible for the projects is the lack of access to financing. First, at the macro-economic level, the international funding of companies is directly affected by the high level of country risk, making it almost impossible to access credit at reasonable rates. Then, most of the participating companies do not have access to the guarantees required by national banks. CAMMESA -a company controlled by the State and responsible for payments of guaranteed prices -presents a balance sheet that does not meet the guarantees required by financial institutions concerning the completion of contracts signed in the field of GENREN project. The lack of confidence has settled despite the existence of guaranteed prices in American dollars with contracts up to 15 years. Thus, only firms with their own funds or with access to internal funding have materialised the renewable projects.

Conclusions

In the Argentinean case, the RET's deployment in the last 20 years was clearly affected by the institutional and economic context -with a big crisis in 2002 and the consequent public control of electricity tariffs. In this high-risk environment, investors have preferred to invest in more conventional and mature technologies or not invest at all. Moreover, we found several problems in the design and implementation of the two first policy instruments (feed-in premium system of 1998 and 2006). The GENREN program was a well-designed economic instrument but it was not accompanied by additional measures to facilitate the funding of projects, i.e. the promotion policies did not respond to the financial failure. More generally, the results show that the malfunctioning of energy markets and also macroeconomic and institutional barriers may prove decisive for the success or failure of promotion policies. The isolated application of economic incentives (like feed-in tariffs, feed-in premium, etc.) may not be enough to promote investment in RETs, especially in developing and emerging countries with a high risk environment. Complementary measures are required: short, medium and long run energy planning, facilitate funding access, facilitate the grid integration of green electricity and education and information diffusion to increase awareness and citizens' participation in RETs projects.

Main references

  -IRENA, 2015. Renewable Energy in Lain America 2015: An Overview of Policies. -Jaffe, A.B., Newell, R.G., Stavins, R.N., 2002. Environmental Policy and Technological Change. Environment and Resource Economics 22, 41-70. doi:10.1023/A:1015519401088 -Kennedy, M., Basu, B., 2013. Overcoming barriers to low carbon technology transfer and deployment: An exploration of the impact of projects in developing and emerging economies. Renewable and Sustainable Energy Reviews 26, 685-693. doi:10.1016/j.rser.2013.05.071 -Maréchal, K., 2007. The economics of climate change and the change of climate in economics. Energy Policy, 35(10), 5181-5194 -Negro, S.O., Alkemade, F., Hekkert, M.P., 2012. Why does renewable energy diffuse so slowly? A review of innovation system problems. Renewable and Sustainable Energy Reviews 16, 3836-3846. doi:10.1016/j.rser.2012.03.043 -Recalde, M. Y., Bouille, D. H., y Girardin, L. O., 2015. Limitaciones para el desarrollo de las energías renovables en Argentina. Problemas del Desarrollo, 46(183), 89-115.


Gautier Laurent

Iterative Thickness Regularization Of Stratigraphic Layers In Discrete Implicit Modeling

Keywords: Structural Modeling, Discrete scheme, Interpolation, Mesh, Norm constraint, Iterative process

Discrete implicit modeling consists in representing structural surfaces as isovalues of threedimensional piece-wise linear scalar fields, which are interpolated from available data points. Data are expressed as local constraints that can enforce the value of the scalar fields as well as their gradients. This paper illustrates some limitations of published discrete implicit methods, related to the difficulty of controlling the norm of scalar field gradient and its evolution over the interpolated domain. It is shown that important artifacts may arise due to the intrinsic dependence between variations in the norm and the direction of the scalar field gradient, from one element to its neighbors. Evidence that these artifacts are related to mesh facet direction with respect to gradient direction are given. The artifacts lead to rapid and uncontrolled variations of thickness that may induce erroneous interpolations. This paper proposes two original approaches to overcome these problems. The first one consists in iteratively adjusting the norm of scalar field gradients in the direction obtained after previous iterations. The second solution consists in optimizing the mesh used by the interpolation. This requires finding appropriate mesh facet orientation with respect to scalar field gradient. These methods demonstrate that the results of discrete implicit surface interpolation can be improved and calls for further development of available interpolation schemes.

Introduction

The representation of geological interfaces as equipotential surfaces of one or several scalar fields, often termed implicit modeling, has become a powerful approach to build geological models [Lajaunie et al., 1997;[START_REF] Cowan | Practical implicit geological modelling[END_REF][START_REF] Chilès | Modelling the geometry of geological units and its uncertainty in 3D from structural data: the potential-field method[END_REF][START_REF] Frank | 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF][START_REF] Calcagno | Geological modelling from field data and geological knowledge Part I. Modelling method coupling 3D potential-field interpolation and geological rules[END_REF][START_REF] Hjelle | A numerical framework for modeling folds in structural geology[END_REF][START_REF] Caumon | 3D implicit stratigraphic model building from remote sensing data on tetrahedral meshes: theory and application to a regional model of La Popa Basin, NE Mexico[END_REF][START_REF] Hillier | Three-dimensional modelling of geological surfaces using generalized interpolation with Radial Basis Functions[END_REF][START_REF] Jessell | Foliation fields and 3D cartography in geology: principles of a method based on potential interpolation[END_REF]. This approach is particularly interesting when representing continuous stratigraphic sequence, because it allows for the representation of a conformable continuum of surfaces at once [START_REF] Jessell | Foliation fields and 3D cartography in geology: principles of a method based on potential interpolation[END_REF]. It can also be useful to represent structural elements, such as foliations and axial surface fields [Lajaunie et al., 1997;[START_REF] Maxelon | A workflow to facilitate threedimensional geometrical modelling of complex poly-deformed geological units[END_REF][START_REF] Massiot | Accounting for axial directions, cleavages and folding style during 3D structural modeling[END_REF][START_REF] Laurent | Folding and poly-deformation modelling in implicit modelling approach[END_REF]. Lastly, implicit surface representation can be used as a way to parameterize geological structures (e.g., fault surface [START_REF] Cherpeau | Method for stochastic inverse modeling of fault geometry and connectivity using flow data[END_REF], fault frame [START_REF] Laurent | A parametric method to model 3D displacements around faults with volumetric vector fields[END_REF], fold frame [START_REF] Laurent | Folding and poly-deformation modelling in implicit modelling approach[END_REF]).

In the context of modeling stratigraphic formations, the gradient of the scalar field can be related to an apparent preserved uncompacted sedimentation rate [START_REF] Mallet | Space-time mathematical framework for sedimentary geology[END_REF][START_REF] Mallet | Elements of mathematical sedimentary geology: the Geochron model[END_REF]. Because thickness changes are key to characterizing basin history and folding style, an important question for implicit methods lies in their ability to accurately capture or predict thickness changes (i.e., to deal with variations of both thickness and norm of the scalar field gradient).

This paper identifies some limitations of the current implicit methods in the management of thickness variations between implicit surfaces, and proposes strategies to address these limitations in discrete methods based on piecewise linear tetrahedra. The basics of discrete implicit modeling are further explained and some undesired effect of the piece-wise linearity are illustrated . Its impact on the thickness of modeled layers is illustrated, demonstrating dramatic artifacts that may be generated when the gradient norm is not sufficiently constrained.

Two approaches are proposed for improving discrete interpolation results for modeling stratigraphic layers. Section 3 presents an adapted gradient norm constraint based on an iterative process. Preliminary results show more limited and better-controlled variations of the thickness of interpolated layers as compared to a classical discrete approach. Section 4 introduces a quantification of the thickness variation with respect to mesh facet orientation. It is shown that using a mesh with optimal orientation is a good way of reducing uncontrolled thickness variations, even without additional norm constraints. Section 5 discusses the advantages and limitations of this approach, which opens new ways to overcome the limitations of published discrete approaches.

2 Implicit Modeling Schemes: Orientation Constraints and Thickness Control

Implicit modeling methods can be divided into two main approaches, depending on the way scalar fields are represented and interpolated from data: polynomial trend based methods, and discrete methods. As gradient and thickness control may be encountered in both methods, they are first briefly presented in this section. This paper then focuses on the discrete approach, for which practical solutions are proposed.

Polynomial Trend Methods: Kriging and RBFs

In implicit methods based on Radial Basis Functions (RBFs) [START_REF] Cowan | Practical implicit geological modelling[END_REF][START_REF] Hillier | Three-dimensional modelling of geological surfaces using generalized interpolation with Radial Basis Functions[END_REF] or Kriging [Lajaunie et al., 1997;[START_REF] Chilès | Modelling the geometry of geological units and its uncertainty in 3D from structural data: the potential-field method[END_REF][START_REF] Calcagno | Geological modelling from field data and geological knowledge Part I. Modelling method coupling 3D potential-field interpolation and geological rules[END_REF], scalar fields are described as a weighted combination of a global polynomial trend and an order-2 stationary residual given by the kriging or RBF weights. Therefore, local thickness variations can be accounted for in the residual, but global thickness changes must be addressed in the polynomial trend.

The instantaneous thickness of strata described by an implicit function ϕ(x, y, z) is given by the inverse of the norm of the gradient. The ability of typical second-order polynomial drift function to approximate general types of thickness variations seems limited. We are not aware of discussions about the impact of this observation in the literature. This suggests that further research may be needed to define trend functions appropriately capturing thickness changes.

Discrete Implicit Modeling and Gradient Norm Variations

In the following, the basics of discrete implicit modeling are briefly presented (Sect. 2.2.1). With this approach, problems of gradient and thickness control come as two main limitations (Sect. 2.2.2): (i) the lack of layer thickness control and (ii) possible rapid variations in gradient norm, referred to as gradient collapse or gradient burst.

Principles of Discrete Implicit Modeling

Discrete implicit modeling represents geological surfaces as isosurfaces of a piece-wise linear scalar field ϕ , interpolated within a set of three-dimensional elements [START_REF] Moyen | 3D-Parameterization of the 3D geological space -The GeoChron model[END_REF][START_REF] Frank | 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF][START_REF] Caumon | 3D implicit stratigraphic model building from remote sensing data on tetrahedral meshes: theory and application to a regional model of La Popa Basin, NE Mexico[END_REF][START_REF] Souche | Volume based modeling-automated construction of complex structural models[END_REF][START_REF] Mallet | Elements of mathematical sedimentary geology: the Geochron model[END_REF][START_REF] Collon-Drouaillet | 3D geomodelling combining implicit surfaces and Voronoi-based remeshing: A case study in the Lorraine Coal Basin (France)[END_REF]. Linear tetrahedron elements are generally used as they provide a simple linear interpolation in three-dimensional space and can be progressively refined [START_REF] Frank | 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF] and adapted to discontinuities such as faults [START_REF] Caumon | 3D implicit stratigraphic model building from remote sensing data on tetrahedral meshes: theory and application to a regional model of La Popa Basin, NE Mexico[END_REF]. Available data and spatial continuity are expressed as a series of linear equations that are solved to find optimal values for the mesh nodes.

This piecewise discrete approach has a number of interesting features: (i) interpolated model can be locally controlled by element size and shape, which makes it possible to easily generate geometric features based both on data and on prior knowledge. For instance, fold hinges may be localized between data points using dip domains [START_REF] Caumon | 3D implicit stratigraphic model building from remote sensing data on tetrahedral meshes: theory and application to a regional model of La Popa Basin, NE Mexico[END_REF] and geometric features may be added with stochastic simulation [START_REF] Caumon | Towards stochastic time-varying geological modeling[END_REF][START_REF] Mallet | Solid earth modeling and geometric uncertainties[END_REF]. (ii) The way each geological constraint is applied can be locally weighted in each element, for example based on the reliability of the information. (iii) Interpolation complexity mainly depends on the number of elements [START_REF] Frank | 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF], which makes it possible to use large numbers of data and control points.

In each tetrahedron, an interpolated scalar field ϕ and its gradient ∇ ϕ are defined by the following linear equations, where ϕ c is a vector holding the values of the corners of the mesh, M and T are matrices defined by the geometry of the tetrahedron [START_REF] Frank | 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF], and • denotes the dot product.

ϕ (x, y, z) = 1, x, y, z • M • ϕ c (1) ∇ ϕ (x, y, z) = T • ϕ c (2)
These two equations are used to transform each data point into a linear equation [START_REF] Frank | 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF]. A regularization term is added to describe how the scalar field should evolve between the data points. This is implemented by minimizing the variation of ∇ ϕ between neighbor elements. A linear equation is produced for each pair of connected elements (i, j). Because neighbors share a facet, their gradient can only vary in the direction n i j normal to this common facet. For this reason, a minimum variation in gradient is enforced by constraining the projection of the gradient onto vector n i j not to vary between elements i and j.

n T i j • T i • ϕ i -T j • ϕ j = 0 (3)
The resulting set of linear equations is gathered into a linear system. This system is solved in the least square sense, which yields a scalar field that balances all data and smoothness constraints.

Classical constraints (Fig. 1(a)), derived from Eq. (1) and Eq. ( 2), can enforce: (i) the value of ϕ at observed location; (ii) the components of ∇ ϕ at observed location; (iii) the variation of ∇ ϕ between given elements. Additional constraints can be developed to improve the control of this scalar field and introduce new information into this basic system [START_REF] Frank | 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF][START_REF] Massiot | Accounting for axial directions, cleavages and folding style during 3D structural modeling[END_REF]. For example, it is possible to only constrain the direction of ∇ ϕ without fully specifying the values of its components. This can be implemented by using the property that the dot product of two orthogonal vectors is null, together with Eq. (2). The following two constraints are then written to enforce the observed direction w at a given location as a direction for ∇ ϕ , using any two non parallel vectors u and v, orthogonal to w.

u T • T • ϕ c = 0 (4) v T • T • ϕ c = 0 (5)
Alternatively, the direction of ∇ ϕ may be controlled by another scalar field [START_REF] Massiot | Accounting for axial directions, cleavages and folding style during 3D structural modeling[END_REF][START_REF] Laurent | A parametric method to model 3D displacements around faults with volumetric vector fields[END_REF], by constraining ∇ ϕ to be orthogonal to the controlling scalar field gradient in each tetrahedron of a given region (Fig. 5(a)).

Intrinsic Limitations of the Gradient Norm Control

The core of the discrete interpolation scheme resides in the regularization term, which minimizes the variation of ∇ ϕ between neighbor elements. Its implementation relies on the continuity of the scalar field, which is ensured by the fact that neighbor elements share the corners of their common facet. This relationship between neighbors also introduces limitations in the interpolation. The difference in gradient direction implicitly constrains the norm of the gradient in the second triangle. To honor the difference of values along the shared edge, the second triangle needs values that grow at a faster rate when the angle increases. The norm of its gradient rapidly increases when the relative angle increases. This simple example illustrates how the gradient norm and direction are intimately linked and represents a good example of gradient burst.

Although this configuration appears over-constrained, it illustrates how each shared edge decreases the degree of freedom of the system. This phenomenon causes the loss of a degree of freedom as soon as a circle of connected triangles is formed. It is also interesting to point out with this example that setting an angle of 45 • in the right hand side triangle would not be possible because the variation of value along the shared edge could not be honored.

In more complex examples, this intrinsic link between gradient norm and direction leads to undesired behaviors. Constraints on the dip of an interpolated stratigraphy results in uncontrolled variations in layer thickness. This situation is illustrated with two synthetic case studies (Figs. 3,5,6,and 9). Figure 3 tests the ability of the interpolator to reconstruct a scalar field which represents a series of constant thickness layers bent around the lower left corner of the picture (Fig. 3(a)). This can be seen as an idealized isopach fold, with an axial surface dipping 45 • to the left and with an interlimb angle of 90 • . Two data points that sample the exact value and gradient of the scalar field around the lower left corner of the picture are considered (Fig. 3(b)). The resulting scalar field is very different from the reference field. Figure 3(c) shows an overlay of the reference scalar field isolines (dotted lines) on top of the interpolated scalar field. The norm of the scalar field gradient appears to be smaller than expected, in particular in the top-left and bottom-right corners (Fig. 3(d)), which implies an increase in layers thickness in these regions. Errors in the orientation of the scalar field are also produced, growing from 0 • along the diagonal to 30 • in the top-left and bottom-right corners. These errors illustrate how the interpolation accommodates the variation of direction in the data by varying both the orientation and the thickness of the interpolated scalar field.

In some configurations, the interpolated scalar field might be an appropriate representation of stratigraphic layering. But in general, geological layers would tend to have a relatively constant thickness or show slow, progressive variations. Stratigraphic orientations may vary more rapidly. In this particular example, only the orientation of the reference scalar field is changing whilst the thickness remains constant. This is not captured by the interpolation and results in erroneous orientation and thickness. When too few elements are used, the results present dramatic thickness overestimation. Figure 4(a) shows that increasing the mesh density reduces this effect. This improvement is limited and further increasing the mesh density does not completely remove thickness over-estimation nor orientation errors in the bottom right and top left corners.

Analysis of Sensitivity to Mesh and Data Density

Figure 4(b) illustrates that this effect becomes even more important when the number or orientation data points increases. Using more data points improves interpolated scalar field direction but increases the thickness over-estimation. As shown by Fig. 2, the interpolator can only account for the variation of orientation by increasing layer thickness.

Increasing the weight of the regularization constraint could limit the variation of thickness. Figure 4(c) shows that increasing the weight of the regularization term (i.e., decreasing the data relative weight) limits the variation of orientation. This yields erroneous orientations without succeeding in correcting the thickness errors.

Finally, the only obvious way to really improve the quality of the result is by considering more value data points, but this is not a practical solution as the density of available data is generally not controlled.

Gradient Collapse with Orientation Constraints

Figure 5 presents a second case study, where an additional scalar field ψ is used for constraining the orientation of ϕ. This is implemented by minimizing the dot product between ∇ϕ and ∇ψ. This type of constraint is used by [START_REF] Laurent | A parametric method to model 3D displacements around faults with volumetric vector fields[END_REF], to build a fault frame orthogonal to a given fault surface, and in [START_REF] Massiot | Accounting for axial directions, cleavages and folding style during 3D structural modeling[END_REF] to constrain the fold axis direction of a folded stratigraphy. This situation, where the density of orientation constraints is relatively high compared to value or gradient norm constraints, can result in unexpected behavior of the interpolated scalar field. In this example, the gradient norm rapidly decreases away from the data point located in the center of the model, which is referred to as a gradient

Discussion of Direction and Thickness Constraints

The situations described in the previous examples seem to be generally avoided in classical use of discrete implicit modeling. This can be attributed to the fact that a higher number of value or norm constraints are available to prevent the interpolated gradient norm from varying too much. However, the results presented in this paper suggest that this problem may still arise, to a lesser extent.

Figure 4 illustrates the difficulty of improving interpolation results by modifying the mesh density or relative weights of constraints. It also shows that extracting more orientation control points from available data would not be a practical solution. In the following sections, it is demonstrated that interpolation results might still be improved. Two avenues are proposed and explored to solve the issues of gradient collapse and norm control: an iterative gradient norm constraint (Sect. 3) and a condition of optimal mesh to minimize the effect of gradient lock (Sect. 4).

Iterative Gradient Norm Constraint

Observations from the previous section suggest that defining an additional constraint capable of controlling the norm of the gradient would help to avoid undesired gradient norm variation and to control layers thickness at a given location. This section proposes a strategy to make such constraints possible by defining a linear norm constraint (Sect. 3.1) and by applying it in an iterative process (Sect. 3.2). Results are presented in Sect. 3.3.

Definition of a Gradient Norm Constraint

It is not trivial to specify the norm of the gradient at a given location without constraining its direction. Where the direction of the gradient is not known, the norm would be defined with the following general equation.

∇ ϕ 2 = ∇ ϕ T • ∇ ϕ = ϕ T c • T T • T • ϕ c = 1 L 2 (6)
In this equation, • refers to the norm, and L refers to the thickness of the layer involved. Equation ( 6) is not linear with respect to ϕ c , so it cannot be directly integrated into the linear system. The resolution of such a system would require quadratic matrix systems, which are more complex to solve [START_REF] Schurbet | Quadratic matrix equations[END_REF].

It is proposed to use a more simplistic approach where the gradient norm constraint can still be expressed as a linear equation. If a unit direction vector u is given as an estimate of the gradient direction, Eq. ( 7) provides an approximation of the gradient norm ∇ ϕ .

∇ ϕ = u T • T • ϕ c = 1 L (7)
This time, Eq. ( 7) is linear with respect to ϕ c , and can be integrated into the linear system.

Iterative Estimation of Gradient Direction Vector

Figure 6: Iterative norm constraint process. A first estimation of the gradient direction is obtained with a classical interpolation. A gradient norm constraint is then iteratively applied based on these directions. weight the constraint is progressively increased. A relatively coarse mesh is used to make the orientation of the scalar field in each element clearly visible.

The main difficulty in using Eq. ( 7) is finding a good estimation of the possible gradient direction for each element where the constraint is to be applied. Ideally, this direction would be derived from the result of the interpolation process, which is not possible as the final solution is not known when defining the constraints. To overcome this hurdle, an iterative process is proposed (Fig. 6), where u is defined as the direction of the gradient obtained at the previous iteration. This process is initialized using a first interpolation without norm constraint.

To avoid propagating erroneous gradient direction throughout the iterations, it is proposed to progressively incorporate this constraint into the system instead of applying it all at once. Three possible strategies are considered: (ii) Applying the constraint in a growing number of randomly selected elements. (iii) Applying the constraint within a growing region around data points.

The following case studies show that this process progressively corrects the direction in which the norm is constrained and attenuates undesired norm variations.

Results of Iterative Gradient Norm Control

Case Study 1: Comparison of Iteration Strategies

The iterative norm correction is applied on case study 1, using the same dataset as in Fig. 3. The process is initialized with the direction of the gradient of the scalar field obtained using a classical interpolation. The gradient norm constraint is then progressively introduced at each iteration. The three proposed strategies have been tested (Fig. 7):

(i) The weight progressively grows from 0.01 to 1.

(ii) The percentage of affected elements progressively grows from 1% to 100%.

(iii) The radius in which elements are constrained goes from 0.1 to 28, which covers the whole model.

One hundred iterations have been run in each case. The norm constraint weight is set to 1 in case B and C.

In each case, the weight of the regularization constraint remains 0.1. For each strategy, the interpolation result is progressively corrected and becomes similar to the reference scalar field (Fig. 3(a)). Small variations are still visible, but the interpolated scalar fields show a good general behavior. Both thickness and layer orientation are honored.

It is difficult to compare the three strategies based on this example because they affect different aspect of the interpolation. It can still be observed that progressively increasing the number of constrained elements (randomly or with a growing radius) seems to provide results that are slightly closer to the reference solution. The second case study shows how the iterative gradient norm control can reduce the problems of gradient collapse illustrated in Fig. 5. In this three-dimensional example, an orthogonality constraint is used, in addition to a data point placed in the center of the model with a value of 0 and gradient of [1, 0, 0]. The scalar field used for the orientation constraint is shown in Fig. 8(a).

Case Study 2: Correcting Gradient Collapse

The scalar field obtained by a classical interpolation (Fig. 8(b)) has two issues: (i) the range of values that are obtained is smaller than expected, because the gradient norm decreases rapidly when moving away from the data point, and (ii) the obtained orientations are not correct as they show a bend in the isosurfaces, particularly visible on the vertical faces of the model.

The gradient constraint is iteratively applied with a growing radius (from 1 to 20, over 20 iterations). Both the problems of norm and direction are progressively fixed (Fig. 8(c)), and the final results correspond to the expected straight isosurfaces honoring the orthogonality constraint.

Case Study 3: Folded Layers

Figure 9 presents a third synthetic example showing a cross-section with four data points defining an anticlinal structure. The norm of the gradient of the control points is 5 for each point. The value of the scalar field is constrained to be 0 for the top two data points and is not constrained for the other two. The result obtained with a classical interpolation shows a sudden decrease of the gradient norm in the region between the data points (Fig. 9(a2)). It is due to the variation of the direction of the gradient set on data points around this area. This variation is accommodated by varying both the direction and the gradient norm of the scalar field.

Figure 9: Gradient norm constraint and norm values. (a1) A scalar field is interpolated from 4 data points with scalar field value and gradient constraint (dots and arrows). The scalar field locally honors the data point direction and gradient norm around the data points, but the gradient norm decreases in between the data points (a2). This results in an increase in the thickness of interpolated layers. (b1 to d1) Iterative gradient norm constraint is applied with an increasing weight. (b2 to d2) This process progressively makes the gradient norm values closer to the data point norm (data gradient norm: 5) When applying a gradient norm correction (Figs.9(b) to (d)), the norm of the gradient of the scalar field gets closer to the control point values. This also has consequences on the resulting geometry of the interpolated fold.

Optimization of Shared Facet Orientation

Section 3 introduces additional constraints that can be used to improve the scalar field interpolation. This section proposes another approach which focuses on the mesh used for the interpolation.

Figure 2 demonstrates that the thickness variation between neighbor elements is controlled by their relative gradient direction. When varying the orientation of their common facet, it appears that the effect of thickness variation can be limited by finding a better orientation. A first step towards a mesh improve-ment would be to quantify this relationship. To that end, two neighbor triangles are considered, where the orientation of the scalar field is known (Fig. 10).

Figure 10: Thickness variation in neighbor elements. (a) Two neighbor triangles t 1 and t 2 containing an interpolated layer (dashed lines). (b) A blow-up around the shared facet which illustrates the layer thickness L and orientations used in Eq. ( 9) It is possible to find an optimal orientation for the shared facet knowing the angle δ between the two gradient directions. In each triangle, the orientation of the shared facet is defined by its angle with the gradient direction, respectively denoted α 1 and α 2 . These angles are linked by the following relationship: α 2 = δα 1 . The apparent thickness L i along the shared facet is obtained by projecting the actual thickness, L 1 and L 2 , onto the facet direction.

L i = L 1 cos α 1 = L 2 cos α 2 (8)
The thickness ratio L 2 /L 1 is then expressed by the following equation.

L 2 /L 1 = cos α 2 cos α 1 = cos (δ -α 1 ) cos α 1 (9)
It appears that L 2 can be the same as L 1 in two situations. A first trivial solution comes for δ equals 0, because cos (-α 1 ) equals cos α 1 . The ratio L 2 /L 1 is then equal to 1.

Another set of solutions is obtained with the following condition.

δ -α 1 = α 1 + 2kπ , k ∈ N (10) ⇐⇒ α 1 = δ/2 + kπ , k ∈ N (11)
These results suggest that another avenue for improving discrete scalar field interpolation would be to use an optimal mesh where each facet honors one of these conditions: (i) the direction of the gradient is the same on both sides of the facet, or (ii) the facet bisects the angle between the gradient directions on both sides.

In very simple cases, it is possible to define the parameters of such a mesh. For example with the data configuration used in Fig. 3 and Fig. 7, the ideal facets would converge towards the bottom left corner of the picture. Figure 11 shows that an interpolation based on such a mesh is able to progressively vary the direction of the scalar field gradient without affecting the norm. This suggests that, in more general cases, a better mesh could result in a better interpolation without unexpected gradient norm variations.

However, it is not trivial to define an algorithm to build a mesh based on these characteristics in the general case. First, it depends on the direction of the interpolated gradient, which is part of the solution of the interpolation. This suggests that an iterative process could be used, as in Sect. 3.

Another limitation comes from the possible contradiction between the ideal position of neighbor facets. Figure 12 presents a slightly more complex mesh with two layers of triangles. Here, only one facet of each triangle can be oriented in the ideal position which lowers the quality of the result. The isolines of the b)). This effect becomes limited when using a higher number of elements because the orientation of the limiting facets become very close to the optimal position (Fig. 12(c)).

An implementation of this criterion could be developed in an optimization process with a cost function based on Eq. ( 9).

Limitations and Discussion

The examples presented above show interesting results but some caveats remain. For the iterative gradient control process, the constraints presented in Sect. 3.1 require that the gradient norm value is known beforehand and remains relatively constant. This would probably not be possible in the general case. To overcome this hurdle, the gradient norm could be estimated from the data and interpolated in the model. Another approach would be to constrain the variation of norm between neighbor elements considering the gradient direction obtained at previous iteration. Adapting the weight of the norm constraint is also a way to control how constant the resulting norm would be. For example, a moderate weight would enforce a globally constant norm while allowing some local variations (e.g., in Fig. 9(b)). This method is also limited by a lack of theoretical proof of convergence. In particular, the system may propagate erroneous orientation generated at the initialization step. For example, in Fig. 13, this approach is applied to the result obtained in Fig. 5, where a strong gradient collapse can be observed. The norm constraint succeeds in limiting the norm variations but propagates erroneous initial orientations. This result is in part due to the complex orientation constraints that are applied (Fig. 5(a)), but also results from the wrong general orientation of the initial interpolation which is kept throughout the iterations. In this case, the norm constraint is not sufficient to correct the gradient direction. This problem suggests that this process may not be robust enough in its current form to be applied in any case study, even if it proves useful for correcting some interpolation results as shown above.

A method to detect and correct erroneous directions used in the norm correction process needs to be developed. A solution could be to use local weighting around the data points, as suggested by the positive results obtained using the growing radius strategy (Fig. 6(c)). This could lead to a solution combining discrete interpolation scheme and radial basis functions (RBFs). Alternatively, a process to independently interpolate gradient directions would be another way to progressively correct these misleading initial orientations.

Conclusion

This paper addresses the limitations of the discrete implicit modeling method related to the norm of the scalar field gradient. The first contribution is to clearly identify and illustrate these limitations with simple examples. They demonstrate that the norm and the direction of the gradient of an interpolated scalar field are closely related. Classical interpolation is not able to independently control these two aspects of scalar field spatial variation.

This paper presents two avenues for improving discrete interpolation schemes: (i) An iterative gradient constraint, with different strategies to progressively apply this constraint. (ii) A mesh optimization criterion for minimizing thickness variation between neighbor elements by finding optimal facet orientation. The preliminary results obtained with both methods are encouraging and suggest that the proposed method may help improve discrete interpolation.

The limitations of the proposed methods are also shown. Its main limitation is its potential difficulty to converge to an appropriate solution in some circumstances (Fig. 13). These limitations could make this approach difficult to apply to more complex case studies in its current form. However, using an iterative scheme proves useful for defining norm constraints and progressively correcting what the classical discrete approach may fail to interpolate in only one first iteration. The mesh considerations presented in this paper should lead to more work on mesh definition for interpolation purposes.

The solutions presented in this paper should not be taken as mature and practical, yet. They clearly demonstrate that it is possible to improve current interpolation schemes available for geomodeling. Two promising directions have been proposed for further research: (i) iterative regularization and (ii) mesh optimization.

Figure 1 :

 1 Figure 1: Discrete two-dimensional interpolation of a scalar field representing stratigraphy. (a) mesh (a set of connected triangles) and data used to define the interpolation. The values imposed on the two isolines are 0 (left hand side line) and 1 (right hand side line). (b) interpolated scalar field

Figure 2 :

 2 Figure2: Gradient norm variation with varying direction constraints. The left hand side triangle has a constrained gradient direction and norm. In the other triangle, only the gradient direction is controlled. An increasing angle is set from 0 • (a) to 44 • (c). As a result, the norm of the gradient in the second triangle increases from (a) to (c), which cannot be controlled in this particular setting Figure2presents a basic two-dimensional example of neighborhood relationship, considering two triangular elements. This example demonstrates how the gradient in both triangles is interdependent because of the shared corners. The left hand side triangle is constrained with a fixed vertical unit gradient, with fixed corner values. The gradient in the other triangle is constrained with an angle from 0 • (Fig.2(a)), to 30 • (Fig. 2(b)), and 44 • (Fig. 2(c)).The difference in gradient direction implicitly constrains the norm of the gradient in the second triangle. To honor the difference of values along the shared edge, the second triangle needs values that grow at a faster rate when the angle increases. The norm of its gradient rapidly increases when the relative angle increases. This simple example illustrates how the gradient norm and direction are intimately linked and represents a good example of gradient burst.Although this configuration appears over-constrained, it illustrates how each shared edge decreases the degree of freedom of the system. This phenomenon causes the loss of a degree of freedom as soon as a circle of connected triangles is formed. It is also interesting to point out with this example that setting an angle of 45 • in the right hand side triangle would not be possible because the variation of value along the shared edge could not be honored.In more complex examples, this intrinsic link between gradient norm and direction leads to undesired behaviors. Constraints on the dip of an interpolated stratigraphy results in uncontrolled variations in layer thickness. This situation is illustrated with two synthetic case studies(Figs. 3, 5, 6, and 9).

Figure 3 :

 3 Figure 3: Case study 1: interpolation of a circular scalar field. (a) The reference scalar field, with the mesh used for the interpolation. (b) Two samples are taken and used as data for the interpolation, which constrain both the value (dot) and the gradient (vector) of the interpolated scalar field. (c) The resulting scalar field (plain lines and shading) differs from the reference field (dotted lines), showing errors in both the norm (d) and the direction of the scalar field (e)

Figure 4

 4 Figure 4 illustrates the sensitivity of the interpolation to: (a) mesh density, (b) density of orientation data, and (c) relative weight of the regularization term and data constraints. None of these three factors suc-

Figure 4 :

 4 Figure 4: Analysis of sensitivity to (a) mesh refinement, (b) orientation data density, and (c) relative weight of regularization term. The reference scalar field is shown for comparison on the right hand side. Its isovalues are also shown on each realization (dotted lines)

Figure 5 :

 5 Figure 5: Case study 2: A case of gradient collapse observed when interpolating using orientation constraints. (a) The scalar field used as a direction constraint. (b) Interpolation result. Other constraints are a data point in the center of the model with value 0 and a gradient vector [1, 0, 0], parallel to the orientation scalar field. The interpolated scalar field shows relatively good orientation and variation above and under the data point, but the gradient norm rapidly decreases in other directions, which also perturbs the gradient direction. (c) A front view of the resulting interpolation

Figure 7 :

 7 Figure7: Application of the iterative gradient norm constraint to case study 1 (Fig.3). Isolines of the reference scalar field are shown (dotted lines) for comparison with the interpolated scalar field (plain lines). The norm constraint is applied: (a) with an increasing weight: 0.02 (a1), 0.1 (a2) and 1 (a3); (b) with an increasing number of elements: 10% (b1), 50% (b2) and 100% (b3); (c) within a growing radius from the data: 1 (c1), 7 (c2) and 28 (c3)

Figure 8 :

 8 Figure 8: Application of the iterative gradient norm constraint to case study 2. (a) A scalar field is used as an orientation constraint as in Fig. 5. (b) The initial scalar field obtained with a value and gradient control point located in the center of the volume, and a constraint of orthogonality to the scalar field shown in (a). (c) The iterative process succeeds in progressively correcting the gradient norm and direction. The constraint is applied in a growing radius around the data point. Results are shown for a radius of 5 (c1), 10 (c2) and 20 (c3)

Figure 11 :

 11 Figure 11: Optimization of the interpolation support. The same dataset as in Fig. 3 and Fig. 7 is used. The isolines of the reference solution are shown as dotted lines. A mesh with respectively 2, 4, 8 and 20 elements is used in examples (a) to (d). Good results are obtained with few elements, without using the gradient norm constraints from Sect. 3

Figure 12 :

 12 Figure 12: Interpolation on a mesh made of two layers of triangles. Here, only one of the facet of each triangle is placed in the optimal orientation. (a) As a result, the interpolated scalar field does not fit the reference isolines (dotted lines). (b and c) Increasing the number of triangles generates facets whose orientation is closer to the theoretical optimum, which improves the interpolation result

Figure 13 :

 13 Figure 13: An attempt to correct gradient collapse in Fig. 5. (a) Initial interpolation: the orientation of the scalar field is correct in the center part (i.e., vertical), but becomes inclined when moving away from the data point. (b to d). Results after a few iterations of the norm correction process with an increasing weight. The orientation is first improved together with the norm (b) but further iterations fail to converge towards the data point orientation (c). Increasing the weight of the norm constraint is not enough to correct the orientation (d)

Acknowledgments

The author would like to gratefully thank anonymous reviewers, as well as Professor Guillaume Caumon and Lachlan Grose, for their very constructive comments and corrections, which helped improve this paper. This work was initiated at Monash University with the Australian Research Council Discovery grant DP110102531 and further developed at GeoRessources (http://ring.georessources. univ-lorraine.fr/) in the framework of the "Investissements d'avenir" Labex RESSOURCES21 (ANR-10-LABX-21).